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

Annotation of OpenXM_contrib2/asir2018/engine/PU.c, Revision 1.2

1.1       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
                     26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
                     27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.2     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2018/engine/PU.c,v 1.1 2018/09/19 05:45:07 noro Exp $
1.1       noro       49: */
                     50: #include "ca.h"
                     51:
                     52: void reorderp(VL nvl,VL ovl,P p,P *pr)
                     53: {
                     54:   DCP dc;
                     55:   P x,m,s,t,c;
                     56:
                     57:   if ( !p )
                     58:     *pr = 0;
                     59:   else if ( NUM(p) )
                     60:     *pr = p;
                     61:   else {
                     62:     MKV(VR(p),x);
                     63:     for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                     64:       reorderp(nvl,ovl,COEF(dc),&c);
                     65:       if ( DEG(dc) ) {
                     66:         pwrp(nvl,x,DEG(dc),&t); mulp(nvl,c,t,&m);
                     67:         addp(nvl,m,s,&t);
                     68:       } else
                     69:         addp(nvl,s,c,&t);
                     70:       s = t;
                     71:     }
                     72:     *pr = s;
                     73:   }
                     74: }
                     75:
                     76: void substp(VL vl,P p,V v0,P p0,P *pr)
                     77: {
                     78:   P x,t,m,c,s,a;
                     79:   DCP dc;
                     80:   Z d;
                     81:
                     82:   if ( !p )
                     83:     *pr = 0;
                     84:   else if ( NUM(p) )
                     85:     *pr = p;
                     86:   else if ( VR(p) != v0 ) {
                     87:     MKV(VR(p),x);
                     88:     for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                     89:       substp(vl,COEF(dc),v0,p0,&t);
                     90:       if ( DEG(dc) ) {
                     91:         pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
                     92:         addp(vl,m,c,&a);
                     93:         c = a;
                     94:       } else {
                     95:         addp(vl,t,c,&a);
                     96:         c = a;
                     97:       }
                     98:     }
                     99:     *pr = c;
                    100:   } else {
                    101:     dc = DC(p);
                    102:     c = COEF(dc);
                    103:     for ( d = DEG(dc), dc = NEXT(dc);
                    104:       dc; d = DEG(dc), dc = NEXT(dc) ) {
                    105:         subz(d,DEG(dc),(Z *)&t); pwrp(vl,p0,(Z)t,&s);
                    106:         mulp(vl,s,c,&m);
                    107:         addp(vl,m,COEF(dc),&c);
                    108:     }
                    109:     if ( d ) {
                    110:       pwrp(vl,p0,d,&t); mulp(vl,t,c,&m);
                    111:       c = m;
                    112:     }
                    113:     *pr = c;
                    114:   }
                    115: }
                    116:
                    117: void substpp(VL vl,P p,V *vvect,P *svect,int nv,P *pr);
                    118:
                    119: void substpp(VL vl,P p,V *vvect,P *svect,int nv,P *pr)
                    120: {
                    121:   P x,t,m,c,s,a,p0,c1;
                    122:   DCP dc;
                    123:   Z d;
                    124:   V v;
                    125:   int i;
                    126:
                    127:   if ( !p )
                    128:     *pr = 0;
                    129:   else if ( NUM(p) )
                    130:     *pr = p;
                    131:   else {
                    132:     v = VR(p);
                    133:     for ( i = 0; i < nv; i++ ) if ( vvect[i] == v ) break;
                    134:     if ( svect[i] && OID(svect[i]) < 0 ) {
                    135:       MKV(VR(p),x);
                    136:       for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                    137:         substpp(vl,COEF(dc),vvect,svect,nv,&t);
                    138:         if ( DEG(dc) ) {
                    139:           pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
                    140:           addp(vl,m,c,&a);
                    141:           c = a;
                    142:         } else {
                    143:           addp(vl,t,c,&a);
                    144:           c = a;
                    145:         }
                    146:       }
                    147:       *pr = c;
                    148:     } else {
                    149:       p0 = svect[i];
                    150:       dc = DC(p);
                    151:       substpp(vl,COEF(dc),vvect,svect,nv,&c);
                    152:       for ( d = DEG(dc), dc = NEXT(dc);
                    153:         dc; d = DEG(dc), dc = NEXT(dc) ) {
                    154:           subz(d,DEG(dc),(Z *)&t); pwrp(vl,p0,(Z)t,&s);
                    155:           mulp(vl,s,c,&m);
                    156:           substpp(vl,COEF(dc),vvect,svect,nv,&c1);
                    157:           addp(vl,m,c1,&c);
                    158:       }
                    159:       if ( d ) {
                    160:         pwrp(vl,p0,d,&t); mulp(vl,t,c,&m);
                    161:         c = m;
                    162:       }
                    163:       *pr = c;
                    164:     }
                    165:   }
                    166: }
                    167:
                    168: void detp(VL vl,P **rmat,int n,P *dp)
                    169: {
                    170:   int i,j,k,l,sgn,nmin,kmin,lmin,ntmp;
                    171:   P mjj,mij,t,s,u,d;
                    172:   P **mat;
                    173:   P *mi,*mj;
                    174:
                    175:   mat = (P **)almat_pointer(n,n);
                    176:   for ( i = 0; i < n; i++ )
                    177:     for ( j = 0; j < n; j++ )
                    178:       mat[i][j] = rmat[i][j];
                    179:   for ( j = 0, d = (P)ONE, sgn = 1; j < n; j++ ) {
                    180:     for ( i = j; (i < n) && !mat[i][j]; i++ );
                    181:     if ( i == n ) {
                    182:       *dp = 0; return;
                    183:     }
                    184:     nmin = nmonop(mat[i][j]);
                    185:     kmin=i; lmin=j;
                    186:     for ( k = j; k < n; k++ )
                    187:       for ( l = j; l < n; l++ )
                    188:         if ( mat[k][l] && ((ntmp=nmonop(mat[k][l])) < nmin) ) {
                    189:           kmin = k; lmin = l; nmin = ntmp;
                    190:         }
                    191:     if ( kmin != j ) {
                    192:       mj = mat[j]; mat[j] = mat[kmin]; mat[kmin] = mj; sgn = -sgn;
                    193:     }
                    194:     if ( lmin != j ) {
                    195:       for ( k = j; k < n; k++ ) {
                    196:         t = mat[k][j]; mat[k][j] = mat[k][lmin]; mat[k][lmin] = t;
                    197:       }
                    198:       sgn = -sgn;
                    199:     }
                    200:     for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
                    201:       for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
                    202:         mulp(vl,mi[k],mjj,&t); mulp(vl,mj[k],mij,&s);
                    203:         subp(vl,t,s,&u); divsp(vl,u,d,&mi[k]);
                    204:       }
                    205:     d = mjj;
                    206:   }
                    207:   if ( sgn < 0 )
                    208:     chsgnp(d,dp);
                    209:   else
                    210:     *dp = d;
                    211: }
                    212:
                    213: void invmatp(VL vl,P **rmat,int n,P ***imatp,P *dnp)
                    214: {
                    215:   int i,j,k,l,n2;
                    216:   P mjj,mij,t,s,u,d;
                    217:   P **mat,**imat;
                    218:   P *mi,*mj,*w;
                    219:
                    220:   n2 = n<<1;
                    221:   mat = (P **)almat_pointer(n,n2);
                    222:   for ( i = 0; i < n; i++ ) {
                    223:     for ( j = 0; j < n; j++ )
                    224:       mat[i][j] = rmat[i][j];
                    225:     mat[i][i+n] = (P)ONE;
                    226:   }
                    227:   for ( j = 0, d = (P)ONE; j < n; j++ ) {
                    228:     for ( i = j; (i < n) && !mat[i][j]; i++ );
                    229:     if ( i == n ) {
                    230:       error("invmatp : input is singular");
                    231:     }
                    232:     for ( k = i; k < n; k++ )
                    233:       if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
                    234:         i = k;
                    235:     if ( j != i ) {
                    236:       mj = mat[j]; mat[j] = mat[i]; mat[i] = mj;
                    237:     }
                    238:     for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
                    239:       for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n2; k++ ) {
                    240:         mulp(vl,mi[k],mjj,&t); mulp(vl,mj[k],mij,&s);
                    241:         subp(vl,t,s,&u); divsp(vl,u,d,&mi[k]);
                    242:       }
                    243:     d = mjj;
                    244:   }
                    245:   /* backward substitution */
                    246:   w = (P *)ALLOCA(n2*sizeof(P));
                    247:   for ( i = n-2; i >= 0; i-- ) {
                    248:     bzero(w,n2*sizeof(P));
                    249:     for ( k = i+1; k < n; k++ )
                    250:       for ( l = k, u = mat[i][l]; l < n2; l++ ) {
                    251:         mulp(vl,mat[k][l],u,&t); addp(vl,w[l],t,&s); w[l] = s;
                    252:       }
                    253:     for ( j = i, u = mat[i][j]; j < n2; j++ ) {
                    254:       mulp(vl,mat[i][j],d,&t); subp(vl,t,w[j],&s);
                    255:       divsp(vl,s,u,&mat[i][j]);
                    256:     }
                    257:   }
                    258:   imat = (P **)almat_pointer(n,n);
                    259:   for ( i = 0; i < n; i++ )
                    260:     for ( j = 0; j < n; j++ )
                    261:       imat[i][j] = mat[i][j+n];
                    262:   *imatp = imat;
                    263:   *dnp = d;
                    264: }
                    265:
                    266: void reordvar(VL vl,V v,VL *nvlp)
                    267: {
                    268:   VL nvl,nvl0;
                    269:
                    270:   for ( NEWVL(nvl0), nvl0->v = v, nvl = nvl0;
                    271:     vl; vl = NEXT(vl) )
                    272:     if ( vl->v == v )
                    273:       continue;
                    274:     else {
                    275:       NEWVL(NEXT(nvl));
                    276:       nvl = NEXT(nvl);
                    277:       nvl->v = vl->v;
                    278:     }
                    279:   NEXT(nvl) = 0;
                    280:   *nvlp = nvl0;
                    281: }
                    282:
                    283: void gcdprsp(VL vl,P p1,P p2,P *pr)
                    284: {
                    285:   P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
                    286:   V v1,v2;
                    287:
                    288:   if ( !p1 )
                    289:     ptozp0(p2,pr);
                    290:   else if ( !p2 )
                    291:     ptozp0(p1,pr);
                    292:   else if ( NUM(p1) || NUM(p2) )
                    293:     *pr = (P)ONE;
                    294:   else {
                    295:     ptozp0(p1,&g1); ptozp0(p2,&g2);
                    296:     if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {
                    297:       gcdcp(vl,g1,&gc1); divsp(vl,g1,gc1,&gp1);
                    298:       gcdcp(vl,g2,&gc2); divsp(vl,g2,gc2,&gp2);
                    299:       gcdprsp(vl,gc1,gc2,&gcr);
                    300:       sprs(vl,v1,gp1,gp2,&g);
                    301:
                    302:       if ( VR(g) == v1 ) {
                    303:         ptozp0(g,&gp);
                    304:         gcdcp(vl,gp,&gc); divsp(vl,gp,gc,&gp1);
                    305:         mulp(vl,gp1,gcr,pr);
                    306:
                    307:       } else
                    308:         *pr = gcr;
                    309:     } else {
                    310:       while ( v1 != vl->v && v2 != vl->v )
                    311:         vl = NEXT(vl);
                    312:       if ( v1 == vl->v ) {
                    313:         gcdcp(vl,g1,&gc1); gcdprsp(vl,gc1,g2,pr);
                    314:       } else {
                    315:         gcdcp(vl,g2,&gc2); gcdprsp(vl,gc2,g1,pr);
                    316:       }
                    317:     }
                    318:   }
                    319: }
                    320:
                    321: void gcdcp(VL vl,P p,P *pr)
                    322: {
                    323:   P g,g1;
                    324:   DCP dc;
                    325:
                    326:   if ( NUM(p) )
                    327:     *pr = (P)ONE;
                    328:   else {
                    329:     dc = DC(p);
                    330:     g = COEF(dc);
                    331:     for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                    332:       gcdprsp(vl,g,COEF(dc),&g1);
                    333:       g = g1;
                    334:     }
                    335:     *pr = g;
                    336:   }
                    337: }
                    338:
                    339: void sprs(VL vl,V v,P p1,P p2,P *pr)
                    340: {
                    341:   P q1,q2,m,m1,m2,x,h,r,g1,g2;
                    342:   int d;
                    343:   Z dq;
                    344:   VL nvl;
                    345:
                    346:   reordvar(vl,v,&nvl);
                    347:   reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    348:
                    349:   if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
                    350:     *pr = 0;
                    351:     return;
                    352:   }
                    353:
                    354:   if ( deg(v,q1) >= deg(v,q2) ) {
                    355:     g1 = q1; g2 = q2;
                    356:   } else {
                    357:     g2 = q1; g1 = q2;
                    358:   }
                    359:
                    360:   for ( h = (P)ONE, x = (P)ONE; ; ) {
                    361:     if ( !deg(v,g2) )
                    362:       break;
                    363:
                    364:     premp(nvl,g1,g2,&r);
                    365:     if ( !r )
                    366:       break;
                    367:
1.2     ! noro      368:     d = deg(v,g1) - deg(v,g2); STOZ(d,dq);
1.1       noro      369:     pwrp(nvl,h,dq,&m); mulp(nvl,m,x,&m1); g1 = g2;
                    370:     divsp(nvl,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */
                    371:     pwrp(nvl,x,dq,&m1); mulp(nvl,m1,h,&m2);
                    372:     divsp(nvl,m2,m,&h);
                    373:   }
                    374:   *pr = g2;
                    375: }
                    376:
                    377: void resultp(VL vl,V v,P p1,P p2,P *pr)
                    378: {
                    379:   P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;
                    380:   int d,d1,d2,j,k;
                    381:   VL nvl;
                    382:   Z dq;
                    383:
                    384:   if ( !p1 || !p2 ) {
                    385:     *pr = 0; return;
                    386:   }
                    387:   reordvar(vl,v,&nvl);
                    388:   reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    389:
                    390:   if ( VR(q1) != v )
                    391:     if ( VR(q2) != v ) {
                    392:       *pr = 0;
                    393:       return;
                    394:     } else {
1.2     ! noro      395:       d = deg(v,q2); STOZ(d,dq);
1.1       noro      396:       pwrp(vl,q1,dq,pr);
                    397:       return;
                    398:     }
                    399:   else if ( VR(q2) != v ) {
1.2     ! noro      400:     d = deg(v,q1); STOZ(d,dq);
1.1       noro      401:     pwrp(vl,q2,dq,pr);
                    402:     return;
                    403:   }
                    404:
                    405:   if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
                    406:     *pr = 0;
                    407:     return;
                    408:   }
                    409:
                    410:   d1 = deg(v,q1); d2 = deg(v,q2);
                    411:   if ( d1 > d2 ) {
                    412:     g1 = q1; g2 = q2; adj = (P)ONE;
                    413:   } else if ( d1 < d2 ) {
                    414:     g2 = q1; g1 = q2;
                    415:     if ( (d1 % 2) && (d2 % 2) ) {
1.2     ! noro      416:       STOZ(-1,dq); adj = (P)dq;
1.1       noro      417:     } else
                    418:       adj = (P)ONE;
                    419:   } else {
                    420:     premp(nvl,q1,q2,&t);
1.2     ! noro      421:     d = deg(v,t); STOZ(d,dq); pwrp(nvl,LC(q2),dq,&adj);
1.1       noro      422:     g1 = q2; g2 = t;
                    423:     if ( d1 % 2 ) {
                    424:       chsgnp(adj,&t); adj = t;
                    425:     }
                    426:   }
                    427:   d1 = deg(v,g1); j = d1 - 1;
                    428:
                    429:   for ( lc = (P)ONE; ; ) {
                    430:     if ( ( k = deg(v,g2) ) < 0 ) {
                    431:       *pr = 0;
                    432:       return;
                    433:     }
                    434:
                    435:     if ( k == j )
                    436:       if ( !k ) {
                    437:         divsp(nvl,g2,adj,pr);
                    438:         return;
                    439:       } else {
                    440:         premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
                    441:         divsp(nvl,r,m,&q);
                    442:         g1 = g2; g2 = q;
                    443:         lc = LC(g1); /* g1 is not const */
                    444:         j = k - 1;
                    445:       }
                    446:     else {
1.2     ! noro      447:       d = j - k; STOZ(d,dq);
1.1       noro      448:       pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
                    449:       mulp(nvl,g2,m,&m1);
                    450:       pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
                    451:       if ( k == 0 ) {
                    452:         divsp(nvl,t,adj,pr);
                    453:         return;
                    454:       } else {
                    455:         premp(nvl,g1,g2,&r);
                    456:         mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
                    457:         divsp(nvl,r,m2,&q);
                    458:         g1 = t; g2 = q;
                    459:         if ( d % 2 ) {
                    460:           chsgnp(g2,&t); g2 = t;
                    461:         }
                    462:         lc = LC(g1); /* g1 is not const */
                    463:         j = k - 1;
                    464:       }
                    465:     }
                    466:   }
                    467: }
                    468:
                    469: void srch2(VL vl,V v,P p1,P p2,P *pr)
                    470: {
                    471:   P q1,q2,m,m1,m2,lc,q,r,t,s,g1,g2,adj;
                    472:   int d,d1,d2,j,k;
                    473:   VL nvl;
                    474:   Z dq;
                    475:
                    476:   reordvar(vl,v,&nvl);
                    477:   reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    478:
                    479:   if ( VR(q1) != v )
                    480:     if ( VR(q2) != v ) {
                    481:       *pr = 0;
                    482:       return;
                    483:     } else {
1.2     ! noro      484:       d = deg(v,q2); STOZ(d,dq);
1.1       noro      485:       pwrp(vl,q1,dq,pr);
                    486:       return;
                    487:     }
                    488:   else if ( VR(q2) != v ) {
1.2     ! noro      489:     d = deg(v,q1); STOZ(d,dq);
1.1       noro      490:     pwrp(vl,q2,dq,pr);
                    491:     return;
                    492:   }
                    493:
                    494:   if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
                    495:     *pr = 0;
                    496:     return;
                    497:   }
                    498:
                    499:   if ( deg(v,q1) >= deg(v,q2) ) {
                    500:     g1 = q1; g2 = q2;
                    501:   } else {
                    502:     g2 = q1; g1 = q2;
                    503:   }
                    504:
                    505:
                    506:   if ( ( d1 = deg(v,g1) ) > ( d2 = deg(v,g2) ) ) {
                    507:     j = d1 - 1;
                    508:     adj = (P)ONE;
                    509:   } else {
                    510:     premp(nvl,g1,g2,&t);
1.2     ! noro      511:     d = deg(v,t); STOZ(d,dq);
1.1       noro      512:     pwrp(nvl,LC(g2),dq,&adj);
                    513:     g1 = g2; g2 = t;
                    514:     j = deg(v,g1) - 1;
                    515:   }
                    516:
                    517:   for ( lc = (P)ONE; ; ) {
                    518:     if ( ( k = deg(v,g2) ) < 0 ) {
                    519:       *pr = 0;
                    520:       return;
                    521:     }
                    522:
                    523:     ptozp(g1,1,(Q *)&t,&s); g1 = s; ptozp(g2,1,(Q *)&t,&s); g2 = s;
                    524:     if ( k == j )
                    525:       if ( !k ) {
                    526:         divsp(nvl,g2,adj,pr);
                    527:         return;
                    528:       } else {
                    529:         premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
                    530:         divsp(nvl,r,m,&q);
                    531:         g1 = g2; g2 = q;
                    532:         lc = LC(g1); /* g1 is not const */
                    533:         j = k - 1;
                    534:       }
                    535:     else {
1.2     ! noro      536:       d = j - k; STOZ(d,dq);
1.1       noro      537:       pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
                    538:       mulp(nvl,g2,m,&m1);
                    539:       pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
                    540:       if ( k == 0 ) {
                    541:         divsp(nvl,t,adj,pr);
                    542:         return;
                    543:       } else {
                    544:         premp(nvl,g1,g2,&r);
                    545:         mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
                    546:         divsp(nvl,r,m2,&q);
                    547:         g1 = t; g2 = q;
                    548:         if ( d % 2 ) {
                    549:           chsgnp(g2,&t); g2 = t;
                    550:         }
                    551:         lc = LC(g1); /* g1 is not const */
                    552:         j = k - 1;
                    553:       }
                    554:     }
                    555:   }
                    556: }
                    557:
                    558: void srcr(VL vl,V v,P p1,P p2,P *pr)
                    559: {
                    560:   P q1,q2,c,c1;
                    561:   P tg,tg1,tg2,resg;
                    562:   int index,mod;
                    563:   Q a,b,f,q,s,u,w;
                    564:   Z n,m,t;
                    565:   VL nvl;
                    566:
                    567:   reordvar(vl,v,&nvl);
                    568:   reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    569:
                    570:   if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
                    571:     *pr = 0;
                    572:     return;
                    573:   }
                    574:   if ( VR(q1) != v ) {
                    575:     pwrp(vl,q1,DEG(DC(q2)),pr);
                    576:     return;
                    577:   }
                    578:   if ( VR(q2) != v ) {
                    579:     pwrp(vl,q2,DEG(DC(q1)),pr);
                    580:     return;
                    581:   }
                    582:   norm1c(q1,&a); norm1c(q2,&b);
                    583:   n = DEG(DC(q1)); m = DEG(DC(q2));
                    584:   pwrq(a,(Q)m,&w); pwrq(b,(Q)n,&s); mulq(w,s,&u);
1.2     ! noro      585:   factorialz(ZTOS(n)+ZTOS(m),&t);
1.1       noro      586:   mulq(u,(Q)t,&s); addq(s,s,&f);
                    587:   for ( index = 0, q = (Q)ONE, c = 0; cmpq(f,q) >= 0; ) {
                    588:     mod = get_lprime(index++);
                    589:     ptomp(mod,LC(q1),&tg);
                    590:     if ( !tg )
                    591:       continue;
                    592:     ptomp(mod,LC(q2),&tg);
                    593:     if ( !tg )
                    594:       continue;
                    595:     ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
                    596:     srchmp(nvl,mod,v,tg1,tg2,&resg);
                    597:     chnremp(nvl,mod,c,(Z)q,resg,&c1); c = c1;
1.2     ! noro      598:     STOZ(mod,t); mulq(q,(Q)t,&s); q = s;
1.1       noro      599:   }
                    600:   *pr = c;
                    601: }
                    602:
                    603: void res_ch_det(VL vl,V v,P p1,P p2,P *pr)
                    604: {
                    605:   P q1,q2,c,c1;
                    606:   P tg,tg1,tg2,resg;
                    607:   int index,mod;
                    608:   Q a,b,f,q,s,u,w;
                    609:   Z n,m,t;
                    610:   VL nvl;
                    611:
                    612:   reordvar(vl,v,&nvl);
                    613:   reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
                    614:
                    615:   if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
                    616:     *pr = 0;
                    617:     return;
                    618:   }
                    619:   if ( VR(q1) != v ) {
                    620:     pwrp(vl,q1,DEG(DC(q2)),pr);
                    621:     return;
                    622:   }
                    623:   if ( VR(q2) != v ) {
                    624:     pwrp(vl,q2,DEG(DC(q1)),pr);
                    625:     return;
                    626:   }
                    627:   norm1c(q1,&a); norm1c(q2,&b);
                    628:   n = DEG(DC(q1)); m = DEG(DC(q2));
                    629:   pwrq(a,(Q)m,&w); pwrq(b,(Q)n,&s); mulq(w,s,&u);
1.2     ! noro      630:   factorialz(ZTOS(n)+ZTOS(m),&t);
1.1       noro      631:   mulq(u,(Q)t,&s); addq(s,s,&f);
                    632:   for ( index = 0, q = (Q)ONE, c = 0; cmpq(f,q) >= 0; ) {
                    633:     mod = get_lprime(index++);
                    634:     ptomp(mod,LC(q1),&tg);
                    635:     if ( !tg )
                    636:       continue;
                    637:     ptomp(mod,LC(q2),&tg);
                    638:     if ( !tg )
                    639:       continue;
                    640:     ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
                    641:     res_detmp(nvl,mod,v,tg1,tg2,&resg);
                    642:     chnremp(nvl,mod,c,(Z)q,resg,&c1); c = c1;
1.2     ! noro      643:     STOZ(mod,t); mulq(q,(Q)t,&s); q = s;
1.1       noro      644:   }
                    645:   *pr = c;
                    646: }
                    647:
                    648: void res_detmp(VL vl,int mod,V v,P p1,P p2,P *dp)
                    649: {
                    650:   int n1,n2,n,sgn;
                    651:   int i,j,k;
                    652:   P mjj,mij,t,s,u,d;
                    653:   P **mat;
                    654:   P *mi,*mj;
                    655:   DCP dc;
                    656:
                    657:   n1 = UDEG(p1); n2 = UDEG(p2); n = n1+n2;
                    658:   mat = (P **)almat_pointer(n,n);
                    659:   for ( dc = DC(p1); dc; dc = NEXT(dc) )
1.2     ! noro      660:     mat[0][n1-ZTOS(DEG(dc))] = COEF(dc);
1.1       noro      661:   for ( i = 1; i < n2; i++ )
                    662:     for ( j = 0; j <= n1; j++ )
                    663:       mat[i][i+j] = mat[0][j];
                    664:   for ( dc = DC(p2); dc; dc = NEXT(dc) )
1.2     ! noro      665:     mat[n2][n2-ZTOS(DEG(dc))] = COEF(dc);
1.1       noro      666:   for ( i = 1; i < n1; i++ )
                    667:     for ( j = 0; j <= n2; j++ )
                    668:       mat[i+n2][i+j] = mat[n2][j];
                    669:   for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
                    670:     for ( i = j; (i < n) && !mat[i][j]; i++ );
                    671:     if ( i == n ) {
                    672:       *dp = 0; return;
                    673:     }
                    674:     for ( k = i; k < n; k++ )
                    675:       if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
                    676:         i = k;
                    677:     if ( j != i ) {
                    678:       mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
                    679:     }
                    680:     for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
                    681:       for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
                    682:         mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
                    683:         submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
                    684:       }
                    685:     d = mjj;
                    686:   }
                    687:   if ( sgn > 0 )
                    688:     *dp = d;
                    689:   else
                    690:     chsgnmp(mod,d,dp);
                    691: }
                    692:
                    693: #if 0
                    694: showmat(VL vl,P **mat,int n)
                    695: {
                    696:   int i,j;
                    697:   P t;
                    698:
                    699:   for ( i = 0; i < n; i++ ) {
                    700:     for ( j = 0; j < n; j++ ) {
                    701:       mptop(mat[i][j],&t); asir_printp(vl,t); fprintf(out," ");
                    702:     }
                    703:     fprintf(out,"\n");
                    704:   }
                    705:   fflush(out);
                    706: }
                    707:
                    708: showmp(VL vl,P p)
                    709: {
                    710:   P t;
                    711:
                    712:   mptop(p,&t); asir_printp(vl,t); fprintf(out,"\n");
                    713: }
                    714: #endif
                    715:
                    716: void premp(VL vl,P p1,P p2,P *pr)
                    717: {
                    718:   P m,m1,m2;
                    719:   P *pw;
                    720:   DCP dc;
                    721:   V v1,v2;
                    722:   register int i,j;
                    723:   int n1,n2,d;
                    724:
                    725:   if ( NUM(p1) )
                    726:     if ( NUM(p2) )
                    727:       *pr = 0;
                    728:     else
                    729:       *pr = p1;
                    730:   else if ( NUM(p2) )
                    731:     *pr = 0;
                    732:   else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
                    733:     n1 = deg(v1,p1); n2 = deg(v1,p2);
                    734:     if ( n1 < n2 ) {
                    735:       *pr = p1;
                    736:       return;
                    737:     }
                    738:     pw = (P *)ALLOCA((n1+1)*sizeof(P));
                    739:     bzero((char *)pw,(int)((n1+1)*sizeof(P)));
                    740:
                    741:     for ( dc = DC(p1); dc; dc = NEXT(dc) )
1.2     ! noro      742:       pw[ZTOS(DEG(dc))] = COEF(dc);
1.1       noro      743:
                    744:     for ( i = n1; i >= n2; i-- ) {
                    745:       if ( pw[i] ) {
                    746:         m = pw[i];
                    747:         for ( j = i; j >= 0; j-- ) {
                    748:           mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
                    749:         }
                    750:
                    751:         for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
                    752:           mulp(vl,COEF(dc),m,&m1);
1.2     ! noro      753:           subp(vl,pw[ZTOS(DEG(dc))+d],m1,&m2);
        !           754:           pw[ZTOS(DEG(dc))+d] = m2;
1.1       noro      755:         }
                    756:       } else
                    757:         for ( j = i; j >= 0; j-- )
                    758:           if ( pw[j] ) {
                    759:             mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
                    760:           }
                    761:     }
                    762:     plisttop(pw,v1,n2-1,pr);
                    763:   } else {
                    764:     while ( v1 != vl->v && v2 != vl->v )
                    765:       vl = NEXT(vl);
                    766:     if ( v1 == vl->v )
                    767:       *pr = 0;
                    768:     else
                    769:       *pr = p1;
                    770:   }
                    771: }
                    772:
                    773: void ptozp0(P p,P *pr)
                    774: {
                    775:   Q c;
                    776:
                    777:   if ( qpcheck((Obj)p) )
                    778:     ptozp(p,1,&c,pr);
                    779:   else
                    780:     *pr = p;
                    781: }
                    782:
                    783: void mindegp(VL vl,P p,VL *mvlp,P *pr)
                    784: {
                    785:   P t;
                    786:   VL nvl,tvl,avl;
                    787:   V v;
                    788:   int n,d;
                    789:
                    790:   clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
                    791:   for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
                    792:     n = getdeg(avl->v,p);
                    793:     if ( n < d ) {
                    794:       v = avl->v; d = n;
                    795:     }
                    796:   }
                    797:   if ( v != nvl->v ) {
                    798:     reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
                    799:     *pr = t; *mvlp = tvl;
                    800:   } else {
                    801:     *pr = p; *mvlp = nvl;
                    802:   }
                    803: }
                    804:
                    805: void maxdegp(VL vl,P p,VL *mvlp,P *pr)
                    806: {
                    807:   P t;
                    808:   VL nvl,tvl,avl;
                    809:   V v;
                    810:   int n,d;
                    811:
                    812:   clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
                    813:   for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
                    814:     n = getdeg(avl->v,p);
                    815:     if ( n > d ) {
                    816:       v = avl->v; d = n;
                    817:     }
                    818:   }
                    819:   if ( v != nvl->v ) {
                    820:     reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
                    821:     *pr = t; *mvlp = tvl;
                    822:   } else {
                    823:     *pr = p; *mvlp = nvl;
                    824:   }
                    825: }
                    826:
                    827: void min_common_vars_in_coefp(VL vl,P p,VL *mvlp,P *pr)
                    828: {
                    829:   P u,p0;
                    830:   VL tvl,cvl,svl,uvl,avl,vl0;
                    831:   int n,n0,d,d0;
                    832:   DCP dc;
                    833:
                    834:   clctv(vl,p,&tvl); vl = tvl;
                    835:   for ( tvl = vl, n0 = 0; tvl; tvl = NEXT(tvl), n0++ );
                    836:   for ( avl = vl; avl; avl = NEXT(avl) ) {
                    837:     if ( VR(p) != avl->v ) {
                    838:       reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
                    839:     } else {
                    840:       uvl = vl; u = p;
                    841:     }
                    842:     for ( cvl = NEXT(uvl), dc = DC(u); dc && cvl; dc = NEXT(dc) ) {
                    843:       clctv(uvl,COEF(dc),&tvl); intersectv(cvl,tvl,&svl); cvl = svl;
                    844:     }
                    845:     if ( !cvl ) {
                    846:       *mvlp = uvl; *pr = u; return;
                    847:     } else {
                    848:       for ( tvl = cvl, n = 0; tvl; tvl = NEXT(tvl), n++ );
                    849:       if ( n < n0 ) {
                    850:         vl0 = uvl; p0 = u; n0 = n; d0 = getdeg(uvl->v,u);
                    851:       } else if ( (n == n0) && ((d = getdeg(uvl->v,u)) < d0) ) {
                    852:           vl0 = uvl; p0 = u; n0 = n; d0 = d;
                    853:       }
                    854:     }
                    855:   }
                    856:   *pr = p0; *mvlp = vl0;
                    857: }
                    858:
                    859: void minlcdegp(VL vl,P p,VL *mvlp,P *pr)
                    860: {
                    861:   P u,p0;
                    862:   VL tvl,uvl,avl,vl0;
                    863:   int d,d0;
                    864:
                    865:   clctv(vl,p,&tvl); vl = tvl;
                    866:   vl0 = vl; p0 = p; d0 = homdeg(COEF(DC(p)));
                    867:   for ( avl = NEXT(vl); avl; avl = NEXT(avl) ) {
                    868:     reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
                    869:     d = homdeg(COEF(DC(u)));
                    870:     if ( d < d0 ) {
                    871:       vl0 = uvl; p0 = u; d0 = d;
                    872:     }
                    873:   }
                    874:   *pr = p0; *mvlp = vl0;
                    875: }
                    876:
                    877: void sort_by_deg(int n,P *p,P *pr)
                    878: {
                    879:   int j,k,d,k0;
                    880:   V v;
                    881:
                    882:   for ( j = 0; j < n; j++ ) {
                    883:     for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
                    884:       k < n; k++ )
                    885:       if ( deg(v,p[k]) < d ) {
                    886:         k0 = k;
                    887:         d = deg(v,p[k]);
                    888:       }
                    889:     pr[j] = p[k0];
                    890:     if ( j != k0 )
                    891:       p[k0] = p[j];
                    892:   }
                    893: }
                    894:
                    895: void sort_by_deg_rev(int n,P *p,P *pr)
                    896: {
                    897:   int j,k,d,k0;
                    898:   V v;
                    899:
                    900:   for ( j = 0; j < n; j++ ) {
                    901:     for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
                    902:       k < n; k++ )
                    903:       if ( deg(v,p[k]) > d ) {
                    904:         k0 = k;
                    905:         d = deg(v,p[k]);
                    906:       }
                    907:     pr[j] = p[k0];
                    908:     if ( j != k0 )
                    909:       p[k0] = p[j];
                    910:   }
                    911: }
                    912:
                    913:
                    914: void getmindeg(V v,P p,Z *dp)
                    915: {
                    916:   Z dt,d;
                    917:   DCP dc;
                    918:
                    919:   if ( !p || NUM(p) )
                    920:     *dp = 0;
                    921:   else if ( v == VR(p) ) {
                    922:     for ( dc = DC(p); NEXT(dc); dc = NEXT(dc) );
                    923:     *dp = DEG(dc);
                    924:   } else {
                    925:     dc = DC(p);
                    926:     getmindeg(v,COEF(dc),&d);
                    927:     for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                    928:       getmindeg(v,COEF(dc),&dt);
                    929:       if ( cmpz(dt,d) < 0 )
                    930:         d = dt;
                    931:     }
                    932:     *dp = d;
                    933:   }
                    934: }
                    935:
                    936: void minchdegp(VL vl,P p,VL *mvlp,P *pr)
                    937: {
                    938:   P t;
                    939:   VL tvl,nvl,avl;
                    940:   int n,d,z;
                    941:   V v;
                    942:
                    943:   if ( NUM(p) ) {
                    944:     *mvlp = vl;
                    945:     *pr = p;
                    946:     return;
                    947:   }
                    948:   clctv(vl,p,&nvl);
                    949:   v = nvl->v;
                    950:   d = getchomdeg(v,p) + getlchomdeg(v,p,&z);
                    951:   for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
                    952:     n = getchomdeg(avl->v,p) + getlchomdeg(avl->v,p,&z);
                    953:     if ( n < d ) {
                    954:       v = avl->v; d = n;
                    955:     }
                    956:   }
                    957:   if ( v != nvl->v ) {
                    958:     reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
                    959:     *pr = t; *mvlp = tvl;
                    960:   } else {
                    961:     *pr = p; *mvlp = nvl;
                    962:   }
                    963: }
                    964:
                    965: int getchomdeg(V v,P p)
                    966: {
                    967:   int m,m1;
                    968:   DCP dc;
                    969:
                    970:   if ( !p || NUM(p) )
                    971:     return ( 0 );
                    972:   else if ( VR(p) == v )
                    973:     return ( dbound(v,p) );
                    974:   else {
                    975:     for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
1.2     ! noro      976:       m1 = getchomdeg(v,COEF(dc))+ZTOS(DEG(dc));
1.1       noro      977:       m = MAX(m,m1);
                    978:     }
                    979:     return ( m );
                    980:   }
                    981: }
                    982:
                    983: int getlchomdeg(V v,P p,int *d)
                    984: {
                    985:   int m0,m1,d0,d1;
                    986:   DCP dc;
                    987:
                    988:   if ( !p || NUM(p) ) {
                    989:     *d = 0;
                    990:     return ( 0 );
                    991:   } else if ( VR(p) == v ) {
1.2     ! noro      992:     *d = ZTOS(DEG(DC(p)));
1.1       noro      993:     return ( homdeg(LC(p)) );
                    994:   } else {
                    995:     for ( dc = DC(p), m0 = 0, d0 = 0; dc; dc = NEXT(dc) ) {
1.2     ! noro      996:       m1 = getlchomdeg(v,COEF(dc),&d1)+ZTOS(DEG(dc));
1.1       noro      997:       if ( d1 > d0 ) {
                    998:         m0 = m1;
                    999:         d0 = d1;
                   1000:       } else if ( d1 == d0 )
                   1001:         m0 = MAX(m0,m1);
                   1002:     }
                   1003:     *d = d0;
                   1004:     return ( m0 );
                   1005:   }
                   1006: }
                   1007:
                   1008: int nmonop(P p)
                   1009: {
                   1010:   int s;
                   1011:   DCP dc;
                   1012:
                   1013:   if ( !p )
                   1014:     return 0;
                   1015:   else
                   1016:     switch ( OID(p) ) {
                   1017:       case O_N:
                   1018:         return 1; break;
                   1019:       case O_P:
                   1020:         for ( dc = DC((P)p), s = 0; dc; dc = NEXT(dc) )
                   1021:           s += nmonop(COEF(dc));
                   1022:         return s; break;
                   1023:       default:
                   1024:         return 0;
                   1025:     }
                   1026: }
                   1027:
                   1028: int qpcheck(Obj p)
                   1029: {
                   1030:   DCP dc;
                   1031:
                   1032:   if ( !p )
                   1033:     return 1;
                   1034:   else
                   1035:     switch ( OID(p) ) {
                   1036:       case O_N:
                   1037:         return RATN((Num)p)?1:0;
                   1038:       case O_P:
                   1039:         for ( dc = DC((P)p); dc; dc = NEXT(dc) )
                   1040:           if ( !qpcheck((Obj)COEF(dc)) )
                   1041:             return 0;
                   1042:         return 1;
                   1043:       default:
                   1044:         return 0;
                   1045:     }
                   1046: }
                   1047:
                   1048: /* check if p is univariate and all coeffs are INT or LM */
                   1049:
                   1050: int uzpcheck(Obj p)
                   1051: {
                   1052:   DCP dc;
                   1053:   P c;
                   1054:
                   1055:   if ( !p )
                   1056:     return 1;
                   1057:   else
                   1058:     switch ( OID(p) ) {
                   1059:       case O_N:
                   1060:         return (RATN((Num)p)&&INT(p))||(NID((Num)p)==N_LM);
                   1061:       case O_P:
                   1062:         for ( dc = DC((P)p); dc; dc = NEXT(dc) ) {
                   1063:           c = COEF(dc);
                   1064:           if ( !NUM(c) || !uzpcheck((Obj)c) )
                   1065:             return 0;
                   1066:         }
                   1067:         return 1;
                   1068:       default:
                   1069:         return 0;
                   1070:     }
                   1071: }
                   1072:
                   1073: int p_mag(P p)
                   1074: {
                   1075:   int s;
                   1076:   DCP dc;
                   1077:
                   1078:   if ( !p )
                   1079:     return 0;
                   1080:   else if ( OID(p) == O_N )
                   1081:     return z_bits((Q)p);
                   1082:   else {
                   1083:     for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) )
                   1084:       s += p_mag(COEF(dc));
                   1085:     return s;
                   1086:   }
                   1087: }
                   1088:
                   1089: int maxblenp(P p)
                   1090: {
                   1091:   int s,t;
                   1092:   DCP dc;
                   1093:
                   1094:   if ( !p )
                   1095:     return 0;
                   1096:   else if ( OID(p) == O_N )
                   1097:     return z_bits((Q)p);
                   1098:   else {
                   1099:     for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
                   1100:       t = maxblenp(COEF(dc));
                   1101:       s = MAX(t,s);
                   1102:     }
                   1103:     return s;
                   1104:   }
                   1105: }

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