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