[BACK]Return to dalg.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Annotation of OpenXM_contrib2/asir2000/engine/dalg.c, Revision 1.13

1.1       noro        1: /*
1.13    ! noro        2:  * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.12 2005/10/12 14:43:36 noro Exp $
1.1       noro        3: */
                      4:
                      5: #include "ca.h"
                      6: #include "base.h"
                      7:
                      8: static NumberField current_numberfield;
                      9: extern struct order_spec *dp_current_spec;
1.2       noro       10: void simpdalg(DAlg da,DAlg *r);
1.8       noro       11: int invdalg(DAlg a,DAlg *c);
1.3       noro       12: void rmcontdalg(DAlg a, DAlg *c);
1.6       noro       13: void algtodalg(Alg a,DAlg *r);
                     14: void dalgtoalg(DAlg da,Alg *r);
1.1       noro       15:
1.5       noro       16: NumberField get_numberfield()
                     17: {
                     18:        return current_numberfield;
                     19: }
                     20:
1.1       noro       21: void setfield_dalg(NODE alist)
                     22: {
                     23:        NumberField nf;
                     24:        VL vl,vl1,vl2;
                     25:        int n,i,dim;
                     26:        Alg *gen;
                     27:        P *defpoly;
                     28:        P p;
                     29:        Q c,iq,two;
                     30:        DP *ps,*mb;
1.3       noro       31:        DP one;
1.1       noro       32:        NODE t,b,b1,b2,hlist,mblist;
                     33:        struct order_spec *current_spec;
                     34:
                     35:        nf = (NumberField)MALLOC(sizeof(struct oNumberField));
                     36:        current_numberfield = nf;
                     37:        vl = 0;
                     38:        for ( t = alist; t; t = NEXT(t) ) {
                     39:                clctalg(BDY((Alg)BDY(t)),&vl1);
                     40:                mergev(ALG,vl,vl1,&vl2); vl = vl2;
                     41:        }
                     42:        for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ );
                     43:        nf->n = n;
                     44:        nf->vl = vl;
                     45:        nf->defpoly = defpoly = (P *)MALLOC(n*sizeof(P));
                     46:        nf->ps = ps = (DP *)MALLOC(n*sizeof(DP));
                     47:        current_spec = dp_current_spec;
                     48:        STOQ(2,two);
                     49:        create_order_spec(0,(Obj)two,&nf->spec);
                     50:        initd(nf->spec);
                     51:        for ( b = hlist = 0, i = 0, vl1 = vl; i < n; vl1 = NEXT(vl1), i++ ) {
                     52:                ptozp(vl1->v->attr,1,&c,&defpoly[i]);
                     53:                ptod(ALG,vl,defpoly[i],&ps[i]);
                     54:                STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
1.3       noro       55:                MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
1.1       noro       56:        }
1.3       noro       57:        ptod(ALG,vl,(P)ONE,&one);
                     58:        MKDAlg(one,ONE,nf->one);
                     59:        nf->ind = b;
                     60:        dp_mbase(hlist,&mblist);
1.1       noro       61:        initd(current_spec);
                     62:        nf->dim = dim = length(mblist);
                     63:        nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
                     64:        for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
1.8       noro       65:                mb[dim-i-1] = (DP)BDY(t);
1.1       noro       66: }
                     67:
1.10      noro       68: void setfield_gb(NODE gb,VL vl,struct order_spec *spec)
                     69: {
                     70:        NumberField nf;
                     71:        VL vl1,vl2;
                     72:        int n,i,dim;
                     73:        Alg *gen;
                     74:        P *defpoly;
                     75:        P p;
                     76:        Q c,iq,two;
                     77:        DP *ps,*mb;
                     78:        DP one;
                     79:        NODE t,b,b1,b2,hlist,mblist;
                     80:        struct order_spec *current_spec;
                     81:
                     82:        nf = (NumberField)MALLOC(sizeof(struct oNumberField));
                     83:        current_numberfield = nf;
                     84:        for ( vl1 = vl, n = 0; vl1; vl1 = NEXT(vl1), n++ );
                     85:        nf->n = n;
                     86:        nf->psn = length(gb);
                     87:        nf->vl = vl;
                     88:        nf->defpoly = defpoly = (P *)MALLOC(nf->psn*sizeof(P));
                     89:        nf->ps = ps = (DP *)MALLOC(nf->psn*sizeof(DP));
                     90:        current_spec = dp_current_spec;
                     91:        nf->spec = spec;
                     92:        initd(nf->spec);
                     93:        for ( b = hlist = 0, i = 0, t = gb; i < nf->psn; t = NEXT(t), i++ ) {
                     94:                ptozp((P)BDY(t),1,&c,&defpoly[i]);
                     95:                ptod(CO,vl,defpoly[i],&ps[i]);
                     96:                STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
                     97:                MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
                     98:        }
                     99:        ptod(ALG,vl,(P)ONE,&one);
                    100:        MKDAlg(one,ONE,nf->one);
                    101:        nf->ind = b;
                    102:        dp_mbase(hlist,&mblist);
                    103:        initd(current_spec);
                    104:        nf->dim = dim = length(mblist);
                    105:        nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
                    106:        for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
                    107:                mb[dim-i-1] = (DP)BDY(t);
                    108: }
                    109:
1.4       noro      110: void qtodalg(Q q,DAlg *r)
                    111: {
                    112:        NumberField nf;
                    113:        Q t;
                    114:        DP nm;
                    115:
                    116:        if ( !(nf=current_numberfield) )
                    117:                error("qtodalg : current_numberfield is not set");
                    118:        if ( !q )
                    119:                *r = 0;
                    120:        else if ( NID(q) == N_DA )
                    121:                *r = (DAlg)q;
                    122:        else if ( NID(q) == N_Q ) {
                    123:                if ( INT(q) ) {
                    124:                        muldc(CO,nf->one->nm,(P)q,&nm);
                    125:                        MKDAlg(nm,ONE,*r);
                    126:                } else {
                    127:                        NTOQ(NM(q),SGN(q),t);
                    128:                        muldc(CO,nf->one->nm,(P)t,&nm);
                    129:                        NTOQ(DN(q),1,t);
                    130:                        MKDAlg(nm,t,*r);
                    131:                }
                    132:        } else
                    133:                error("qtodalg : invalid argument");
1.6       noro      134: }
                    135:
                    136: void obj_algtodalg(Obj obj,Obj *r)
                    137: {
                    138:        DAlg d;
                    139:        DCP dc,dcr0,dcr;
                    140:        P c,p;
                    141:        Obj t;
                    142:        Obj nm,dn;
                    143:        NODE b,s,s0;
                    144:        R rat;
                    145:        VECT v;
                    146:        MAT mat;
                    147:        LIST list;
                    148:        pointer *a;
                    149:        pointer **m;
                    150:        int len,row,col,i,j,l;
                    151:
                    152:        if ( !obj ) {
                    153:                *r = 0;
                    154:                return;
                    155:        }
                    156:        switch ( OID(obj) ) {
                    157:                case O_N:
                    158:                        algtodalg((Alg)obj,&d); *r = (Obj)d;
                    159:                        break;
                    160:                case O_P:
                    161:                        for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {
                    162:                                obj_algtodalg((Obj)COEF(dc),&t);
                    163:                                if ( t ) {
                    164:                                        NEXTDC(dcr0,dcr);
                    165:                                        COEF(dcr) = (P)t;
                    166:                                        DEG(dcr) = DEG(dc);
                    167:                                }
                    168:                        }
                    169:                        if ( dcr0 ) {
                    170:                                MKP(VR((P)obj),dcr0,p);
                    171:                                *r = (Obj)p;
                    172:                        } else
                    173:                                *r = 0;
                    174:                        break;
                    175:                case O_R:
                    176:                        obj_algtodalg((Obj)NM((R)obj),&nm);
                    177:                        obj_algtodalg((Obj)DN((R)obj),&dn);
                    178:                        if ( !dn )
                    179:                                error("obj_algtodalg : division by 0");
                    180:                        if ( !nm )
                    181:                                *r = 0;
                    182:                        else {
                    183:                                MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;
                    184:                        }
                    185:                        break;
                    186:                case O_LIST:
                    187:                        s0 = 0;
                    188:                        for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {
                    189:                                NEXTNODE(s0,s);
                    190:                                obj_algtodalg((Obj)BDY(b),&t);
                    191:                                BDY(s) = (pointer)t;
                    192:                        }
                    193:                        NEXT(s) = 0;
                    194:                        MKLIST(list,s0);
                    195:                        *r = (Obj)list;
                    196:                        break;
                    197:                case O_VECT:
                    198:                        l = ((VECT)obj)->len;
                    199:                        a = BDY((VECT)obj);
                    200:                        MKVECT(v,l);
                    201:                        for ( i = 0; i < l; i++ ) {
                    202:                                obj_algtodalg((Obj)a[i],&t);
                    203:                                BDY(v)[i] = (pointer)t;
                    204:                        }
                    205:                        *r = (Obj)v;
                    206:                        break;
                    207:                case O_MAT:
                    208:                        row = ((MAT)obj)->row; col = ((MAT)obj)->col;
                    209:                        m = BDY((MAT)obj);
                    210:                        MKMAT(mat,row,col);
                    211:                        for ( i = 0; i < row; i++ )
                    212:                                for ( j = 0; j < col; j++ ) {
                    213:                                        obj_algtodalg((Obj)m[i][j],&t);
                    214:                                        BDY(mat)[i][j] = (pointer)t;
                    215:                                }
                    216:                        *r = (Obj)mat;
                    217:                        break;
                    218:                default:
                    219:                        *r = obj;
                    220:                        break;
                    221:        }
                    222: }
                    223:
                    224: void obj_dalgtoalg(Obj obj,Obj *r)
                    225: {
                    226:        Alg d;
                    227:        DCP dc,dcr0,dcr;
                    228:        P c,p;
                    229:        Obj t;
                    230:        Obj nm,dn;
                    231:        NODE b,s,s0;
                    232:        R rat;
                    233:        VECT v;
                    234:        MAT mat;
                    235:        LIST list;
                    236:        pointer *a;
                    237:        pointer **m;
                    238:        int len,row,col,i,j,l;
                    239:
                    240:        if ( !obj ) {
                    241:                *r = 0;
                    242:                return;
                    243:        }
                    244:        switch ( OID(obj) ) {
                    245:                case O_N:
                    246:                        dalgtoalg((DAlg)obj,&d); *r = (Obj)d;
                    247:                        break;
                    248:                case O_P:
                    249:                        for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {
                    250:                                obj_dalgtoalg((Obj)COEF(dc),&t);
                    251:                                if ( t ) {
                    252:                                        NEXTDC(dcr0,dcr);
                    253:                                        COEF(dcr) = (P)t;
                    254:                                        DEG(dcr) = DEG(dc);
                    255:                                }
                    256:                        }
                    257:                        if ( dcr0 ) {
                    258:                                MKP(VR((P)obj),dcr0,p);
                    259:                                *r = (Obj)p;
                    260:                        } else
                    261:                                *r = 0;
                    262:                        break;
                    263:                case O_R:
                    264:                        obj_dalgtoalg((Obj)NM((R)obj),&nm);
                    265:                        obj_dalgtoalg((Obj)DN((R)obj),&dn);
                    266:                        if ( !dn )
                    267:                                error("obj_dalgtoalg : division by 0");
                    268:                        if ( !nm )
                    269:                                *r = 0;
                    270:                        else {
                    271:                                MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;
                    272:                        }
                    273:                        break;
                    274:                case O_LIST:
                    275:                        s0 = 0;
                    276:                        for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {
                    277:                                NEXTNODE(s0,s);
                    278:                                obj_dalgtoalg((Obj)BDY(b),&t);
                    279:                                BDY(s) = (pointer)t;
                    280:                        }
                    281:                        NEXT(s) = 0;
                    282:                        MKLIST(list,s0);
                    283:                        *r = (Obj)list;
                    284:                        break;
                    285:                case O_VECT:
                    286:                        l = ((VECT)obj)->len;
                    287:                        a = BDY((VECT)obj);
                    288:                        MKVECT(v,l);
                    289:                        for ( i = 0; i < l; i++ ) {
                    290:                                obj_dalgtoalg((Obj)a[i],&t);
                    291:                                BDY(v)[i] = (pointer)t;
                    292:                        }
                    293:                        *r = (Obj)v;
                    294:                        break;
                    295:                case O_MAT:
                    296:                        row = ((MAT)obj)->row; col = ((MAT)obj)->col;
                    297:                        m = BDY((MAT)obj);
                    298:                        MKMAT(mat,row,col);
                    299:                        for ( i = 0; i < row; i++ )
                    300:                                for ( j = 0; j < col; j++ ) {
                    301:                                        obj_dalgtoalg((Obj)m[i][j],&t);
                    302:                                        BDY(mat)[i][j] = (pointer)t;
                    303:                                }
                    304:                        *r = (Obj)mat;
                    305:                        break;
                    306:                default:
                    307:                        *r = obj;
                    308:                        break;
                    309:        }
1.4       noro      310: }
                    311:
1.1       noro      312: void algtodalg(Alg a,DAlg *r)
                    313: {
                    314:        P ap,p,p1;
1.3       noro      315:        Q c,c1,d1,dn,nm;
1.1       noro      316:        DP dp;
                    317:        DAlg da;
                    318:        NumberField nf;
                    319:        struct order_spec *current_spec;
1.3       noro      320:        VL vl,tvl,svl;
                    321:        V v;
1.1       noro      322:
                    323:        if ( !(nf=current_numberfield) )
                    324:                error("algtodalg : current_numberfield is not set");
1.3       noro      325:        if ( !a ) {
                    326:                *r = 0;
                    327:                return;
                    328:        }
                    329:        switch (NID((Num)a) ) {
                    330:                case N_Q:
                    331:                        c = (Q)a;
                    332:                        if ( INT(c) ) {
                    333:                                muldc(CO,nf->one->nm,(P)c,&dp);
                    334:                                MKDAlg(dp,ONE,*r);
                    335:                        } else {
                    336:                                NTOQ(NM(c),SGN(c),c1);
                    337:                                NTOQ(DN(c),1,d1);
                    338:                                muldc(CO,nf->one->nm,(P)c1,&dp);
1.9       noro      339:                                MKDAlg(dp,d1,*r);
1.3       noro      340:                        }
                    341:                        break;
                    342:                case N_A:
                    343:                        ap = (P)BDY(a);
                    344:                        ptozp(ap,1,&c,&p);
                    345:                        if ( INT(c) ) {
                    346:                                p = ap;
                    347:                                dn = ONE;
                    348:                        } else {
                    349:                                NTOQ(NM(c),SGN(c),nm);
                    350:                                NTOQ(DN(c),1,dn);
                    351:                                mulpq(p,(P)nm,&p1); p = p1;
                    352:                        }
                    353:                        current_spec = dp_current_spec; initd(nf->spec);
                    354:                        get_vars(p,&vl);
                    355:                        for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
                    356:                                v = tvl->v;
                    357:                                for ( svl = nf->vl; svl; svl = NEXT(svl) )
                    358:                                        if ( v == svl->v )
                    359:                                                break;
                    360:                                if ( !svl )
                    361:                                        error("algtodalg : incompatible numberfield");
                    362:                        }
                    363:                        ptod(ALG,nf->vl,p,&dp);
                    364:                        MKDAlg(dp,dn,da);
                    365:                        simpdalg(da,r);
                    366:                        break;
                    367:                default:
                    368:                        error("algtodalg : invalid argument");
                    369:                        break;
1.1       noro      370:        }
                    371: }
                    372:
1.2       noro      373: void dalgtoalg(DAlg da,Alg *r)
1.1       noro      374: {
1.2       noro      375:        NumberField nf;
                    376:        P p,p1;
                    377:        Q inv;
                    378:
                    379:        if ( !(nf=current_numberfield) )
1.3       noro      380:                error("dalgtoalg : current_numberfield is not set");
1.2       noro      381:        dtop(ALG,nf->vl,da->nm,&p);
                    382:        invq(da->dn,&inv);
                    383:        mulpq(p,(P)inv,&p1);
                    384:        MKAlg(p1,*r);
1.1       noro      385: }
                    386:
                    387: void simpdalg(DAlg da,DAlg *r)
                    388: {
1.2       noro      389:        NumberField nf;
                    390:        DP nm;
1.3       noro      391:        DAlg d;
1.2       noro      392:        Q dn,dn1;
1.5       noro      393:        struct order_spec *current_spec;
1.2       noro      394:
                    395:        if ( !(nf=current_numberfield) )
1.3       noro      396:                error("simpdalg : current_numberfield is not set");
                    397:        if ( !da ) {
                    398:                *r = 0;
                    399:                return;
                    400:        }
1.5       noro      401:        current_spec = dp_current_spec; initd(nf->spec);
1.2       noro      402:        dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);
1.10      noro      403:        if ( !nm ) *r = 0;
                    404:        else {
                    405:                initd(current_spec);
                    406:                mulq(da->dn,dn,&dn1);
                    407:                MKDAlg(nm,dn1,d);
                    408:                rmcontdalg(d,r);
                    409:        }
