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