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