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>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to up_lm.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine