[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.2

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:  *
1.2     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2018/engine/up_lm.c,v 1.1 2018/09/19 05:45:07 noro Exp $
1.1       noro       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");
1.2     ! noro       95:     STOZ(mod,m1); mulz(m,m1,&m2); m = m2;
1.1       noro       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");
1.2     ! noro      140:     STOZ(mod,m1); mulz(m,m1,&m2); m = m2;
1.1       noro      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");
1.2     ! noro      188:     STOZ(mod,m1); mulz(m,m1,&m2); m = m2;
1.1       noro      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>