Return to A.c CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / engine

Annotation of OpenXM_contrib2/asir2018/engine/A.c, Revision 1.1

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:  *
        !            48:  * $OpenXM$
        !            49: */
        !            50: #include "ca.h"
        !            51: #include "parse.h"
        !            52:
        !            53: int get_lprime();
        !            54:
        !            55: void pdiva(VL vl,P p0,P p1,P p2,P *pr)
        !            56: {
        !            57:   Z m,d;
        !            58:   P t,q,r;
        !            59:
        !            60:   if ( VR(p2) == VR(p0) ) {
        !            61:     *pr = p1;
        !            62:     return;
        !            63:   }
        !            64:   if ( VR(p1) == VR(p0) ) {
        !            65:     error("pdiva : ???");
        !            66:     return;
        !            67:   }
        !            68:   subz(DEG(DC(p1)),DEG(DC(p2)),&m); addz(m,ONE,&d);
        !            69:   pwrz((Z)COEF(DC(p2)),d,&m); mulp(vl,p1,(P)m,&t);
        !            70:   divsrp(vl,t,p2,&q,&r); remsdcp(vl,q,p0,pr);
        !            71: }
        !            72:
        !            73: void rema(VL vl,P p0,P p1,P p2,P *pr)
        !            74: {
        !            75:   P m,m1,m2,g,lc,x,t;
        !            76:   Z q;
        !            77:   V v;
        !            78:   int d1,d2;
        !            79:
        !            80:   if ( !p1 || !p2 || NUM(p2) || ( VR(p2) == VR(p0) ) )
        !            81:     *pr = 0;
        !            82:   else if ( NUM(p1) )
        !            83:     *pr = p1;
        !            84:   else if ( ( v = VR(p1) ) == VR(p0) )
        !            85:     remsdcp(vl,p1,p0,pr);
        !            86:   else if ( ( d1 = deg(v,p1) ) < ( d2 = deg(v,p2) ) )
        !            87:     *pr = p1;
        !            88:   else if ( d1 == d2 ) {
        !            89:     mulp(vl,p1,LC(p2),&m); mulp(vl,p2,LC(p1),&m1);
        !            90:     subp(vl,m,m1,&m2);
        !            91:     remsdcp(vl,m2,p0,pr);
        !            92:   } else {
        !            93:     g = p1; lc = LC(p2); MKV(v,x);
        !            94:     while ( ( d1 = deg(v,g) ) >= d2 ) {
        !            95:       mulp(vl,g,lc,&m); mulp(vl,p2,LC(g),&m1);
        !            96:       STOQ(d1-d2,q); pwrp(vl,x,q,&t);
        !            97:       mulp(vl,m1,t,&m2); subp(vl,m,m2,&m1);
        !            98:       remsdcp(vl,m1,p0,&g);
        !            99:     }
        !           100:     *pr = g;
        !           101:   }
        !           102: }
        !           103:
        !           104: void gcda(VL vl,P p0,P p1,P p2,P *pr)
        !           105: {
        !           106:   P g1,g2,r,t,c;
        !           107:
        !           108:   remsdcp(vl,p1,p0,&t); pcp(vl,t,&g1,&c); remsdcp(vl,p2,p0,&g2);
        !           109:   while ( 1 ) {
        !           110:     if ( NUM(g2) || (VR(g2) == VR(p0)) ) {
        !           111:       *pr = (P)ONE;
        !           112:       return;
        !           113:     }
        !           114:     pcp(vl,g2,&t,&c); pmonic(vl,p0,t,&g2); rema(vl,p0,g1,g2,&r);
        !           115:     if ( !r ) {
        !           116:       pmonic(vl,p0,g2,pr);
        !           117:       return;
        !           118:     }
        !           119:     g1 = g2; g2 = r;
        !           120:   }
        !           121: }
        !           122:
        !           123: void sqa(VL vl,P p0,P p,DCP *dcp)
        !           124: {
        !           125:   V v;
        !           126:   int i,j,n;
        !           127:   P f,f1,f2,gcd,flat,q,r,t;
        !           128:   P *a;
        !           129:   int *b;
        !           130:   DCP dc,dc0;
        !           131:
        !           132:   if ( ( v = VR(p) ) == VR(p0) ) {
        !           133:     NEWDC(dc); *dcp = dc; COEF(dc) = (P)ONE; DEG(dc) = ONE;
        !           134:     NEXT(dc) = 0;
        !           135:   } else {
        !           136:     n = deg(v,p); W_CALLOC(n,P,a); W_CALLOC(n,int,b);
        !           137:     for ( i = 0, f = p; ;i++ ) {
        !           138:       if ( VR(f) != v )
        !           139:         break;
        !           140:       remsdcp(vl,f,p0,&f1); diffp(vl,f1,VR(f1),&f2);
        !           141:       while ( 1 ) {
        !           142:         pmonic(vl,p0,f2,&t); f2 = t; rema(vl,p0,f1,f2,&r);
        !           143:         if ( !r ) {
        !           144:           pmonic(vl,p0,f2,&gcd); pqra(vl,p0,f,gcd,&flat,&r);
        !           145:           break;
        !           146:         }
        !           147:         f1 = f2; f2 = r;
        !           148:       }
        !           149:       if ( VR(gcd) != v ) {
        !           150:         a[i] = f; b[i] = 1;
        !           151:         break;
        !           152:       }
        !           153:       for ( j = 1, f = gcd; ; j++ ) {
        !           154:         pqra(vl,p0,f,flat,&q,&r);
        !           155:         if ( r )
        !           156:           break;
        !           157:         else
        !           158:           f = q;
        !           159:       }
        !           160:       a[i] = flat; b[i] = j;
        !           161:     }
        !           162:     for ( i = 0, j = 0, dc0 = 0; a[i]; i++ ) {
        !           163:       NEXTDC(dc0,dc); j += b[i]; STOQ(j,DEG(dc));
        !           164:       if ( a[i+1] )
        !           165:         pqra(vl,p0,a[i],a[i+1],&COEF(dc),&t);
        !           166:       else
        !           167:         COEF(dc) = a[i];
        !           168:     }
        !           169:     NEXT(dc) = 0; *dcp = dc0;
        !           170:   }
        !           171: }
        !           172:
        !           173: void pqra(VL vl,P p0,P p1,P p2,P *p,P *pr)
        !           174: {
        !           175:   V v,v1;
        !           176:   P  m,m1,q,r;
        !           177:   Z tq;
        !           178:   int n;
        !           179:
        !           180:   if ( ( v = VR(p1) ) != ( v1 = VR(p2) ) ) {
        !           181:     while ( ( VR(vl) != v ) && ( VR(vl) != v1 ) )
        !           182:       vl = NEXT(vl);
        !           183:     if ( VR(vl) == v ) {
        !           184:       *p = p1; *pr = 0;
        !           185:     } else {
        !           186:       *p = 0; *pr = p1;
        !           187:     }
        !           188:   } else if ( (n = deg(v,p1) - deg(v,p2)) < 0 ) {
        !           189:     *pr = p1; *p = 0;
        !           190:   } else {
        !           191:     n++; STOQ(n,tq);
        !           192:     pwrp(vl,LC(p2),tq,&m); mulp(vl,m,p1,&m1); divsrp(vl,m1,p2,&q,&r);
        !           193:     remsdcp(vl,q,p0,p); remsdcp(vl,r,p0,pr);
        !           194:   }
        !           195: }
        !           196:
        !           197: void pmonic(VL vl, P p0, P p, P *pr)
        !           198: {
        !           199:   P d,m,pp,c,t;
        !           200:
        !           201:   if ( NUM(p) || ( VR(p) == VR(p0) ) )
        !           202:     *pr = (P)ONE;
        !           203:   else if ( NUM(COEF(DC(p))) )
        !           204:     *pr = p;
        !           205:   else {
        !           206:     ptozp(COEF(DC(p)),1,(Q *)&c,&pp); pinva(p0,pp,&d);
        !           207:     mulp(vl,p,d,&m); remsdcp(vl,m,p0,&t); ptozp(t,1,(Q *)&c,pr);
        !           208:   }
        !           209: }
        !           210:
        !           211: void remsdcp(VL vl, P p, P p0, P *pr)
        !           212: {
        !           213:   P q,r,c;
        !           214:   DCP dc,dcr,dcr0;
        !           215:
        !           216:   if ( !p )
        !           217:     *pr = 0;
        !           218:   else if ( NUM(p) )
        !           219:     *pr = (P)ONE;
        !           220:   else if ( VR(p) == VR(p0) ) {
        !           221:     divsrp(vl,p,p0,&q,&r); ptozp(r,1,(Q *)&c,pr);
        !           222:   } else {
        !           223:     for ( dcr = dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
        !           224:       divsrp(vl,COEF(dc),p0,&q,&r);
        !           225:       if ( r ) {
        !           226:         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = r;
        !           227:       }
        !           228:     }
        !           229:     if ( !dcr0 )
        !           230:       *pr = 0;
        !           231:     else {
        !           232:       NEXT(dcr) = 0; MKP(VR(p),dcr0,r); ptozp(r,1,(Q *)&c,pr);
        !           233:     }
        !           234:   }
        !           235: }
        !           236:
        !           237: void pinva(P p0, P p, P *pr)
        !           238: {
        !           239:   V v;
        !           240:   Z w,q,q1;
        !           241:   Q f;
        !           242:   Q t,a,b,s,ts;
        !           243:   P c,c1;
        !           244:   P tg,th,tr;
        !           245:   P z,zz,zzz;
        !           246:   UM wg,wh,ws,wt,wm;
        !           247:   int n,m,modres,mod,index,i;
        !           248:
        !           249:   v = VR(p); n=deg(v,p); m=deg(v,p0);
        !           250:   norm(p,&a); norm(p0,&b);
        !           251:   STOQ(m,w); pwrq(a,(Q)w,&t); STOQ(n,w); pwrq(b,(Q)w,&s);
        !           252:   mulq(t,s,&ts);
        !           253:   factorialz(n+m,&w); mulq(ts,(Q)w,&s); addq(s,s,&f);
        !           254:   wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
        !           255:   wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
        !           256:   wm = W_UMALLOC(m+n);
        !           257:   for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
        !           258:     mod = get_lprime(index++);
        !           259:     if ( !remqi((Q)LC(p),mod) || !remqi((Q)LC(p0),mod) )
        !           260:       continue;
        !           261:     ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr);
        !           262:     if ( !tr )
        !           263:       continue;
        !           264:     else
        !           265:       modres = CONT((MQ)tr);
        !           266:     mptoum(tg,wg); mptoum(th,wh);
        !           267:     eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
        !           268:     DEG(wm) = 0; COEF(wm)[0] = modres;
        !           269:     mulum(mod,ws,wm,wt);
        !           270:     for ( i = DEG(wt); i >= 0; i-- )
        !           271:       if ( ( COEF(wt)[i] * 2 ) > mod )
        !           272:         COEF(wt)[i] -= mod;
        !           273:     chnrem(mod,v,c,q,wt,&c1,&q1);
        !           274:     if ( !ucmpp(c,c1) ) {
        !           275:       mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);
        !           276:       if ( NUM(zzz) ) {
        !           277:         q = q1; c = c1; break;
        !           278:       }
        !           279:     }
        !           280:     q = q1; c = c1;
        !           281:     if ( cmpq(f,(Q)q) < 0 )
        !           282:       break;
        !           283:   }
        !           284:   ptozp(c,1,&s,pr);
        !           285: }
        !           286:
        !           287: /*
        !           288:  * compute g s.t. h*p0+g*p = res(v,p0,p)
        !           289:  */
        !           290:
        !           291: void inva_chrem(P p0, P p, NODE *pr)
        !           292: {
        !           293:   V v;
        !           294:   Z w,q,q1;
        !           295:   Q f,ts;
        !           296:   Q u,t,a,b,s,dn;
        !           297:   P c,c0,c1;
        !           298:   P tg,th,tr;
        !           299:   P z,zz,zzz;
        !           300:   UM wg,wh,ws,wt,wm;
        !           301:   int n,m,modres,mod,index,i;
        !           302:
        !           303:   v = VR(p); n=deg(v,p); m=deg(v,p0);
        !           304:   norm(p,&a); norm(p0,&b);
        !           305:   STOQ(m,w); pwrq(a,(Q)w,&t); STOQ(n,w); pwrq(b,(Q)w,&s);
        !           306:   mulq(t,s,&ts);
        !           307:   factorialz(n+m,&w); mulq(ts,(Q)w,&s); addq(s,s,&f);
        !           308:   wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
        !           309:   wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
        !           310:   wm = W_UMALLOC(m+n);
        !           311:   for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
        !           312:     mod = get_lprime(index++);
        !           313:     if ( !remqi((Q)LC(p),mod) || !remqi((Q)LC(p0),mod) )
        !           314:       continue;
        !           315:     ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr);
        !           316:     if ( !tr )
        !           317:       continue;
        !           318:     else
        !           319:       modres = CONT((MQ)tr);
        !           320:     mptoum(tg,wg); mptoum(th,wh);
        !           321:     eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
        !           322:     DEG(wm) = 0; COEF(wm)[0] = modres;
        !           323:     mulum(mod,ws,wm,wt);
        !           324:     for ( i = DEG(wt); i >= 0; i-- )
        !           325:       if ( ( COEF(wt)[i] * 2 ) > mod )
        !           326:         COEF(wt)[i] -= mod;
        !           327:     chnrem(mod,v,c,q,wt,&c1,&q1);
        !           328:     if ( !ucmpp(c,c1) ) {
        !           329:       mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);
        !           330:       if ( NUM(zzz) ) {
        !           331:         q = (Z)zzz; break;
        !           332:       }
        !           333:     }
        !           334:     q = q1; c = c1;
        !           335:     if ( cmpq(f,(Q)q) < 0 ) {
        !           336:       mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);
        !           337:       q = (Z)zzz;
        !           338:       break;
        !           339:     }
        !           340:   }
        !           341:   /* c*p = q mod p0 */
        !           342:   /* c = s*c0 */
        !           343:   ptozp(c,1,&s,&c0);
        !           344:   /* c/q = c0/(q/s) */
        !           345:   divq((Q)q,s,&dn);
        !           346:   *pr = (NODE)mknode(2,c0,dn);
        !           347: }