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

Annotation of OpenXM_contrib2/asir2000/engine/up_lm.c, Revision 1.12

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.12    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/up_lm.c,v 1.11 2015/08/14 13:51:55 fujimoto Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include <math.h>
                     52:
                     53: extern GC_dont_gc;
                     54: extern struct oEGT eg_chrem,eg_fore,eg_back;
                     55: extern int debug_up;
                     56: extern int up_lazy;
1.4       noro       57: extern int current_ff;
1.1       noro       58:
1.6       noro       59: void fft_mulup_lm(UP n1,UP n2,UP *nr)
1.1       noro       60: {
1.12    ! noro       61:   ModNum *f1,*f2,*w,*fr;
        !            62:   ModNum *frarray[1024];
        !            63:   int modarray[1024];
        !            64:   int frarray_index = 0;
        !            65:   N m,m1,m2,lm_mod;
        !            66:   int d1,d2,dmin,i,mod,root,d,cond,bound;
        !            67:   UP r;
        !            68:
        !            69:   if ( !n1 || !n2 ) {
        !            70:     *nr = 0; return;
        !            71:   }
        !            72:   d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
        !            73:   if ( !d1 || !d2 ) {
        !            74:     mulup(n1,n2,nr); return;
        !            75:   }
        !            76:   getmod_lm(&lm_mod);
        !            77:   if ( !lm_mod )
        !            78:     error("fft_mulup_lm : current_mod_lm is not set");
        !            79:   m = ONEN;
        !            80:   bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;
        !            81:   f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !            82:   f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !            83:   w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
        !            84:   for ( i = 0; i < NPrimes; i++ ) {
        !            85:     FFT_primes(i,&mod,&root,&d);
        !            86:     if ( (1<<d) < d1+d2+1 )
        !            87:       continue;
        !            88:     modarray[frarray_index] = mod;
        !            89:     frarray[frarray_index++] = fr
        !            90:       = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !            91:     uptofmarray(mod,n1,f1);
        !            92:     uptofmarray(mod,n2,f2);
        !            93:     cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
        !            94:     if ( cond )
        !            95:       error("fft_mulup : error in FFT_pol_product");
        !            96:     STON(mod,m1); muln(m,m1,&m2); m = m2;
        !            97:     if ( n_bits(m) > bound ) {
        !            98: /*      GC_dont_gc = 1; */
        !            99:       crup_lm(frarray,d1+d2,modarray,frarray_index,m,lm_mod,&r);
        !           100:       uptolmup(r,nr);
        !           101:       GC_dont_gc = 0;
        !           102:       return;
        !           103:     }
        !           104:   }
        !           105:   error("fft_mulup : FFT_primes exhausted");
1.1       noro      106: }
                    107:
1.6       noro      108: void fft_squareup_lm(UP n1,UP *nr)
1.1       noro      109: {
1.12    ! noro      110:   ModNum *f1,*w,*fr;
        !           111:   ModNum *frarray[1024];
        !           112:   int modarray[1024];
        !           113:   int frarray_index = 0;
        !           114:   N m,m1,m2,lm_mod;
        !           115:   int d1,dmin,i,mod,root,d,cond,bound;
        !           116:   UP r;
        !           117:
        !           118:   if ( !n1 ) {
        !           119:     *nr = 0; return;
        !           120:   }
        !           121:   d1 = n1->d; dmin = d1;
        !           122:   if ( !d1 ) {
        !           123:     mulup(n1,n1,nr); return;
        !           124:   }
        !           125:   getmod_lm(&lm_mod);
        !           126:   if ( !lm_mod )
        !           127:     error("fft_squareup_lm : current_mod_lm is not set");
        !           128:   m = ONEN;
        !           129:   bound = 2*maxblenup(n1)+int_bits(d1)+2;
        !           130:   f1 = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));
        !           131:   w = (ModNum *)ALLOCA(6*(1<<int_bits(2*d1+1))*sizeof(ModNum));
        !           132:   for ( i = 0; i < NPrimes; i++ ) {
        !           133:     FFT_primes(i,&mod,&root,&d);
        !           134:     if ( (1<<d) < 2*d1+1 )
        !           135:       continue;
        !           136:     modarray[frarray_index] = mod;
        !           137:     frarray[frarray_index++] = fr
        !           138:       = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));
        !           139:     uptofmarray(mod,n1,f1);
        !           140:     cond = FFT_pol_square(d1,f1,fr,i,w);
        !           141:     if ( cond )
        !           142:       error("fft_mulup : error in FFT_pol_product");
        !           143:     STON(mod,m1); muln(m,m1,&m2); m = m2;
        !           144:     if ( n_bits(m) > bound ) {
        !           145: /*      GC_dont_gc = 1; */
        !           146:       crup_lm(frarray,2*d1,modarray,frarray_index,m,lm_mod,&r);
        !           147:       uptolmup(r,nr);
        !           148:       GC_dont_gc = 0;
        !           149:       return;
        !           150:     }
        !           151:   }
        !           152:   error("fft_squareup : FFT_primes exhausted");
