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