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

Annotation of OpenXM_contrib2/asir2000/engine/P.c, Revision 1.10

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.10    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/P.c,v 1.9 2004/03/11 09:52:56 noro Exp $
1.2       noro       49: */
1.1       noro       50: #ifndef FBASE
                     51: #define FBASE
                     52: #endif
                     53:
                     54: #include "b.h"
                     55: #include "ca.h"
                     56:
1.10    ! noro       57: #include "poly.c"
1.1       noro       58:
                     59: void ptozp(p,sgn,c,pr)
                     60: P p;
                     61: int sgn;
                     62: Q *c;
                     63: P *pr;
                     64: {
                     65:        N nm,dn;
                     66:
                     67:        if ( !p ) {
                     68:                *c = 0; *pr = 0;
                     69:        } else {
                     70:                lgp(p,&nm,&dn);
                     71:                if ( UNIN(dn) )
                     72:                        NTOQ(nm,sgn,*c);
                     73:                else
                     74:                        NDTOQ(nm,dn,sgn,*c);
                     75:                divcp(p,*c,pr);
                     76:        }
                     77: }
                     78:
                     79: void lgp(p,g,l)
                     80: P p;
                     81: N *g,*l;
                     82: {
                     83:        N g1,g2,l1,l2,l3,l4,tmp;
                     84:        DCP dc;
                     85:
                     86:        if ( NUM(p) ) {
                     87:                *g = NM((Q)p);
                     88:                if ( INT((Q)p) )
                     89:                        *l = ONEN;
                     90:                else
                     91:                        *l = DN((Q)p);
                     92:        } else {
                     93:                dc = DC(p); lgp(COEF(dc),g,l);
                     94:                for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                     95:                        lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
                     96:                        gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
                     97:                }
                     98:        }
                     99: }
                    100:
                    101: void divcp(f,q,rp)
                    102: P f;
                    103: Q q;
                    104: P *rp;
                    105: {
                    106:        DCP dc,dcr,dcr0;
                    107:
                    108:        if ( !f )
                    109:                *rp = 0;
                    110:        else if ( NUM(f) )
                    111:                DIVNUM(f,q,rp);
                    112:        else {
                    113:                for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
                    114:                        NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
                    115:                        divcp(COEF(dc),q,&COEF(dcr));
                    116:                }
                    117:                NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
                    118:        }
                    119: }
                    120:
                    121: void diffp(vl,p,v,r)
                    122: VL vl;
                    123: P p;
                    124: V v;
                    125: P *r;
                    126: {
                    127:        P t;
                    128:        DCP dc,dcr,dcr0;
                    129:
                    130:        if ( !p || NUM(p) )
                    131:                *r = 0;
                    132:        else {
                    133:                if ( v == VR(p) ) {
                    134:                        for ( dc = DC(p), dcr0 = 0;
                    135:                                dc && DEG(dc); dc = NEXT(dc) ) {
                    136:                                MULPQ(COEF(dc),(P)DEG(dc),&t);
                    137:                                if ( t ) {
                    138:                                        NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
                    139:                                        COEF(dcr) = t;
                    140:                                }
                    141:                        }
                    142:                        if ( !dcr0 )
                    143:                                *r = 0;
                    144:                        else {
                    145:                                NEXT(dcr) = 0; MKP(v,dcr0,*r);
                    146:                        }
                    147:                } else {
                    148:                        for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
                    149:                                diffp(vl,COEF(dc),v,&t);
1.8       ohara     150:                                if ( t ) {
                    151:                                        NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
                    152:                                }
                    153:                        }
                    154:                        if ( !dcr0 )
                    155:                                *r = 0;
                    156:                        else {
                    157:                                NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
                    158:                        }
                    159:                }
                    160:        }
                    161: }
                    162:
                    163: /* Euler derivation */
                    164: void ediffp(vl,p,v,r)
                    165: VL vl;
                    166: P p;
                    167: V v;
                    168: P *r;
                    169: {
                    170:        P t;
                    171:        DCP dc,dcr,dcr0;
                    172:
                    173:        if ( !p || NUM(p) )
                    174:                *r = 0;
                    175:        else {
                    176:                if ( v == VR(p) ) {
                    177:                        for ( dc = DC(p), dcr0 = 0;
                    178:                                  dc && DEG(dc); dc = NEXT(dc) ) {
                    179:                                MULPQ(COEF(dc),(P)DEG(dc),&t);
                    180:                                if ( t ) {
                    181:                                        NEXTDC(dcr0,dcr);
                    182:                                        DEG(dcr) = DEG(dc);
                    183:                                        COEF(dcr) = t;
                    184:                                }
                    185:                        }
                    186:                        if ( !dcr0 )
                    187:                                *r = 0;
                    188:                        else {
                    189:                                NEXT(dcr) = 0; MKP(v,dcr0,*r);
                    190:                        }
                    191:                } else {
                    192:                        for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
                    193:                                ediffp(vl,COEF(dc),v,&t);
1.1       noro      194:                                if ( t ) {
                    195:                                        NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
                    196:                                }
                    197:                        }
                    198:                        if ( !dcr0 )
                    199:                                *r = 0;
                    200:                        else {
                    201:                                NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
                    202:                        }
                    203:                }
                    204:        }
                    205: }
                    206:
                    207: void coefp(p,d,pr)
                    208: P p;
                    209: int d;
                    210: P *pr;
                    211: {
                    212:        DCP dc;
                    213:        int sgn;
                    214:        Q dq;
                    215:
                    216:        if ( NUM(p) )
                    217:                if ( d == 0 )
                    218:                        *pr = p;
                    219:                else
                    220:                        *pr = 0;
                    221:        else {
                    222:                for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
                    223:                        if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
                    224:                                continue;
                    225:                        else if ( sgn == 0 ) {
                    226:                                *pr = COEF(dc);
                    227:                                return;
                    228:                        } else {
                    229:                                *pr = 0;
                    230:                                break;
                    231:                        }
                    232:                *pr = 0;
                    233:        }
                    234: }
                    235:
                    236: int compp(vl,p1,p2)
                    237: VL vl;
                    238: P p1,p2;
                    239: {
                    240:        DCP dc1,dc2;
                    241:        V v1,v2;
                    242:
                    243:        if ( !p1 )
                    244:                return p2 ? -1 : 0;
                    245:        else if ( !p2 )
                    246:                return 1;
                    247:        else if ( NUM(p1) )
                    248:                return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
                    249:        else if ( NUM(p2) )
                    250:                return 1;
                    251:        if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
                    252:                for ( dc1 = DC(p1), dc2 = DC(p2);
                    253:                        dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
                    254:                        switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
                    255:                                case 1:
                    256:                                        return 1; break;
                    257:                                case -1:
                    258:                                        return -1; break;
                    259:                                default:
                    260:                                        switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
                    261:                                                case 1:
                    262:                                                        return 1; break;
                    263:                                                case -1:
                    264:                                                        return -1; break;
                    265:                                                default:
                    266:                                                        break;
                    267:                                        }
                    268:                                        break;
                    269:                        }
                    270:                return dc1 ? 1 : (dc2 ? -1 : 0);
                    271:        } else {
                    272:                for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
                    273:                return v1 == VR(vl) ? 1 : -1;
                    274:        }
