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

Annotation of OpenXM_contrib2/asir2000/engine/PU.c, Revision 1.1.1.1

1.1       noro        1: /* $OpenXM: OpenXM/src/asir99/engine/PU.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
                      2: #include "ca.h"
                      3:
                      4: void reorderp(nvl,ovl,p,pr)
                      5: VL nvl,ovl;
                      6: P p;
                      7: P *pr;
                      8: {
                      9:        DCP dc;
                     10:        P x,m,s,t,c;
                     11:
                     12:        if ( !p )
                     13:                *pr = 0;
                     14:        else if ( NUM(p) )
                     15:                *pr = p;
                     16:        else {
                     17:                MKV(VR(p),x);
                     18:                for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                     19:                        reorderp(nvl,ovl,COEF(dc),&c);
                     20:                        if ( DEG(dc) ) {
                     21:                                pwrp(nvl,x,DEG(dc),&t); mulp(nvl,c,t,&m);
                     22:                                addp(nvl,m,s,&t);
                     23:                        } else
                     24:                                addp(nvl,s,c,&t);
                     25:                        s = t;
                     26:                }
                     27:                *pr = s;
                     28:        }
                     29: }
                     30:
                     31: void substp(vl,p,v0,p0,pr)
                     32: VL vl;
                     33: V v0;
                     34: P p,p0;
                     35: P *pr;
                     36: {
                     37:        P x,t,m,c,s,a;
                     38:        DCP dc;
                     39:        Q d;
                     40:
                     41:        if ( !p )
                     42:                *pr = 0;
                     43:        else if ( NUM(p) )
                     44:                *pr = p;
                     45:        else if ( VR(p) != v0 ) {
                     46:                MKV(VR(p),x);
                     47:                for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                     48:                        substp(vl,COEF(dc),v0,p0,&t);
                     49:                        if ( DEG(dc) ) {
                     50:                                pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
                     51:                                addp(vl,m,c,&a);
                     52:                                c = a;
                     53:                        } else {
                     54:                                addp(vl,t,c,&a);
                     55:                                c = a;
                     56:                        }
                     57:                }
                     58:                *pr = c;
                     59:        } else {
                     60:                dc = DC(p);
                     61:                c = COEF(dc);
                     62:                for ( d = DEG(dc), dc = NEXT(dc);
                     63:                        dc; d = DEG(dc), dc = NEXT(dc) ) {
                     64:                                subq(d,DEG(dc),(Q *)&t); pwrp(vl,p0,(Q)t,&s);
                     65:                                mulp(vl,s,c,&m);
                     66:                                addp(vl,m,COEF(dc),&c);
                     67:                }
                     68:                if ( d ) {
                     69:                        pwrp(vl,p0,d,&t); mulp(vl,t,c,&m);
                     70:                        c = m;
                     71:                }
                     72:                *pr = c;
                     73:        }
                     74: }
                     75:
                     76: void detp(vl,rmat,n,dp)
                     77: VL vl;
                     78: P **rmat;
                     79: int n;
                     80: P *dp;
                     81: {
                     82:        int i,j,k,sgn;
                     83:        P mjj,mij,t,s,u,d;
                     84:        P **mat;
                     85:        P *mi,*mj;
                     86:
                     87:        mat = (P **)almat_pointer(n,n);
                     88:        for ( i = 0; i < n; i++ )
                     89:                for ( j = 0; j < n; j++ )
                     90:                        mat[i][j] = rmat[i][j];
                     91:        for ( j = 0, d = (P)ONE, sgn = 1; j < n; j++ ) {
                     92:                for ( i = j; (i < n) && !mat[i][j]; i++ );
                     93:                if ( i == n ) {
                     94:                        *dp = 0; return;
                     95:                }
                     96:                for ( k = i; k < n; k++ )
                     97:                        if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
                     98:                                i = k;
                     99:                if ( j != i ) {
                    100:                        mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
                    101:                }
                    102:                for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
                    103:                        for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
                    104:                                mulp(vl,mi[k],mjj,&t); mulp(vl,mj[k],mij,&s);
                    105:                                subp(vl,t,s,&u); divsp(vl,u,d,&mi[k]);
                    106:                        }
                    107:                d = mjj;
                    108:        }
                    109:        if ( sgn < 0 )
                    110:                chsgnp(d,dp);
                    111:        else
                    112:                *dp = d;
                    113: }
                    114:
                    115: void reordvar(vl,v,nvlp)
                    116: VL vl;
                    117: V v;
                    118: VL *nvlp;
                    119: {
                    120:        VL nvl,nvl0;
                    121:
                    122:        for ( NEWVL(nvl0), nvl0->v = v, nvl = nvl0;
                    123:                vl; vl = NEXT(vl) )
                    124:                if ( vl->v == v )
                    125:                        continue;
                    126:                else {
                    127:                        NEWVL(NEXT(nvl));
                    128:                        nvl = NEXT(nvl);
                    129:                        nvl->v = vl->v;
                    130:                }
                    131:        NEXT(nvl) = 0;
                    132:        *nvlp = nvl0;
                    133: }
                    134:
                    135: void gcdprsp(vl,p1,p2,pr)
                    136: VL vl;
                    137: P p1,p2,*pr;
                    138: {
                    139:        P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
                    140:        V v1,v2;
                    141:
                    142:        if ( !p1 )
                    143:                ptozp0(p2,pr);
                    144:        else if ( !p2 )
                    145:                ptozp0(p1,pr);
                    146:        else if ( NUM(p1) || NUM(p2) )
                    147:                *pr = (P)ONE;
                    148:        else {
                    149:                ptozp0(p1,&g1); ptozp0(p2,&g2);
                    150:                if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {
                    151:                        gcdcp(vl,g1,&gc1); divsp(vl,g1,gc1,&gp1);
                    152:                        gcdcp(vl,g2,&gc2); divsp(vl,g2,gc2,&gp2);
                    153:                        gcdprsp(vl,gc1,gc2,&gcr);
                    154:                        sprs(vl,v1,gp1,gp2,&g);
                    155:
                    156:                        if ( VR(g) == v1 ) {
                    157:                                ptozp0(g,&gp);
                    158:                                gcdcp(vl,gp,&gc); divsp(vl,gp,gc,&gp1);
                    159:                                mulp(vl,gp1,gcr,pr);
                    160:
                    161:                        } else
                    162:                                *pr = gcr;
                    163:                } else {
                    164:                        while ( v1 != vl->v && v2 != vl->v )
                    165:                                vl = NEXT(vl);
                    166:                        if ( v1 == vl->v ) {
                    167:                                gcdcp(vl,g1,&gc1); gcdprsp(vl,gc1,g2,pr);
                    168:                        } else {
                    169:                                gcdcp(vl,g2,&gc2); gcdprsp(vl,gc2,g1,pr);
                    170:                        }
                    171:                }
                    172:        }
                    173: }
                    174:
                    175: void gcdcp(vl,p,pr)
                    176: VL vl;
                    177: P p,*pr;
                    178: {
                    179:        P g,g1;
                    180:        DCP dc;
                    181:
                    182:        if ( NUM(p) )
                    183:                *pr = (P)ONE;
                    184:        else {
                    185:                dc = DC(p);
                    186:                g = COEF(dc);
                    187:                for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                    188:                        gcdprsp(vl,g,COEF(dc),&g1);
                    189:                        g = g1;
                    190:                }
                    191:                *pr = g;
                    192:        }
                    193: }
                    194:
                    195: void sprs(vl,v,p1,p2,pr)
                    196: VL vl;
                    197: V v;
                    198: P p1,p2,*pr;
                    199: {
                    200:        P q1,q2,m,m1,m2,x,h,r,g1,g2;
                    201:        int d;
                    202:        Q dq;
                    203:        VL nvl;
                    204:
                    205:        reordvar(vl,v,&nvl);
                    206:        reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    207:
                    208:        if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
                    209:                *pr = 0;
                    210:                return;
                    211:        }
                    212:
                    213:        if ( deg(v,q1) >= deg(v,q2) ) {
                    214:                g1 = q1; g2 = q2;
                    215:        } else {
                    216:                g2 = q1; g1 = q2;
                    217:        }
                    218:
                    219:        for ( h = (P)ONE, x = (P)ONE; ; ) {
                    220:                if ( !deg(v,g2) )
                    221:                        break;
                    222:
                    223:                premp(nvl,g1,g2,&r);
                    224:                if ( !r )
                    225:                        break;
                    226:
                    227:                d = deg(v,g1) - deg(v,g2); STOQ(d,dq);
                    228:                pwrp(nvl,h,dq,&m); mulp(nvl,m,x,&m1); g1 = g2;
                    229:                divsp(nvl,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */
                    230:                pwrp(nvl,x,dq,&m1); mulp(nvl,m1,h,&m2);
                    231:                divsp(nvl,m2,m,&h);
                    232:        }
                    233:        *pr = g2;
                    234: }
                    235:
                    236: void resultp(vl,v,p1,p2,pr)
                    237: VL vl;
                    238: V v;
                    239: P p1,p2,*pr;
                    240: {
                    241:        P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;
                    242:        int d,d1,d2,j,k;
                    243:        VL nvl;
                    244:        Q dq;
                    245:
                    246:        if ( !p1 || !p2 ) {
                    247:                *pr = 0; return;
                    248:        }
                    249:        reordvar(vl,v,&nvl);
                    250:        reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    251:
                    252:        if ( VR(q1) != v )
                    253:                if ( VR(q2) != v ) {
                    254:                        *pr = 0;
                    255:                        return;
                    256:                } else {
                    257:                        d = deg(v,q2); STOQ(d,dq);
                    258:                        pwrp(vl,q1,dq,pr);
                    259:                        return;
                    260:                }
                    261:        else if ( VR(q2) != v ) {
                    262:                d = deg(v,q1); STOQ(d,dq);
                    263:                pwrp(vl,q2,dq,pr);
                    264:                return;
                    265:        }
                    266:
                    267:        if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
                    268:                *pr = 0;
                    269:                return;
                    270:        }
                    271:
                    272:        d1 = deg(v,q1); d2 = deg(v,q2);
                    273:        if ( d1 > d2 ) {
                    274:                g1 = q1; g2 = q2; adj = (P)ONE;
                    275:        } else if ( d1 < d2 ) {
                    276:                g2 = q1; g1 = q2;
                    277:                if ( (d1 % 2) && (d2 % 2) ) {
                    278:                        STOQ(-1,dq); adj = (P)dq;
                    279:                } else
                    280:                        adj = (P)ONE;
                    281:        } else {
                    282:                premp(nvl,q1,q2,&t);
                    283:                d = deg(v,t); STOQ(d,dq); pwrp(nvl,LC(q2),dq,&adj);
                    284:                g1 = q2; g2 = t;
                    285:                if ( d1 % 2 ) {
                    286:                        chsgnp(adj,&t); adj = t;
                    287:                }
                    288:        }
                    289:        d1 = deg(v,g1); j = d1 - 1;
                    290:
                    291:        for ( lc = (P)ONE; ; ) {
                    292:                if ( ( k = deg(v,g2) ) < 0 ) {
                    293:                        *pr = 0;
                    294:                        return;
                    295:                }
                    296:
                    297:                if ( k == j )
                    298:                        if ( !k ) {
                    299:                                divsp(nvl,g2,adj,pr);
                    300:                                return;
                    301:                        } else {
                    302:                                premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
                    303:                                divsp(nvl,r,m,&q);
                    304:                                g1 = g2; g2 = q;
                    305:                                lc = LC(g1); /* g1 is not const */
                    306:                                j = k - 1;
                    307:                        }
                    308:                else {
                    309:                        d = j - k; STOQ(d,dq);
                    310:                        pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
                    311:                        mulp(nvl,g2,m,&m1);
                    312:                        pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
                    313:                        if ( k == 0 ) {
                    314:                                divsp(nvl,t,adj,pr);
                    315:                                return;
                    316:                        } else {
                    317:                                premp(nvl,g1,g2,&r);
                    318:                                mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
                    319:                                divsp(nvl,r,m2,&q);
                    320:                                g1 = t; g2 = q;
                    321:                                if ( d % 2 ) {
                    322:                                        chsgnp(g2,&t); g2 = t;
                    323:                                }
                    324:                                lc = LC(g1); /* g1 is not const */
                    325:                                j = k - 1;
                    326:                        }
                    327:                }
                    328:        }
                    329: }
                    330:
                    331: void srch2(vl,v,p1,p2,pr)
                    332: VL vl;
                    333: V v;
                    334: P p1,p2,*pr;
                    335: {
                    336:        P q1,q2,m,m1,m2,lc,q,r,t,s,g1,g2,adj;
                    337:        int d,d1,d2,j,k;
                    338:        VL nvl;
                    339:        Q dq;
                    340:
                    341:        reordvar(vl,v,&nvl);
                    342:        reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    343:
                    344:        if ( VR(q1) != v )
                    345:                if ( VR(q2) != v ) {
                    346:                        *pr = 0;
                    347:                        return;
                    348:                } else {
                    349:                        d = deg(v,q2); STOQ(d,dq);
                    350:                        pwrp(vl,q1,dq,pr);
                    351:                        return;
                    352:                }
                    353:        else if ( VR(q2) != v ) {
                    354:                d = deg(v,q1); STOQ(d,dq);
                    355:                pwrp(vl,q2,dq,pr);
                    356:                return;
                    357:        }
                    358:
                    359:        if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
                    360:                *pr = 0;
                    361:                return;
                    362:        }
                    363:
                    364:        if ( deg(v,q1) >= deg(v,q2) ) {
                    365:                g1 = q1; g2 = q2;
                    366:        } else {
                    367:                g2 = q1; g1 = q2;
                    368:        }
                    369:
                    370:
                    371:        if ( ( d1 = deg(v,g1) ) > ( d2 = deg(v,g2) ) ) {
                    372:                j = d1 - 1;
                    373:                adj = (P)ONE;
                    374:        } else {
                    375:                premp(nvl,g1,g2,&t);
                    376:                d = deg(v,t); STOQ(d,dq);
                    377:                pwrp(nvl,LC(g2),dq,&adj);
                    378:                g1 = g2; g2 = t;
                    379:                j = deg(v,g1) - 1;
                    380:        }
                    381:
                    382:        for ( lc = (P)ONE; ; ) {
                    383:                if ( ( k = deg(v,g2) ) < 0 ) {
                    384:                        *pr = 0;
                    385:                        return;
                    386:                }
                    387:
                    388:                ptozp(g1,1,(Q *)&t,&s); g1 = s; ptozp(g2,1,(Q *)&t,&s); g2 = s;
                    389:                if ( k == j )
                    390:                        if ( !k ) {
                    391:                                divsp(nvl,g2,adj,pr);
                    392:                                return;
                    393:                        } else {
                    394:                                premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
                    395:                                divsp(nvl,r,m,&q);
                    396:                                g1 = g2; g2 = q;
                    397:                                lc = LC(g1); /* g1 is not const */
                    398:                                j = k - 1;
                    399:                        }
                    400:                else {
                    401:                        d = j - k; STOQ(d,dq);
                    402:                        pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
                    403:                        mulp(nvl,g2,m,&m1);
                    404:                        pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
                    405:                        if ( k == 0 ) {
                    406:                                divsp(nvl,t,adj,pr);
                    407:                                return;
                    408:                        } else {
                    409:                                premp(nvl,g1,g2,&r);
                    410:                                mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
                    411:                                divsp(nvl,r,m2,&q);
                    412:                                g1 = t; g2 = q;
                    413:                                if ( d % 2 ) {
                    414:                                        chsgnp(g2,&t); g2 = t;
                    415:                                }
                    416:                                lc = LC(g1); /* g1 is not const */
                    417:                                j = k - 1;
                    418:                        }
                    419:                }
                    420:        }
                    421: }
                    422:
                    423: void srcr(vl,v,p1,p2,pr)
                    424: VL vl;
                    425: V v;
                    426: P p1,p2,*pr;
                    427: {
                    428:        P q1,q2,c,c1;
                    429:        P tg,tg1,tg2,resg;
                    430:        int index,mod;
                    431:        Q a,b,f,q,t,s,u,n,m;
                    432:        VL nvl;
                    433:
                    434:        reordvar(vl,v,&nvl);
                    435:        reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    436:
                    437:        if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
                    438:                *pr = 0;
                    439:                return;
                    440:        }
                    441:        if ( VR(q1) != v ) {
                    442:                pwrp(vl,q1,DEG(DC(q2)),pr);
                    443:                return;
                    444:        }
                    445:        if ( VR(q2) != v ) {
                    446:                pwrp(vl,q2,DEG(DC(q1)),pr);
                    447:                return;
                    448:        }
                    449:        norm1c(q1,&a); norm1c(q2,&b);
                    450:        n = DEG(DC(q1)); m = DEG(DC(q2));
                    451:        pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
                    452:        factorial(QTOS(n)+QTOS(m),&t);
                    453:        mulq(u,t,&s); addq(s,s,&f);
                    454:        for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
                    455:                mod = lprime[index++];
                    456:                if ( !mod )
                    457:                        error("sqfrum : lprime[] exhausted.");
                    458:                ptomp(mod,LC(q1),&tg);
                    459:                if ( !tg )
                    460:                        continue;
                    461:                ptomp(mod,LC(q2),&tg);
                    462:                if ( !tg )
                    463:                        continue;
                    464:                ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
                    465:                srchmp(nvl,mod,v,tg1,tg2,&resg);
                    466:                chnremp(nvl,mod,c,q,resg,&c1); c = c1;
                    467:                STOQ(mod,t); mulq(q,t,&s); q = s;
                    468:        }
                    469:        *pr = c;
                    470: }
                    471:
                    472: void res_ch_det(vl,v,p1,p2,pr)
                    473: VL vl;
                    474: V v;
                    475: P p1,p2,*pr;
                    476: {
                    477:        P q1,q2,c,c1;
                    478:        P tg,tg1,tg2,resg;
                    479:        int index,mod;
                    480:        Q a,b,f,q,t,s,u,n,m;
                    481:        VL nvl;
                    482:
                    483:        reordvar(vl,v,&nvl);
                    484:        reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    485:
                    486:        if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
                    487:                *pr = 0;
                    488:                return;
                    489:        }
                    490:        if ( VR(q1) != v ) {
                    491:                pwrp(vl,q1,DEG(DC(q2)),pr);
                    492:                return;
                    493:        }
                    494:        if ( VR(q2) != v ) {
                    495:                pwrp(vl,q2,DEG(DC(q1)),pr);
                    496:                return;
                    497:        }
                    498:        norm1c(q1,&a); norm1c(q2,&b);
                    499:        n = DEG(DC(q1)); m = DEG(DC(q2));
                    500:        pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
                    501:        factorial(QTOS(n)+QTOS(m),&t);
                    502:        mulq(u,t,&s); addq(s,s,&f);
                    503:        for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
                    504:                mod = lprime[index++];
                    505:                if ( !mod )
                    506:                        error("sqfrum : lprime[] exhausted.");
                    507:                ptomp(mod,LC(q1),&tg);
                    508:                if ( !tg )
                    509:                        continue;
                    510:                ptomp(mod,LC(q2),&tg);
                    511:                if ( !tg )
                    512:                        continue;
                    513:                ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
                    514:                res_detmp(nvl,mod,v,tg1,tg2,&resg);
                    515:                chnremp(nvl,mod,c,q,resg,&c1); c = c1;
                    516:                STOQ(mod,t); mulq(q,t,&s); q = s;
                    517:        }
                    518:        *pr = c;
                    519: }
                    520:
                    521: void res_detmp(vl,mod,v,p1,p2,dp)
                    522: VL vl;
                    523: int mod;
                    524: V v;
                    525: P p1,p2;
                    526: P *dp;
                    527: {
                    528:        int n1,n2,n,sgn;
                    529:        int i,j,k;
                    530:        P mjj,mij,t,s,u,d;
                    531:        P **mat;
                    532:        P *mi,*mj;
                    533:        DCP dc;
                    534:
                    535:        n1 = UDEG(p1); n2 = UDEG(p2); n = n1+n2;
                    536:        mat = (P **)almat_pointer(n,n);
                    537:        for ( dc = DC(p1); dc; dc = NEXT(dc) )
                    538:                mat[0][n1-QTOS(DEG(dc))] = COEF(dc);
                    539:        for ( i = 1; i < n2; i++ )
                    540:                for ( j = 0; j <= n1; j++ )
                    541:                        mat[i][i+j] = mat[0][j];
                    542:        for ( dc = DC(p2); dc; dc = NEXT(dc) )
                    543:                mat[n2][n2-QTOS(DEG(dc))] = COEF(dc);
                    544:        for ( i = 1; i < n1; i++ )
                    545:                for ( j = 0; j <= n2; j++ )
                    546:                        mat[i+n2][i+j] = mat[n2][j];
                    547:        for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
                    548:                for ( i = j; (i < n) && !mat[i][j]; i++ );
                    549:                if ( i == n ) {
                    550:                        *dp = 0; return;
                    551:                }
                    552:                for ( k = i; k < n; k++ )
                    553:                        if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
                    554:                                i = k;
                    555:                if ( j != i ) {
                    556:                        mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
                    557:                }
                    558:                for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
                    559:                        for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
                    560:                                mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
                    561:                                submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
                    562:                        }
                    563:                d = mjj;
                    564:        }
                    565:        if ( sgn > 0 )
                    566:                *dp = d;
                    567:        else
                    568:                chsgnmp(mod,d,dp);
                    569: }
                    570:
                    571: #if 0
                    572: showmat(vl,mat,n)
                    573: VL vl;
                    574: P **mat;
                    575: int n;
                    576: {
                    577:        int i,j;
                    578:        P t;
                    579:
                    580:        for ( i = 0; i < n; i++ ) {
                    581:                for ( j = 0; j < n; j++ ) {
                    582:                        mptop(mat[i][j],&t); printp(vl,t); fprintf(out," ");
                    583:                }
                    584:                fprintf(out,"\n");
                    585:        }
                    586:        fflush(out);
                    587: }
                    588:
                    589: showmp(vl,p)
                    590: VL vl;
                    591: P p;
                    592: {
                    593:        P t;
                    594:
                    595:        mptop(p,&t); printp(vl,t); fprintf(out,"\n");
                    596: }
                    597: #endif
                    598:
                    599: void premp(vl,p1,p2,pr)
                    600: VL vl;
                    601: P p1,p2,*pr;
                    602: {
                    603:        P m,m1,m2;
                    604:        P *pw;
                    605:        DCP dc;
                    606:        V v1,v2;
                    607:        register int i,j;
                    608:        int n1,n2,d;
                    609:
                    610:        if ( NUM(p1) )
                    611:                if ( NUM(p2) )
                    612:                        *pr = 0;
                    613:                else
                    614:                        *pr = p1;
                    615:        else if ( NUM(p2) )
                    616:                *pr = 0;
                    617:        else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
                    618:                n1 = deg(v1,p1); n2 = deg(v1,p2);
                    619:                pw = (P *)ALLOCA((n1+1)*sizeof(P));
                    620:                bzero((char *)pw,(int)((n1+1)*sizeof(P)));
                    621:
                    622:                for ( dc = DC(p1); dc; dc = NEXT(dc) )
                    623:                        pw[QTOS(DEG(dc))] = COEF(dc);
                    624:
                    625:                for ( i = n1; i >= n2; i-- ) {
                    626:                        if ( pw[i] ) {
                    627:                                m = pw[i];
                    628:                                for ( j = i; j >= 0; j-- ) {
                    629:                                        mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
                    630:                                }
                    631:
                    632:                                for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
                    633:                                        mulp(vl,COEF(dc),m,&m1);
                    634:                                        subp(vl,pw[QTOS(DEG(dc))+d],m1,&m2);
                    635:                                        pw[QTOS(DEG(dc))+d] = m2;
                    636:                                }
                    637:                        } else
                    638:                                for ( j = i; j >= 0; j-- )
                    639:                                        if ( pw[j] ) {
                    640:                                                mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
                    641:                                        }
                    642:                }
                    643:                plisttop(pw,v1,n2-1,pr);
                    644:        } else {
                    645:                while ( v1 != vl->v && v2 != vl->v )
                    646:                        vl = NEXT(vl);
                    647:                if ( v1 == vl->v )
                    648:                        *pr = 0;
                    649:                else
                    650:                        *pr = p1;
                    651:        }
                    652: }
                    653:
                    654: void ptozp0(p,pr)
                    655: P p;
                    656: P *pr;
                    657: {
                    658:        Q c;
                    659:
                    660:        ptozp(p,1,&c,pr);
                    661: }
                    662:
                    663: void mindegp(vl,p,mvlp,pr)
                    664: VL vl,*mvlp;
                    665: P p,*pr;
                    666: {
                    667:        P t;
                    668:        VL nvl,tvl,avl;
                    669:        V v;
                    670:        int n,d;
                    671:
                    672:        clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
                    673:        for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
                    674:                n = getdeg(avl->v,p);
                    675:                if ( n < d ) {
                    676:                        v = avl->v; d = n;
                    677:                }
                    678:        }
                    679:        if ( v != nvl->v ) {
                    680:                reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
                    681:                *pr = t; *mvlp = tvl;
                    682:        } else {
                    683:                *pr = p; *mvlp = nvl;
                    684:        }
                    685: }
                    686:
                    687: void maxdegp(vl,p,mvlp,pr)
                    688: VL vl,*mvlp;
                    689: P p,*pr;
                    690: {
                    691:        P t;
                    692:        VL nvl,tvl,avl;
                    693:        V v;
                    694:        int n,d;
                    695:
                    696:        clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
                    697:        for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
                    698:                n = getdeg(avl->v,p);
                    699:                if ( n > d ) {
                    700:                        v = avl->v; d = n;
                    701:                }
                    702:        }
                    703:        if ( v != nvl->v ) {
                    704:                reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
                    705:                *pr = t; *mvlp = tvl;
                    706:        } else {
                    707:                *pr = p; *mvlp = nvl;
                    708:        }
                    709: }
                    710:
                    711: void min_common_vars_in_coefp(vl,p,mvlp,pr)
                    712: VL vl,*mvlp;
                    713: P p,*pr;
                    714: {
                    715:        P u,p0;
                    716:        VL tvl,cvl,svl,uvl,avl,vl0;
                    717:        int n,n0,d,d0;
                    718:        DCP dc;
                    719:
                    720:        clctv(vl,p,&tvl); vl = tvl;
                    721:        for ( tvl = vl, n0 = 0; tvl; tvl = NEXT(tvl), n0++ );
                    722:        for ( avl = vl; avl; avl = NEXT(avl) ) {
                    723:                if ( VR(p) != avl->v ) {
                    724:                        reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
                    725:                } else {
                    726:                        uvl = vl; u = p;
                    727:                }
                    728:                for ( cvl = NEXT(uvl), dc = DC(u); dc && cvl; dc = NEXT(dc) ) {
                    729:                        clctv(uvl,COEF(dc),&tvl); intersectv(cvl,tvl,&svl); cvl = svl;
                    730:                }
                    731:                if ( !cvl ) {
                    732:                        *mvlp = uvl; *pr = u; return;
                    733:                } else {
                    734:                        for ( tvl = cvl, n = 0; tvl; tvl = NEXT(tvl), n++ );
                    735:                        if ( n < n0 ) {
                    736:                                vl0 = uvl; p0 = u; n0 = n; d0 = getdeg(uvl->v,u);
                    737:                        } else if ( (n == n0) && ((d = getdeg(uvl->v,u)) < d0) ) {
                    738:                                        vl0 = uvl; p0 = u; n0 = n; d0 = d;
                    739:                        }
                    740:                }
                    741:        }
                    742:        *pr = p0; *mvlp = vl0;
                    743: }
                    744:
                    745: void minlcdegp(vl,p,mvlp,pr)
                    746: VL vl,*mvlp;
                    747: P p,*pr;
                    748: {
                    749:        P u,p0;
                    750:        VL tvl,uvl,avl,vl0;
                    751:        int d,d0;
                    752:
                    753:        clctv(vl,p,&tvl); vl = tvl;
                    754:        vl0 = vl; p0 = p; d0 = homdeg(COEF(DC(p)));
                    755:        for ( avl = NEXT(vl); avl; avl = NEXT(avl) ) {
                    756:                reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
                    757:                d = homdeg(COEF(DC(u)));
                    758:                if ( d < d0 ) {
                    759:                        vl0 = uvl; p0 = u; d0 = d;
                    760:                }
                    761:        }
                    762:        *pr = p0; *mvlp = vl0;
                    763: }
                    764:
                    765: void sort_by_deg(n,p,pr)
                    766: int n;
                    767: P *p,*pr;
                    768: {
                    769:        int j,k,d,k0;
                    770:        V v;
                    771:
                    772:        for ( j = 0; j < n; j++ ) {
                    773:                for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
                    774:                        k < n; k++ )
                    775:                        if ( deg(v,p[k]) < d ) {
                    776:                                k0 = k;
                    777:                                d = deg(v,p[k]);
                    778:                        }
                    779:                pr[j] = p[k0];
                    780:                if ( j != k0 )
                    781:                        p[k0] = p[j];
                    782:        }
                    783: }
                    784:
                    785: void sort_by_deg_rev(n,p,pr)
                    786: int n;
                    787: P *p,*pr;
                    788: {
                    789:        int j,k,d,k0;
                    790:        V v;
                    791:
                    792:        for ( j = 0; j < n; j++ ) {
                    793:                for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
                    794:                        k < n; k++ )
                    795:                        if ( deg(v,p[k]) > d ) {
                    796:                                k0 = k;
                    797:                                d = deg(v,p[k]);
                    798:                        }
                    799:                pr[j] = p[k0];
                    800:                if ( j != k0 )
                    801:                        p[k0] = p[j];
                    802:        }
                    803: }
                    804:
                    805:
                    806: void getmindeg(v,p,dp)
                    807: V v;
                    808: P p;
                    809: Q *dp;
                    810: {
                    811:        Q dt,d;
                    812:        DCP dc;
                    813:
                    814:        if ( !p || NUM(p) )
                    815:                *dp = 0;
                    816:        else if ( v == VR(p) ) {
                    817:                for ( dc = DC(p); NEXT(dc); dc = NEXT(dc) );
                    818:                *dp = DEG(dc);
                    819:        } else {
                    820:                dc = DC(p);
                    821:                getmindeg(v,COEF(dc),&d);
                    822:                for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                    823:                        getmindeg(v,COEF(dc),&dt);
                    824:                        if ( cmpq(dt,d) < 0 )
                    825:                                d = dt;
                    826:                }
                    827:                *dp = d;
                    828:        }
                    829: }
                    830:
                    831: void minchdegp(vl,p,mvlp,pr)
                    832: VL vl,*mvlp;
                    833: P p,*pr;
                    834: {
                    835:        P t;
                    836:        VL tvl,nvl,avl;
                    837:        int n,d,z;
                    838:        V v;
                    839:
                    840:        if ( NUM(p) ) {
                    841:                *mvlp = vl;
                    842:                *pr = p;
                    843:                return;
                    844:        }
                    845:        clctv(vl,p,&nvl);
                    846:        v = nvl->v;
                    847:        d = getchomdeg(v,p) + getlchomdeg(v,p,&z);
                    848:        for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
                    849:                n = getchomdeg(avl->v,p) + getlchomdeg(avl->v,p,&z);
                    850:                if ( n < d ) {
                    851:                        v = avl->v; d = n;
                    852:                }
                    853:        }
                    854:        if ( v != nvl->v ) {
                    855:                reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
                    856:                *pr = t; *mvlp = tvl;
                    857:        } else {
                    858:                *pr = p; *mvlp = nvl;
                    859:        }
                    860: }
                    861:
                    862: int getchomdeg(v,p)
                    863: V v;
                    864: P p;
                    865: {
                    866:        int m,m1;
                    867:        DCP dc;
                    868:
                    869:        if ( !p || NUM(p) )
                    870:                return ( 0 );
                    871:        else if ( VR(p) == v )
                    872:                return ( dbound(v,p) );
                    873:        else {
                    874:                for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
                    875:                        m1 = getchomdeg(v,COEF(dc))+QTOS(DEG(dc));
                    876:                        m = MAX(m,m1);
                    877:                }
                    878:                return ( m );
                    879:        }
                    880: }
                    881:
                    882: int getlchomdeg(v,p,d)
                    883: V v;
                    884: P p;
                    885: int *d;
                    886: {
                    887:        int m0,m1,d0,d1;
                    888:        DCP dc;
                    889:
                    890:        if ( !p || NUM(p) ) {
                    891:                *d = 0;
                    892:                return ( 0 );
                    893:        } else if ( VR(p) == v ) {
                    894:                *d = QTOS(DEG(DC(p)));
                    895:                return ( homdeg(LC(p)) );
                    896:        } else {
                    897:                for ( dc = DC(p), m0 = 0, d0 = 0; dc; dc = NEXT(dc) ) {
                    898:                        m1 = getlchomdeg(v,COEF(dc),&d1)+QTOS(DEG(dc));
                    899:                        if ( d1 > d0 ) {
                    900:                                m0 = m1;
                    901:                                d0 = d1;
                    902:                        } else if ( d1 == d0 )
                    903:                                m0 = MAX(m0,m1);
                    904:                }
                    905:                *d = d0;
                    906:                return ( m0 );
                    907:        }
                    908: }
                    909:
                    910: int nmonop(p)
                    911: P p;
                    912: {
                    913:        int s;
                    914:        DCP dc;
                    915:
                    916:        if ( !p )
                    917:                return 0;
                    918:        else
                    919:                switch ( OID(p) ) {
                    920:                        case O_N:
                    921:                                return 1; break;
                    922:                        case O_P:
                    923:                                for ( dc = DC((P)p), s = 0; dc; dc = NEXT(dc) )
                    924:                                        s += nmonop(COEF(dc));
                    925:                                return s; break;
                    926:                        default:
                    927:                                return 0;
                    928:                }
                    929: }
                    930:
                    931: int qpcheck(p)
                    932: Obj p;
                    933: {
                    934:        DCP dc;
                    935:
                    936:        if ( !p )
                    937:                return 1;
                    938:        else
                    939:                switch ( OID(p) ) {
                    940:                        case O_N:
                    941:                                return RATN((Num)p)?1:0;
                    942:                        case O_P:
                    943:                                for ( dc = DC((P)p); dc; dc = NEXT(dc) )
                    944:                                        if ( !qpcheck((Obj)COEF(dc)) )
                    945:                                                return 0;
                    946:                                return 1;
                    947:                        default:
                    948:                                return 0;
                    949:                }
                    950: }
                    951:
                    952: /* check if p is univariate and all coeffs are INT or LM */
                    953:
                    954: int uzpcheck(p)
                    955: Obj p;
                    956: {
                    957:        DCP dc;
                    958:        P c;
                    959:
                    960:        if ( !p )
                    961:                return 1;
                    962:        else
                    963:                switch ( OID(p) ) {
                    964:                        case O_N:
                    965:                                return (RATN((Num)p)&&INT(p))||(NID((Num)p)==N_LM);
                    966:                        case O_P:
                    967:                                for ( dc = DC((P)p); dc; dc = NEXT(dc) ) {
                    968:                                        c = COEF(dc);
                    969:                                        if ( !NUM(c) || !uzpcheck(c) )
                    970:                                                return 0;
                    971:                                }
                    972:                                return 1;
                    973:                        default:
                    974:                                return 0;
                    975:                }
                    976: }
                    977:
                    978: int p_mag(p)
                    979: P p;
                    980: {
                    981:        int s;
                    982:        DCP dc;
                    983:
                    984:        if ( !p )
                    985:                return 0;
                    986:        else if ( OID(p) == O_N )
                    987:                return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
                    988:        else {
                    989:                for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) )
                    990:                        s += p_mag(COEF(dc));
                    991:                return s;
                    992:        }
                    993: }

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