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

Annotation of OpenXM_contrib2/asir2000/engine/PDM.c, Revision 1.3

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

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