1.1       noro      410: }
                    411:
                    412: void adddalg(DAlg a,DAlg b,DAlg *c)
                    413: {
1.3       noro      414:        NumberField nf;
                    415:        Q dna,dnb,a1,b1,dn,g;
                    416:        N an,bn,gn;
1.4       noro      417:        DAlg t;
1.3       noro      418:        DP ta,tb,nm;
                    419:        struct order_spec *current_spec;
                    420:
                    421:        if ( !(nf=current_numberfield) )
                    422:                error("adddalg : current_numberfield is not set");
                    423:        if ( !a )
                    424:                *c = b;
                    425:        else if ( !b )
                    426:                *c = a;
                    427:        else {
1.4       noro      428:                qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3       noro      429:                dna = a->dn;
                    430:                dnb = b->dn;
                    431:                gcdn(NM(dna),NM(dnb),&gn);
                    432:                divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
                    433:                NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
                    434:                /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */
                    435:                muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
                    436:                current_spec = dp_current_spec; initd(nf->spec);
                    437:                addd(CO,ta,tb,&nm);
                    438:                initd(current_spec);
                    439:                if ( !nm )
                    440:                        *c = 0;
                    441:                else {
                    442:                        mulq(dna,b1,&dn);
                    443:                        MKDAlg(nm,dn,*c);
                    444:                }
                    445:        }
1.1       noro      446: }
                    447:
                    448: void subdalg(DAlg a,DAlg b,DAlg *c)
                    449: {
1.3       noro      450:        NumberField nf;
                    451:        Q dna,dnb,a1,b1,dn,g;
                    452:        N an,bn,gn;
                    453:        DP ta,tb,nm;
1.4       noro      454:        DAlg t;
1.3       noro      455:        struct order_spec *current_spec;
                    456:
                    457:        if ( !(nf=current_numberfield) )
                    458:                error("subdalg : current_numberfield is not set");
                    459:        if ( !a )
                    460:                *c = b;
                    461:        else if ( !b )
                    462:                *c = a;
                    463:        else {
1.4       noro      464:                qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3       noro      465:                dna = a->dn;
                    466:                dnb = b->dn;
                    467:                gcdn(NM(dna),NM(dnb),&gn);
                    468:                divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
                    469:                NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
                    470:                /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */
                    471:                muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
                    472:                current_spec = dp_current_spec; initd(nf->spec);
                    473:                subd(CO,ta,tb,&nm);
                    474:                initd(current_spec);
                    475:                if ( !nm )
                    476:                        *c = 0;
                    477:                else {
                    478:                        mulq(dna,b1,&dn);
                    479:                        MKDAlg(nm,dn,*c);
                    480:                }
                    481:        }
1.1       noro      482: }
                    483:
1.2       noro      484: void muldalg(DAlg a,DAlg b,DAlg *c)
                    485: {
1.3       noro      486:        NumberField nf;
                    487:        DP nm;
                    488:        Q dn;
                    489:        DAlg t;
                    490:        struct order_spec *current_spec;
                    491:
                    492:        if ( !(nf=current_numberfield) )
                    493:                error("muldalg : current_numberfield is not set");
                    494:        if ( !a || !b )
                    495:                *c = 0;
                    496:        else {
1.4       noro      497:                qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.3       noro      498:                current_spec = dp_current_spec; initd(nf->spec);
                    499:                muld(CO,a->nm,b->nm,&nm);
                    500:                initd(current_spec);
                    501:                mulq(a->dn,b->dn,&dn);
                    502:                MKDAlg(nm,dn,t);
                    503:                simpdalg(t,c);
                    504:        }
1.2       noro      505: }
                    506:
                    507:
1.1       noro      508: void divdalg(DAlg a,DAlg b,DAlg *c)
                    509: {
1.4       noro      510:        DAlg inv,t;
1.8       noro      511:        int ret;
1.3       noro      512:
1.1       noro      513:        if ( !current_numberfield )
1.3       noro      514:                error("divdalg : current_numberfield is not set");
                    515:        if ( !b )
                    516:                error("divdalg : division by 0");
                    517:        if ( !a )
                    518:                c = 0;
                    519:        else {
1.4       noro      520:                qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
1.8       noro      521:                ret = invdalg(b,&inv);
                    522:                if ( !ret ) {
                    523:                        error("divdalg : the denominator is not invertible");
                    524:                }
1.3       noro      525:                muldalg(a,inv,c);
                    526:        }
                    527: }
                    528:
                    529: void rmcontdalg(DAlg a, DAlg *r)
                    530: {
                    531:        DP u,u1;
                    532:        Q cont,c,d;
                    533:        N gn,cn,dn;
                    534:
                    535:        if ( !a )
                    536:                *r = a;
                    537:        else {
                    538:                dp_ptozp(a->nm,&u);
                    539:                divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);
                    540:                gcdn(NM(cont),NM(a->dn),&gn);
                    541:                divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
                    542:                divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
                    543:                muldc(CO,u,(P)c,&u1);
                    544:                MKDAlg(u1,d,*r);
                    545:        }
1.1       noro      546: }
                    547:
1.8       noro      548: int invdalg(DAlg a,DAlg *c)
1.1       noro      549: {
1.3       noro      550:        NumberField nf;
1.8       noro      551:        int dim,n,i,j,k,l;
1.3       noro      552:        DP *mb;
                    553:        DP m,d,u;
                    554:        N ln,gn,qn;
                    555:        DAlg *simp;
                    556:        DAlg t,a0,r;
1.7       noro      557:        Q dn,dnsol,mul,nmc,dn1;
1.3       noro      558:        MAT mobj,sol;
                    559:        Q **mat,**solmat;
                    560:        MP mp0,mp;
                    561:        int *rinfo,*cinfo;
1.8       noro      562:        int rank,nparam;
                    563:        NODE nd0,nd,ndt;
1.3       noro      564:        struct order_spec *current_spec;
1.5       noro      565:        struct oEGT eg0,eg1;
                    566:        extern struct oEGT eg_le;
1.3       noro      567:
                    568:        if ( !(nf=current_numberfield) )
                    569:                error("invdalg : current_numberfield is not set");
                    570:        if ( !a )
                    571:                error("invdalg : division by 0");
1.4       noro      572:        else if ( NID(a) == N_Q ) {
                    573:                invq((Q)a,&dn); *c = (DAlg)dn;
                    574:                return;
                    575:        }
1.3       noro      576:        dim = nf->dim;
                    577:        mb = nf->mb;
                    578:        n = nf->n;
                    579:        ln = ONEN;
1.7       noro      580:        dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
                    581:        MKDAlg(u,ONE,a0);
1.3       noro      582:        simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
                    583:        current_spec = dp_current_spec; initd(nf->spec);
1.8       noro      584:        for ( i = 0; i < dim; i++ ) {
1.3       noro      585:                m = mb[i];
1.8       noro      586:                for ( j = i-1; j >= 0; j-- )
1.3       noro      587:                        if ( dp_redble(m,mb[j]) )
                    588:                                break;
1.8       noro      589:                if ( j >= 0 ) {
1.3       noro      590:                        dp_subd(m,mb[j],&d);
                    591:                        muld(CO,d,simp[j]->nm,&u);
                    592:                        MKDAlg(u,simp[j]->dn,t);
                    593:                        simpdalg(t,&simp[i]);
                    594:                } else {
                    595:                        MKDAlg(m,ONE,t);
                    596:                        muldalg(t,a0,&simp[i]);
                    597:                }
                    598:                gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
                    599:                muln(NM(simp[i]->dn),qn,&ln);
                    600:        }
                    601:        initd(current_spec);
                    602:        NTOQ(ln,1,dn);
                    603:        MKMAT(mobj,dim,dim+1);
                    604:        mat = (Q **)BDY(mobj);
1.8       noro      605:        mulq(dn,a->dn,&mat[0][dim]);
1.3       noro      606:        for ( j = 0; j < dim; j++ ) {
                    607:                divq(dn,simp[j]->dn,&mul);
1.8       noro      608:                for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
1.3       noro      609:                        if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
                    610:                                mulq(mul,(Q)mp->c,&mat[i][j]);
                    611:                                mp = NEXT(mp);
                    612:                        }
                    613:        }
