[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.9

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

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