Annotation of OpenXM_contrib2/asir2018/engine/up_lm.c, Revision 1.3
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.3 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/engine/up_lm.c,v 1.2 2018/09/28 08:20:28 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;
1.3 ! noro 61: ModNum *frarray[NPrimes];
! 62: int modarray[NPrimes];
1.1 noro 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;
1.3 ! noro 108: ModNum *frarray[NPrimes];
! 109: int modarray[NPrimes];
1.1 noro 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;
1.3 ! noro 153: ModNum *frarray[NPrimes];
! 154: int modarray[NPrimes];
1.1 noro 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: }
1.3 ! noro 208: for ( i = d; (i >= 0) && (w->c[i] == 0); i-- );
1.1 noro 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;
1.3 ! noro 278: UP x,y,t,invf,s,s1;
1.1 noro 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