1.5       noro      614:        get_eg(&eg0);
1.8       noro      615:        rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
1.5       noro      616:        get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
1.8       noro      617:        if ( cinfo[0] == dim ) {
                    618:                /* the input is invertible */
                    619:                solmat = (Q **)BDY(sol);
                    620:                for ( i = dim-1, mp0 = 0; i >= 0; i-- )
                    621:                        if ( solmat[i][0] ) {
                    622:                                NEXTMP(mp0,mp);
                    623:                                mp->c = (P)solmat[i][0];
                    624:                                mp->dl = BDY(mb[i])->dl;
                    625:                        }
                    626:                NEXT(mp) = 0; MKDP(n,mp0,u);
                    627:                mulq(dnsol,nmc,&dn1);
                    628:                MKDAlg(u,dn1,r);
                    629:                rmcontdalg(r,c);
                    630:                return 1;
                    631:        } else
                    632:                return 0;
                    633: }
                    634:
                    635: NODE inv_or_split_dalg(DAlg a,DAlg *c)
                    636: {
                    637:        NumberField nf;
                    638:        int dim,n,i,j,k,l;
                    639:        DP *mb;
                    640:        DP m,d,u;
                    641:        N ln,gn,qn;
                    642:        DAlg *simp;
                    643:        DAlg t,a0,r;
                    644:        Q dn,dnsol,mul,nmc,dn1;
                    645:        MAT mobj,sol;
                    646:        Q **mat,**solmat;
                    647:        MP mp0,mp;
                    648:        int *rinfo,*cinfo;
                    649:        int rank,nparam;
                    650:        NODE nd0,nd,ndt;
                    651:        struct order_spec *current_spec;
                    652:        struct oEGT eg0,eg1;
                    653:        extern struct oEGT eg_le;
1.12      noro      654:        extern int DP_Print;
1.8       noro      655:
                    656:        if ( !(nf=current_numberfield) )
                    657:                error("invdalg : current_numberfield is not set");
                    658:        if ( !a )
                    659:                error("invdalg : division by 0");
                    660:        else if ( NID(a) == N_Q ) {
                    661:                invq((Q)a,&dn); *c = (DAlg)dn;
                    662:                return;
                    663:        }
                    664:        dim = nf->dim;
                    665:        mb = nf->mb;
                    666:        n = nf->n;
                    667:        ln = ONEN;
                    668:        dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
                    669:        MKDAlg(u,ONE,a0);
                    670:        simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
                    671:        current_spec = dp_current_spec; initd(nf->spec);
                    672:        for ( i = 0; i < dim; i++ ) {
1.12      noro      673:                if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
1.8       noro      674:                m = mb[i];
                    675:                for ( j = i-1; j >= 0; j-- )
                    676:                        if ( dp_redble(m,mb[j]) )
                    677:                                break;
                    678:                if ( j >= 0 ) {
                    679:                        dp_subd(m,mb[j],&d);
1.11      noro      680:                        if ( simp[j] ) {
                    681:                                muld(CO,d,simp[j]->nm,&u);
                    682:                                MKDAlg(u,simp[j]->dn,t);
                    683:                                simpdalg(t,&simp[i]);
                    684:                        } else
                    685:                                simp[i] = 0;
1.8       noro      686:                } else {
                    687:                        MKDAlg(m,ONE,t);
                    688:                        muldalg(t,a0,&simp[i]);
                    689:                }
1.11      noro      690:                if ( simp[i] ) {
                    691:                        gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
                    692:                        muln(NM(simp[i]->dn),qn,&ln);
                    693:                }
1.8       noro      694:        }
                    695:        initd(current_spec);
                    696:        NTOQ(ln,1,dn);
                    697:        MKMAT(mobj,dim,dim+1);
                    698:        mat = (Q **)BDY(mobj);
                    699:        mulq(dn,a->dn,&mat[0][dim]);
                    700:        for ( j = 0; j < dim; j++ ) {
1.11      noro      701:                if ( simp[j] ) {
                    702:                        divq(dn,simp[j]->dn,&mul);
                    703:                        for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
                    704:                                if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
                    705:                                        mulq(mul,(Q)mp->c,&mat[i][j]);
                    706:                                        mp = NEXT(mp);
                    707:                                }
                    708:                }
1.8       noro      709:        }
                    710:        get_eg(&eg0);
1.13    ! noro      711:        rank = generic_gauss_elim_hensel_dalg(mobj,&sol,&dnsol,&rinfo,&cinfo);
1.8       noro      712:        get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
                    713:        if ( cinfo[0] == dim ) {
                    714:                /* the input is invertible */
                    715:                solmat = (Q **)BDY(sol);
                    716:                for ( i = dim-1, mp0 = 0; i >= 0; i-- )
                    717:                        if ( solmat[i][0] ) {
                    718:                                NEXTMP(mp0,mp);
                    719:                                mp->c = (P)solmat[i][0];
                    720:                                mp->dl = BDY(mb[i])->dl;
                    721:                        }
                    722:                NEXT(mp) = 0; MKDP(n,mp0,u);
                    723:                mulq(dnsol,nmc,&dn1);
                    724:                MKDAlg(u,dn1,r);
                    725:                rmcontdalg(r,c);
                    726:                return 0;
                    727:        } else {
                    728:                /* the input is not invertible */
1.13    ! noro      729:                nparam = sol->col;
1.8       noro      730:                solmat = (Q **)BDY(sol);
                    731:                nd0 = 0;
                    732:                for ( k = 0; k < nparam; k++ ) {
1.13    ! noro      733:                        /* construct a new basis element */
1.8       noro      734:                        m = mb[cinfo[k]];
                    735:                        mp0 = 0;
1.3       noro      736:                        NEXTMP(mp0,mp);
1.8       noro      737:                        chsgnq(dnsol,&dn1); mp->c = (P)dn1;
                    738:                        mp->dl = BDY(m)->dl;
                    739:                        /* skip the last parameter */
                    740:                        for ( l = rank-2; l >= 0; l-- ) {
                    741:                                if ( solmat[l][k] ) {
                    742:                                        NEXTMP(mp0,mp);
                    743:                                        mp->c = (P)solmat[l][k];
                    744:                                        mp->dl = BDY(mb[rinfo[l]])->dl;
                    745:                                }
                    746:                        }
                    747:                        NEXT(mp) = 0; MKDP(n,mp0,u);
                    748:                        NEXTNODE(nd0,nd);
                    749:                        BDY(nd) = (pointer)u;
                    750:                        NEXT(nd) = 0;
1.3       noro      751:                }
1.8       noro      752:                NEXT(nd) = 0;
                    753:                return nd0;
                    754:        }