1.1       noro      153: }
                    154:
1.6       noro      155: void trunc_fft_mulup_lm(UP n1,UP n2,int dbd,UP *nr)
1.1       noro      156: {
1.12    ! noro      157:   ModNum *f1,*f2,*fr,*w;
        !           158:   ModNum *frarray[1024];
        !           159:   int modarray[1024];
        !           160:   int frarray_index = 0;
        !           161:   N m,m1,m2,lm_mod;
        !           162:   int d1,d2,dmin,i,mod,root,d,cond,bound;
        !           163:   UP r;
        !           164:
        !           165:   if ( !n1 || !n2 ) {
        !           166:     *nr = 0; return;
        !           167:   }
        !           168:   d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
        !           169:   if ( !d1 || !d2 ) {
        !           170:     tmulup(n1,n2,dbd,nr); return;
        !           171:   }
        !           172:   getmod_lm(&lm_mod);
        !           173:   if ( !lm_mod )
        !           174:     error("trunc_fft_mulup_lm : current_mod_lm is not set");
        !           175:   m = ONEN;
        !           176:   bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;
        !           177:   f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !           178:   f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !           179:   w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
        !           180:   for ( i = 0; i < NPrimes; i++ ) {
        !           181:     FFT_primes(i,&mod,&root,&d);
        !           182:     if ( (1<<d) < d1+d2+1 )
        !           183:       continue;
        !           184:
        !           185:     modarray[frarray_index] = mod;
        !           186:     frarray[frarray_index++] = fr
        !           187:       = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !           188:     uptofmarray(mod,n1,f1);
        !           189:     uptofmarray(mod,n2,f2);
        !           190:     cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
        !           191:     if ( cond )
        !           192:       error("fft_mulup : error in FFT_pol_product");
        !           193:     STON(mod,m1); muln(m,m1,&m2); m = m2;
        !           194:     if ( n_bits(m) > bound ) {
        !           195: /*      GC_dont_gc = 1; */
        !           196:       crup_lm(frarray,MIN(dbd-1,d1+d2),modarray,frarray_index,m,lm_mod,&r);
        !           197:       uptolmup(r,nr);
        !           198:       GC_dont_gc = 0;
        !           199:       return;
        !           200:     }
        !           201:   }
        !           202:   error("trunc_fft_mulup : FFT_primes exhausted");
1.1       noro      203: }
                    204:
1.6       noro      205: void crup_lm(ModNum **f,int d,int *mod,int index,N m,N lm_mod,UP *r)
1.1       noro      206: {
1.12    ! noro      207:   double *k;
        !           208:   double c2;
        !           209:   unsigned int *w;
        !           210:   unsigned int zi,au,al;
        !           211:   UL a;
        !           212:   int i,j,l,len;
        !           213:   UP s,s1,u;
        !           214:   struct oUP c;
        !           215:   N t,ci,mi,qn;
        !           216:   unsigned int **sum;
        !           217:   unsigned int *sum_b;
        !           218:   Q q;
        !           219:
        !           220:   if ( !lm_mod )
        !           221:     error("crup_lm : current_mod_lm is not set");
        !           222:   k = (double *)ALLOCA((d+1)*sizeof(double));
        !           223:   for ( j = 0; j <= d; j++ )
        !           224:     k[j] = 0.5;
        !           225:   up_lazy = 1;
        !           226:   sum = (unsigned int **)ALLOCA((d+1)*sizeof(unsigned int *));
        !           227:   len = (int_bits(index)+n_bits(lm_mod)+31)/32+1;
        !           228:   w = (unsigned int *)ALLOCA((len+1)*sizeof(unsigned int));
        !           229:   sum_b = (unsigned int *)MALLOC_ATOMIC((d+1)*len*sizeof(unsigned int));
        !           230:   for ( j = 0; j <= d; j++ ) {
        !           231:     sum[j] = sum_b+len*j;
        !           232:     bzero((char *)sum[j],len*sizeof(unsigned int));
        !           233:   }
        !           234:   for ( i = 0, s = 0; i < index; i++ ) {
        !           235:     divin(m,mod[i],&ci);
        !           236:     zi = (unsigned int)invm((unsigned int)rem(ci,mod[i]),mod[i]);
        !           237:
        !           238:     STON(zi,t); muln(ci,t,&mi);
        !           239:     divn(mi,lm_mod,&qn,&t);
        !           240:     if ( t )
        !           241:       for ( j = 0; j <= d; j++ ) {
        !           242:         bzero((char *)w,(len+1)*sizeof(unsigned int));
        !           243:         muln_1(BD(t),PL(t),(unsigned int)f[i][j],w);
        !           244:         for ( l = PL(t); l >= 0 && !w[l]; l--);
        !           245:         l++;
        !           246:         if ( l )
        !           247:           addarray_to(w,l,sum[j],len);
        !           248:       }
        !           249:     c2 = (double)zi/(double)mod[i];
        !           250:     for ( j = 0; j <= d; j++ )
        !           251:       k[j] += c2*f[i][j];
        !           252:   }
        !           253:   uiarraytoup(sum,len,d,&s);
        !           254:   GCFREE(sum_b);
1.1       noro      255:
1.12    ! noro      256:   u = UPALLOC(d);
        !           257:   for ( j = 0; j <= d; j++ ) {
1.1       noro      258: #if 1
1.12    ! noro      259:     a = (UL)floor(k[j]);
1.11      fujimoto  260: #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
1.12    ! noro      261:     au = ((unsigned int *)&a)[1];
        !           262:     al = ((unsigned int *)&a)[0];
1.1       noro      263: #else
1.12    ! noro      264:     al = ((unsigned int *)&a)[1];
        !           265:     au = ((unsigned int *)&a)[0];
1.1       noro      266: #endif
1.12    ! noro      267:     if ( au ) {
        !           268:       NEWQ(q); SGN(q) = 1; NM(q)=NALLOC(2); DN(q)=0;
        !           269:       PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
        !           270:     } else
        !           271:       UTOQ(al,q);
1.1       noro      272: #else
1.12    ! noro      273:     al = (int)floor(k[j]); STOQ(al,q);
1.1       noro      274: #endif
1.12    ! noro      275:      u->c[j] = (Num)q;
        !           276:   }
        !           277:   for ( j = d; j >= 0 && !u->c[j]; j-- );
        !           278:   if ( j < 0 )
        !           279:     u = 0;
        !           280:   else
        !           281:     u->d = j;
        !           282:   divn(m,lm_mod,&qn,&t); NTOQ(t,-1,q);
        !           283:   c.d = 0; c.c[0] = (Num)q;
        !           284:   mulup(u,&c,&s1);
        !           285:   up_lazy = 0;
1.1       noro      286:
1.12    ! noro      287:   addup(s,s1,r);
1.1       noro      288: }
                    289:
1.6       noro      290: void fft_rembymulup_special_lm(UP n1,UP n2,UP inv2,UP *nr)
1.1       noro      291: {
1.12    ! noro      292:   int d1,d2,d;
        !           293:   UP r1,t,s,q,u;
1.1       noro      294:
1.12    ! noro      295:   if ( !n2 )
        !           296:     error("rembymulup : division by 0");
        !           297:   else if ( !n1 || !n2->d )
        !           298:     *nr = 0;
        !           299:   else if ( (d1 = n1->d) < (d2 = n2->d) )
        !           300:     *nr = n1;
        !           301:   else {
        !           302:     d = d1-d2;
        !           303:     reverseup(n1,d1,&t); truncup(t,d+1,&r1);
        !           304:     trunc_fft_mulup_lm(r1,inv2,d+1,&t);
        !           305:     truncup(t,d+1,&s);
        !           306:     reverseup(s,d,&q);
        !           307:     trunc_fft_mulup_lm(q,n2,d2,&t); truncup(t,d2,&u);
        !           308:     truncup(n1,d2,&s);
        !           309:     subup(s,u,nr);
        !           310:   }
1.1       noro      311: }
                    312:
1.6       noro      313: void uptolmup(UP n,UP *nr)
1.1       noro      314: {
1.12    ! noro      315:   int i,d;
        !           316:   Q *c;
        !           317:   LM *cr;
        !           318:   UP r;
        !           319:
        !           320:   if ( !n )
        !           321:     *nr = 0;
        !           322:   else {
        !           323:     d = n->d; c = (Q *)n->c;
        !           324:     *nr = r = UPALLOC(d); cr = (LM *)r->c;
        !           325:     for ( i = 0; i <= d; i++ )
        !           326:       qtolm(c[i],&cr[i]);
        !           327:     for ( i = d; i >= 0 && !cr[i]; i-- );
        !           328:     if ( i < 0 )
        !           329:       *nr = 0;
        !           330:     else
        !           331:       r->d = i;
        !           332:   }
1.1       noro      333: }
                    334:
1.6       noro      335: void save_up(UP obj,char *name)
1.1       noro      336: {
1.12    ! noro      337:   P p;
        !           338:   Obj ret;
        !           339:   NODE n0,n1;
        !           340:   STRING s;
        !           341:   void Pbsave();
        !           342:
        !           343:   uptop(obj,&p);
        !           344:   MKSTR(s,name);
        !           345:   MKNODE(n1,s,0); MKNODE(n0,p,n1);
        !           346:   Pbsave(n0,&ret);
1.1       noro      347: }
                    348:
1.6       noro      349: void hybrid_powermodup(UP f,UP *xp)
1.1       noro      350: {
1.12    ! noro      351:   N n;
        !           352:   UP x,y,t,invf,s;
        !           353:   int k;
        !           354:   LM lm;
        !           355:
        !           356:   getmod_lm(&n);
        !           357:   if ( !n )
        !           358:     error("hybrid_powermodup : current_mod_lm is not set");
        !           359:   MKLM(ONEN,lm);
        !           360:   x = UPALLOC(1); x->d = 1; x->c[1] = (Num)lm;
        !           361:   y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
        !           362:
        !           363:   reverseup(f,f->d,&t);
        !           364:   invmodup(t,f->d,&s); uptolmup(s,&invf);
        !           365:   for ( k = n_bits(n)-1; k >= 0; k-- ) {
        !           366:     hybrid_squareup(FF_GFP,y,&t);
        !           367:     hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);
        !           368:     y = s;
        !           369:     if ( n->b[k/32] & (1<<(k%32)) ) {
        !           370:       mulup(y,x,&t);
        !           371:       remup(t,f,&s);
        !           372:       y = s;
        !           373:     }
        !           374:   }
        !           375:   *xp = y;
1.1       noro      376: }
                    377:
1.6       noro      378: void powermodup(UP f,UP *xp)
1.1       noro      379: {
1.12    ! noro      380:   N n;
        !           381:   UP x,y,t,invf,s;
        !           382:   int k;
        !           383:   Num c;
        !           384:
        !           385:   if ( !current_ff )
        !           386:     error("powermodup : current_ff is not set");
        !           387:   field_order_ff(&n);
        !           388:   one_ff(&c);
        !           389:   x = UPALLOC(1); x->d = 1; x->c[1] = c;
        !           390:   y = UPALLOC(0); y->d = 0; y->c[0] = c;
        !           391:
        !           392:   reverseup(f,f->d,&t);
        !           393:   invmodup(t,f->d,&s);
        !           394:   switch ( current_ff ) {
        !           395:     case FF_GFP:
        !           396:     case FF_GFPN:
        !           397:       uptolmup(s,&invf);
        !           398:       break;
        !           399:     case FF_GFS:
        !           400:     case FF_GFSN:
        !           401:       invf = s; /* XXX */
        !           402:       break;
        !           403:     default:
        !           404:       error("powermodup : not implemented yet");
        !           405:   }
        !           406:   for ( k = n_bits(n)-1; k >= 0; k-- ) {
        !           407:     ksquareup(y,&t);
        !           408:     rembymulup_special(t,f,invf,&s);
        !           409:     y = s;
        !           410:     if ( n->b[k/32] & (1<<(k%32)) ) {
        !           411:       mulup(y,x,&t);
        !           412:       remup(t,f,&s);
        !           413:       y = s;
        !           414:     }
        !           415:   }
        !           416:   *xp = y;
1.1       noro      417: }
                    418:
                    419: /* g^d mod f */
                    420:
1.6       noro      421: void hybrid_generic_powermodup(UP g,UP f,Q d,UP *xp)
1.1       noro      422: {
1.12    ! noro      423:   N e;
        !           424:   UP x,y,t,invf,s;
        !           425:   int k;
        !           426:   LM lm;
        !           427:
        !           428:   e = NM(d);
        !           429:   MKLM(ONEN,lm);
        !           430:   y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
        !           431:   remup(g,f,&x);
        !           432:   if ( !x ) {
        !           433:     *xp = !d ? y : 0;
        !           434:     return;
        !           435:   } else if ( !x->d ) {
        !           436:     pwrup(x,d,xp);
        !           437:     return;
        !           438:   }
        !           439:   reverseup(f,f->d,&t);
        !           440:   invmodup(t,f->d,&invf);
        !           441:   for ( k = n_bits(e)-1; k >= 0; k-- ) {
        !           442:     hybrid_squareup(FF_GFP,y,&t);
        !           443:     hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);
        !           444:     y = s;
        !           445:     if ( e->b[k/32] & (1<<(k%32)) ) {
        !           446:       mulup(y,x,&t);
        !           447:       remup(t,f,&s);
        !           448:       y = s;
        !           449:     }
        !           450:   }
        !           451:   *xp = y;
1.1       noro      452: }
                    453:
1.6       noro      454: void generic_powermodup(UP g,UP f,Q d,UP *xp)
1.1       noro      455: {
1.12    ! noro      456:   N e;
        !           457:   UP x,y,t,invf,s;
        !           458:   int k;
        !           459:   Num c;
        !           460:
        !           461:   e = NM(d);
        !           462:   one_ff(&c);
        !           463:   y = UPALLOC(0); y->d = 0; y->c[0] = c;
        !           464:   remup(g,f,&x);
        !           465:   if ( !x ) {
        !           466:     *xp = !d ? y : 0;
        !           467:     return;
        !           468:   } else if ( !x->d ) {
        !           469:     pwrup(x,d,xp);
        !           470:     return;
        !           471:   }
        !           472:   reverseup(f,f->d,&t);
        !           473:   invmodup(t,f->d,&invf);
        !           474:   for ( k = n_bits(e)-1; k >= 0; k-- ) {
        !           475:     ksquareup(y,&t);
        !           476:     rembymulup_special(t,f,invf,&s);
        !           477:     y = s;
        !           478:     if ( e->b[k/32] & (1<<(k%32)) ) {
        !           479:       mulup(y,x,&t);
        !           480:       remup(t,f,&s);
        !           481:       y = s;
        !           482:     }
        !           483:   }
        !           484:   *xp = y;
1.1       noro      485: }
                    486:
1.6       noro      487: void hybrid_powertabup(UP f,UP xp,UP *tab)
1.1       noro      488: {
1.12    ! noro      489:   UP y,t,invf;
        !           490:   int i,d;
        !           491:   LM lm;
        !           492:
        !           493:   d = f->d;
        !           494:   MKLM(ONEN,lm);
        !           495:   y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
        !           496:   tab[0] = y;
        !           497:   tab[1] = xp;
        !           498:
        !           499:   reverseup(f,f->d,&t);
        !           500:   invmodup(t,f->d,&invf);
        !           501:
        !           502:   for ( i = 2; i < d; i++ ) {
        !           503:     if ( debug_up ){
        !           504:       fprintf(stderr,".");
        !           505:     }
        !           506:     hybrid_mulup(FF_GFP,tab[i-1],xp,&t);
        !           507:     hybrid_rembymulup_special(FF_GFP,t,f,invf,&tab[i]);
        !           508:   }
1.1       noro      509: }
                    510:
1.6       noro      511: void powertabup(UP f,UP xp,UP *tab)
1.1       noro      512: {
1.12    ! noro      513:   UP y,t,invf;
        !           514:   int i,d;
        !           515:   Num c;
        !           516:
        !           517:   d = f->d;
        !           518:   one_ff(&c);
        !           519:   y = UPALLOC(0); y->d = 0; y->c[0] = c;
        !           520:   tab[0] = y;
        !           521:   tab[1] = xp;
        !           522:
        !           523:   reverseup(f,f->d,&t);
        !           524:   invmodup(t,f->d,&invf);
        !           525:
        !           526:   for ( i = 2; i < d; i++ ) {
        !           527:     if ( debug_up ){
        !           528:       fprintf(stderr,".");
        !           529:     }
        !           530:     kmulup(tab[i-1],xp,&t);
        !           531:     rembymulup_special(t,f,invf,&tab[i]);
        !           532:   }
1.1       noro      533: }

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