![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Q.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 Q.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.4 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.5 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.4 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.9 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.8 2002/08/08 10:57:01 noro Exp $
1.4 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
52: #include "inline.h"
53:
1.7 noro 54: void addq(Q n1,Q n2,Q *nr)
1.1 noro 55: {
56: N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;
57: int sgn;
58: unsigned t,t1;
59:
60: if ( !n1 )
61: *nr = n2;
62: else if ( !n2 )
63: *nr = n1;
64: else if ( INT(n1) )
65: if ( INT(n2) ) {
66: nm1 = NM(n1); nm2 = NM(n2);
67: if ( SGN(n1) == SGN(n2) ) {
68: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
69: t1 = BD(nm1)[0]; t = t1+BD(nm2)[0];
70: if ( t < t1 ) {
71: m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1;
72: } else
73: UTON(t,m);
74: } else
75: addn(NM(n1),NM(n2),&m);
76: NTOQ(m,SGN(n1),*nr);
77: } else {
78: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
79: t1 = BD(nm1)[0]; t = t1-BD(nm2)[0];
80: if ( !t )
81: sgn = 0;
82: else if ( t > t1 ) {
83: sgn = -1; t = -((int)t); UTON(t,m);
84: } else {
85: sgn = 1; UTON(t,m);
86: }
87: } else
88: sgn = subn(NM(n1),NM(n2),&m);
89: if ( sgn ) {
90: if ( SGN(n1) == sgn )
91: NTOQ(m,1,*nr);
92: else
93: NTOQ(m,-1,*nr);
94: } else
95: *nr = 0;
96: }
97: } else {
98: kmuln(NM(n1),DN(n2),&m);
99: if ( SGN(n1) == SGN(n2) ) {
100: sgn = SGN(n1); addn(m,NM(n2),&nm);
101: } else
102: sgn = SGN(n1)*subn(m,NM(n2),&nm);
103: if ( sgn ) {
104: dn = DN(n2); NDTOQ(nm,dn,sgn,*nr);
105: } else
106: *nr = 0;
107: }
108: else if ( INT(n2) ) {
109: kmuln(NM(n2),DN(n1),&m);
110: if ( SGN(n1) == SGN(n2) ) {
111: sgn = SGN(n1); addn(m,NM(n1),&nm);
112: } else
113: sgn = SGN(n1)*subn(NM(n1),m,&nm);
114: if ( sgn ) {
115: dn = DN(n1); NDTOQ(nm,dn,sgn,*nr);
116: } else
117: *nr = 0;
118: } else {
119: gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2);
120: kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2);
121: if ( SGN(n1) == SGN(n2) ) {
122: sgn = SGN(n1); addn(nm1,nm2,&nm3);
123: } else
124: sgn = SGN(n1)*subn(nm1,nm2,&nm3);
125: if ( sgn ) {
126: gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0);
127: kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn);
128: if ( UNIN(dn) )
129: NTOQ(nm,sgn,*nr);
130: else
131: NDTOQ(nm,dn,sgn,*nr);
132: } else
133: *nr = 0;
134: }
135: }
136:
1.7 noro 137: void subq(Q n1,Q n2,Q *nr)
1.1 noro 138: {
139: Q m;
140:
141: if ( !n1 )
142: if ( !n2 )
143: *nr = 0;
144: else {
145: DUPQ(n2,*nr); SGN(*nr) = -SGN(n2);
146: }
147: else if ( !n2 )
148: *nr = n1;
149: else if ( n1 == n2 )
150: *nr = 0;
151: else {
152: DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr);
153: }
154: }
155:
1.7 noro 156: void mulq(Q n1,Q n2,Q *nr)
1.1 noro 157: {
158: N nm,nm1,nm2,dn,dn1,dn2,g;
159: int sgn;
160: unsigned u,l;
161:
162: if ( !n1 || !n2 ) *nr = 0;
163: else if ( INT(n1) )
164: if ( INT(n2) ) {
165: nm1 = NM(n1); nm2 = NM(n2);
166: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
167: DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l)
168: if ( u ) {
169: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u;
170: } else {
171: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l;
172: }
173: } else
174: kmuln(nm1,nm2,&nm);
175: sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr);
176: } else {
177: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn);
178: kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2);
179: if ( UNIN(dn) )
180: NTOQ(nm,sgn,*nr);
181: else
182: NDTOQ(nm,dn,sgn,*nr);
183: }
184: else if ( INT(n2) ) {
185: gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn);
186: kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2);
187: if ( UNIN(dn) )
188: NTOQ(nm,sgn,*nr);
189: else
190: NDTOQ(nm,dn,sgn,*nr);
191: } else {
192: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2);
193: gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2);
194: kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2);
195: if ( UNIN(dn) )
196: NTOQ(nm,sgn,*nr);
197: else
198: NDTOQ(nm,dn,sgn,*nr);
199: }
200: }
201:
1.7 noro 202: void divq(Q n1,Q n2,Q *nq)
1.1 noro 203: {
204: Q m;
205:
206: if ( !n2 ) {
207: error("division by 0");
208: *nq = 0;
209: return;
210: } else if ( !n1 )
211: *nq = 0;
212: else if ( n1 == n2 )
213: *nq = ONE;
214: else {
215: invq(n2,&m); mulq(n1,m,nq);
1.6 noro 216: }
217: }
218:
1.7 noro 219: void divsq(Q n1,Q n2,Q *nq)
1.6 noro 220: {
221: int sgn;
222: N tn;
223:
224: if ( !n2 ) {
225: error("division by 0");
226: *nq = 0;
227: return;
228: } else if ( !n1 )
229: *nq = 0;
230: else if ( n1 == n2 )
231: *nq = ONE;
232: else {
233: divsn(NM(n1),NM(n2),&tn);
234: sgn = SGN(n1)*SGN(n2);
235: NTOQ(tn,sgn,*nq);
1.1 noro 236: }
237: }
238:
1.7 noro 239: void invq(Q n,Q *nr)
1.1 noro 240: {
241: N nm,dn;
242:
243: if ( INT(n) )
244: if ( UNIN(NM(n)) )
245: if ( SGN(n) > 0 )
246: *nr = ONE;
247: else {
248: nm = ONEN; NTOQ(nm,SGN(n),*nr);
249: }
250: else {
251: nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
252: }
253: else if ( UNIN(NM(n)) ) {
254: nm = DN(n); NTOQ(nm,SGN(n),*nr);
255: } else {
256: nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
257: }
258: }
259:
1.7 noro 260: void chsgnq(Q n,Q *nr)
1.1 noro 261: {
262: Q t;
263:
264: if ( !n )
265: *nr = 0;
266: else {
267: DUPQ(n,t); SGN(t) = -SGN(t); *nr = t;
268: }
269: }
270:
1.7 noro 271: void pwrq(Q n1,Q n,Q *nr)
1.1 noro 272: {
273: N nm,dn;
274: int sgn;
275:
276: if ( !n || UNIQ(n1) )
277: *nr = ONE;
278: else if ( !n1 )
279: *nr = 0;
280: else if ( !INT(n) ) {
281: error("can't calculate fractional power."); *nr = 0;
282: } else if ( MUNIQ(n1) )
283: *nr = BD(NM(n))[0]%2 ? n1 : ONE;
284: else if ( PL(NM(n)) > 1 ) {
285: error("exponent too big."); *nr = 0;
286: } else {
287: sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1);
288: pwrn(NM(n1),BD(NM(n))[0],&nm);
289: if ( INT(n1) ) {
290: if ( UNIN(nm) )
291: if ( sgn == 1 )
292: *nr = ONE;
293: else
294: NTOQ(ONEN,sgn,*nr);
295: else if ( SGN(n) > 0 )
296: NTOQ(nm,sgn,*nr);
297: else
298: NDTOQ(ONEN,nm,sgn,*nr);
299: } else {
300: pwrn(DN(n1),BD(NM(n))[0],&dn);
301: if ( SGN(n) > 0 )
302: NDTOQ(nm,dn,sgn,*nr);
303: else if ( UNIN(nm) )
304: NTOQ(dn,sgn,*nr);
305: else
306: NDTOQ(dn,nm,sgn,*nr);
307: }
308: }
309: }
310:
1.7 noro 311: int cmpq(Q q1,Q q2)
1.1 noro 312: {
313: int sgn;
314: N t,s;
315:
316: if ( !q1 )
317: if ( !q2 )
318: return ( 0 );
319: else
320: return ( -SGN(q2) );
321: else if ( !q2 )
322: return ( SGN(q1) );
323: else if ( SGN(q1) != SGN(q2) )
324: return ( SGN(q1) );
325: else if ( INT(q1) )
326: if ( INT(q2) ) {
327: sgn = cmpn(NM(q1),NM(q2));
328: if ( !sgn )
329: return ( 0 );
330: else
331: return ( SGN(q1)==sgn?1:-1 );
332: } else {
333: kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2));
334: return ( SGN(q1) * sgn );
335: }
336: else if ( INT(q2) ) {
337: kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t);
338: return ( SGN(q1) * sgn );
339: } else {
340: kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s);
341: return ( SGN(q1) * sgn );
342: }
343: }
344:
1.7 noro 345: void remq(Q n,Q m,Q *nr)
1.1 noro 346: {
347: N q,r,s;
348:
349: if ( !n ) {
350: *nr = 0;
351: return;
352: }
353: divn(NM(n),NM(m),&q,&r);
354: if ( !r )
355: *nr = 0;
356: else if ( SGN(n) > 0 )
357: NTOQ(r,1,*nr);
358: else {
359: subn(NM(m),r,&s); NTOQ(s,1,*nr);
360: }
361: }
362:
1.2 noro 363: /* t = [nC0 nC1 ... nCn] */
364:
1.7 noro 365: void mkbc(int n,Q *t)
1.1 noro 366: {
367: int i;
368: N c,d;
369:
370: for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
371: c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;
372: kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
373: }
374: for ( ; i <= n; i++ )
375: t[i] = t[n-i];
1.2 noro 376: }
377:
378: /*
379: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
380: *
381: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
382: * where W(k,l,i) = i! * kCi * lCi
383: */
384:
1.7 noro 385: void mkwc(int k,int l,Q *t)
1.2 noro 386: {
387: int i,n,up,low;
388: N nm,d,c;
389:
390: n = MIN(k,l);
391: for ( t[0] = ONE, i = 1; i <= n; i++ ) {
392: DM(k-i+1,l-i+1,up,low);
393: if ( up ) {
394: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
395: } else {
396: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
397: }
398: kmuln(NM(t[i-1]),nm,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
1.3 noro 399: }
400: }
401:
402: /* mod m table */
403: /* XXX : should be optimized */
404:
1.7 noro 405: void mkwcm(int k,int l,int m,int *t)
1.3 noro 406: {
407: int i,n;
408: Q *s;
409:
410: n = MIN(k,l);
411: s = (Q *)ALLOCA((n+1)*sizeof(Q));
412: mkwc(k,l,s);
413: for ( i = 0; i <= n; i++ ) {
414: t[i] = rem(NM(s[i]),m);
1.2 noro 415: }
1.1 noro 416: }
417:
1.8 noro 418: #if 0
1.7 noro 419: void factorial(int n,Q *r)
1.1 noro 420: {
421: Q t,iq,s;
422: unsigned int i,m,m0;
423: unsigned int c;
424:
425: if ( !n )
426: *r = ONE;
427: else if ( n < 0 )
428: *r = 0;
429: else {
430: for ( i = 1, t = ONE; (int)i <= n; ) {
431: for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
432: m0 = m; DM(m0,i,c,m)
433: }
434: if ( c ) {
435: m = m0; i--;
436: }
437: UTOQ(m,iq); mulq(t,iq,&s); t = s;
438: }
439: *r = t;
440: }
1.8 noro 441: }
442: #endif
443:
444: void partial_factorial(int s,int e,N *r);
445:
446: void factorial(int n,Q *r)
447: {
448: N nm;
449:
450: if ( !n )
451: *r = ONE;
452: else if ( n < 0 )
453: *r = 0;
454: else {
455: partial_factorial(1,n,&nm);
456: NTOQ(nm,1,*r);
457: }
458: }
459:
460: /* s*(s+1)*...*e */
461: void partial_factorial(int s,int e,N *r)
462: {
463: unsigned int len,i,m,m0,c;
464: unsigned int *p,*p1,*p2;
465: N nm,r1,r2;
466:
467: /* XXX */
468: if ( e-s > 2000 ) {
469: m = (e+s)/2;
470: partial_factorial(s,m,&r1);
471: partial_factorial(m+1,e,&r2);
472: kmuln(r1,r2,r);
473: return;
474: }
475: /* find i s.t. 2^(i-1) < m <= 2^i */
476: for ( i = 0, m = 1; m < e; m <<=1, i++ );
477: /* XXX estimate of word length of the result */
478: len = (i*(e-s+1)+1+31)/32;
1.9 ! noro 479: p = ALLOCA((len+1)*sizeof(int));
! 480: p1 = ALLOCA((len+1)*sizeof(int));
1.8 noro 481: p[0] = s;
482: len = 1;
483: for ( i = s+1; (int)i <= e; ) {
484: for ( m0 = m = 1, c = 0; ((int)i <= e) && !c; i++ ) {
485: m0 = m; DM(m0,i,c,m)
486: }
487: if ( c ) {
488: m = m0; i--;
489: }
490: bzero(p1,(len+1)*sizeof(int));
491: muln_1(p,len,m,p1);
492: if ( p1[len] )
493: len++;
494: p2 = p; p = p1; p1 = p2;
495: }
496: *r = nm = NALLOC(len);
497: bcopy(p,BD(nm),sizeof(int)*len);
498: PL(nm) = len;
1.1 noro 499: }
500:
1.7 noro 501: void invl(Q a,Q mod,Q *ar)
1.1 noro 502: {
503: Q f1,f2,a1,a2,q,m,s;
504: N qn,rn;
505:
506: a1 = ONE; a2 = 0;
507: DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
508:
509: while ( 1 ) {
510: divn(NM(f1),NM(f2),&qn,&rn);
511: if ( !qn )
512: q = 0;
513: else
514: NTOQ(qn,1,q);
515:
516: f1 = f2;
517: if ( !rn )
518: break;
519: else
520: NTOQ(rn,1,f2);
521:
522: mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
523: if ( !s )
524: a2 = 0;
525: else
526: remq(s,mod,&a2);
527: }
528: if ( SGN(a) < 0 )
529: chsgnq(a2,&m);
530: else
531: m = a2;
532:
533: if ( SGN(m) >= 0 )
534: *ar = m;
535: else
536: addq(m,mod,ar);
537: }
538:
539: int kara_mag=100;
540:
1.7 noro 541: void kmuln(N n1,N n2,N *nr)
1.1 noro 542: {
543: N n,t,s,m,carry;
544: int d,d1,d2,len,i,l;
545: int *r,*r0;
546:
547: if ( !n1 || !n2 ) {
548: *nr = 0; return;
549: }
550: d1 = PL(n1); d2 = PL(n2);
551: if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
552: muln(n1,n2,nr); return;
553: }
554: if ( d1 < d2 ) {
555: n = n1; n1 = n2; n2 = n;
556: d = d1; d1 = d2; d2 = d;
557: }
558: if ( d2 > (d1+1)/2 ) {
559: kmulnmain(n1,n2,nr); return;
560: }
561: d = (d1/d2)+((d1%d2)!=0);
562: len = (d+1)*d2;
563: r0 = (int *)ALLOCA(len*sizeof(int));
564: bzero((char *)r0,len*sizeof(int));
565: for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
566: extractn(n1,i*d2,d2,&m);
567: if ( m ) {
568: kmuln(m,n2,&t);
569: addn(t,carry,&s);
570: copyn(s,d2,r);
571: extractn(s,d2,d2,&carry);
572: } else {
573: copyn(carry,d2,r);
574: carry = 0;
575: }
576: }
577: copyn(carry,d2,r);
578: for ( l = len - 1; !r0[l]; l-- );
579: l++;
580: *nr = t = NALLOC(l); PL(t) = l;
581: bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
582: }
583:
1.7 noro 584: void extractn(N n,int index,int len,N *nr)
1.1 noro 585: {
586: unsigned int *m;
587: int l;
588: N t;
589:
590: if ( !n ) {
591: *nr = 0; return;
592: }
593: m = BD(n)+index;
594: if ( (l = PL(n)-index) >= len ) {
595: for ( l = len - 1; (l >= 0) && !m[l]; l-- );
596: l++;
597: }
598: if ( l <= 0 ) {
599: *nr = 0; return;
600: } else {
601: *nr = t = NALLOC(l); PL(t) = l;
602: bcopy((char *)m,(char *)BD(t),l*sizeof(int));
603: }
604: }
605:
1.7 noro 606: void copyn(N n,int len,int *p)
1.1 noro 607: {
608: if ( n )
609: bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
610: }
611:
1.7 noro 612: void dupn(N n,N p)
1.1 noro 613: {
614: if ( n )
615: bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
616: }
617:
1.7 noro 618: void kmulnmain(N n1,N n2,N *nr)
1.1 noro 619: {
620: int d1,d2,h,sgn,sgn1,sgn2,len;
621: N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
622:
623: d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
624: extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
625: extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
626: kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
627: len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
628: bzero((char *)BD(t1),len*sizeof(int));
629: if ( lo )
630: bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
631: if ( hi )
632: bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
633:
634: addn(hi,lo,&mid1);
635: sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
636: kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
637: if ( sgn > 0 )
638: addn(mid1,mid2,&mid);
639: else
640: subn(mid1,mid2,&mid);
641: if ( mid ) {
642: len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
643: bzero((char *)BD(t2),len*sizeof(int));
644: bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
645: addn(t1,t2,nr);
646: } else
647: *nr = t1;
648: }