1.1       noro      755: }
                    756:
                    757: void chsgndalg(DAlg a,DAlg *c)
                    758: {
1.3       noro      759:        DP nm;
1.4       noro      760:        Q t;
1.3       noro      761:
                    762:        if ( !a ) *c = 0;
1.4       noro      763:        else if ( NID(a) == N_Q ) {
                    764:                chsgnq((Q)a,&t); *c = (DAlg)t;
                    765:        } else {
1.3       noro      766:                chsgnd(a->nm,&nm);
                    767:                MKDAlg(nm,a->dn,*c);
                    768:        }
1.1       noro      769: }
                    770:
1.3       noro      771: void pwrdalg(DAlg a,Q e,DAlg *c)
1.1       noro      772: {
1.3       noro      773:        NumberField nf;
                    774:        DAlg t,z,y;
1.4       noro      775:        Q q;
1.3       noro      776:        N en,qn;
                    777:        int r;
1.8       noro      778:        int ret;
1.3       noro      779:
                    780:        if ( !(nf=current_numberfield) )
                    781:                error("pwrdalg : current_numberfield is not set");
1.4       noro      782:        if ( !a )
                    783:                *c = !e ? (DAlg)ONE : 0;
                    784:        else if ( NID(a) == N_Q ) {
                    785:                pwrq((Q)a,e,&q); *c = (DAlg)q;
                    786:        } else if ( !e )
1.3       noro      787:                *c = nf->one;
                    788:        else if ( UNIQ(e) )
                    789:                *c = a;
                    790:        else {
                    791:                if ( SGN(e) < 0 ) {
1.8       noro      792:                        ret = invdalg(a,&t);
                    793:                        if ( !ret )
                    794:                                error("pwrdalg : the denominator is not invertible");
                    795:                        a = t;
1.3       noro      796:                }
                    797:                en = NM(e);
                    798:                y = nf->one;
                    799:                z = a;
                    800:                while ( 1 ) {
                    801:                        r = divin(en,2,&qn); en = qn;
                    802:                        if ( r ) {
                    803:                                muldalg(z,y,&t); y = t;
                    804:                                if ( !en ) {
                    805:                                        *c = y;
                    806:                                        return;
                    807:                                }
                    808:                        }
                    809:                        muldalg(z,z,&t); z = t;
                    810:                }
                    811:        }
1.1       noro      812: }
                    813:
1.3       noro      814: int cmpdalg(DAlg a,DAlg b)
1.1       noro      815: {
1.3       noro      816:        DAlg c;
                    817:
                    818:        subdalg(a,b,&c);
                    819:        if ( !c ) return 0;
                    820:        else
                    821:                return SGN((Q)BDY(c->nm)->c);
1.1       noro      822: }
1.9       noro      823:
                    824: /* convert da to a univariate poly; return the position of variable */
                    825:
                    826: int dalgtoup(DAlg da,P *up,Q *dn)
                    827: {
                    828:        int nv,i,hi,current_d;
                    829:        DCP dc0,dc;
                    830:        MP h,mp0,mp,t;
                    831:        DL hd,d;
                    832:        DP c;
                    833:        DAlg cc;
                    834:        P v;
                    835:
                    836:        nv = da->nm->nv;
                    837:        h = BDY(da->nm);
                    838:        *dn = da->dn;
                    839:        hd = h->dl;
                    840:        for ( i = 0; i < nv; i++ )
                    841:                if ( hd->d[i] ) break;
                    842:        hi = i;
                    843:        current_d = hd->d[i];
                    844:        dc0 = 0;
                    845:        mp0 = 0;
                    846:        for ( t = h; t; t = NEXT(t) ) {
                    847:                NEWDL(d,nv);
                    848:                for ( i = 0; i <= hi; i++ ) d->d[i] = 0;
                    849:                for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i];
                    850:                d->td = t->dl->td - t->dl->d[hi];
                    851:                if ( t->dl->d[hi] != current_d ) {
                    852:                        NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
                    853:                        NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
                    854:                        current_d = t->dl->d[hi];
                    855:                        mp0 = 0;
                    856:                }
                    857:                NEXTMP(mp0,mp);
                    858:                mp->c = t->c; mp->dl = d;
                    859:        }
                    860:        NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
                    861:        NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
                    862:        NEXT(dc) = 0;
                    863:        makevar("x",&v);
                    864:        MKP(VR(v),dc0,*up);
                    865:        return hi;
                    866: }
                    867:

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>