1.9       noro      275: }
                    276:
                    277: int equalp(vl,p1,p2)
                    278: VL vl;
                    279: P p1,p2;
                    280: {
                    281:        DCP dc1,dc2;
                    282:        V v1,v2;
                    283:
                    284:        if ( !p1 ) {
                    285:                if ( !p2 ) return 1;
                    286:                else return 0;
                    287:        }
                    288:        /* p1 != 0 */
                    289:        if ( !p2 ) return 0;
                    290:
                    291:        /* p1 != 0, p2 != 0 */
                    292:        if ( NUM(p1) ) {
                    293:                if ( !NUM(p2) ) return 0;
                    294:                /* p1 and p2 are numbers */
                    295:                if ( NID((Num)p1) != NID((Num)p2) ) return 0;
                    296:                if ( !(*cmpnumt[NID(p1),NID(p1)])(p1,p2) ) return 1;
                    297:                return 0;
                    298:        }
                    299:        if ( VR(p1) != VR(p2) ) return 0;
                    300:        for ( dc1 = DC(p1), dc2 = DC(p2);
                    301:                dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) ) {
                    302:                if ( cmpq(DEG(dc1),DEG(dc2)) ) return 0;
                    303:                else if ( !equalp(vl,COEF(dc1),COEF(dc2)) ) return 0;
                    304:        }
                    305:        if ( dc1 || dc2 ) return 0;
                    306:        else return 1;
1.1       noro      307: }
                    308:
                    309: void csump(vl,p,c)
                    310: VL vl;
                    311: P p;
                    312: Q *c;
                    313: {
                    314:        DCP dc;
                    315:        Q s,s1,s2;
                    316:
                    317:        if ( !p || NUM(p) )
                    318:                *c = (Q)p;
                    319:        else {
                    320:                for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
                    321:                        csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
                    322:                }
                    323:                *c = s;
                    324:        }
                    325: }
                    326:
                    327: void degp(v,p,d)
                    328: V v;
                    329: P p;
                    330: Q *d;
                    331: {
                    332:        DCP dc;
                    333:        Q m,m1;
                    334:
                    335:        if ( !p || NUM(p) )
                    336:                *d = 0;
                    337:        else if ( v == VR(p) )
                    338:                *d = DEG(DC(p));
                    339:        else {
                    340:                for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
                    341:                        degp(v,COEF(dc),&m1);
                    342:                        if ( cmpq(m,m1) < 0 )
                    343:                                m = m1;
                    344:                }
                    345:                *d = m;
1.6       noro      346:        }
                    347: }
                    348:
                    349: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr);
                    350: void mulpq_trunc(P p,Q q,VN vn,P *pr);
                    351: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr);
                    352:
                    353: void mulp_trunc(VL vl,P p1,P p2,VN vn,P *pr)
                    354: {
                    355:        register DCP dc,dct,dcr,dcr0;
                    356:        V v1,v2;
                    357:        P t,s,u;
                    358:        int n1,n2,i,d,d1;
                    359:
                    360:        if ( !p1 || !p2 ) *pr = 0;
                    361:        else if ( NUM(p1) )
                    362:                mulpq_trunc(p2,(Q)p1,vn,pr);
                    363:        else if ( NUM(p2) )
                    364:                mulpq_trunc(p1,(Q)p2,vn,pr);
                    365:        else if ( ( v1 = VR(p1) ) ==  ( v2 = VR(p2) ) ) {
                    366:                for ( ; vn->v && vn->v != v1; vn++ )
                    367:                        if ( vn->n ) {
                    368:                                /* p1,p2 do not contain vn->v */
                    369:                                *pr = 0;
                    370:                                return;
                    371:                        }
                    372:                if ( !vn->v )
                    373:                        error("mulp_trunc : invalid vn");
                    374:                d = vn->n;
                    375:                for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
                    376:                        for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
                    377:                                d1 = QTOS(DEG(dct))+QTOS(DEG(dc));
                    378:                                if ( d1 >= d ) {
                    379:                                        mulp_trunc(vl,COEF(dct),COEF(dc),vn+1,&t);
                    380:                                        if ( t ) {
                    381:                                                NEXTDC(dcr0,dcr);
                    382:                                                STOQ(d1,DEG(dcr));
                    383:                                                COEF(dcr) = t;
                    384:                                        }
                    385:                                }
                    386:                        }
                    387:                        if ( dcr0 ) {
                    388:                                NEXT(dcr) = 0; MKP(v1,dcr0,t);
                    389:                                addp(vl,s,t,&u); s = u; t = u = 0;
                    390:                        }
                    391:                }
                    392:                *pr = s;
                    393:        } else {
                    394:                while ( v1 != VR(vl) && v2 != VR(vl) )
                    395:                        vl = NEXT(vl);
                    396:                if ( v1 == VR(vl) )
                    397:                        mulpc_trunc(vl,p1,p2,vn,pr);
                    398:                else
                    399:                        mulpc_trunc(vl,p2,p1,vn,pr);
                    400:        }
                    401: }
                    402:
                    403: void mulpq_trunc(P p,Q q,VN vn,P *pr)
                    404: {
                    405:        DCP dc,dcr,dcr0;
                    406:        P t;
                    407:        int i,d;
                    408:        V v;
                    409:
                    410:        if (!p || !q)
                    411:                *pr = 0;
                    412:        else if ( NUM(p) ) {
                    413:                for ( ; vn->v; vn++ )
                    414:                        if ( vn->n ) {
                    415:                                *pr = 0;
                    416:                                return;
                    417:                        }
                    418:                MULNUM(p,q,pr);
                    419:        } else {
                    420:                v = VR(p);
                    421:                for ( ; vn->v && vn->v != v; vn++ ) {
                    422:                        if ( vn->n ) {
                    423:                                /* p does not contain vn->v */
                    424:                                *pr = 0;
                    425:                                return;
                    426:                        }
                    427:                }
                    428:                if ( !vn->v )
                    429:                        error("mulpq_trunc : invalid vn");
                    430:                d = vn->n;
                    431:                for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
                    432:                        mulpq_trunc(COEF(dc),q,vn+1,&t);
                    433:                        if ( t ) {
                    434:                                NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
                    435:                        }
                    436:                }
                    437:                if ( dcr0 ) {
                    438:                        NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
                    439:                } else
                    440:                        *pr = 0;
                    441:        }
                    442: }
                    443:
                    444: void mulpc_trunc(VL vl,P p,P c,VN vn,P *pr)
                    445: {
                    446:        DCP dc,dcr,dcr0;
                    447:        P t;
                    448:        V v;
                    449:        int i,d;
                    450:
                    451:        if ( NUM(c) )
                    452:                mulpq_trunc(p,(Q)c,vn,pr);
                    453:        else {
                    454:                v = VR(p);
                    455:                for ( ; vn->v && vn->v != v; vn++ )
                    456:                        if ( vn->n ) {
                    457:                                /* p,c do not contain vn->v */
                    458:                                *pr = 0;
                    459:                                return;
                    460:                        }
                    461:                if ( !vn->v )
                    462:                        error("mulpc_trunc : invalid vn");
                    463:                d = vn->n;
                    464:                for ( dcr0 = 0, dc = DC(p); dc && QTOS(DEG(dc)) >= d; dc = NEXT(dc) ) {
                    465:                        mulp_trunc(vl,COEF(dc),c,vn+1,&t);
                    466:                        if ( t ) {
                    467:                                NEXTDC(dcr0,dcr); COEF(dcr) = t; DEG(dcr) = DEG(dc);
                    468:                        }
                    469:                }
                    470:                if ( dcr0 ) {
                    471:                        NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
                    472:                } else
                    473:                        *pr = 0;
1.7       noro      474:        }
                    475: }
                    476:
                    477: void quop_trunc(VL vl,P p1,P p2,VN vn,P *pr)
                    478: {
                    479:        DCP dc,dcq0,dcq;
                    480:        P t,s,m,lc2,qt;
                    481:        V v1,v2;
                    482:        Q d2;
                    483:        VN vn1;
                    484:
                    485:        if ( !p1 )
                    486:                *pr = 0;
                    487:        else if ( NUM(p2) )
                    488:                divsp(vl,p1,p2,pr);
                    489:        else if ( (v1 = VR(p1)) != (v2 = VR(p2)) ) {
                    490:                for ( dcq0 = 0, dc = DC(p1); dc; dc = NEXT(dc) ) {
                    491:                        NEXTDC(dcq0,dcq);
                    492:                        DEG(dcq) = DEG(dc);
                    493:                        quop_trunc(vl,COEF(dc),p2,vn,&COEF(dcq));
                    494:                }
                    495:                NEXT(dcq) = 0;
                    496:                MKP(v1,dcq0,*pr);
                    497:        } else {
                    498:                d2 = DEG(DC(p2));
                    499:                lc2 = COEF(DC(p2));
                    500:                t = p1;
                    501:                dcq0 = 0;
                    502:                /* vn1 = degree list of LC(p2) */
                    503:                for ( vn1 = vn; vn1->v != v1; vn1++ );
                    504:                vn1++;
                    505:                while ( t ) {
                    506:                        dc = DC(t);
                    507:                        NEXTDC(dcq0,dcq);
                    508:                        subq(DEG(dc),d2,&DEG(dcq));
                    509:                        quop_trunc(vl,COEF(dc),lc2,vn1,&COEF(dcq));
                    510:                        NEXT(dcq) = 0;
                    511:                        MKP(v1,dcq,qt);
                    512:                        mulp_trunc(vl,p2,qt,vn,&m);
                    513:                        subp(vl,t,m,&s); t = s;
                    514:                }
                    515:                NEXT(dcq) = 0;
                    516:                MKP(v1,dcq0,*pr);
1.1       noro      517:        }
                    518: }

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