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

Annotation of OpenXM_contrib2/asir2000/engine/up_gfpn.c, Revision 1.9

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

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