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

Annotation of OpenXM_contrib2/asir2018/engine/up_lm.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: extern struct oEGT eg_chrem,eg_fore,eg_back;
        !            54: extern int debug_up;
        !            55: extern int up_lazy;
        !            56: extern int current_ff;
        !            57:
        !            58: void fft_mulup_lm(UP n1,UP n2,UP *nr)
        !            59: {
        !            60:   ModNum *f1,*f2,*w,*fr;
        !            61:   ModNum *frarray[1024];
        !            62:   int modarray[1024];
        !            63:   int frarray_index = 0;
        !            64:   Z m,m1,m2,lm_mod;
        !            65:   int d1,d2,dmin,i,mod,root,d,cond,bound;
        !            66:   UP r;
        !            67:
        !            68:   if ( !n1 || !n2 ) {
        !            69:     *nr = 0; return;
        !            70:   }
        !            71:   d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
        !            72:   if ( !d1 || !d2 ) {
        !            73:     mulup(n1,n2,nr); return;
        !            74:   }
        !            75:   getmod_lm(&lm_mod);
        !            76:   if ( !lm_mod )
        !            77:     error("fft_mulup_lm : current_mod_lm is not set");
        !            78:   m = ONE;
        !            79:   bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;
        !            80:   f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !            81:   f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !            82:   w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
        !            83:   for ( i = 0; i < NPrimes; i++ ) {
        !            84:     FFT_primes(i,&mod,&root,&d);
        !            85:     if ( (1<<d) < d1+d2+1 )
        !            86:       continue;
        !            87:     modarray[frarray_index] = mod;
        !            88:     frarray[frarray_index++] = fr
        !            89:       = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !            90:     uptofmarray(mod,n1,f1);
        !            91:     uptofmarray(mod,n2,f2);
        !            92:     cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
        !            93:     if ( cond )
        !            94:       error("fft_mulup : error in FFT_pol_product");
        !            95:     STOQ(mod,m1); mulz(m,m1,&m2); m = m2;
        !            96:     if ( z_bits((Q)m) > bound ) {
        !            97:       crup_lm(frarray,d1+d2,modarray,frarray_index,m,lm_mod,&r);
        !            98:       uptolmup(r,nr);
        !            99:       return;
        !           100:     }
        !           101:   }
        !           102:   error("fft_mulup : FFT_primes exhausted");
        !           103: }
        !           104:
        !           105: void fft_squareup_lm(UP n1,UP *nr)
        !           106: {
        !           107:   ModNum *f1,*w,*fr;
        !           108:   ModNum *frarray[1024];
        !           109:   int modarray[1024];
        !           110:   int frarray_index = 0;
        !           111:   Z m,m1,m2,lm_mod;
        !           112:   int d1,dmin,i,mod,root,d,cond,bound;
        !           113:   UP r;
        !           114:
        !           115:   if ( !n1 ) {
        !           116:     *nr = 0; return;
        !           117:   }
        !           118:   d1 = n1->d; dmin = d1;
        !           119:   if ( !d1 ) {
        !           120:     mulup(n1,n1,nr); return;
        !           121:   }
        !           122:   getmod_lm(&lm_mod);
        !           123:   if ( !lm_mod )
        !           124:     error("fft_squareup_lm : current_mod_lm is not set");
        !           125:   m = ONE;
        !           126:   bound = 2*maxblenup(n1)+int_bits(d1)+2;
        !           127:   f1 = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));
        !           128:   w = (ModNum *)ALLOCA(6*(1<<int_bits(2*d1+1))*sizeof(ModNum));
        !           129:   for ( i = 0; i < NPrimes; i++ ) {
        !           130:     FFT_primes(i,&mod,&root,&d);
        !           131:     if ( (1<<d) < 2*d1+1 )
        !           132:       continue;
        !           133:     modarray[frarray_index] = mod;
        !           134:     frarray[frarray_index++] = fr
        !           135:       = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));
        !           136:     uptofmarray(mod,n1,f1);
        !           137:     cond = FFT_pol_square(d1,f1,fr,i,w);
        !           138:     if ( cond )
        !           139:       error("fft_mulup : error in FFT_pol_product");
        !           140:     STOQ(mod,m1); mulz(m,m1,&m2); m = m2;
        !           141:     if ( z_bits((Q)m) > bound ) {
        !           142:       crup_lm(frarray,2*d1,modarray,frarray_index,m,lm_mod,&r);
        !           143:       uptolmup(r,nr);
        !           144:       return;
        !           145:     }
        !           146:   }
        !           147:   error("fft_squareup : FFT_primes exhausted");
        !           148: }
        !           149:
        !           150: void trunc_fft_mulup_lm(UP n1,UP n2,int dbd,UP *nr)
        !           151: {
        !           152:   ModNum *f1,*f2,*fr,*w;
        !           153:   ModNum *frarray[1024];
        !           154:   int modarray[1024];
        !           155:   int frarray_index = 0;
        !           156:   Z m,m1,m2,lm_mod;
        !           157:   int d1,d2,dmin,i,mod,root,d,cond,bound;
        !           158:   UP r;
        !           159:
        !           160:   if ( !n1 || !n2 ) {
        !           161:     *nr = 0; return;
        !           162:   }
        !           163:   d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
        !           164:   if ( !d1 || !d2 ) {
        !           165:     tmulup(n1,n2,dbd,nr); return;
        !           166:   }
        !           167:   getmod_lm(&lm_mod);
        !           168:   if ( !lm_mod )
        !           169:     error("trunc_fft_mulup_lm : current_mod_lm is not set");
        !           170:   m = ONE;
        !           171:   bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;
        !           172:   f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !           173:   f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !           174:   w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
        !           175:   for ( i = 0; i < NPrimes; i++ ) {
        !           176:     FFT_primes(i,&mod,&root,&d);
        !           177:     if ( (1<<d) < d1+d2+1 )
        !           178:       continue;
        !           179:
        !           180:     modarray[frarray_index] = mod;
        !           181:     frarray[frarray_index++] = fr
        !           182:       = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
        !           183:     uptofmarray(mod,n1,f1);
        !           184:     uptofmarray(mod,n2,f2);
        !           185:     cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
        !           186:     if ( cond )
        !           187:       error("fft_mulup : error in FFT_pol_product");
        !           188:     STOQ(mod,m1); mulz(m,m1,&m2); m = m2;
        !           189:     if ( z_bits((Q)m) > bound ) {
        !           190:       crup_lm(frarray,MIN(dbd-1,d1+d2),modarray,frarray_index,m,lm_mod,&r);
        !           191:       uptolmup(r,nr);
        !           192:       return;
        !           193:     }
        !           194:   }
        !           195:   error("trunc_fft_mulup : FFT_primes exhausted");
        !           196: }
        !           197:
        !           198: void crup_lm(ModNum **f,int d,int *mod,int index,Z m,Z lm_mod,UP *r)
        !           199: {
        !           200:   UP w;
        !           201:   Z t;
        !           202:   int i;
        !           203:
        !           204:   crup(f,d,mod,index,m,&w);
        !           205:   for ( i = 0; i <= d; i++ ) {
        !           206:     remz((Z)w->c[i],lm_mod,&t); w->c[i] = (Num)t;
        !           207:   }
        !           208:   for ( i = d; (i >= 0) && (w->c[i] != 0); i-- );
        !           209:   if ( i < 0 ) *r = 0;
        !           210:   else {
        !           211:     w->d = i;
        !           212:     *r = w;
        !           213:   }
        !           214: }
        !           215:
        !           216: void fft_rembymulup_special_lm(UP n1,UP n2,UP inv2,UP *nr)
        !           217: {
        !           218:   int d1,d2,d;
        !           219:   UP r1,t,s,q,u;
        !           220:
        !           221:   if ( !n2 )
        !           222:     error("rembymulup : division by 0");
        !           223:   else if ( !n1 || !n2->d )
        !           224:     *nr = 0;
        !           225:   else if ( (d1 = n1->d) < (d2 = n2->d) )
        !           226:     *nr = n1;
        !           227:   else {
        !           228:     d = d1-d2;
        !           229:     reverseup(n1,d1,&t); truncup(t,d+1,&r1);
        !           230:     trunc_fft_mulup_lm(r1,inv2,d+1,&t);
        !           231:     truncup(t,d+1,&s);
        !           232:     reverseup(s,d,&q);
        !           233:     trunc_fft_mulup_lm(q,n2,d2,&t); truncup(t,d2,&u);
        !           234:     truncup(n1,d2,&s);
        !           235:     subup(s,u,nr);
        !           236:   }
        !           237: }
        !           238:
        !           239: void uptolmup(UP n,UP *nr)
        !           240: {
        !           241:   int i,d;
        !           242:   Q *c;
        !           243:   LM *cr;
        !           244:   UP r;
        !           245:
        !           246:   if ( !n )
        !           247:     *nr = 0;
        !           248:   else {
        !           249:     d = n->d; c = (Q *)n->c;
        !           250:     *nr = r = UPALLOC(d); cr = (LM *)r->c;
        !           251:     for ( i = 0; i <= d; i++ )
        !           252:       qtolm(c[i],&cr[i]);
        !           253:     for ( i = d; i >= 0 && !cr[i]; i-- );
        !           254:     if ( i < 0 )
        !           255:       *nr = 0;
        !           256:     else
        !           257:       r->d = i;
        !           258:   }
        !           259: }
        !           260:
        !           261: void save_up(UP obj,char *name)
        !           262: {
        !           263:   P p;
        !           264:   Obj ret;
        !           265:   NODE n0,n1;
        !           266:   STRING s;
        !           267:   void Pbsave();
        !           268:
        !           269:   uptop(obj,&p);
        !           270:   MKSTR(s,name);
        !           271:   MKNODE(n1,s,0); MKNODE(n0,p,n1);
        !           272:   Pbsave(n0,&ret);
        !           273: }
        !           274:
        !           275: void hybrid_powermodup(UP f,UP *xp)
        !           276: {
        !           277:   Z n;
        !           278:   UP x,y,t,invf,s;
        !           279:   int k;
        !           280:   LM lm;
        !           281:
        !           282:   getmod_lm(&n);
        !           283:   if ( !n )
        !           284:     error("hybrid_powermodup : current_mod_lm is not set");
        !           285:   lm = ONELM;
        !           286:   x = UPALLOC(1); x->d = 1; x->c[1] = (Num)lm;
        !           287:   y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
        !           288:
        !           289:   reverseup(f,f->d,&t);
        !           290:   invmodup(t,f->d,&s); uptolmup(s,&invf);
        !           291:   for ( k = z_bits((Q)n)-1; k >= 0; k-- ) {
        !           292:     hybrid_squareup(FF_GFP,y,&t);
        !           293:     hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);
        !           294:     y = s;
        !           295:     if ( tstbitz(n,k) ) {
        !           296:       mulup(y,x,&t);
        !           297:       remup(t,f,&s);
        !           298:       y = s;
        !           299:     }
        !           300:   }
        !           301:   *xp = y;
        !           302: }
        !           303:
        !           304: void powermodup(UP f,UP *xp)
        !           305: {
        !           306:   Z n;
        !           307:   UP x,y,t,invf,s;
        !           308:   int k;
        !           309:   Num c;
        !           310:
        !           311:   if ( !current_ff )
        !           312:     error("powermodup : current_ff is not set");
        !           313:   field_order_ff(&n);
        !           314:   one_ff(&c);
        !           315:   x = UPALLOC(1); x->d = 1; x->c[1] = c;
        !           316:   y = UPALLOC(0); y->d = 0; y->c[0] = c;
        !           317:
        !           318:   reverseup(f,f->d,&t);
        !           319:   invmodup(t,f->d,&s);
        !           320:   switch ( current_ff ) {
        !           321:     case FF_GFP:
        !           322:     case FF_GFPN:
        !           323:       uptolmup(s,&invf);
        !           324:       break;
        !           325:     case FF_GFS:
        !           326:     case FF_GFSN:
        !           327:       invf = s; /* XXX */
        !           328:       break;
        !           329:     default:
        !           330:       error("powermodup : not implemented yet");
        !           331:   }
        !           332:   for ( k = z_bits((Q)n)-1; k >= 0; k-- ) {
        !           333:     ksquareup(y,&t);
        !           334:     rembymulup_special(t,f,invf,&s);
        !           335:     y = s;
        !           336:     if ( tstbitz(n,k) ) {
        !           337:       mulup(y,x,&t);
        !           338:       remup(t,f,&s);
        !           339:       y = s;
        !           340:     }
        !           341:   }
        !           342:   *xp = y;
        !           343: }
        !           344:
        !           345: /* g^d mod f */
        !           346:
        !           347: void hybrid_generic_powermodup(UP g,UP f,Z e,UP *xp)
        !           348: {
        !           349:   UP x,y,t,invf,s;
        !           350:   int k;
        !           351:   LM lm;
        !           352:
        !           353:   lm = ONELM;
        !           354:   y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
        !           355:   remup(g,f,&x);
        !           356:   if ( !x ) {
        !           357:     *xp = !e ? y : 0;
        !           358:     return;
        !           359:   } else if ( !x->d ) {
        !           360:     pwrup(x,e,xp);
        !           361:     return;
        !           362:   }
        !           363:   reverseup(f,f->d,&t);
        !           364:   invmodup(t,f->d,&invf);
        !           365:   for ( k = z_bits((Q)e)-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 ( tstbitz(e,k) ) {
        !           370:       mulup(y,x,&t);
        !           371:       remup(t,f,&s);
        !           372:       y = s;
        !           373:     }
        !           374:   }
        !           375:   *xp = y;
        !           376: }
        !           377:
        !           378: void generic_powermodup(UP g,UP f,Z e,UP *xp)
        !           379: {
        !           380:   UP x,y,t,invf,s;
        !           381:   int k;
        !           382:   Num c;
        !           383:
        !           384:   one_ff(&c);
        !           385:   y = UPALLOC(0); y->d = 0; y->c[0] = c;
        !           386:   remup(g,f,&x);
        !           387:   if ( !x ) {
        !           388:     *xp = !e ? y : 0;
        !           389:     return;
        !           390:   } else if ( !x->d ) {
        !           391:     pwrup(x,e,xp);
        !           392:     return;
        !           393:   }
        !           394:   reverseup(f,f->d,&t);
        !           395:   invmodup(t,f->d,&invf);
        !           396:   for ( k = z_bits((Q)e)-1; k >= 0; k-- ) {
        !           397:     ksquareup(y,&t);
        !           398:     rembymulup_special(t,f,invf,&s);
        !           399:     y = s;
        !           400:     if ( tstbitz(e,k) ) {
        !           401:       mulup(y,x,&t);
        !           402:       remup(t,f,&s);
        !           403:       y = s;
        !           404:     }
        !           405:   }
        !           406:   *xp = y;
        !           407: }
        !           408:
        !           409: void hybrid_powertabup(UP f,UP xp,UP *tab)
        !           410: {
        !           411:   UP y,t,invf;
        !           412:   int i,d;
        !           413:   LM lm;
        !           414:
        !           415:   d = f->d;
        !           416:   lm = ONELM;
        !           417:   y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
        !           418:   tab[0] = y;
        !           419:   tab[1] = xp;
        !           420:
        !           421:   reverseup(f,f->d,&t);
        !           422:   invmodup(t,f->d,&invf);
        !           423:
        !           424:   for ( i = 2; i < d; i++ ) {
        !           425:     if ( debug_up ){
        !           426:       fprintf(stderr,".");
        !           427:     }
        !           428:     hybrid_mulup(FF_GFP,tab[i-1],xp,&t);
        !           429:     hybrid_rembymulup_special(FF_GFP,t,f,invf,&tab[i]);
        !           430:   }
        !           431: }
        !           432:
        !           433: void powertabup(UP f,UP xp,UP *tab)
        !           434: {
        !           435:   UP y,t,invf;
        !           436:   int i,d;
        !           437:   Num c;
        !           438:
        !           439:   d = f->d;
        !           440:   one_ff(&c);
        !           441:   y = UPALLOC(0); y->d = 0; y->c[0] = c;
        !           442:   tab[0] = y;
        !           443:   tab[1] = xp;
        !           444:
        !           445:   reverseup(f,f->d,&t);
        !           446:   invmodup(t,f->d,&invf);
        !           447:
        !           448:   for ( i = 2; i < d; i++ ) {
        !           449:     if ( debug_up ){
        !           450:       fprintf(stderr,".");
        !           451:     }
        !           452:     kmulup(tab[i-1],xp,&t);
        !           453:     rembymulup_special(t,f,invf,&tab[i]);
        !           454:   }
        !           455: }

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