![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to up.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.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.3 ! 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@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: OpenXM_contrib2/asir2000/engine/up.c,v 1.2 1999/12/27 04:16:31 noro Exp $
! 49: */
1.1 noro 50: #include "ca.h"
51: #include <math.h>
52:
53: struct oEGT eg_chrem,eg_fore,eg_back;
54:
55: /*
56: int up_kara_mag=15;
57: int up_tkara_mag=15;
58: */
59: /*
60: int up_kara_mag=30;
61: int up_tkara_mag=20;
62: */
63: int up_kara_mag=20;
64: int up_tkara_mag=15;
65:
66: #if defined(sparc)
67: int up_fft_mag=50;
68: #else
69: int up_fft_mag=80;
70: #endif
71:
72: int up_rem_mag=1;
73: int debug_up;
74: int up_lazy;
75: extern int lm_lazy;
76: extern int current_ff;
77: extern int GC_dont_gc;
78:
79: void monicup(a,b)
80: UP a;
81: UP *b;
82: {
83: UP w;
84:
85: if ( !a )
86: *b = 0;
87: else {
88: w = W_UPALLOC(0); w->d = 0;
89: divnum(0,(Num)ONE,(Num)a->c[a->d],(Num *)&w->c[0]);
90: mulup(a,w,b);
91: }
92: }
93:
94: void simpup(a,b)
95: UP a;
96: UP *b;
97: {
98: int i,d;
99: UP c;
100:
101: if ( !a )
102: *b = 0;
103: else {
104: d = a->d;
105: c = UPALLOC(d);
106:
107: for ( i = 0; i <= d; i++ )
108: simpnum(a->c[i],&c->c[i]);
109: for ( i = d; i >= 0 && !c->c[i]; i-- );
110: if ( i < 0 )
111: *b = 0;
112: else {
113: c->d = i;
114: *b = c;
115: }
116: }
117: }
118:
119: void simpnum(a,b)
120: Num a;
121: Num *b;
122: {
123: LM lm;
124: GF2N gf;
125: GFPN gfpn;
126:
127: if ( !a )
128: *b = 0;
129: else
130: switch ( NID(a) ) {
131: case N_LM:
132: simplm((LM)a,&lm); *b = (Num)lm; break;
133: case N_GF2N:
134: simpgf2n((GF2N)a,&gf); *b = (Num)gf; break;
135: case N_GFPN:
136: simpgfpn((GFPN)a,&gfpn); *b = (Num)gfpn; break;
137: default:
138: *b = a; break;
139: }
140: }
141:
142: void uremp(p1,p2,rp)
143: P p1,p2;
144: P *rp;
145: {
146: UP n1,n2,r;
147:
148: if ( !p1 || NUM(p1) )
149: *rp = p1;
150: else {
151: ptoup(p1,&n1); ptoup(p2,&n2);
152: remup(n1,n2,&r);
153: uptop(r,rp);
154: }
155: }
156:
157: void ugcdp(p1,p2,rp)
158: P p1,p2;
159: P *rp;
160: {
161: UP n1,n2,r;
162:
163: ptoup(p1,&n1); ptoup(p2,&n2);
164: gcdup(n1,n2,&r);
165: uptop(r,rp);
166: }
167:
168: void reversep(p1,d,rp)
169: P p1;
170: Q d;
171: P *rp;
172: {
173: UP n1,r;
174:
175: if ( !p1 )
176: *rp = 0;
177: else {
178: ptoup(p1,&n1);
179: reverseup(n1,QTOS(d),&r);
180: uptop(r,rp);
181: }
182: }
183:
184: void invmodp(p1,d,rp)
185: P p1;
186: Q d;
187: P *rp;
188: {
189: UP n1,r;
190:
191: if ( !p1 )
192: error("invmodp : invalid input");
193: else {
194: ptoup(p1,&n1);
195: invmodup(n1,QTOS(d),&r);
196: uptop(r,rp);
197: }
198: }
199:
200: void addup(n1,n2,nr)
201: UP n1,n2;
202: UP *nr;
203: {
204: UP r,t;
205: int i,d1,d2;
206: Num *c,*c1,*c2;
207:
208: if ( !n1 ) {
209: *nr = n2;
210: } else if ( !n2 ) {
211: *nr = n1;
212: } else {
213: if ( n2->d > n1->d ) {
214: t = n1; n1 = n2; n2 = t;
215: }
216: d1 = n1->d; d2 = n2->d;
217: *nr = r = UPALLOC(d1);
218: c1 = n1->c; c2 = n2->c; c = r->c;
219: for ( i = 0; i <= d2; i++ )
220: addnum(0,*c1++,*c2++,c++);
221: for ( ; i <= d1; i++ )
222: *c++ = *c1++;
223: for ( i = d1, c = r->c; i >=0 && !c[i]; i-- );
224: if ( i < 0 )
225: *nr = 0;
226: else
227: r->d = i;
228: }
229: }
230:
231: void subup(n1,n2,nr)
232: UP n1,n2;
233: UP *nr;
234: {
235: UP r;
236: int i,d1,d2,d;
237: Num *c,*c1,*c2;
238:
239: if ( !n1 ) {
240: chsgnup(n2,nr);
241: } else if ( !n2 ) {
242: *nr = n1;
243: } else {
244: d1 = n1->d; d2 = n2->d; d = MAX(d1,d2);
245: *nr = r = UPALLOC(d);
246: c1 = n1->c; c2 = n2->c; c = r->c;
247: if ( d1 >= d2 ) {
248: for ( i = 0; i <= d2; i++ )
249: subnum(0,*c1++,*c2++,c++);
250: for ( ; i <= d1; i++ )
251: *c++ = *c1++;
252: } else {
253: for ( i = 0; i <= d1; i++ )
254: subnum(0,*c1++,*c2++,c++);
255: for ( ; i <= d2; i++ )
256: chsgnnum(*c2++,c++);
257: }
258: for ( i = d, c = r->c; i >=0 && !c[i]; i-- );
259: if ( i < 0 )
260: *nr = 0;
261: else
262: r->d = i;
263: }
264: }
265:
266: void chsgnup(n1,nr)
267: UP n1;
268: UP *nr;
269: {
270: UP r;
271: int d1,i;
272: Num *c1,*c;
273:
274: if ( !n1 ) {
275: *nr = 0;
276: } else {
277: d1 = n1->d;
278: *nr = r = UPALLOC(d1); r->d = d1;
279: c1 = n1->c; c = r->c;
280: for ( i = 0; i <= d1; i++ )
281: chsgnnum(*c1++,c++);
282: }
283: }
284:
285: void hybrid_mulup(ff,n1,n2,nr)
286: int ff;
287: UP n1,n2;
288: UP *nr;
289: {
290: if ( !n1 || !n2 )
291: *nr = 0;
292: else if ( MAX(n1->d,n2->d) < up_fft_mag )
293: kmulup(n1,n2,nr);
294: else
295: switch ( ff ) {
296: case FF_GFP:
297: fft_mulup_lm(n1,n2,nr); break;
298: case FF_GF2N:
299: kmulup(n1,n2,nr); break;
300: default:
301: fft_mulup(n1,n2,nr); break;
302: }
303: }
304:
305: void hybrid_squareup(ff,n1,nr)
306: int ff;
307: UP n1;
308: UP *nr;
309: {
310: if ( !n1 )
311: *nr = 0;
312: else if ( n1->d < up_fft_mag )
313: ksquareup(n1,nr);
314: else
315: switch ( ff ) {
316: case FF_GFP:
317: fft_squareup_lm(n1,nr); break;
318: case FF_GF2N:
319: ksquareup(n1,nr); break;
320: default:
321: fft_squareup(n1,nr); break;
322: }
323: }
324:
325: void hybrid_tmulup(ff,n1,n2,d,nr)
326: int ff;
327: UP n1,n2;
328: int d;
329: UP *nr;
330: {
331: if ( !n1 || !n2 )
332: *nr = 0;
333: else if ( MAX(n1->d,n2->d) < up_fft_mag )
334: tkmulup(n1,n2,d,nr);
335: else
336: switch ( ff ) {
337: case FF_GFP:
338: trunc_fft_mulup_lm(n1,n2,d,nr); break;
339: case FF_GF2N:
340: tkmulup(n1,n2,d,nr); break;
341: default:
342: trunc_fft_mulup(n1,n2,d,nr); break;
343: }
344: }
345:
346: void mulup(n1,n2,nr)
347: UP n1,n2;
348: UP *nr;
349: {
350: UP r;
351: Num *pc1,*pc,*c1,*c2,*c;
352: Num mul,t,s;
353: int i,j,d1,d2;
354:
355: if ( !n1 || !n2 )
356: *nr = 0;
357: else {
358: d1 = n1->d; d2 = n2->d;
359: *nr = r = UPALLOC(d1+d2); r->d = d1+d2;
360: c1 = n1->c; c2 = n2->c; c = r->c;
361: lm_lazy = 1;
362: for ( i = 0; i <= d2; i++, c++ )
363: if ( mul = *c2++ )
364: for ( j = 0, pc1 = c1, pc = c; j <= d1; j++ ) {
365: mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;
366: }
367: lm_lazy = 0;
368: if ( !up_lazy )
369: for ( i = 0, c = r->c; i <= r->d; i++ ) {
370: simpnum(c[i],&t); c[i] = t;
371: }
372: }
373: }
374:
375: /* nr = c*n1 */
376:
377: void mulcup(c,n1,nr)
378: Num c;
379: UP n1;
380: UP *nr;
381: {
382: int d;
383: UP r;
384: Num t;
385: Num *c1,*cr;
386: int i;
387:
388: if ( !c || !n1 )
389: *nr = 0;
390: else {
391: d = n1->d;
392: *nr = r = UPALLOC(d); r->d = d;
393: c1 = n1->c; cr = r->c;
394: lm_lazy = 1;
395: for ( i = 0; i <= d; i++ )
396: mulnum(0,c1[i],c,&cr[i]);
397: lm_lazy = 0;
398: if ( !up_lazy )
399: for ( i = 0; i <= d; i++ ) {
400: simpnum(cr[i],&t); cr[i] = t;
401: }
402: }
403: }
404:
405: void tmulup(n1,n2,d,nr)
406: UP n1,n2;
407: int d;
408: UP *nr;
409: {
410: UP r;
411: Num *pc1,*pc,*c1,*c2,*c;
412: Num mul,t,s;
413: int i,j,d1,d2,dr;
414:
415: if ( !n1 || !n2 )
416: *nr = 0;
417: else {
418: d1 = n1->d; d2 = n2->d;
419: dr = MAX(d-1,d1+d2);
420: *nr = r = UPALLOC(dr);
421: c1 = n1->c; c2 = n2->c; c = r->c;
422: lm_lazy = 1;
423: for ( i = 0; i <= d2; i++, c++ )
424: if ( mul = *c2++ )
425: for ( j = 0, pc1 = c1, pc = c; (j <= d1) && (i+j < d); j++ ) {
426: mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;
427: }
428: lm_lazy = 0;
429: c = r->c;
430: if ( !up_lazy )
431: for ( i = 0; i < d; i++ ) {
432: simpnum(c[i],&t); c[i] = t;
433: }
434: for ( i = d-1; i >= 0 && !c[i]; i-- );
435: if ( i < 0 )
436: *nr = 0;
437: else
438: r->d = i;
439: }
440: }
441:
442: void squareup(n1,nr)
443: UP n1;
444: UP *nr;
445: {
446: UP r;
447: Num *c1,*c;
448: Num t,s,u;
449: int i,j,h,d1,d;
450:
451: if ( !n1 )
452: *nr = 0;
453: else if ( !n1->d ) {
454: *nr = r = UPALLOC(0); r->d = 0;
455: mulnum(0,n1->c[0],n1->c[0],&r->c[0]);
456: } else {
457: d1 = n1->d;
458: d = 2*d1;
459: *nr = r = UPALLOC(2*d); r->d = d;
460: c1 = n1->c; c = r->c;
461: lm_lazy = 1;
462: for ( i = 0; i <= d1; i++ )
463: mulnum(0,c1[i],c1[i],&c[2*i]);
464: for ( j = 1; j < d; j++ ) {
465: h = (j-1)/2;
466: for ( i = MAX(j-d1,0), s = 0; i <= h; i++ ) {
467: mulnum(0,c1[i],c1[j-i],&t); addnum(0,t,s,&u); s = u;
468: }
469: addnum(0,s,s,&t); addnum(0,t,c[j],&u); c[j] = u;
470: }
471: lm_lazy = 0;
472: if ( !up_lazy )
473: for ( i = 0, c = r->c; i <= d; i++ ) {
474: simpnum(c[i],&t); c[i] = t;
475: }
476: }
477: }
478:
479: void remup(n1,n2,nr)
480: UP n1,n2;
481: UP *nr;
482: {
483: UP w,r;
484:
485: if ( !n2 )
486: error("remup : division by 0");
487: else if ( !n1 || !n2->d )
488: *nr = 0;
489: else if ( n1->d < n2->d )
490: *nr = n1;
491: else {
492: w = W_UPALLOC(n1->d); copyup(n1,w);
493: remup_destructive(w,n2);
494: if ( w->d < 0 )
495: *nr = 0;
496: else {
497: *nr = r = UPALLOC(w->d); copyup(w,r);
498: }
499: }
500: }
501:
502: void remup_destructive(n1,n2)
503: UP n1,n2;
504: {
505: Num *c1,*c2;
506: int d1,d2,i,j;
507: Num m,t,s,mhc;
508:
509: c1 = n1->c; c2 = n2->c;
510: d1 = n1->d; d2 = n2->d;
511: chsgnnum(c2[d2],&mhc);
512: for ( i = d1; i >= d2; i-- ) {
513: simpnum(c1[i],&t); c1[i] = t;
514: if ( !c1[i] )
515: continue;
516: divnum(0,c1[i],mhc,&m);
517: lm_lazy = 1;
518: for ( j = d2; j >= 0; j-- ) {
519: mulnum(0,c2[j],m,&t); addnum(0,c1[i+j-d2],t,&s);
520: c1[i+j-d2] = s;
521: }
522: lm_lazy = 0;
523: }
524: for ( i = 0; i < d2; i++ ) {
525: simpnum(c1[i],&t); c1[i] = t;
526: }
527: for ( i = d2-1; i >= 0 && !c1[i]; i-- );
528: n1->d = i;
529: }
530:
531: void qrup(n1,n2,nq,nr)
532: UP n1,n2;
533: UP *nq,*nr;
534: {
535: UP w,r,q;
536: struct oUP t;
537: Num m;
538: int d;
539:
540: if ( !n2 )
541: error("qrup : division by 0");
542: else if ( !n1 ) {
543: *nq = 0; *nr = 0;
544: } else if ( !n2->d )
545: if ( !n1->d ) {
546: divnum(0,n1->c[0],n2->c[0],(Num *)nq); *nr = 0;
547: } else {
548: divnum(0,(Num)ONE,n2->c[0],&m);
549: t.d = 0; t.c[0] = m;
550: mulup(n1,&t,nq); *nr = 0;
551: }
552: else if ( n1->d < n2->d ) {
553: *nq = 0; *nr = n1;
554: } else {
555: w = W_UPALLOC(n1->d); copyup(n1,w);
556: qrup_destructive(w,n2);
557: d = n1->d-n2->d;
558: *nq = q = UPALLOC(d); q->d = d;
559: bcopy(w->c+n2->d,q->c,(d+1)*sizeof(Num));
560: if ( w->d < 0 )
561: *nr = 0;
562: else {
563: *nr = r = UPALLOC(w->d); copyup(w,r);
564: }
565: }
566: }
567:
568: void qrup_destructive(n1,n2)
569: UP n1,n2;
570: {
571: Num *c1,*c2;
572: int d1,d2,i,j;
573: Num q,t,s,hc;
574:
575: c1 = n1->c; c2 = n2->c;
576: d1 = n1->d; d2 = n2->d;
577: hc = c2[d2];
578: for ( i = d1; i >= d2; i-- ) {
579: simpnum(c1[i],&t); c1[i] = t;
580: if ( c1[i] ) {
581: divnum(0,c1[i],hc,&q);
582: lm_lazy = 1;
583: for ( j = d2; j >= 0; j-- ) {
584: mulnum(0,c2[j],q,&t); subnum(0,c1[i+j-d2],t,&s);
585: c1[i+j-d2] = s;
586: }
587: lm_lazy = 0;
588: } else
589: q = 0;
590: c1[i] = q;
591: }
592: for ( i = 0; i < d2; i++ ) {
593: simpnum(c1[i],&t); c1[i] = t;
594: }
595: for ( i = d2-1; i >= 0 && !c1[i]; i-- );
596: n1->d = i;
597: }
598:
599: void gcdup(n1,n2,nr)
600: UP n1,n2;
601: UP *nr;
602: {
603: UP w1,w2,t,r;
604: int d1,d2;
605:
606: if ( !n1 )
607: *nr = n2;
608: else if ( !n2 )
609: *nr = n1;
610: else if ( !n1->d || !n2->d ) {
611: *nr = r = UPALLOC(0); r->d = 0;
612: divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]);
613: } else {
614: d1 = n1->d; d2 = n2->d;
615: if ( d2 > d1 ) {
616: w1 = W_UPALLOC(d2); copyup(n2,w1);
617: w2 = W_UPALLOC(d1); copyup(n1,w2);
618: } else {
619: w1 = W_UPALLOC(d1); copyup(n1,w1);
620: w2 = W_UPALLOC(d2); copyup(n2,w2);
621: }
622: d1 = w1->d; d2 = w2->d;
623: while ( 1 ) {
624: remup_destructive(w1,w2);
625: if ( w1->d < 0 ) {
626: *nr = r = UPALLOC(w2->d); copyup(w2,r); return;
627: } else if ( !w1->d ) {
628: *nr = r = UPALLOC(0); r->d = 0;
629: divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]); return;
630: } else {
631: t = w1; w1 = w2; w2 = t;
632: }
633: }
634: }
635: }
636:
637: /* compute r s.t. a*r = 1 mod m */
638:
639: void extended_gcdup(a,m,r)
640: UP a,m;
641: UP *r;
642: {
643: UP one,g1,g2,a1,a2,a3,b1,b2,b3,inv,quo,rem,t;
644: Num i;
645:
646: one = UPALLOC(0); one->d = 0; one->c[0] = (Num)ONELM;
647: g1 = m; g2 = a;
648: a1 = one; a2 = 0;
649: b1 = 0; b2 = one;
650: /* a2*m+b2*a = g2 always holds */
651: while ( g2->d ) {
652: qrup(g1,g2,&quo,&rem); g1 = g2; g2 = rem;
653: mulup(a2,quo,&t); subup(a1,t,&a3); a1 = a2; a2 = a3;
654: mulup(b2,quo,&t); subup(b1,t,&b3); b1 = b2; b2 = b3;
655: }
656: /* now a2*m+b2*a = g2 (const) */
657: divnum(0,(Num)ONE,g2->c[0],&i);
658: inv = UPALLOC(0); inv->d = 0; inv->c[0] = i;
659: mulup(b2,inv,r);
660: }
661:
662: void reverseup(n1,d,nr)
663: UP n1;
664: int d;
665: UP *nr;
666: {
667: Num *c1,*c;
668: int i,d1;
669: UP r;
670:
671: if ( !n1 )
672: *nr = 0;
673: else if ( d < (d1 = n1->d) )
674: error("reverseup : invalid input");
675: else {
676: *nr = r = UPALLOC(d);
677: c = r->c;
678: bzero((char *)c,(d+1)*sizeof(Num));
679: c1 = n1->c;
680: for ( i = 0; i <= d1; i++ )
681: c[d-i] = c1[i];
682: for ( i = d; i >= 0 && !c[i]; i-- );
683: r->d = i;
684: if ( i < 0 )
685: *nr = 0;
686: }
687: }
688:
689: void invmodup(n1,d,nr)
690: UP n1;
691: int d;
692: UP *nr;
693: {
694: UP r;
695: Num s,t,u,hinv;
696: int i,j,dmin;
697: Num *w,*c,*cr;
698:
699: if ( !n1 || !n1->c[0] )
700: error("invmodup : invalid input");
701: else if ( !n1->d ) {
702: *nr = r = UPALLOC(0); r->d = 0;
703: divnum(0,(Num)ONE,n1->c[0],&r->c[0]);
704: } else {
705: *nr = r = UPALLOC(d);
706: w = (Num *)ALLOCA((d+1)*sizeof(Q));
707: bzero((char *)w,(d+1)*sizeof(Q));
708: dmin = MIN(d,n1->d);
709: divnum(0,(Num)ONE,n1->c[0],&hinv);
710: for ( i = 0, c = n1->c; i <= dmin; i++ )
711: mulnum(0,c[i],hinv,&w[i]);
712: cr = r->c;
713: cr[0] = w[0];
714: for ( i = 1; i <= d; i++ ) {
715: for ( j = 1, s = w[i]; j < i; j++ ) {
716: mulnum(0,cr[j],w[i-j],&t);
717: addnum(0,s,t,&u);
718: s = u;
719: }
720: chsgnnum(s,&cr[i]);
721: }
722: for ( i = 0; i <= d; i++ ) {
723: mulnum(0,cr[i],hinv,&t); cr[i] = t;
724: }
725: for ( i = d; i >= 0 && !cr[i]; i-- );
726: r->d = i;
727: }
728: }
729:
730: void pwrup(n,e,nr)
731: UP n;
732: Q e;
733: UP *nr;
734: {
735: UP y,z,t;
736: N en,qn;
737: int r;
738:
739: y = UPALLOC(0); y->d = 0; y->c[0] = (Num)ONE;
740: z = n;
741: if ( !e )
742: *nr = y;
743: else if ( UNIQ(e) )
744: *nr = n;
745: else {
746: en = NM(e);
747: while ( 1 ) {
748: r = divin(en,2,&qn); en = qn;
749: if ( r ) {
750: mulup(z,y,&t); y = t;
751: if ( !en ) {
752: *nr = y;
753: return;
754: }
755: }
756: mulup(z,z,&t); z = t;
757: }
758: }
759: }
760:
761: int compup(n1,n2)
762: UP n1,n2;
763: {
764: int i,r;
765:
766: if ( !n1 )
767: return n2 ? -1 : 0;
768: else if ( !n2 )
769: return 1;
770: else if ( n1->d > n2->d )
771: return 1;
772: else if ( n1->d < n2->d )
773: return -1;
774: else {
775: for ( i = n1->d; i >= 0; i-- ) {
776: r = compnum(CO,n1->c[i],n2->c[i]);
777: if ( r )
778: return r;
779: }
780: return 0;
781: }
782: }
783:
784: void kmulp(vl,n1,n2,nr)
785: VL vl;
786: P n1,n2;
787: P *nr;
788: {
789: UP b1,b2,br;
790:
791: if ( !n1 || !n2 )
792: *nr = 0;
793: else if ( OID(n1) == O_N || OID(n2) == O_N )
794: mulp(vl,n1,n2,nr);
795: else {
796: ptoup(n1,&b1); ptoup(n2,&b2);
797: kmulup(b1,b2,&br);
798: uptop(br,nr);
799: }
800: }
801:
802: void ksquarep(vl,n1,nr)
803: VL vl;
804: P n1;
805: P *nr;
806: {
807: UP b1,br;
808:
809: if ( !n1 )
810: *nr = 0;
811: else if ( OID(n1) == O_N )
812: mulp(vl,n1,n1,nr);
813: else {
814: ptoup(n1,&b1);
815: ksquareup(b1,&br);
816: uptop(br,nr);
817: }
818: }
819:
820: void kmulup(n1,n2,nr)
821: UP n1,n2,*nr;
822: {
823: UP n,t,s,m,carry;
824: int d,d1,d2,len,i,l;
825: pointer *r,*r0;
826:
827: if ( !n1 || !n2 ) {
828: *nr = 0; return;
829: }
830: d1 = DEG(n1)+1; d2 = DEG(n2)+1;
831: if ( (d1 < up_kara_mag) && (d2 < up_kara_mag) )
832: mulup(n1,n2,nr);
833: else
834: kmulupmain(n1,n2,nr);
835: }
836:
837: void ksquareup(n1,nr)
838: UP n1,*nr;
839: {
840: int d1;
841: extern Q TWO;
842:
843: if ( !n1 ) {
844: *nr = 0; return;
845: }
846: d1 = DEG(n1)+1;
847: if ( (d1 < up_kara_mag) )
848: squareup(n1,nr);
849: else
850: ksquareupmain(n1,nr);
851: }
852:
853: void copyup(n1,n2)
854: UP n1,n2;
855: {
856: n2->d = n1->d;
857: bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(Q));
858: }
859:
860: void c_copyup(n,len,p)
861: UP n;
862: int len;
863: pointer *p;
864: {
865: if ( n )
866: bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(pointer));
867: }
868:
869: void kmulupmain(n1,n2,nr)
870: UP n1,n2,*nr;
871: {
872: int d1,d2,h,len;
873: UP n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
874:
875: d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (MAX(d1,d2)+1)/2;
876: decompup(n1,h,&n1lo,&n1hi);
877: decompup(n2,h,&n2lo,&n2hi);
878: kmulup(n1lo,n2lo,&lo);
879: kmulup(n1hi,n2hi,&hi);
880: shiftup(hi,2*h,&t2);
881: addup(lo,t2,&t1);
882: addup(hi,lo,&mid1);
883: subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2); kmulup(s1,s2,&mid2);
884: addup(mid1,mid2,&mid);
885: shiftup(mid,h,&t2);
886: addup(t1,t2,nr);
887: }
888:
889: void ksquareupmain(n1,nr)
890: UP n1,*nr;
891: {
892: int d1,h,len;
893: UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
894:
895: d1 = DEG(n1)+1; h = (d1+1)/2;
896: decompup(n1,h,&n1lo,&n1hi);
897: ksquareup(n1hi,&hi); ksquareup(n1lo,&lo);
898: shiftup(hi,2*h,&t2);
899: addup(lo,t2,&t1);
900: addup(hi,lo,&mid1);
901: subup(n1hi,n1lo,&s1); ksquareup(s1,&mid2);
902: subup(mid1,mid2,&mid);
903: shiftup(mid,h,&t2);
904: addup(t1,t2,nr);
905: }
906:
907: void rembymulup(n1,n2,nr)
908: UP n1,n2;
909: UP *nr;
910: {
911: int d1,d2,d;
912: UP r1,r2,inv2,t,s,q;
913:
914: if ( !n2 )
915: error("rembymulup : division by 0");
916: else if ( !n1 || !n2->d )
917: *nr = 0;
918: else if ( (d1 = n1->d) < (d2 = n2->d) )
919: *nr = n1;
920: else {
921: d = d1-d2;
922: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
923: reverseup(n2,d2,&r2);
924: invmodup(r2,d,&inv2);
925: kmulup(r1,inv2,&t);
926: truncup(t,d+1,&s);
927: reverseup(s,d,&q);
928: kmulup(q,n2,&t); subup(n1,t,nr);
929: }
930: }
931:
932: /*
933: deg n2 = d
934: deg n1 <= 2*d-1
935: inv2 = inverse of reversep(n2) mod x^d
936: */
937:
938: void hybrid_rembymulup_special(ff,n1,n2,inv2,nr)
939: int ff;
940: UP n1,n2,inv2;
941: UP *nr;
942: {
943: int d1,d2,d;
944: UP r1,t,s,q,u;
945:
946: if ( !n2 )
947: error("hybrid_rembymulup : division by 0");
948: else if ( !n1 || !n2->d )
949: *nr = 0;
950: else if ( (d1 = n1->d) < (d2 = n2->d) )
951: *nr = n1;
952: else {
953: d = d1-d2;
954: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
955: hybrid_tmulup(ff,r1,inv2,d+1,&t);
956: truncup(t,d+1,&s);
957: reverseup(s,d,&q);
958:
959: hybrid_tmulup(ff,q,n2,d2,&t);
960: truncup(n1,d2,&s);
961: subup(s,t,nr);
962: }
963: }
964:
965: void rembymulup_special(n1,n2,inv2,nr)
966: UP n1,n2,inv2;
967: UP *nr;
968: {
969: int d1,d2,d;
970: UP r1,t,s,q,u;
971:
972: if ( !n2 )
973: error("rembymulup : division by 0");
974: else if ( !n1 || !n2->d )
975: *nr = 0;
976: else if ( (d1 = n1->d) < (d2 = n2->d) )
977: *nr = n1;
978: else {
979: d = d1-d2;
980: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
981: tkmulup(r1,inv2,d+1,&t);
982: truncup(t,d+1,&s);
983: reverseup(s,d,&q);
984:
985: tkmulup(q,n2,d2,&t);
986: truncup(n1,d2,&s);
987: subup(s,t,nr);
988: }
989: }
990:
991: /* *nr = n1*n2 mod x^d */
992:
993: void tkmulup(n1,n2,d,nr)
994: UP n1,n2;
995: int d;
996: UP *nr;
997: {
998: int m;
999: UP n1l,n1h,n2l,n2h,l,h,t,s,u,afo;
1000:
1001: if ( d < 0 )
1002: error("tkmulup : invalid argument");
1003: else if ( d == 0 )
1004: *nr = 0;
1005: else {
1006: truncup(n1,d,&t); n1 = t;
1007: truncup(n2,d,&t); n2 = t;
1008: if ( !n1 || !n2 )
1009: *nr = 0;
1010: else if ( n1->d < up_tkara_mag || n2->d < up_tkara_mag )
1011: tmulup(n1,n2,d,nr);
1012: else {
1013: m = (d+1)/2;
1014: decompup(n1,m,&n1l,&n1h);
1015: decompup(n2,m,&n2l,&n2h);
1016: kmulup(n1l,n2l,&l);
1017: tkmulup(n1l,n2h,d-m,&t);
1018: tkmulup(n2l,n1h,d-m,&s);
1019: addup(t,s,&u);
1020: shiftup(u,m,&h);
1021: addup(l,h,&t);
1022: truncup(t,d,nr);
1023: }
1024: }
1025: }
1026:
1027: /* n->n*x^d */
1028:
1029: void shiftup(n,d,nr)
1030: UP n;
1031: int d;
1032: UP *nr;
1033: {
1034: int dr;
1035: UP r;
1036:
1037: if ( !n )
1038: *nr = 0;
1039: else {
1040: dr = n->d+d;
1041: *nr = r = UPALLOC(dr); r->d = dr;
1042: bzero(r->c,(dr+1)*sizeof(Num));
1043: bcopy(n->c,r->c+d,(n->d+1)*sizeof(Num));
1044: }
1045: }
1046:
1047: void fft_rembymulup_special(n1,n2,inv2,nr)
1048: UP n1,n2,inv2;
1049: UP *nr;
1050: {
1051: int d1,d2,d;
1052: UP r1,t,s,q,u;
1053:
1054: if ( !n2 )
1055: error("rembymulup : division by 0");
1056: else if ( !n1 || !n2->d )
1057: *nr = 0;
1058: else if ( (d1 = n1->d) < (d2 = n2->d) )
1059: *nr = n1;
1060: else {
1061: d = d1-d2;
1062: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
1063: trunc_fft_mulup(r1,inv2,d+1,&t);
1064: truncup(t,d+1,&s);
1065: reverseup(s,d,&q);
1066: trunc_fft_mulup(q,n2,d2,&t); truncup(t,d2,&u);
1067: truncup(n1,d2,&s);
1068: subup(s,u,nr);
1069: }
1070: }
1071:
1072: void set_degreeup(n,d)
1073: UP n;
1074: int d;
1075: {
1076: int i;
1077:
1078: for ( i = d; i >= 0 && !n->c[i]; i-- );
1079: n->d = i;
1080: }
1081:
1082: /* n -> n0 + x^d n1 */
1083:
1084: void decompup(n,d,n0,n1)
1085: UP n;
1086: int d;
1087: UP *n0,*n1;
1088: {
1089: int dn;
1090: UP r0,r1;
1091:
1092: if ( !n || d > n->d ) {
1093: *n0 = n; *n1 = 0;
1094: } else if ( d < 0 )
1095: error("decompup : invalid argument");
1096: else if ( d == 0 ) {
1097: *n0 = 0; *n1 = n;
1098: } else {
1099: dn = n->d;
1100: *n0 = r0 = UPALLOC(d-1);
1101: *n1 = r1 = UPALLOC(dn-d);
1102: bcopy(n->c,r0->c,d*sizeof(Num));
1103: bcopy(n->c+d,r1->c,(dn-d+1)*sizeof(Num));
1104: set_degreeup(r0,d-1);
1105: if ( r0->d < 0 )
1106: *n0 = 0;
1107: set_degreeup(r1,dn-d);
1108: if ( r1->d < 0 )
1109: *n1 = 0;
1110: }
1111: }
1112:
1113: /* n -> n mod x^d */
1114:
1115: void truncup(n1,d,nr)
1116: UP n1;
1117: int d;
1118: UP *nr;
1119: {
1120: int i;
1121: UP r;
1122:
1123: if ( !n1 || d > n1->d )
1124: *nr = n1;
1125: else if ( d < 0 )
1126: error("truncup : invalid argument");
1127: else if ( d == 0 )
1128: *nr = 0;
1129: else {
1130: *nr = r = UPALLOC(d-1);
1131: bcopy(n1->c,r->c,d*sizeof(Num));
1132: for ( i = d-1; i >= 0 && !r->c[i]; i-- );
1133: if ( i < 0 )
1134: *nr = 0;
1135: else
1136: r->d = i;
1137: }
1138: }
1139:
1140: int int_bits(t)
1141: int t;
1142: {
1143: int k;
1144:
1145: for ( k = 0; t; t>>=1, k++);
1146: return k;
1147: }
1148:
1149: /* n is assumed to be LM or integer coefficient */
1150:
1151: int maxblenup(n)
1152: UP n;
1153: {
1154: int m,r,i,d;
1155: Num *c;
1156:
1157: if ( !n )
1158: return 0;
1159: d = n->d; c = (Num *)n->c;
1160: for ( r = 0, i = 0; i <= d; i++ )
1161: if ( c[i] ) {
1162: switch ( NID(c[i]) ) {
1163: case N_LM:
1164: m = n_bits(((LM)c[i])->body);
1165: break;
1166: case N_Q:
1167: m = n_bits(((Q)c[i])->nm);
1168: break;
1169: default:
1170: error("maxblenup : invalid coefficient");
1171: }
1172: r = MAX(m,r);
1173: }
1174: return r;
1175: }
1176:
1177: void uptofmarray(mod,n,f)
1178: int mod;
1179: UP n;
1180: ModNum *f;
1181: {
1182: int d,i;
1183: unsigned int r;
1184: Num *c;
1185:
1186: if ( !n )
1187: return;
1188: else {
1189: d = n->d; c = (Num *)n->c;
1190: for ( i = 0; i <= d; i++ ) {
1191: if ( c[i] ) {
1192: switch ( NID(c[i]) ) {
1193: case N_LM:
1194: f[i] = (ModNum)rem(((LM)c[i])->body,mod);
1195: break;
1196: case N_Q:
1197: r = rem(NM((Q)c[i]),mod);
1198: if ( r && (SGN((Q)c[i])<0) )
1199: r = mod-r;
1200: f[i] = (ModNum)r;
1201: break;
1202: default:
1203: error("uptofmarray : invalid coefficient");
1204: }
1205: } else
1206: f[i] = 0;
1207: }
1208: }
1209: }
1210:
1211: void fmarraytoup(f,d,nr)
1212: ModNum *f;
1213: int d;
1214: UP *nr;
1215: {
1216: int i;
1217: Q *c;
1218: UP r;
1219:
1220: if ( d < 0 ) {
1221: *nr = 0;
1222: } else {
1223: *nr = r = UPALLOC(d); c = (Q *)r->c;
1224: for ( i = 0; i <= d; i++ ) {
1225: UTOQ((unsigned int)f[i],c[i]);
1226: }
1227: for ( i = d; i >= 0 && !c[i]; i-- );
1228: if ( i < 0 )
1229: *nr = 0;
1230: else
1231: r->d = i;
1232: }
1233: }
1234:
1235: /* f[i]: an array of length n */
1236:
1237: void uiarraytoup(f,n,d,nr)
1238: unsigned int **f;
1239: int n,d;
1240: UP *nr;
1241: {
1242: int i,j;
1243: unsigned int *fi;
1244: N ci;
1245: Q *c;
1246: UP r;
1247:
1248: if ( d < 0 ) {
1249: *nr = 0;
1250: } else {
1251: *nr = r = UPALLOC(d); c = (Q *)r->c;
1252: for ( i = 0; i <= d; i++ ) {
1253: fi = f[i];
1254: for ( j = n-1; j >= 0 && !fi[j]; j-- );
1255: j++;
1256: if ( j ) {
1257: ci = NALLOC(j);
1258: PL(ci) = j;
1259: bcopy(fi,BD(ci),j*sizeof(unsigned int));
1260: NTOQ(ci,1,c[i]);
1261: } else
1262: c[i] = 0;
1263: }
1264: for ( i = d; i >= 0 && !c[i]; i-- );
1265: if ( i < 0 )
1266: *nr = 0;
1267: else
1268: r->d = i;
1269: }
1270: }
1271:
1272: void adj_coefup(n,m,m2,nr)
1273: UP n;
1274: N m,m2;
1275: UP *nr;
1276: {
1277: int d;
1278: Q *c,*cr;
1279: Q mq;
1280: int i;
1281: UP r;
1282:
1283: if ( !n )
1284: *nr = 0;
1285: else {
1286: d = n->d; c = (Q *)n->c;
1287: *nr = r = UPALLOC(d); cr = (Q *)r->c;
1288: NTOQ(m,1,mq);
1289: for ( i = 0; i <= d; i++ ) {
1290: if ( !c[i] )
1291: continue;
1292: if ( cmpn(NM(c[i]),m2) > 0 )
1293: subq(c[i],mq,&cr[i]);
1294: else
1295: cr[i] = c[i];
1296: }
1297: for ( i = d; i >= 0 && !cr[i]; i-- );
1298: if ( i < 0 )
1299: *nr = 0;
1300: else
1301: r->d = i;
1302: }
1303: }
1304:
1305: /* n is assumed to have positive integer coefficients. */
1306:
1307: void remcup(n,mod,nr)
1308: UP n;
1309: N mod;
1310: UP *nr;
1311: {
1312: int i,d;
1313: Q *c,*cr;
1314: UP r;
1315: N t;
1316:
1317: if ( !n )
1318: *nr = 0;
1319: else {
1320: d = n->d; c = (Q *)n->c;
1321: *nr = r = UPALLOC(d); cr = (Q *)r->c;
1322: for ( i = 0; i <= d; i++ )
1323: if ( c[i] ) {
1324: remn(NM(c[i]),mod,&t);
1325: if ( t )
1326: NTOQ(t,1,cr[i]);
1327: else
1328: cr[i] = 0;
1329: } else
1330: cr[i] = 0;
1331: for ( i = d; i >= 0 && !cr[i]; i-- );
1332: if ( i < 0 )
1333: *nr = 0;
1334: else
1335: r->d = i;
1336: }
1337: }
1338:
1339: void fft_mulup_main(UP,UP,int,UP *);
1340:
1341: void fft_mulup(n1,n2,nr)
1342: UP n1,n2;
1343: UP *nr;
1344: {
1345: int d1,d2,d,b1,b2,h;
1346: UP n1lo,n1hi,n2lo,n2hi,lo,hi,t1,t2,mid1,mid2,mid,s1,s2;
1347:
1348: if ( !n1 || !n2 )
1349: *nr = 0;
1350: else {
1351: d1 = n1->d; d2 = n2->d;
1352: if ( MAX(d1,d2) < up_fft_mag )
1353: kmulup(n1,n2,nr);
1354: else {
1355: b1 = maxblenup(n1); b2 = maxblenup(n2);
1356: if ( fft_available(d1,b1,d2,b2) )
1357: fft_mulup_main(n1,n2,0,nr);
1358: else {
1359: d = MAX(d1,d2)+1; h = (d+1)/2;
1360: decompup(n1,h,&n1lo,&n1hi);
1361: decompup(n2,h,&n2lo,&n2hi);
1362: fft_mulup(n1lo,n2lo,&lo);
1363: fft_mulup(n1hi,n2hi,&hi);
1364: shiftup(hi,2*h,&t2);
1365: addup(lo,t2,&t1);
1366: addup(hi,lo,&mid1);
1367: subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2);
1368: fft_mulup(s1,s2,&mid2);
1369: addup(mid1,mid2,&mid);
1370: shiftup(mid,h,&t2);
1371: addup(t1,t2,nr);
1372: }
1373: }
1374: }
1375: }
1376:
1377: void trunc_fft_mulup(n1,n2,dbd,nr)
1378: UP n1,n2;
1379: int dbd;
1380: UP *nr;
1381: {
1382: int d1,d2,b1,b2,m;
1383: UP n1l,n1h,n2l,n2h,l,h,t,s,u;
1384:
1385: if ( !n1 || !n2 )
1386: *nr = 0;
1387: else if ( dbd < 0 )
1388: error("trunc_fft_mulup : invalid argument");
1389: else if ( dbd == 0 )
1390: *nr = 0;
1391: else {
1392: truncup(n1,dbd,&t); n1 = t;
1393: truncup(n2,dbd,&t); n2 = t;
1394: d1 = n1->d; d2 = n2->d;
1395: if ( MAX(d1,d2) < up_fft_mag )
1396: tkmulup(n1,n2,dbd,nr);
1397: else {
1398: b1 = maxblenup(n1); b2 = maxblenup(n2);
1399: if ( fft_available(d1,b1,d2,b2) )
1400: fft_mulup_main(n1,n2,dbd,nr);
1401: else {
1402: m = (dbd+1)/2;
1403: decompup(n1,m,&n1l,&n1h);
1404: decompup(n2,m,&n2l,&n2h);
1405: fft_mulup(n1l,n2l,&l);
1406: trunc_fft_mulup(n1l,n2h,dbd-m,&t);
1407: trunc_fft_mulup(n2l,n1h,dbd-m,&s);
1408: addup(t,s,&u);
1409: shiftup(u,m,&h);
1410: addup(l,h,&t);
1411: truncup(t,dbd,nr);
1412: }
1413: }
1414: }
1415: }
1416:
1417: void fft_squareup(n1,nr)
1418: UP n1;
1419: UP *nr;
1420: {
1421: int d1,d,h,len,b1;
1422: UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
1423:
1424: if ( !n1 )
1425: *nr = 0;
1426: else {
1427: d1 = n1->d;
1428: if ( d1 < up_fft_mag )
1429: ksquareup(n1,nr);
1430: else {
1431: b1 = maxblenup(n1);
1432: if ( fft_available(d1,b1,d1,b1) )
1433: fft_mulup_main(n1,n1,0,nr);
1434: else {
1435: d = d1+1; h = (d1+1)/2;
1436: decompup(n1,h,&n1lo,&n1hi);
1437: fft_squareup(n1hi,&hi);
1438: fft_squareup(n1lo,&lo);
1439: shiftup(hi,2*h,&t2);
1440: addup(lo,t2,&t1);
1441: addup(hi,lo,&mid1);
1442: subup(n1hi,n1lo,&s1);
1443: fft_squareup(s1,&mid2);
1444: subup(mid1,mid2,&mid);
1445: shiftup(mid,h,&t2);
1446: addup(t1,t2,nr);
1447: }
1448: }
1449: }
1450: }
1451:
1452: /*
1453: * dbd == 0 => n1 * n2
1454: * dbd > 0 => n1 * n2 mod x^dbd
1455: * n1 == n2 => squaring
1456: */
1457:
1458: void fft_mulup_main(n1,n2,dbd,nr)
1459: UP n1,n2;
1460: UP *nr;
1461: {
1462: ModNum *f1,*f2,*w,*fr;
1463: ModNum **frarray,**fa;
1464: unsigned int *modarray,*ma;
1465: int modarray_len;
1466: int frarray_index = 0;
1467: N m,m1,m2;
1468: int d1,d2,dmin,i,mod,root,d,cond,bound;
1469: UP r;
1470:
1471: if ( !n1 || !n2 ) {
1472: *nr = 0; return;
1473: }
1474: d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
1475: if ( !d1 || !d2 ) {
1476: mulup(n1,n2,nr); return;
1477: }
1478: m = ONEN;
1479: bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;
1480: f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
1481: if ( n1 == n2 )
1482: f2 = 0;
1483: else
1484: f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
1485: w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
1486: modarray_len = 1024;
1487: modarray = (unsigned int *)MALLOC_ATOMIC(modarray_len*sizeof(unsigned int));
1488: frarray = (ModNum **)MALLOC(modarray_len*sizeof(ModNum *));
1489: for ( i = 0; i < NPrimes; i++ ) {
1490: FFT_primes(i,&mod,&root,&d);
1491: if ( (1<<d) < d1+d2+1 )
1492: continue;
1493: if ( frarray_index == modarray_len ) {
1494: ma = (unsigned int *)MALLOC_ATOMIC(2*modarray_len*sizeof(unsigned int));
1495: bcopy(modarray,ma,modarray_len*sizeof(unsigned int));
1496: modarray = ma;
1497: fa = (ModNum **)MALLOC(2*modarray_len*sizeof(ModNum *));
1498: bcopy(frarray,fa,modarray_len*sizeof(ModNum *));
1499: frarray = fa;
1500: modarray_len *= 2;
1501: }
1502: modarray[frarray_index] = mod;
1503: frarray[frarray_index++] = fr
1504: = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
1505: uptofmarray(mod,n1,f1);
1506: if ( !f2 )
1507: cond = FFT_pol_square(d1,f1,fr,i,w);
1508: else {
1509: uptofmarray(mod,n2,f2);
1510: cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
1511: }
1512: if ( cond )
1513: error("fft_mulup : error in FFT_pol_product");
1514: STON(mod,m1); muln(m,m1,&m2); m = m2;
1515: if ( n_bits(m) > bound ) {
1516: if ( !dbd )
1517: dbd = d1+d2+1;
1518: crup(frarray,MIN(d1+d2,dbd-1),modarray,frarray_index,m,&r);
1519: divin(m,2,&m2);
1520: adj_coefup(r,m,m2,nr);
1521: return;
1522: }
1523: }
1524: error("fft_mulup : FFT_primes exhausted");
1525: }
1526: #if 0
1527: /* inefficient version */
1528:
1529: void crup(f,d,mod,index,m,r)
1530: ModNum **f;
1531: int d;
1532: int *mod;
1533: int index;
1534: N m;
1535: UP *r;
1536: {
1537: N *cof,*c;
1538: int *inv;
1539: struct oUP w;
1540: int i;
1541: UP s,s1,s2;
1542: N t;
1543: UP *g;
1544: Q q;
1545: struct oEGT eg0,eg1;
1546:
1547: get_eg(&eg0);
1548: w.d = 0;
1549: inv = (int *)ALLOCA(index*sizeof(int));
1550: cof = (N *)ALLOCA(index*sizeof(N));
1551: c = (N *)ALLOCA(index*sizeof(N));
1552: g = (UP *)ALLOCA(index*sizeof(UP));
1553: up_lazy = 1;
1554: for ( i = 0, s = 0; i < index; i++ ) {
1555: divin(m,mod[i],&cof[i]);
1556: inv[i] = invm(rem(cof[i],mod[i]),mod[i]);
1557: STON(inv[i],t);
1558: muln(cof[i],t,&c[i]);
1559: NTOQ(c[i],1,q); w.c[0] = (Num)q;
1560: fmarraytoup(f[i],d,&g[i]);
1561: mulup(g[i],&w,&s1);
1562: addup(s,s1,&s2);
1563: s = s2;
1564: }
1565: up_lazy = 0;
1566: remcup(s,m,r);
1567: get_eg(&eg1);
1568: add_eg(&eg_chrem,&eg0,&eg1);
1569: }
1570:
1571: #else
1572: /* improved version */
1573:
1574: void crup(f,d,mod,index,m,r)
1575: ModNum **f;
1576: int d;
1577: int *mod;
1578: int index;
1579: N m;
1580: UP *r;
1581: {
1582: N cof,c,t,w,w1;
1583: struct oN fc;
1584: N *s;
1585: UP u;
1586: Q q;
1587: int inv,i,j,rlen;
1588:
1589: rlen = PL(m)+10; /* XXX */
1590: PL(&fc) = 1;
1591: c = NALLOC(rlen);
1592: w = NALLOC(rlen);
1593: w1 = NALLOC(rlen);
1594: s = (N *)MALLOC((d+1)*sizeof(N));
1595: for ( i = 0; i <= d; i++ ) {
1596: s[i] = NALLOC(rlen);
1597: PL(s[i]) = 0;
1598: bzero(BD(s[i]),rlen*sizeof(unsigned int));
1599: }
1600: for ( i = 0; i < index; i++ ) {
1601: divin(m,mod[i],&cof);
1602: inv = invm(rem(cof,mod[i]),mod[i]); STON(inv,t);
1603: _muln(cof,t,c);
1604: /* s += c*f[i] */
1605: for ( j = 0; j <= d; j++ )
1606: if ( f[i][j] ) {
1607: BD(&fc)[0] = f[i][j];
1608: _muln(c,&fc,w);
1609: _addn(s[j],w,w1);
1610: dupn(w1,s[j]);
1611: }
1612: }
1613: for ( i = d; i >= 0; i-- )
1614: if ( s[i] )
1615: break;
1616: if ( i < 0 )
1617: *r = 0;
1618: else {
1619: u = UPALLOC(i);
1620: DEG(u) = i;
1621: for ( j = 0; j <= i; j++ ) {
1622: if ( PL(s[j]) )
1623: NTOQ(s[j],1,q);
1624: else
1625: q = 0;
1626: COEF(u)[j] = (Num)q;
1627: }
1628: remcup(u,m,r);
1629: }
1630: }
1631: #endif
1632:
1.2 noro 1633: /*
1634: * dbd == 0 => n1 * n2
1635: * dbd > 0 => n1 * n2 mod x^dbd
1636: * n1 == n2 => squaring
1637: * return: n1*n2 mod Primes[modind[0]]*.prime...*Primes[modind[nmod-1]].prime
1638: */
1639:
1640: void fft_mulup_specialmod_main(n1,n2,dbd,modind,nmod,nr)
1641: UP n1,n2;
1642: int dbd;
1643: int *modind;
1644: int nmod;
1645: UP *nr;
1646: {
1647: ModNum *f1,*f2,*w,*fr;
1648: ModNum **frarray,**fa;
1649: N m,m1,m2;
1650: unsigned int *modarray;
1651: int d1,d2,dmin,i,mod,root,d,cond,bound;
1652: UP r;
1653:
1654: if ( !n1 || !n2 ) {
1655: *nr = 0; return;
1656: }
1657: d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
1658: if ( !d1 || !d2 ) {
1659: mulup(n1,n2,nr); return;
1660: }
1661: m = ONEN;
1662: bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;
1663: f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
1664: if ( n1 == n2 )
1665: f2 = 0;
1666: else
1667: f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
1668: w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
1669: frarray = (ModNum **)MALLOC(nmod*sizeof(ModNum *));
1670: modarray = (unsigned int *)MALLOC_ATOMIC(nmod*sizeof(unsigned int *));
1671:
1672: for ( i = 0; i < nmod; i++ ) {
1673: FFT_primes(modind[i],&modarray[i],&root,&d);
1674: if ( (1<<d) < d1+d2+1 )
1675: error("fft_mulup_specialmod_main : invalid modulus");
1676: frarray[i] = fr
1677: = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
1678: uptofmarray(modarray[i],n1,f1);
1679: if ( !f2 )
1680: cond = FFT_pol_square(d1,f1,fr,modind[i],w);
1681: else {
1682: uptofmarray(modarray[i],n2,f2);
1683: cond = FFT_pol_product(d1,f1,d2,f2,fr,modind[i],w);
1684: }
1685: if ( cond )
1686: error("fft_mulup_specialmod_main : error in FFT_pol_product");
1687: STON(modarray[i],m1); muln(m,m1,&m2); m = m2;
1688: }
1689: if ( !dbd )
1690: dbd = d1+d2+1;
1691: crup(frarray,MIN(d1+d2,dbd-1),modarray,nmod,m,nr);
1692: }