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

Annotation of OpenXM_contrib2/asir2018/engine/PD.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/PD.c,v 1.1 2018/09/19 05:45:07 noro Exp $
1.1       noro       49: */
                     50: #ifndef FBASE
                     51: #define FBASE
                     52: #endif
                     53:
                     54: #include "b.h"
                     55: #include "ca.h"
                     56:
                     57: #include "polydiv.c"
                     58:
                     59: void plisttop(P *f,V v,int n,P *gp)
                     60: {
                     61:   int i;
                     62:   DCP dc,dc0;
                     63:
                     64:   for ( i = n; (i >= 0) && !f[i]; i-- );
                     65:   if ( i < 0 )
                     66:     *gp = 0;
                     67:   else if ( i == 0 )
                     68:     *gp = f[0];
                     69:   else {
                     70:     for ( dc0 = 0; i >= 0; i-- ) {
                     71:       if ( !f[i] )
                     72:         continue;
                     73:       NEXTDC(dc0,dc);
                     74:       if ( i )
1.2     ! noro       75:         STOZ(i,DEG(dc));
1.1       noro       76:       else
                     77:         DEG(dc) = 0;
                     78:       COEF(dc) = f[i];
                     79:     }
                     80:     NEXT(dc) = 0; MKP(v,dc0,*gp);
                     81:   }
                     82: }
                     83:
                     84: /* for multivariate polynomials over fields */
                     85:
                     86: int divtp(VL vl,P p1,P p2,P *q)
                     87: {
                     88:   register int i,j;
                     89:   register DCP dc1,dc2,dc;
                     90:   P m,m1,s,dvr,t;
                     91:   P *pq,*pr,*pd;
                     92:   V v1,v2;
                     93:   Z deg1,deg2;
                     94:   int d1,d2,sgn;
                     95:
                     96:   if ( !p1 ) {
                     97:     *q = 0;
                     98:     return 1;
                     99:   } else if ( NUM(p2) ) {
                    100:     divsp(vl,p1,p2,q);
                    101:     return 1;
                    102:   } else if ( NUM(p1) ) {
                    103:     *q = 0;
                    104:     return 0;
                    105:   } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
                    106:     dc1 = DC(p1); dc2 = DC(p2);
                    107:     deg1 = DEG(dc1); deg2 = DEG(dc2);
                    108:     sgn = cmpz(deg1,deg2);
                    109:     if ( sgn == 0 )
                    110:       if ( !divtp(vl,COEF(dc1),COEF(dc2),&m) ) {
                    111:         *q = 0;
                    112:         return 0;
                    113:       } else {
                    114:         mulp(vl,p2,m,&m1); subp(vl,p1,m1,&s);
                    115:         if ( !s ) {
                    116:           *q = m;
                    117:           return 1;
                    118:         } else {
                    119:           *q = 0;
                    120:           return 0;
                    121:         }
                    122:       }
                    123:     else if ( sgn < 0 ) {
                    124:       *q = 0;
                    125:       return 0;
                    126:     } else {
                    127:       if ( !smallz(deg1) ) {
                    128:         error("divtp : invalid input");
                    129:         *q = 0;
                    130:         return ( 0 );
                    131:       }
1.2     ! noro      132:       d1 = ZTOS(deg1); d2 = ZTOS(deg2);
1.1       noro      133:       W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
                    134:       for ( dc = dc1; dc; dc = NEXT(dc) )
1.2     ! noro      135:         pr[ZTOS(DEG(dc))] = COEF(dc);
1.1       noro      136:       for ( dc = dc2; dc; dc = NEXT(dc) )
1.2     ! noro      137:         pd[ZTOS(DEG(dc))] = COEF(dc);
1.1       noro      138:       for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- )
                    139:         if ( !pr[i+d2] )
                    140:           continue;
                    141:         else if ( !divtp(vl,pr[i+d2],dvr,&m) ) {
                    142:           *q = 0;
                    143:           return 0;
                    144:         } else {
                    145:           pq[i] = m;
                    146:           for ( j = d2; j >= 0; j-- ) {
                    147:             mulp(vl,pq[i],pd[j],&m);
                    148:             subp(vl,pr[i + j],m,&s); pr[i + j] = s;
                    149:           }
                    150:         }
                    151:       plisttop(pq,v1,d1 - d2,&m); plisttop(pr,v1,d1 - 1,&t);
                    152:       if ( t ) {
                    153:         *q = 0;
                    154:         return 0;
                    155:       } else {
                    156:         *q = m;
                    157:         return 1;
                    158:       }
                    159:     }
                    160:   } else {
                    161:     for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
                    162:     if ( v2 == vl->v ) {
                    163:       *q = 0;
                    164:       return 0;
                    165:     } else
                    166:       return divtdcp(vl,p1,p2,q);
                    167:   }
                    168: }
                    169:
                    170: int divtdcp(VL vl,P p1,P p2,P *q)
                    171: {
                    172:
                    173:   P m;
                    174:   register DCP dc,dcr,dcr0;
                    175:
                    176:   for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) )
                    177:     if ( !divtp(vl,COEF(dc),p2,&m) ) {
                    178:       *q = 0;
                    179:       return 0;
                    180:     } else {
                    181:       NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
                    182:     }
                    183:   MKP(VR(p1),dcr0,*q);
                    184:   return 1;
                    185: }
                    186:
                    187: int divtpz(VL vl,P p1,P p2,P *q)
                    188: {
                    189:   register int i,j;
                    190:   register DCP dc1,dc2,dc;
                    191:   P m,m1,s,dvr,t;
                    192:   P *pq,*pr,*pd;
                    193:   V v1,v2;
                    194:   Z deg1,deg2;
                    195:   int d1,d2,sgn;
                    196:
                    197:   if ( !p1 ) {
                    198:     *q = 0;
                    199:     return ( 1 );
                    200:   } else if ( NUM(p2) )
                    201:     if ( NUM(p1) ) {
                    202:       divq((Q)p1,(Q)p2,(Q *)&s);
                    203:       if ( INT((Q)s) ) {
                    204:         *q = s;
                    205:         return ( 1 );
                    206:       } else {
                    207:         *q = 0;
                    208:         return ( 0 );
                    209:       }
                    210:     } else
                    211:       return ( divtdcpz(vl,p1,p2,q) );
                    212:   else if ( NUM(p1) ) {
                    213:     *q = 0;
                    214:     return ( 0 );
                    215:   } else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
                    216:     Q csum1,csum2;
                    217:
                    218:     csump(vl,p1,&csum1); csump(vl,p2,&csum2);
                    219:     if ( csum2 && !divtpz(vl,(P)csum1,(P)csum2,&t) ) {
                    220:       *q = 0;
                    221:       return 0;
                    222:     }
                    223:     dc1 = DC(p1); dc2 = DC(p2);
                    224:     deg1 = DEG(dc1); deg2 = DEG(dc2);
                    225:     sgn = cmpz(deg1,deg2);
                    226:     if ( sgn == 0 )
                    227:       if ( !divtpz(vl,COEF(dc1),COEF(dc2),&m) ) {
                    228:         *q = 0;
                    229:         return ( 0 );
                    230:       } else {
                    231:         mulp(vl,p2,m,&m1); subp(vl,p1,m1,&s);
                    232:         if ( !s ) {
                    233:           *q = m;
                    234:           return ( 1 );
                    235:         } else {
                    236:           *q = 0;
                    237:           return ( 0 );
                    238:         }
                    239:       }
                    240:     else if ( sgn < 0 ) {
                    241:       *q = 0;
                    242:       return ( 0 );
                    243:     } else {
                    244:       if ( !smallz(deg1) ) {
                    245:         error("divtpz : invalid input");
                    246:         *q = 0;
                    247:         return ( 0 );
                    248:       }
1.2     ! noro      249:       d1 = ZTOS(deg1); d2 = ZTOS(deg2);
1.1       noro      250:       W_CALLOC(d1-d2,P,pq); W_CALLOC(d1,P,pr); W_CALLOC(d2,P,pd);
                    251:       for ( dc = dc1; dc; dc = NEXT(dc) )
1.2     ! noro      252:         pr[ZTOS(DEG(dc))] = COEF(dc);
1.1       noro      253:       for ( dc = dc2; dc; dc = NEXT(dc) )
1.2     ! noro      254:         pd[ZTOS(DEG(dc))] = COEF(dc);
1.1       noro      255:       for ( dvr = COEF(dc2), i = d1 - d2; i >= 0; i-- )
                    256:         if ( !pr[i+d2] )
                    257:           continue;
                    258:         else if ( !divtpz(vl,pr[i+d2],dvr,&m) ) {
                    259:           *q = 0;
                    260:           return ( 0 );
                    261:         } else {
                    262:           pq[i] = m;
                    263:           for ( j = d2; j >= 0; j-- ) {
                    264:             mulp(vl,pq[i],pd[j],&m);
                    265:             subp(vl,pr[i + j],m,&s); pr[i + j] = s;
                    266:           }
                    267:         }
                    268:       plisttop(pq,v1,d1 - d2,&m); plisttop(pr,v1,d1 - 1,&t);
                    269:       if ( t ) {
                    270:         *q = 0;
                    271:         return ( 0 );
                    272:       } else {
                    273:         *q = m;
                    274:         return ( 1 );
                    275:       }
                    276:     }
                    277:   } else {
                    278:     for ( ; (v1 != vl->v) && (v2 != vl->v); vl = NEXT(vl) );
                    279:     if ( v2 == vl->v ) {
                    280:       *q = 0;
                    281:       return ( 0 );
                    282:     } else
                    283:       return ( divtdcpz(vl,p1,p2,q) ) ;
                    284:   }
                    285: }
                    286:
                    287: int divtdcpz(VL vl,P p1,P p2,P *q)
                    288: {
                    289:
                    290:   P m;
                    291:   register DCP dc,dcr,dcr0;
                    292:
                    293:   for ( dc = DC(p1), dcr0 = 0; dc; dc = NEXT(dc) )
                    294:     if ( !divtpz(vl,COEF(dc),p2,&m) ) {
                    295:       *q = 0;
                    296:       return ( 0 );
                    297:     } else {
                    298:       NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = m; NEXT(dcr) = 0;
                    299:     }
                    300:   MKP(VR(p1),dcr0,*q);
                    301:   return ( 1 );
                    302: }
                    303:
                    304: void udivpz(P f1,P f2,P *fqp,P *frp)
                    305: {
                    306:   int n1,n2,i,j,sgn;
                    307:   Z *pq,*pr,*pd,d,m,s;
                    308:   DCP dc;
                    309:   Z qn,rn,q;
                    310:
                    311:   if ( !f2 )
                    312:     error("udivpz: division by 0");
                    313:   else if ( !f1 ) {
                    314:     *fqp = *frp = 0;
                    315:     return;
                    316:   } else if ( NUM(f1) ) {
                    317:     if ( NUM(f2) ) {
                    318:       divqrz((Z)f1,(Z)f2,&qn,&rn);
                    319:       if ( rn ) {
                    320:         *fqp = *frp = 0;
                    321:       } else {
                    322:         *fqp = (P)qn; *frp = 0;
                    323:       }
                    324:       return;
                    325:     } else {
                    326:       *fqp = 0; *frp = f1;
                    327:       return;
                    328:     }
                    329:   } else if ( NUM(f2) ) {
                    330:     n1 = UDEG(f1); W_CALLOC(n1,Z,pq);
                    331:     for ( dc = DC(f1); dc; dc = NEXT(dc) ) {
                    332:       divqrz((Z)COEF(dc),(Z)f2,&qn,&rn);
                    333:       if ( rn ) {
                    334:         *fqp = *frp = 0;
                    335:         return;
                    336:       } else {
1.2     ! noro      337:         pq[ZTOS(DEG(dc))] = qn;
1.1       noro      338:       }
                    339:     }
                    340:     plisttop((P *)pq,VR(f1),n1,fqp);
                    341:     *frp = 0;
                    342:     return;
                    343:   }
                    344:   n1  = UDEG(f1); n2  = UDEG(f2);
                    345:   if ( n1 < n2 ) {
                    346:     *fqp = NULL; *frp = f1;
                    347:     return;
                    348:   }
                    349:   W_CALLOC(n1-n2,Z,pq); W_CALLOC(n1,Z,pr); W_CALLOC(n2,Z,pd);
                    350:   for ( dc = DC(f1); dc; dc = NEXT(dc) )
1.2     ! noro      351:     pr[ZTOS(DEG(dc))] = (Z)COEF(dc);
1.1       noro      352:   for ( dc = DC(f2); dc; dc = NEXT(dc) )
1.2     ! noro      353:     pd[ZTOS(DEG(dc))] = (Z)COEF(dc);
1.1       noro      354:   for ( d = (Z)UCOEF(f2), i = n1 - n2; i >= 0; i-- ) {
                    355:     if ( !pr[i+n2] )
                    356:       continue;
                    357:     divqrz(pr[i+n2],d,&qn,&rn);
                    358:     if ( rn ) {
                    359:       *fqp = *frp = 0;
                    360:       return;
                    361:     }
                    362:     pq[i] = qn;
                    363:     for ( j = n2; j >= 0; j-- ) {
                    364:       mulz(pq[i],pd[j],&m); subz(pr[i+j],m,&s); pr[i+j] = s;
                    365:     }
                    366:   }
                    367:   plisttop((P *)pq,VR(f1),n1-n2,fqp); plisttop((P *)pr,VR(f1),n2-1,frp);
                    368: }
                    369:
                    370: /* YYY */
                    371:
                    372: void udivpwm(Z mod,P p1,P p2,P *q,P *r)
                    373: {
                    374:   P s,t,u,tq,tr;
                    375:
                    376:   invz((Z)UCOEF(p2),mod,(Z *)&t); mulpq(p2,t,&s); cmp(mod,s,&u);
                    377:   udivpzwm(mod,p1,u,&tq,&tr);
                    378:   cmp(mod,tr,r); mulpq(tq,t,&s); cmp(mod,s,q);
                    379: }
                    380:
                    381: void udivpzwm(Z mod,P f1,P f2,P *fqp,P *frp)
                    382: {
                    383:   int n1,n2,i,j;
                    384:   Z *pq,*pr,*pd,d,m,s;
                    385:   DCP dc;
                    386:   Z qn,rn,q;
                    387:
                    388:   if ( !f2 )
                    389:     error("udivpz: division by 0");
                    390:   else if ( !f1 ) {
                    391:     *fqp = *frp = 0;
                    392:     return;
                    393:   } else if ( NUM(f1) ) {
                    394:     if ( NUM(f2) ) {
                    395:       divqrz((Z)f1,(Z)f2,&qn,&rn);
                    396:       if ( rn ) {
                    397:         *fqp = *frp = 0;
                    398:       } else {
                    399:         *fqp = (P)qn; *frp = 0;
                    400:       }
                    401:       return;
                    402:     } else {
                    403:       *fqp = 0; *frp = f1;
                    404:       return;
                    405:     }
                    406:   } else if ( NUM(f2) ) {
                    407:     n1 = UDEG(f1); W_CALLOC(n1,Z,pq);
                    408:     for ( dc = DC(f1); dc; dc = NEXT(dc) ) {
                    409:       divqrz((Z)COEF(dc),(Z)f2,&qn,&rn);
                    410:       if ( rn ) {
                    411:         *fqp = *frp = 0;
                    412:         return;
                    413:       } else {
1.2     ! noro      414:         pq[ZTOS(DEG(dc))] = qn;
1.1       noro      415:       }
                    416:     }
                    417:     plisttop((P *)pq,VR(f1),n1,fqp);
                    418:     *frp = 0;
                    419:     return;
                    420:   }
                    421:   n1  = UDEG(f1); n2  = UDEG(f2);
                    422:   if ( n1 < n2 ) {
                    423:     *fqp = NULL; *frp = f1;
                    424:     return;
                    425:   }
                    426:   W_CALLOC(n1-n2,Z,pq); W_CALLOC(n1,Z,pr); W_CALLOC(n2,Z,pd);
                    427:   for ( dc = DC(f1); dc; dc = NEXT(dc) )
1.2     ! noro      428:     pr[ZTOS(DEG(dc))] = (Z)COEF(dc);
1.1       noro      429:   for ( dc = DC(f2); dc; dc = NEXT(dc) )
1.2     ! noro      430:     pd[ZTOS(DEG(dc))] = (Z)COEF(dc);
1.1       noro      431:   for ( d = (Z)UCOEF(f2), i = n1 - n2; i >= 0; i-- ) {
                    432:     if ( !pr[i+n2] )
                    433:       continue;
                    434:     divqrz(pr[i+n2],d,&qn,&rn);
                    435:     if ( rn ) {
                    436:       *fqp = *frp = 0;
                    437:       return;
                    438:     } else
                    439:       pq[i] = qn;
                    440:     for ( j = n2; j >= 0; j-- ) {
                    441:       mulz(pq[i],pd[j],&m); remz(m,mod,&s);
                    442:       subz(pr[i+j],s,&m); remz(m,mod,&s); pr[i+j] = s;
                    443:     }
                    444:   }
                    445:   plisttop((P *)pq,VR(f1),n1-n2,fqp); plisttop((P *)pr,VR(f1),n2-1,frp);
                    446: }

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