Return to polydiv.c CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / engine

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

1.1       noro        1: #ifdef MODULAR
                      2: void divsrmp(VL vl,int mod,P p1,P p2,P *q,P *r)
                      3: #else
                      4: void divsrp(VL vl,P p1,P p2,P *q,P *r)
                      5: #endif
                      6: {
                      7:   register int i,j;
                      8:   register DCP dc1,dc2,dc;
                      9:   P m,s,dvr;
                     10:   P *pq,*pr,*pd;
                     11:   V v1,v2;
                     12:   Z deg1,deg2;
                     13:   int d1,d2,sgn;
                     14:
                     15:   if ( !p1 ) {
                     16:     *q = 0; *r = 0;
                     17:   } else if ( NUM(p2) )
                     18:     if ( NUM(p1) ) {
                     19:       DIVNUM(p1,p2,q); *r = 0;
                     20:     } else {
                     21:       DIVSDCP(vl,p1,p2,q); *r = 0;
                     22:     }
                     23:   else if ( NUM(p1) ) {
                     24:     *q = 0; *r = p1;
                     25:   } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
                     26:     dc1 = DC(p1); dc2 = DC(p2);
                     27:     deg1 = DEG(dc1); deg2 = DEG(dc2);
                     28:     sgn = cmpz(deg1,deg2);
                     29:     if ( sgn == 0 ) {
                     30:       DIVSP(vl,COEF(dc1),COEF(dc2),q);
                     31:       MULP(vl,p2,*q,&m); SUBP(vl,p1,m,r);
                     32:     } else if ( sgn < 0 ) {
                     33:       *q = 0; *r = p1;
                     34:     } else {
                     35:       if ( !smallz(deg1)  )
                     36:         error("divsrp : invalid input");
1.2     ! noro       37:       d1 = ZTOS(deg1); d2 = ZTOS(deg2);
1.1       noro       38:       W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
                     39:       for ( dc = dc1; dc; dc = NEXT(dc) )
1.2     ! noro       40:         pr[ZTOS(DEG(dc))] = COEF(dc);
1.1       noro       41:       for ( dc = dc2; dc; dc = NEXT(dc) )
1.2     ! noro       42:         pd[ZTOS(DEG(dc))] = COEF(dc);
1.1       noro       43:       for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- ) {
                     44:         if ( !pr[i+d2] )
                     45:           continue;
                     46:         DIVSP(vl,pr[i+d2],dvr,&pq[i]);
                     47:         for ( j = d2; j >= 0; j-- ) {
                     48:           MULP(vl,pq[i],pd[j],&m);
                     49:           SUBP(vl,pr[i + j],m,&s); pr[i + j] = s;
                     50:         }
                     51:       }
                     52:       plisttop(pq,v1,d1 - d2,q); plisttop(pr,v1,d1 - 1,r);
                     53:     }
                     54:   } else {
                     55:     for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
                     56:     if ( v2 == vl->v ) {
                     57:       *q = 0; *r = p1;
                     58:     } else
                     59:       DIVSRDCP(vl,p1,p2,q,r);
                     60:   }
                     61: }
                     62:
                     63: #ifdef MODULAR
                     64: void divsrdcmp(VL vl,int mod,P p1,P p2,P *q,P *r)
                     65: #else
                     66: void divsrdcp(VL vl,P p1,P p2,P *q,P *r)
                     67: #endif
                     68: {
                     69:
                     70:   P qc,rc;
                     71:   DCP dc,dcq,dcq0,dcr,dcr0;
                     72:
                     73:   for ( dc = DC(p1), dcq0 = 0, dcr0 = 0; dc; dc = NEXT(dc) ) {
                     74:     DIVSRP(vl,COEF(dc),p2,&qc,&rc);
                     75:     if ( qc ) {
                     76:       NEXTDC(dcq0,dcq); DEG(dcq) = DEG(dc); COEF(dcq) = qc;
                     77:     }
                     78:     if ( rc ) {
                     79:       NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = rc;
                     80:     }
                     81:   }
                     82:   if ( dcq0 ) {
                     83:     NEXT(dcq) = 0; MKP(VR(p1),dcq0,*q);
                     84:   } else
                     85:     *q = 0;
                     86:   if ( dcr0 ) {
                     87:     NEXT(dcr) = 0; MKP(VR(p1),dcr0,*r);
                     88:   } else
                     89:     *r = 0;
                     90: }
                     91:
                     92: #ifdef MODULAR
                     93: void divsmp(VL vl,int mod,P p1,P p2,P *q)
                     94: #else
                     95: void divsp(VL vl,P p1,P p2,P *q)
                     96: #endif
                     97: {
                     98:   P t;
                     99:
                    100:   DIVSRP(vl,p1,p2,q,&t);
                    101:   if ( t )
                    102:     error("divsp: cannot happen");
                    103: }
                    104:
                    105: #ifdef MODULAR
                    106: void divsdcmp(VL vl,int mod,P p1,P p2,P *q)
                    107: #else
                    108: void divsdcp(VL vl,P p1,P p2,P *q)
                    109: #endif
                    110: {
                    111:
                    112:   P m;
                    113:   register DCP dc,dcr,dcr0;
                    114:
                    115:   for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) ) {
                    116:     DIVSP(vl,COEF(dc),p2,&m);
                    117:     NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
                    118:   }
                    119:   MKP(VR(p1),dcr0,*q);
                    120: }
                    121: