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

Annotation of OpenXM_contrib2/asir2018/engine/up.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 <math.h>
        !            52:
        !            53: struct oEGT eg_chrem,eg_fore,eg_back;
        !            54:
        !            55: /*
        !            56: int up_kara_mag=15;
        !            57: int up_tkara_mag=15;
        !            58: */
        !            59: /*
        !            60: int up_kara_mag=30;
        !            61: int up_tkara_mag=20;
        !            62: */
        !            63: int up_kara_mag=20;
        !            64: int up_tkara_mag=15;
        !            65:
        !            66: #if defined(sparc)
        !            67: int up_fft_mag=50;
        !            68: #else
        !            69: int up_fft_mag=80;
        !            70: #endif
        !            71:
        !            72: int up_rem_mag=1;
        !            73: int debug_up;
        !            74: int up_lazy;
        !            75: extern int lm_lazy;
        !            76: extern int current_ff;
        !            77: extern int GC_dont_gc;
        !            78:
        !            79: void monicup(UP a,UP *b)
        !            80: {
        !            81:   UP w;
        !            82:
        !            83:   if ( !a )
        !            84:     *b = 0;
        !            85:   else {
        !            86:     w = W_UPALLOC(0); w->d = 0;
        !            87:     divnum(0,(Num)ONE,(Num)a->c[a->d],(Num *)&w->c[0]);
        !            88:     mulup(a,w,b);
        !            89:   }
        !            90: }
        !            91:
        !            92: void simpup(UP a,UP *b)
        !            93: {
        !            94:   int i,d;
        !            95:   UP c;
        !            96:
        !            97:   if ( !a )
        !            98:     *b = 0;
        !            99:   else {
        !           100:     d = a->d;
        !           101:     c = UPALLOC(d);
        !           102:
        !           103:     for ( i = 0; i <= d; i++ )
        !           104:       simpnum(a->c[i],&c->c[i]);
        !           105:     for ( i = d; i >= 0 && !c->c[i]; i-- );
        !           106:     if ( i < 0 )
        !           107:       *b = 0;
        !           108:     else {
        !           109:       c->d = i;
        !           110:       *b = c;
        !           111:     }
        !           112:   }
        !           113: }
        !           114:
        !           115: void simpnum(Num a,Num *b)
        !           116: {
        !           117:   LM lm;
        !           118:   GF2N gf;
        !           119:   GFPN gfpn;
        !           120:
        !           121:   if ( !a )
        !           122:     *b = 0;
        !           123:   else
        !           124:     switch ( NID(a) ) {
        !           125:       case N_LM:
        !           126:         simplm((LM)a,&lm); *b = (Num)lm; break;
        !           127:       case N_GF2N:
        !           128:         simpgf2n((GF2N)a,&gf); *b = (Num)gf; break;
        !           129:       case N_GFPN:
        !           130:         simpgfpn((GFPN)a,&gfpn); *b = (Num)gfpn; break;
        !           131:       default:
        !           132:         *b = a; break;
        !           133:     }
        !           134: }
        !           135:
        !           136: void uremp(P p1,P p2,P *rp)
        !           137: {
        !           138:   UP n1,n2,r;
        !           139:
        !           140:   if ( !p1 || NUM(p1) )
        !           141:     *rp = p1;
        !           142:   else {
        !           143:     ptoup(p1,&n1); ptoup(p2,&n2);
        !           144:     remup(n1,n2,&r);
        !           145:     uptop(r,rp);
        !           146:   }
        !           147: }
        !           148:
        !           149: void ugcdp(P p1,P p2,P *rp)
        !           150: {
        !           151:   UP n1,n2,r;
        !           152:
        !           153:   ptoup(p1,&n1); ptoup(p2,&n2);
        !           154:   gcdup(n1,n2,&r);
        !           155:   uptop(r,rp);
        !           156: }
        !           157:
        !           158: void reversep(P p1,Q d,P *rp)
        !           159: {
        !           160:   UP n1,r;
        !           161:
        !           162:   if ( !p1 )
        !           163:     *rp = 0;
        !           164:   else {
        !           165:     ptoup(p1,&n1);
        !           166:     reverseup(n1,QTOS(d),&r);
        !           167:     uptop(r,rp);
        !           168:   }
        !           169: }
        !           170:
        !           171: void invmodp(P p1,Q d,P *rp)
        !           172: {
        !           173:   UP n1,r;
        !           174:
        !           175:   if ( !p1 )
        !           176:     error("invmodp : invalid input");
        !           177:   else {
        !           178:     ptoup(p1,&n1);
        !           179:     invmodup(n1,QTOS(d),&r);
        !           180:     uptop(r,rp);
        !           181:   }
        !           182: }
        !           183:
        !           184: void addup(UP n1,UP n2,UP *nr)
        !           185: {
        !           186:   UP r,t;
        !           187:   int i,d1,d2;
        !           188:   Num *c,*c1,*c2;
        !           189:
        !           190:   if ( !n1 ) {
        !           191:     *nr = n2;
        !           192:   } else if ( !n2 ) {
        !           193:     *nr = n1;
        !           194:   } else {
        !           195:     if ( n2->d > n1->d ) {
        !           196:       t = n1; n1 = n2; n2 = t;
        !           197:     }
        !           198:     d1 = n1->d; d2 = n2->d;
        !           199:     *nr = r = UPALLOC(d1);
        !           200:     c1 = n1->c; c2 = n2->c; c = r->c;
        !           201:     for ( i = 0; i <= d2; i++ )
        !           202:       addnum(0,*c1++,*c2++,c++);
        !           203:     for ( ; i <= d1; i++ )
        !           204:       *c++ = *c1++;
        !           205:     for ( i = d1, c = r->c; i >=0 && !c[i]; i-- );
        !           206:     if ( i < 0 )
        !           207:       *nr = 0;
        !           208:     else
        !           209:       r->d = i;
        !           210:   }
        !           211: }
        !           212:
        !           213: void subup(UP n1,UP n2,UP *nr)
        !           214: {
        !           215:   UP r;
        !           216:   int i,d1,d2,d;
        !           217:   Num *c,*c1,*c2;
        !           218:
        !           219:   if ( !n1 ) {
        !           220:     chsgnup(n2,nr);
        !           221:   } else if ( !n2 ) {
        !           222:     *nr = n1;
        !           223:   } else {
        !           224:     d1 = n1->d; d2 = n2->d; d = MAX(d1,d2);
        !           225:     *nr = r = UPALLOC(d);
        !           226:     c1 = n1->c; c2 = n2->c; c = r->c;
        !           227:     if ( d1 >= d2 ) {
        !           228:       for ( i = 0; i <= d2; i++ )
        !           229:         subnum(0,*c1++,*c2++,c++);
        !           230:       for ( ; i <= d1; i++ )
        !           231:         *c++ = *c1++;
        !           232:     } else {
        !           233:       for ( i = 0; i <= d1; i++ )
        !           234:         subnum(0,*c1++,*c2++,c++);
        !           235:       for ( ; i <= d2; i++ )
        !           236:         chsgnnum(*c2++,c++);
        !           237:     }
        !           238:     for ( i = d, c = r->c; i >=0 && !c[i]; i-- );
        !           239:     if ( i < 0 )
        !           240:       *nr = 0;
        !           241:     else
        !           242:       r->d = i;
        !           243:   }
        !           244: }
        !           245:
        !           246: void chsgnup(UP n1,UP *nr)
        !           247: {
        !           248:   UP r;
        !           249:   int d1,i;
        !           250:   Num *c1,*c;
        !           251:
        !           252:   if ( !n1 ) {
        !           253:     *nr = 0;
        !           254:   } else {
        !           255:     d1 = n1->d;
        !           256:     *nr = r = UPALLOC(d1); r->d = d1;
        !           257:     c1 = n1->c; c = r->c;
        !           258:     for ( i = 0; i <= d1; i++ )
        !           259:       chsgnnum(*c1++,c++);
        !           260:   }
        !           261: }
        !           262:
        !           263: void hybrid_mulup(int ff,UP n1,UP n2,UP *nr)
        !           264: {
        !           265:   if ( !n1 || !n2 )
        !           266:     *nr = 0;
        !           267:   else if ( MAX(n1->d,n2->d) < up_fft_mag )
        !           268:     kmulup(n1,n2,nr);
        !           269:   else
        !           270:     switch ( ff ) {
        !           271:       case FF_GFP:
        !           272:         fft_mulup_lm(n1,n2,nr); break;
        !           273:       case FF_GF2N:
        !           274:         kmulup(n1,n2,nr); break;
        !           275:       default:
        !           276:         fft_mulup(n1,n2,nr); break;
        !           277:     }
        !           278: }
        !           279:
        !           280: void hybrid_squareup(int ff,UP n1,UP *nr)
        !           281: {
        !           282:   if ( !n1 )
        !           283:     *nr = 0;
        !           284:   else if ( n1->d < up_fft_mag )
        !           285:     ksquareup(n1,nr);
        !           286:   else
        !           287:     switch ( ff ) {
        !           288:       case FF_GFP:
        !           289:         fft_squareup_lm(n1,nr); break;
        !           290:       case FF_GF2N:
        !           291:         ksquareup(n1,nr); break;
        !           292:       default:
        !           293:         fft_squareup(n1,nr); break;
        !           294:     }
        !           295: }
        !           296:
        !           297: void hybrid_tmulup(int ff,UP n1,UP n2,int d,UP *nr)
        !           298: {
        !           299:   if ( !n1 || !n2 )
        !           300:     *nr = 0;
        !           301:   else if ( MAX(n1->d,n2->d) < up_fft_mag )
        !           302:     tkmulup(n1,n2,d,nr);
        !           303:   else
        !           304:     switch ( ff ) {
        !           305:       case FF_GFP:
        !           306:         trunc_fft_mulup_lm(n1,n2,d,nr); break;
        !           307:       case FF_GF2N:
        !           308:         tkmulup(n1,n2,d,nr); break;
        !           309:       default:
        !           310:         trunc_fft_mulup(n1,n2,d,nr); break;
        !           311:     }
        !           312: }
        !           313:
        !           314: void mulup(UP n1,UP n2,UP *nr)
        !           315: {
        !           316:   UP r;
        !           317:   Num *pc1,*pc,*c1,*c2,*c;
        !           318:   Num mul,t,s;
        !           319:   int i,j,d1,d2;
        !           320:
        !           321:   if ( !n1 || !n2 )
        !           322:     *nr = 0;
        !           323:   else {
        !           324:     d1 = n1->d; d2 = n2->d;
        !           325:     *nr = r = UPALLOC(d1+d2); r->d = d1+d2;
        !           326:     c1 = n1->c; c2 = n2->c; c = r->c;
        !           327:     lm_lazy = 1;
        !           328:     for ( i = 0; i <= d2; i++, c++ )
        !           329:       if ( (mul = *c2++) != 0 )
        !           330:         for ( j = 0, pc1 = c1, pc = c; j <= d1; j++ ) {
        !           331:           mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;
        !           332:       }
        !           333:     lm_lazy = 0;
        !           334:     if ( !up_lazy )
        !           335:       for ( i = 0, c = r->c; i <= r->d; i++ ) {
        !           336:         simpnum(c[i],&t); c[i] = t;
        !           337:       }
        !           338:   }
        !           339: }
        !           340:
        !           341: /* nr = c*n1 */
        !           342:
        !           343: void mulcup(Num c,UP n1,UP *nr)
        !           344: {
        !           345:   int d;
        !           346:   UP r;
        !           347:   Num t;
        !           348:   Num *c1,*cr;
        !           349:   int i;
        !           350:
        !           351:   if ( !c || !n1 )
        !           352:     *nr = 0;
        !           353:   else {
        !           354:     d = n1->d;
        !           355:     *nr = r = UPALLOC(d); r->d = d;
        !           356:     c1 = n1->c; cr = r->c;
        !           357:     lm_lazy = 1;
        !           358:     for ( i = 0; i <= d; i++ )
        !           359:       mulnum(0,c1[i],c,&cr[i]);
        !           360:     lm_lazy = 0;
        !           361:     if ( !up_lazy )
        !           362:       for ( i = 0; i <= d; i++ ) {
        !           363:         simpnum(cr[i],&t); cr[i] = t;
        !           364:       }
        !           365:   }
        !           366: }
        !           367:
        !           368: void tmulup(UP n1,UP n2,int d,UP *nr)
        !           369: {
        !           370:   UP r;
        !           371:   Num *pc1,*pc,*c1,*c2,*c;
        !           372:   Num mul,t,s;
        !           373:   int i,j,d1,d2,dr;
        !           374:
        !           375:   if ( !n1 || !n2 )
        !           376:     *nr = 0;
        !           377:   else {
        !           378:     d1 = n1->d; d2 = n2->d;
        !           379:     dr = MAX(d-1,d1+d2);
        !           380:     *nr = r = UPALLOC(dr);
        !           381:     c1 = n1->c; c2 = n2->c; c = r->c;
        !           382:     lm_lazy = 1;
        !           383:     for ( i = 0; i <= d2; i++, c++ )
        !           384:       if ( ( mul = *c2++ ) != 0 )
        !           385:         for ( j = 0, pc1 = c1, pc = c; (j <= d1) && (i+j < d); j++ ) {
        !           386:           mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;
        !           387:       }
        !           388:     lm_lazy = 0;
        !           389:     c = r->c;
        !           390:     if ( !up_lazy )
        !           391:       for ( i = 0; i < d; i++ ) {
        !           392:         simpnum(c[i],&t); c[i] = t;
        !           393:       }
        !           394:     for ( i = d-1; i >= 0 && !c[i]; i-- );
        !           395:     if ( i < 0 )
        !           396:       *nr = 0;
        !           397:     else
        !           398:       r->d = i;
        !           399:   }
        !           400: }
        !           401:
        !           402: void squareup(UP n1,UP *nr)
        !           403: {
        !           404:   UP r;
        !           405:   Num *c1,*c;
        !           406:   Num t,s,u;
        !           407:   int i,j,h,d1,d;
        !           408:
        !           409:   if ( !n1 )
        !           410:     *nr = 0;
        !           411:   else if ( !n1->d ) {
        !           412:     *nr = r = UPALLOC(0); r->d = 0;
        !           413:     mulnum(0,n1->c[0],n1->c[0],&r->c[0]);
        !           414:   } else {
        !           415:     d1 = n1->d;
        !           416:     d = 2*d1;
        !           417:     *nr = r = UPALLOC(2*d); r->d = d;
        !           418:     c1 = n1->c; c = r->c;
        !           419:     lm_lazy = 1;
        !           420:     for ( i = 0; i <= d1; i++ )
        !           421:       mulnum(0,c1[i],c1[i],&c[2*i]);
        !           422:     for ( j = 1; j < d; j++ ) {
        !           423:       h = (j-1)/2;
        !           424:       for ( i = MAX(j-d1,0), s = 0; i <= h; i++ ) {
        !           425:         mulnum(0,c1[i],c1[j-i],&t); addnum(0,t,s,&u); s = u;
        !           426:       }
        !           427:       addnum(0,s,s,&t); addnum(0,t,c[j],&u); c[j] = u;
        !           428:     }
        !           429:     lm_lazy = 0;
        !           430:     if ( !up_lazy )
        !           431:       for ( i = 0, c = r->c; i <= d; i++ ) {
        !           432:         simpnum(c[i],&t); c[i] = t;
        !           433:       }
        !           434:   }
        !           435: }
        !           436:
        !           437: void remup(UP n1,UP n2,UP *nr)
        !           438: {
        !           439:   UP w,r;
        !           440:
        !           441:   if ( !n2 )
        !           442:     error("remup : division by 0");
        !           443:   else if ( !n1 || !n2->d )
        !           444:     *nr = 0;
        !           445:   else if ( n1->d < n2->d )
        !           446:     *nr = n1;
        !           447:   else {
        !           448:     w = W_UPALLOC(n1->d); copyup(n1,w);
        !           449:     remup_destructive(w,n2);
        !           450:     if ( w->d < 0 )
        !           451:       *nr = 0;
        !           452:     else {
        !           453:       *nr = r = UPALLOC(w->d); copyup(w,r);
        !           454:     }
        !           455:   }
        !           456: }
        !           457:
        !           458: void remup_destructive(UP n1,UP n2)
        !           459: {
        !           460:   Num *c1,*c2;
        !           461:   int d1,d2,i,j;
        !           462:   Num m,t,s,mhc;
        !           463:
        !           464:   c1 = n1->c; c2 = n2->c;
        !           465:   d1 = n1->d; d2 = n2->d;
        !           466:   chsgnnum(c2[d2],&mhc);
        !           467:   for ( i = d1; i >= d2; i-- ) {
        !           468:     simpnum(c1[i],&t); c1[i] = t;
        !           469:     if ( !c1[i] )
        !           470:       continue;
        !           471:     divnum(0,c1[i],mhc,&m);
        !           472:     lm_lazy = 1;
        !           473:     for ( j = d2; j >= 0; j-- ) {
        !           474:       mulnum(0,c2[j],m,&t); addnum(0,c1[i+j-d2],t,&s);
        !           475:       c1[i+j-d2] = s;
        !           476:     }
        !           477:     lm_lazy = 0;
        !           478:   }
        !           479:   for ( i = 0; i < d2; i++ ) {
        !           480:     simpnum(c1[i],&t); c1[i] = t;
        !           481:   }
        !           482:   for ( i = d2-1; i >= 0 && !c1[i]; i-- );
        !           483:   n1->d = i;
        !           484: }
        !           485:
        !           486: void qrup(UP n1,UP n2,UP *nq,UP *nr)
        !           487: {
        !           488:   UP w,r,q;
        !           489:   struct oUP t;
        !           490:   Num m;
        !           491:   int d;
        !           492:
        !           493:   if ( !n2 )
        !           494:     error("qrup : division by 0");
        !           495:   else if ( !n1 ) {
        !           496:     *nq = 0; *nr = 0;
        !           497:   } else if ( !n2->d )
        !           498:     if ( !n1->d ) {
        !           499:       divnum(0,n1->c[0],n2->c[0],(Num *)nq); *nr = 0;
        !           500:     } else {
        !           501:       divnum(0,(Num)ONE,n2->c[0],&m);
        !           502:       t.d = 0; t.c[0] = m;
        !           503:       mulup(n1,&t,nq); *nr = 0;
        !           504:     }
        !           505:   else if ( n1->d < n2->d ) {
        !           506:     *nq = 0; *nr = n1;
        !           507:   } else {
        !           508:     w = W_UPALLOC(n1->d); copyup(n1,w);
        !           509:     qrup_destructive(w,n2);
        !           510:     d = n1->d-n2->d;
        !           511:     *nq = q = UPALLOC(d); q->d = d;
        !           512:     bcopy(w->c+n2->d,q->c,(d+1)*sizeof(Num));
        !           513:     if ( w->d < 0 )
        !           514:       *nr = 0;
        !           515:     else {
        !           516:       *nr = r = UPALLOC(w->d); copyup(w,r);
        !           517:     }
        !           518:   }
        !           519: }
        !           520:
        !           521: void qrup_destructive(UP n1,UP n2)
        !           522: {
        !           523:   Num *c1,*c2;
        !           524:   int d1,d2,i,j;
        !           525:   Num q,t,s,hc;
        !           526:
        !           527:   c1 = n1->c; c2 = n2->c;
        !           528:   d1 = n1->d; d2 = n2->d;
        !           529:   hc = c2[d2];
        !           530:   for ( i = d1; i >= d2; i-- ) {
        !           531:     simpnum(c1[i],&t); c1[i] = t;
        !           532:     if ( c1[i] ) {
        !           533:       divnum(0,c1[i],hc,&q);
        !           534:       lm_lazy = 1;
        !           535:       for ( j = d2; j >= 0; j-- ) {
        !           536:         mulnum(0,c2[j],q,&t); subnum(0,c1[i+j-d2],t,&s);
        !           537:         c1[i+j-d2] = s;
        !           538:       }
        !           539:       lm_lazy = 0;
        !           540:     } else
        !           541:       q = 0;
        !           542:     c1[i] = q;
        !           543:   }
        !           544:   for ( i = 0; i < d2; i++ ) {
        !           545:     simpnum(c1[i],&t); c1[i] = t;
        !           546:   }
        !           547:   for ( i = d2-1; i >= 0 && !c1[i]; i-- );
        !           548:   n1->d = i;
        !           549: }
        !           550:
        !           551: void gcdup(UP n1,UP n2,UP *nr)
        !           552: {
        !           553:   UP w1,w2,t,r;
        !           554:   int d1,d2;
        !           555:
        !           556:   if ( !n1 )
        !           557:     *nr = n2;
        !           558:   else if ( !n2 )
        !           559:     *nr = n1;
        !           560:   else if ( !n1->d || !n2->d ) {
        !           561:     *nr = r = UPALLOC(0); r->d = 0;
        !           562:     divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]);
        !           563:   } else {
        !           564:     d1 = n1->d; d2 = n2->d;
        !           565:     if ( d2 > d1 ) {
        !           566:       w1 = W_UPALLOC(d2); copyup(n2,w1);
        !           567:       w2 = W_UPALLOC(d1); copyup(n1,w2);
        !           568:     } else {
        !           569:       w1 = W_UPALLOC(d1); copyup(n1,w1);
        !           570:       w2 = W_UPALLOC(d2); copyup(n2,w2);
        !           571:     }
        !           572:     d1 = w1->d; d2 = w2->d;
        !           573:     while ( 1 ) {
        !           574:       remup_destructive(w1,w2);
        !           575:       if ( w1->d < 0 ) {
        !           576:         *nr = r = UPALLOC(w2->d); copyup(w2,r); return;
        !           577:       } else if ( !w1->d ) {
        !           578:         *nr = r = UPALLOC(0); r->d = 0;
        !           579:         divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]); return;
        !           580:       } else {
        !           581:         t = w1; w1 = w2; w2 = t;
        !           582:       }
        !           583:     }
        !           584:   }
        !           585: }
        !           586:
        !           587: /* compute r s.t. a*r = 1 mod m */
        !           588:
        !           589: void extended_gcdup(UP a,UP m,UP *r)
        !           590: {
        !           591:   UP one,g1,g2,a1,a2,a3,b1,b2,b3,inv,quo,rem,t;
        !           592:   Num i;
        !           593:
        !           594:   one = UPALLOC(0); one->d = 0; one->c[0] = (Num)ONELM;
        !           595:   g1 = m; g2 = a;
        !           596:   a1 = one; a2 = 0;
        !           597:   b1 = 0; b2 = one;
        !           598:   /* a2*m+b2*a = g2 always holds */
        !           599:   while ( g2->d ) {
        !           600:     qrup(g1,g2,&quo,&rem); g1 = g2; g2 = rem;
        !           601:     mulup(a2,quo,&t); subup(a1,t,&a3); a1 = a2; a2 = a3;
        !           602:     mulup(b2,quo,&t); subup(b1,t,&b3); b1 = b2; b2 = b3;
        !           603:   }
        !           604:   /* now a2*m+b2*a = g2 (const) */
        !           605:   divnum(0,(Num)ONE,g2->c[0],&i);
        !           606:   inv = UPALLOC(0); inv->d = 0; inv->c[0] = i;
        !           607:   mulup(b2,inv,r);
        !           608: }
        !           609:
        !           610: void reverseup(UP n1,int d,UP *nr)
        !           611: {
        !           612:   Num *c1,*c;
        !           613:   int i,d1;
        !           614:   UP r;
        !           615:
        !           616:   if ( !n1 )
        !           617:     *nr = 0;
        !           618:   else if ( d < (d1 = n1->d) )
        !           619:     error("reverseup : invalid input");
        !           620:   else {
        !           621:     *nr = r = UPALLOC(d);
        !           622:     c = r->c;
        !           623:     bzero((char *)c,(d+1)*sizeof(Num));
        !           624:     c1 = n1->c;
        !           625:     for ( i = 0; i <= d1; i++ )
        !           626:       c[d-i] = c1[i];
        !           627:     for ( i = d; i >= 0 && !c[i]; i-- );
        !           628:     r->d = i;
        !           629:     if ( i < 0 )
        !           630:       *nr = 0;
        !           631:   }
        !           632: }
        !           633:
        !           634: void invmodup(UP n1,int d,UP *nr)
        !           635: {
        !           636:   UP r;
        !           637:   Num s,t,u,hinv;
        !           638:   int i,j,dmin;
        !           639:   Num *w,*c,*cr;
        !           640:
        !           641:   if ( !n1 || !n1->c[0] )
        !           642:     error("invmodup : invalid input");
        !           643:   else if ( !n1->d ) {
        !           644:     *nr = r = UPALLOC(0); r->d = 0;
        !           645:     divnum(0,(Num)ONE,n1->c[0],&r->c[0]);
        !           646:   } else {
        !           647:     *nr = r = UPALLOC(d);
        !           648:     w = (Num *)ALLOCA((d+1)*sizeof(Q));
        !           649:     bzero((char *)w,(d+1)*sizeof(Q));
        !           650:     dmin = MIN(d,n1->d);
        !           651:     divnum(0,(Num)ONE,n1->c[0],&hinv);
        !           652:     for ( i = 0, c = n1->c; i <= dmin; i++ )
        !           653:       mulnum(0,c[i],hinv,&w[i]);
        !           654:     cr = r->c;
        !           655:     cr[0] = w[0];
        !           656:     for ( i = 1; i <= d; i++ ) {
        !           657:       for ( j = 1, s = w[i]; j < i; j++ ) {
        !           658:         mulnum(0,cr[j],w[i-j],&t);
        !           659:         addnum(0,s,t,&u);
        !           660:         s = u;
        !           661:       }
        !           662:       chsgnnum(s,&cr[i]);
        !           663:     }
        !           664:     for ( i = 0; i <= d; i++ ) {
        !           665:       mulnum(0,cr[i],hinv,&t); cr[i] = t;
        !           666:     }
        !           667:     for ( i = d; i >= 0 && !cr[i]; i-- );
        !           668:     r->d = i;
        !           669:   }
        !           670: }
        !           671:
        !           672: void pwrup(UP n,Z e,UP *nr)
        !           673: {
        !           674:   UP y,z,t;
        !           675:   Z q,r,two;
        !           676:
        !           677:   y = UPALLOC(0); y->d = 0; y->c[0] = (Num)ONE;
        !           678:   z = n;
        !           679:   if ( !e )
        !           680:     *nr = y;
        !           681:   else if ( UNIQ(e) )
        !           682:     *nr = n;
        !           683:   else {
        !           684:     STOQ(2,two);
        !           685:     while ( 1 ) {
        !           686:       divqrz(e,two,&q,&r); e = q;
        !           687:       if ( r ) {
        !           688:         mulup(z,y,&t); y = t;
        !           689:         if ( !e ) {
        !           690:           *nr = y;
        !           691:           return;
        !           692:         }
        !           693:       }
        !           694:       mulup(z,z,&t); z = t;
        !           695:     }
        !           696:   }
        !           697: }
        !           698:
        !           699: int compup(UP n1,UP n2)
        !           700: {
        !           701:   int i,r;
        !           702:
        !           703:   if ( !n1 )
        !           704:     return n2 ? -1 : 0;
        !           705:   else if ( !n2 )
        !           706:     return 1;
        !           707:   else if ( n1->d > n2->d )
        !           708:     return 1;
        !           709:   else if ( n1->d < n2->d )
        !           710:     return -1;
        !           711:   else {
        !           712:     for ( i = n1->d; i >= 0; i-- ) {
        !           713:       r = compnum(CO,n1->c[i],n2->c[i]);
        !           714:       if ( r )
        !           715:         return r;
        !           716:     }
        !           717:     return 0;
        !           718:   }
        !           719: }
        !           720:
        !           721: void kmulp(VL vl,P n1,P n2,P *nr)
        !           722: {
        !           723:   UP b1,b2,br;
        !           724:
        !           725:   if ( !n1 || !n2 )
        !           726:     *nr = 0;
        !           727:   else if ( OID(n1) == O_N || OID(n2) == O_N )
        !           728:     mulp(vl,n1,n2,nr);
        !           729:   else {
        !           730:     ptoup(n1,&b1); ptoup(n2,&b2);
        !           731:     kmulup(b1,b2,&br);
        !           732:     uptop(br,nr);
        !           733:   }
        !           734: }
        !           735:
        !           736: void ksquarep(VL vl,P n1,P *nr)
        !           737: {
        !           738:   UP b1,br;
        !           739:
        !           740:   if ( !n1 )
        !           741:     *nr = 0;
        !           742:   else if ( OID(n1) == O_N )
        !           743:     mulp(vl,n1,n1,nr);
        !           744:   else {
        !           745:     ptoup(n1,&b1);
        !           746:     ksquareup(b1,&br);
        !           747:     uptop(br,nr);
        !           748:   }
        !           749: }
        !           750:
        !           751: void kmulup(UP n1,UP n2,UP *nr)
        !           752: {
        !           753:   int d1,d2;
        !           754:
        !           755:   if ( !n1 || !n2 ) {
        !           756:     *nr = 0; return;
        !           757:   }
        !           758:   d1 = DEG(n1)+1; d2 = DEG(n2)+1;
        !           759:   if ( (d1 < up_kara_mag) && (d2 < up_kara_mag) )
        !           760:     mulup(n1,n2,nr);
        !           761:   else
        !           762:     kmulupmain(n1,n2,nr);
        !           763: }
        !           764:
        !           765: void ksquareup(UP n1,UP *nr)
        !           766: {
        !           767:   int d1;
        !           768:
        !           769:   if ( !n1 ) {
        !           770:     *nr = 0; return;
        !           771:   }
        !           772:   d1 = DEG(n1)+1;
        !           773:   if ( (d1 < up_kara_mag) )
        !           774:     squareup(n1,nr);
        !           775:   else
        !           776:     ksquareupmain(n1,nr);
        !           777: }
        !           778:
        !           779: void copyup(UP n1,UP n2)
        !           780: {
        !           781:   n2->d = n1->d;
        !           782:   bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(Q));
        !           783: }
        !           784:
        !           785: void c_copyup(UP n,int len,pointer *p)
        !           786: {
        !           787:   if ( n )
        !           788:     bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(pointer));
        !           789: }
        !           790:
        !           791: void kmulupmain(UP n1,UP n2,UP *nr)
        !           792: {
        !           793:   int d1,d2,h;
        !           794:   UP n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
        !           795:
        !           796:   d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (MAX(d1,d2)+1)/2;
        !           797:   decompup(n1,h,&n1lo,&n1hi);
        !           798:   decompup(n2,h,&n2lo,&n2hi);
        !           799:   kmulup(n1lo,n2lo,&lo);
        !           800:   kmulup(n1hi,n2hi,&hi);
        !           801:   shiftup(hi,2*h,&t2);
        !           802:   addup(lo,t2,&t1);
        !           803:   addup(hi,lo,&mid1);
        !           804:   subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2); kmulup(s1,s2,&mid2);
        !           805:   addup(mid1,mid2,&mid);
        !           806:   shiftup(mid,h,&t2);
        !           807:   addup(t1,t2,nr);
        !           808: }
        !           809:
        !           810: void ksquareupmain(UP n1,UP *nr)
        !           811: {
        !           812:   int d1,h;
        !           813:   UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
        !           814:
        !           815:   d1 = DEG(n1)+1; h = (d1+1)/2;
        !           816:   decompup(n1,h,&n1lo,&n1hi);
        !           817:   ksquareup(n1hi,&hi); ksquareup(n1lo,&lo);
        !           818:   shiftup(hi,2*h,&t2);
        !           819:   addup(lo,t2,&t1);
        !           820:   addup(hi,lo,&mid1);
        !           821:   subup(n1hi,n1lo,&s1); ksquareup(s1,&mid2);
        !           822:   subup(mid1,mid2,&mid);
        !           823:   shiftup(mid,h,&t2);
        !           824:   addup(t1,t2,nr);
        !           825: }
        !           826:
        !           827: void rembymulup(UP n1,UP n2,UP *nr)
        !           828: {
        !           829:   int d1,d2,d;
        !           830:   UP r1,r2,inv2,t,s,q;
        !           831:
        !           832:   if ( !n2 )
        !           833:     error("rembymulup : division by 0");
        !           834:   else if ( !n1 || !n2->d )
        !           835:     *nr = 0;
        !           836:   else if ( (d1 = n1->d) < (d2 = n2->d) )
        !           837:     *nr = n1;
        !           838:   else {
        !           839:     d = d1-d2;
        !           840:     reverseup(n1,d1,&t); truncup(t,d+1,&r1);
        !           841:     reverseup(n2,d2,&r2);
        !           842:     invmodup(r2,d,&inv2);
        !           843:     kmulup(r1,inv2,&t);
        !           844:     truncup(t,d+1,&s);
        !           845:     reverseup(s,d,&q);
        !           846:     kmulup(q,n2,&t); subup(n1,t,nr);
        !           847:   }
        !           848: }
        !           849:
        !           850: /*
        !           851:   deg n2 = d
        !           852:   deg n1 <= 2*d-1
        !           853:   inv2 = inverse of reversep(n2) mod x^d
        !           854: */
        !           855:
        !           856: void hybrid_rembymulup_special(int ff,UP n1,UP n2,UP inv2,UP *nr)
        !           857: {
        !           858:   int d1,d2,d;
        !           859:   UP r1,t,s,q;
        !           860:
        !           861:   if ( !n2 )
        !           862:     error("hybrid_rembymulup : division by 0");
        !           863:   else if ( !n1 || !n2->d )
        !           864:     *nr = 0;
        !           865:   else if ( (d1 = n1->d) < (d2 = n2->d) )
        !           866:     *nr = n1;
        !           867:   else {
        !           868:     d = d1-d2;
        !           869:     reverseup(n1,d1,&t); truncup(t,d+1,&r1);
        !           870:     hybrid_tmulup(ff,r1,inv2,d+1,&t);
        !           871:     truncup(t,d+1,&s);
        !           872:     reverseup(s,d,&q);
        !           873:
        !           874:     hybrid_tmulup(ff,q,n2,d2,&t);
        !           875:     truncup(n1,d2,&s);
        !           876:     subup(s,t,nr);
        !           877:   }
        !           878: }
        !           879:
        !           880: void rembymulup_special(UP n1,UP n2,UP inv2,UP *nr)
        !           881: {
        !           882:   int d1,d2,d;
        !           883:   UP r1,t,s,q;
        !           884:
        !           885:   if ( !n2 )
        !           886:     error("rembymulup : division by 0");
        !           887:   else if ( !n1 || !n2->d )
        !           888:     *nr = 0;
        !           889:   else if ( (d1 = n1->d) < (d2 = n2->d) )
        !           890:     *nr = n1;
        !           891:   else {
        !           892:     d = d1-d2;
        !           893:     reverseup(n1,d1,&t); truncup(t,d+1,&r1);
        !           894:     tkmulup(r1,inv2,d+1,&t);
        !           895:     truncup(t,d+1,&s);
        !           896:     reverseup(s,d,&q);
        !           897:
        !           898:     tkmulup(q,n2,d2,&t);
        !           899:     truncup(n1,d2,&s);
        !           900:     subup(s,t,nr);
        !           901:   }
        !           902: }
        !           903:
        !           904: /* *nr = n1*n2 mod x^d */
        !           905:
        !           906: void tkmulup(UP n1,UP n2,int d,UP *nr)
        !           907: {
        !           908:   int m;
        !           909:   UP n1l,n1h,n2l,n2h,l,h,t,s,u;
        !           910:
        !           911:   if ( d < 0 )
        !           912:     error("tkmulup : invalid argument");
        !           913:   else if ( d == 0 )
        !           914:     *nr = 0;
        !           915:   else {
        !           916:     truncup(n1,d,&t); n1 = t;
        !           917:     truncup(n2,d,&t); n2 = t;
        !           918:     if ( !n1 || !n2 )
        !           919:       *nr = 0;
        !           920:     else if ( n1->d < up_tkara_mag || n2->d < up_tkara_mag )
        !           921:       tmulup(n1,n2,d,nr);
        !           922:     else {
        !           923:       m = (d+1)/2;
        !           924:       decompup(n1,m,&n1l,&n1h);
        !           925:       decompup(n2,m,&n2l,&n2h);
        !           926:       kmulup(n1l,n2l,&l);
        !           927:       tkmulup(n1l,n2h,d-m,&t);
        !           928:       tkmulup(n2l,n1h,d-m,&s);
        !           929:       addup(t,s,&u);
        !           930:       shiftup(u,m,&h);
        !           931:       addup(l,h,&t);
        !           932:       truncup(t,d,nr);
        !           933:     }
        !           934:   }
        !           935: }
        !           936:
        !           937: /* n->n*x^d */
        !           938:
        !           939: void shiftup(UP n,int d,UP *nr)
        !           940: {
        !           941:   int dr;
        !           942:   UP r;
        !           943:
        !           944:   if ( !n )
        !           945:     *nr = 0;
        !           946:   else {
        !           947:     dr = n->d+d;
        !           948:     *nr = r = UPALLOC(dr); r->d = dr;
        !           949:     bzero(r->c,(dr+1)*sizeof(Num));
        !           950:     bcopy(n->c,r->c+d,(n->d+1)*sizeof(Num));
        !           951:   }
        !           952: }
        !           953:
        !           954: void fft_rembymulup_special(UP n1,UP n2,UP inv2,UP *nr)
        !           955: {
        !           956:   int d1,d2,d;
        !           957:   UP r1,t,s,q,u;
        !           958:
        !           959:   if ( !n2 )
        !           960:     error("rembymulup : division by 0");
        !           961:   else if ( !n1 || !n2->d )
        !           962:     *nr = 0;
        !           963:   else if ( (d1 = n1->d) < (d2 = n2->d) )
        !           964:     *nr = n1;
        !           965:   else {
        !           966:     d = d1-d2;
        !           967:     reverseup(n1,d1,&t); truncup(t,d+1,&r1);
        !           968:     trunc_fft_mulup(r1,inv2,d+1,&t);
        !           969:     truncup(t,d+1,&s);
        !           970:     reverseup(s,d,&q);
        !           971:     trunc_fft_mulup(q,n2,d2,&t); truncup(t,d2,&u);
        !           972:     truncup(n1,d2,&s);
        !           973:     subup(s,u,nr);
        !           974:   }
        !           975: }
        !           976:
        !           977: void set_degreeup(UP n,int d)
        !           978: {
        !           979:   int i;
        !           980:
        !           981:   for ( i = d; i >= 0 && !n->c[i]; i-- );
        !           982:   n->d = i;
        !           983: }
        !           984:
        !           985: /* n -> n0 + x^d n1 */
        !           986:
        !           987: void decompup(UP n,int d,UP *n0,UP *n1)
        !           988: {
        !           989:   int dn;
        !           990:   UP r0,r1;
        !           991:
        !           992:   if ( !n || d > n->d ) {
        !           993:     *n0 = n; *n1 = 0;
        !           994:   } else if ( d < 0 )
        !           995:     error("decompup : invalid argument");
        !           996:   else if ( d == 0 ) {
        !           997:     *n0 = 0; *n1 = n;
        !           998:   } else {
        !           999:     dn = n->d;
        !          1000:     *n0 = r0 = UPALLOC(d-1);
        !          1001:     *n1 = r1 = UPALLOC(dn-d);
        !          1002:     bcopy(n->c,r0->c,d*sizeof(Num));
        !          1003:     bcopy(n->c+d,r1->c,(dn-d+1)*sizeof(Num));
        !          1004:     set_degreeup(r0,d-1);
        !          1005:     if ( r0->d < 0 )
        !          1006:       *n0 = 0;
        !          1007:     set_degreeup(r1,dn-d);
        !          1008:     if ( r1->d < 0 )
        !          1009:       *n1 = 0;
        !          1010:   }
        !          1011: }
        !          1012:
        !          1013: /* n -> n mod x^d */
        !          1014:
        !          1015: void truncup(UP n1,int d,UP *nr)
        !          1016: {
        !          1017:   int i;
        !          1018:   UP r;
        !          1019:
        !          1020:   if ( !n1 || d > n1->d )
        !          1021:     *nr = n1;
        !          1022:   else if ( d < 0 )
        !          1023:     error("truncup : invalid argument");
        !          1024:   else if ( d == 0 )
        !          1025:     *nr = 0;
        !          1026:   else {
        !          1027:      *nr = r = UPALLOC(d-1);
        !          1028:     bcopy(n1->c,r->c,d*sizeof(Num));
        !          1029:     for ( i = d-1; i >= 0 && !r->c[i]; i-- );
        !          1030:     if ( i < 0 )
        !          1031:       *nr = 0;
        !          1032:     else
        !          1033:       r->d = i;
        !          1034:   }
        !          1035: }
        !          1036:
        !          1037: int int_bits(int t)
        !          1038: {
        !          1039:   int k;
        !          1040:
        !          1041:   for ( k = 0; t; t>>=1, k++);
        !          1042:   return k;
        !          1043: }
        !          1044:
        !          1045: /* n is assumed to be LM or integer coefficient */
        !          1046:
        !          1047: int maxblenup(UP n)
        !          1048: {
        !          1049:   int m,r,i,d;
        !          1050:   Num *c;
        !          1051:
        !          1052:   if ( !n )
        !          1053:     return 0;
        !          1054:   d = n->d; c = (Num *)n->c;
        !          1055:   for ( r = 0, i = 0; i <= d; i++ )
        !          1056:     if ( c[i] ) {
        !          1057:       switch ( NID(c[i]) ) {
        !          1058:         case N_LM:
        !          1059:           m = mpz_sizeinbase(BDY((LM)c[i]),2);
        !          1060:           break;
        !          1061:         case N_Q:
        !          1062:           m = mpz_sizeinbase(BDY((Z)c[i]),2);
        !          1063:           break;
        !          1064:         default:
        !          1065:           error("maxblenup : invalid coefficient");
        !          1066:       }
        !          1067:       r = MAX(m,r);
        !          1068:     }
        !          1069:   return r;
        !          1070: }
        !          1071:
        !          1072: /* YYY */
        !          1073:
        !          1074: void uptofmarray(int mod,UP n,ModNum *f)
        !          1075: {
        !          1076:   int d,i;
        !          1077:   unsigned int r;
        !          1078:   Num *c;
        !          1079:   Z z;
        !          1080:
        !          1081:   if ( !n )
        !          1082:     return;
        !          1083:   else {
        !          1084:     d = n->d; c = (Num *)n->c;
        !          1085:     for ( i = 0; i <= d; i++ ) {
        !          1086:       if ( c[i] ) {
        !          1087:         switch ( NID(c[i]) ) {
        !          1088:           case N_LM:
        !          1089:             MPZTOZ(BDY((LM)c[i]),z);
        !          1090:             f[i] = remqi((Q)z,mod);
        !          1091:             break;
        !          1092:           case N_Q:
        !          1093:             f[i] = remqi((Q)c[i],mod);
        !          1094:             break;
        !          1095:           default:
        !          1096:             error("uptofmarray : invalid coefficient");
        !          1097:         }
        !          1098:       } else
        !          1099:         f[i] = 0;
        !          1100:     }
        !          1101:   }
        !          1102: }
        !          1103:
        !          1104: void fmarraytoup(ModNum *f,int d,UP *nr)
        !          1105: {
        !          1106:   int i;
        !          1107:   Z *c;
        !          1108:   UP r;
        !          1109:
        !          1110:   if ( d < 0 ) {
        !          1111:     *nr = 0;
        !          1112:   } else {
        !          1113:     *nr = r = UPALLOC(d); c = (Z *)r->c;
        !          1114:     for ( i = 0; i <= d; i++ ) {
        !          1115:       UTOQ((unsigned int)f[i],c[i]);
        !          1116:     }
        !          1117:     for ( i = d; i >= 0 && !c[i]; i-- );
        !          1118:     if ( i < 0 )
        !          1119:       *nr = 0;
        !          1120:     else
        !          1121:       r->d = i;
        !          1122:   }
        !          1123: }
        !          1124:
        !          1125: /* f[i]: an array of length n */
        !          1126:
        !          1127: void uiarraytoup(unsigned int **f,int n,int d,UP *nr)
        !          1128: {
        !          1129:   int i,j;
        !          1130:   Z *c;
        !          1131:   UP r;
        !          1132:
        !          1133:   if ( d < 0 ) {
        !          1134:     *nr = 0;
        !          1135:   } else {
        !          1136:     *nr = r = UPALLOC(d); c = (Z *)r->c;
        !          1137:     for ( i = 0; i <= d; i++ )
        !          1138:       intarraytoz(n,f[i],&c[i]);
        !          1139:     for ( i = d; i >= 0 && !c[i]; i-- );
        !          1140:     if ( i < 0 )
        !          1141:       *nr = 0;
        !          1142:     else
        !          1143:       r->d = i;
        !          1144:   }
        !          1145: }
        !          1146:
        !          1147: void adj_coefup(UP n,Z m,Z m2,UP *nr)
        !          1148: {
        !          1149:   int d;
        !          1150:   Z *c,*cr;
        !          1151:   Z mq;
        !          1152:   int i;
        !          1153:   UP r;
        !          1154:
        !          1155:   if ( !n )
        !          1156:     *nr = 0;
        !          1157:   else {
        !          1158:     d = n->d; c = (Z *)n->c;
        !          1159:     *nr = r = UPALLOC(d); cr = (Z *)r->c;
        !          1160:     for ( i = 0; i <= d; i++ ) {
        !          1161:       if ( !c[i] )
        !          1162:         continue;
        !          1163:       if ( cmpz(c[i],m2) > 0 )
        !          1164:         subz(c[i],m,&cr[i]);
        !          1165:       else
        !          1166:         cr[i] = c[i];
        !          1167:     }
        !          1168:     for ( i = d; i >= 0 && !cr[i]; i-- );
        !          1169:     if ( i < 0 )
        !          1170:       *nr = 0;
        !          1171:     else
        !          1172:       r->d = i;
        !          1173:   }
        !          1174: }
        !          1175:
        !          1176: /* n is assumed to have positive integer coefficients. */
        !          1177:
        !          1178: void remcup(UP n,Z mod,UP *nr)
        !          1179: {
        !          1180:   int i,d;
        !          1181:   Z *c,*cr;
        !          1182:   UP r;
        !          1183:   Z t;
        !          1184:
        !          1185:   if ( !n )
        !          1186:     *nr = 0;
        !          1187:   else {
        !          1188:     d = n->d; c = (Z *)n->c;
        !          1189:     *nr = r = UPALLOC(d); cr = (Z *)r->c;
        !          1190:     for ( i = 0; i <= d; i++ )
        !          1191:       remz(c[i],mod,&cr[i]);
        !          1192:     for ( i = d; i >= 0 && !cr[i]; i-- );
        !          1193:     if ( i < 0 )
        !          1194:       *nr = 0;
        !          1195:     else
        !          1196:       r->d = i;
        !          1197:   }
        !          1198: }
        !          1199:
        !          1200: void fft_mulup_main(UP,UP,int,UP *);
        !          1201:
        !          1202: void fft_mulup(UP n1,UP n2,UP *nr)
        !          1203: {
        !          1204:   int d1,d2,d,b1,b2,h;
        !          1205:   UP n1lo,n1hi,n2lo,n2hi,lo,hi,t1,t2,mid1,mid2,mid,s1,s2;
        !          1206:
        !          1207:   if ( !n1 || !n2 )
        !          1208:     *nr = 0;
        !          1209:   else {
        !          1210:     d1 = n1->d; d2 = n2->d;
        !          1211:     if ( MAX(d1,d2) < up_fft_mag )
        !          1212:       kmulup(n1,n2,nr);
        !          1213:     else {
        !          1214:       b1 = maxblenup(n1); b2 = maxblenup(n2);
        !          1215:       if ( fft_available(d1,b1,d2,b2) )
        !          1216:         fft_mulup_main(n1,n2,0,nr);
        !          1217:       else {
        !          1218:         d = MAX(d1,d2)+1; h = (d+1)/2;
        !          1219:         decompup(n1,h,&n1lo,&n1hi);
        !          1220:         decompup(n2,h,&n2lo,&n2hi);
        !          1221:         fft_mulup(n1lo,n2lo,&lo);
        !          1222:         fft_mulup(n1hi,n2hi,&hi);
        !          1223:         shiftup(hi,2*h,&t2);
        !          1224:         addup(lo,t2,&t1);
        !          1225:         addup(hi,lo,&mid1);
        !          1226:         subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2);
        !          1227:         fft_mulup(s1,s2,&mid2);
        !          1228:         addup(mid1,mid2,&mid);
        !          1229:         shiftup(mid,h,&t2);
        !          1230:         addup(t1,t2,nr);
        !          1231:       }
        !          1232:     }
        !          1233:   }
        !          1234: }
        !          1235:
        !          1236: void trunc_fft_mulup(UP n1,UP n2,int dbd,UP *nr)
        !          1237: {
        !          1238:   int d1,d2,b1,b2,m;
        !          1239:   UP n1l,n1h,n2l,n2h,l,h,t,s,u;
        !          1240:
        !          1241:   if ( !n1 || !n2 )
        !          1242:     *nr = 0;
        !          1243:   else if ( dbd < 0 )
        !          1244:     error("trunc_fft_mulup : invalid argument");
        !          1245:   else if ( dbd == 0 )
        !          1246:     *nr = 0;
        !          1247:   else {
        !          1248:     truncup(n1,dbd,&t); n1 = t;
        !          1249:     truncup(n2,dbd,&t); n2 = t;
        !          1250:     d1 = n1->d; d2 = n2->d;
        !          1251:     if ( MAX(d1,d2) < up_fft_mag )
        !          1252:       tkmulup(n1,n2,dbd,nr);
        !          1253:     else {
        !          1254:       b1 = maxblenup(n1); b2 = maxblenup(n2);
        !          1255:       if ( fft_available(d1,b1,d2,b2) )
        !          1256:         fft_mulup_main(n1,n2,dbd,nr);
        !          1257:       else {
        !          1258:         m = (dbd+1)/2;
        !          1259:         decompup(n1,m,&n1l,&n1h);
        !          1260:         decompup(n2,m,&n2l,&n2h);
        !          1261:         fft_mulup(n1l,n2l,&l);
        !          1262:         trunc_fft_mulup(n1l,n2h,dbd-m,&t);
        !          1263:         trunc_fft_mulup(n2l,n1h,dbd-m,&s);
        !          1264:         addup(t,s,&u);
        !          1265:         shiftup(u,m,&h);
        !          1266:         addup(l,h,&t);
        !          1267:         truncup(t,dbd,nr);
        !          1268:       }
        !          1269:     }
        !          1270:   }
        !          1271: }
        !          1272:
        !          1273: void fft_squareup(UP n1,UP *nr)
        !          1274: {
        !          1275:   int d1,d,h,b1;
        !          1276:   UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
        !          1277:
        !          1278:   if ( !n1 )
        !          1279:     *nr = 0;
        !          1280:   else {
        !          1281:     d1 = n1->d;
        !          1282:     if ( d1 < up_fft_mag )
        !          1283:       ksquareup(n1,nr);
        !          1284:     else {
        !          1285:       b1 = maxblenup(n1);
        !          1286:       if ( fft_available(d1,b1,d1,b1) )
        !          1287:         fft_mulup_main(n1,n1,0,nr);
        !          1288:       else {
        !          1289:         d = d1+1; h = (d1+1)/2;
        !          1290:         decompup(n1,h,&n1lo,&n1hi);
        !          1291:         fft_squareup(n1hi,&hi);
        !          1292:         fft_squareup(n1lo,&lo);
        !          1293:         shiftup(hi,2*h,&t2);
        !          1294:         addup(lo,t2,&t1);
        !          1295:         addup(hi,lo,&mid1);
        !          1296:         subup(n1hi,n1lo,&s1);
        !          1297:         fft_squareup(s1,&mid2);
        !          1298:         subup(mid1,mid2,&mid);
        !          1299:         shiftup(mid,h,&t2);
        !          1300:         addup(t1,t2,nr);
        !          1301:       }
        !          1302:     }
        !          1303:   }
        !          1304: }
        !          1305:
        !          1306: /*
        !          1307:  * dbd == 0 => n1 * n2
        !          1308:  * dbd > 0  => n1 * n2 mod x^dbd
        !          1309:  * n1 == n2 => squaring
        !          1310:  */
        !          1311:
        !          1312: void fft_mulup_main(UP n1,UP n2,int dbd,UP *nr)
        !          1313: {
        !          1314:   ModNum *f1,*f2,*w,*fr;
        !          1315:   ModNum **frarray,**fa;
        !          1316:   int *modarray,*ma;
        !          1317:   int modarray_len;
        !          1318:   int frarray_index = 0;
        !          1319:   Z m,m1,m2,two,rem;
        !          1320:   int d1,d2,dmin,i,mod,root,d,cond,bound;
        !          1321:   UP r;
        !          1322:
        !          1323:   if ( !n1 || !n2 ) {
        !          1324:     *nr = 0; return;
        !          1325:   }
        !          1326:   d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
        !          1327:   if ( !d1 || !d2 ) {
        !          1328:     mulup(n1,n2,nr); return;
        !          1329:   }
        !          1330:   m = ONE;
        !          1331:   bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;
        !          1332:   f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
        !          1333:   if ( n1 == n2 )
        !          1334:     f2 = 0;
        !          1335:   else
        !          1336:     f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
        !          1337:   w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
        !          1338:   modarray_len = 1024;
        !          1339:   modarray = (int *)MALLOC_ATOMIC(modarray_len*sizeof(unsigned int));
        !          1340:   frarray = (ModNum **)MALLOC(modarray_len*sizeof(ModNum *));
        !          1341:   for ( i = 0; i < NPrimes; i++ ) {
        !          1342:     FFT_primes(i,&mod,&root,&d);
        !          1343:     if ( (1<<d) < d1+d2+1 )
        !          1344:       continue;
        !          1345:     if ( frarray_index == modarray_len ) {
        !          1346:       ma = (int *)MALLOC_ATOMIC(2*modarray_len*sizeof(unsigned int));
        !          1347:       bcopy(modarray,ma,modarray_len*sizeof(unsigned int));
        !          1348:       modarray = ma;
        !          1349:       fa = (ModNum **)MALLOC(2*modarray_len*sizeof(ModNum *));
        !          1350:       bcopy(frarray,fa,modarray_len*sizeof(ModNum *));
        !          1351:       frarray = fa;
        !          1352:       modarray_len *= 2;
        !          1353:     }
        !          1354:     modarray[frarray_index] = mod;
        !          1355:     frarray[frarray_index++] = fr
        !          1356:       = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
        !          1357:     uptofmarray(mod,n1,f1);
        !          1358:     if ( !f2 )
        !          1359:       cond = FFT_pol_square(d1,f1,fr,i,w);
        !          1360:     else {
        !          1361:       uptofmarray(mod,n2,f2);
        !          1362:       cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
        !          1363:     }
        !          1364:     if ( cond )
        !          1365:       error("fft_mulup : error in FFT_pol_product");
        !          1366:     STOQ(mod,m1); mulz(m,m1,&m2); m = m2;
        !          1367:     if ( z_bits((Q)m) > bound ) {
        !          1368:       if ( !dbd )
        !          1369:         dbd = d1+d2+1;
        !          1370:       crup(frarray,MIN(d1+d2,dbd-1),modarray,frarray_index,m,&r);
        !          1371:       STOQ(2,two);
        !          1372:       divqrz(m,two,&m2,&rem);
        !          1373:       adj_coefup(r,m,m2,nr);
        !          1374:       return;
        !          1375:     }
        !          1376:   }
        !          1377:   error("fft_mulup : FFT_primes exhausted");
        !          1378: }
        !          1379:
        !          1380: /* improved version */
        !          1381:
        !          1382: void crup(ModNum **f,int d,int *mod,int index,Z m,UP *r)
        !          1383: {
        !          1384:   mpz_t cof,c,rem;
        !          1385:   mpz_t *s;
        !          1386:   UP u;
        !          1387:   int inv,i,j,t;
        !          1388:
        !          1389:   mpz_init(c); mpz_init(cof); mpz_init(rem);
        !          1390:   s = (mpz_t *)MALLOC((d+1)*sizeof(mpz_t));
        !          1391:   for ( j = 0; j <= d; j++ )
        !          1392:     mpz_init_set_ui(s[j],0);
        !          1393:   for ( i = 0; i < index; i++ ) {
        !          1394:     mpz_fdiv_q_ui(cof,BDY(m),mod[i]);
        !          1395:     t = mpz_fdiv_r_ui(rem,cof,mod[i]);
        !          1396:     inv = invm(t,mod[i]);
        !          1397:     mpz_mul_ui(c,cof,inv);
        !          1398:     /* s += c*f[i] */
        !          1399:     for ( j = 0; j <= d; j++ )
        !          1400:       if ( f[i][j] )
        !          1401:         mpz_addmul_ui(s[j],c,f[i][j]);
        !          1402:   }
        !          1403:   for ( i = d; i >= 0; i-- )
        !          1404:     if ( s[i] )
        !          1405:       break;
        !          1406:   if ( i < 0 )
        !          1407:     *r = 0;
        !          1408:   else {
        !          1409:     u = UPALLOC(i);
        !          1410:     DEG(u) = i;
        !          1411:     for ( j = 0; j <= i; j++ )
        !          1412:       COEF(u)[j] = (Num)s[j];
        !          1413:     remcup(u,m,r);
        !          1414:   }
        !          1415: }
        !          1416:
        !          1417: /*
        !          1418:  * dbd == 0 => n1 * n2
        !          1419:  * dbd > 0  => n1 * n2 mod x^dbd
        !          1420:  * n1 == n2 => squaring
        !          1421:  * return: n1*n2 mod Primes[modind[0]]*.prime...*Primes[modind[nmod-1]].prime
        !          1422:  */
        !          1423:
        !          1424: void fft_mulup_specialmod_main(UP n1,UP n2,int dbd,int *modind,int nmod,UP *nr)
        !          1425: {
        !          1426:   ModNum *f1,*f2,*w,*fr;
        !          1427:   ModNum **frarray;
        !          1428:   Z m,m1,m2;
        !          1429:   int *modarray;
        !          1430:   int d1,d2,dmin,i,root,d,cond,bound;
        !          1431:
        !          1432:   if ( !n1 || !n2 ) {
        !          1433:     *nr = 0; return;
        !          1434:   }
        !          1435:   d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
        !          1436:   if ( !d1 || !d2 ) {
        !          1437:     mulup(n1,n2,nr); return;
        !          1438:   }
        !          1439:   m = ONE;
        !          1440:   bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;
        !          1441:   f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
        !          1442:   if ( n1 == n2 )
        !          1443:     f2 = 0;
        !          1444:   else
        !          1445:     f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
        !          1446:   w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
        !          1447:   frarray = (ModNum **)MALLOC(nmod*sizeof(ModNum *));
        !          1448:   modarray = (int *)MALLOC_ATOMIC(nmod*sizeof(unsigned int *));
        !          1449:
        !          1450:   for ( i = 0; i < nmod; i++ ) {
        !          1451:     FFT_primes(modind[i],&modarray[i],&root,&d);
        !          1452:       if ( (1<<d) < d1+d2+1 )
        !          1453:         error("fft_mulup_specialmod_main : invalid modulus");
        !          1454:     frarray[i] = fr
        !          1455:       = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
        !          1456:     uptofmarray(modarray[i],n1,f1);
        !          1457:     if ( !f2 )
        !          1458:       cond = FFT_pol_square(d1,f1,fr,modind[i],w);
        !          1459:     else {
        !          1460:       uptofmarray(modarray[i],n2,f2);
        !          1461:       cond = FFT_pol_product(d1,f1,d2,f2,fr,modind[i],w);
        !          1462:     }
        !          1463:     if ( cond )
        !          1464:       error("fft_mulup_specialmod_main : error in FFT_pol_product");
        !          1465:     STOQ(modarray[i],m1); mulz(m,m1,&m2); m = m2;
        !          1466:   }
        !          1467:   if ( !dbd )
        !          1468:     dbd = d1+d2+1;
        !          1469:   crup(frarray,MIN(d1+d2,dbd-1),modarray,nmod,m,nr);
        !          1470: }

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