![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Z.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 Z.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: #include "ca.h"
1.5 noro 2: #include "base.h"
1.1 noro 3: #include "inline.h"
4:
1.5 noro 5: inline void _addz(Z n1,Z n2,Z nr);
6: inline void _subz(Z n1,Z n2,Z nr);
7: inline void _mulz(Z n1,Z n2,Z nr);
8: inline int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
9: inline int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
1.3 noro 10:
11: /* immediate int -> Z */
1.5 noro 12: #define UTOZ(c,n) (n)=(!((unsigned int)(c))?0:(((unsigned int)(c))<=IMM_MAX?((Z)((((unsigned int)(c))<<1)|1)):utoz((unsigned int)(c))))
13: #define STOZ(c,n) (n)=(!((int)(c))?0:(((int)(c))>=IMM_MIN&&((int)(c))<=IMM_MAX?((Z)((((int)(c))<<1)|1)):stoz((int)(c))))
1.3 noro 14: /* immediate Z ? */
15: #define IS_IMM(c) (((unsigned int)c)&1)
16: /* Z can be conver to immediate ? */
1.5 noro 17: #define IS_SZ(n) (((SL(n) == 1)||(SL(n)==-1))&&BD(n)[0]<=IMM_MAX)
1.3 noro 18: /* Z -> immediate Z */
1.5 noro 19: #define SZTOZ(n,z) (z)=(Z)(SL(n)<0?(((-BD(n)[0])<<1)|1):(((BD(n)[0])<<1)|1))
1.3 noro 20: /* Z -> immediate int */
1.5 noro 21: #define ZTOS(c) (((int)(c))>>1)
22:
23: int uniz(Z a)
24: {
25: if ( IS_IMM(a) && ZTOS(a) == 1 ) return 1;
26: else return 0;
27: }
28:
29: int cmpz(Z a,Z b)
30: {
31: int *ma,*mb;
32: int sa,sb,da,db,ca,cb,i;
33:
34: if ( !a )
35: return -sgnz(b);
36: else if ( !b )
37: return sgnz(a);
38: else {
39: sa = sgnz(a); sb = sgnz(b);
40: if ( sa > sb ) return 1;
41: else if ( sa < sb ) return -1;
42: else if ( IS_IMM(a) )
43: if ( IS_IMM(b) ) {
44: ca = ZTOS(a); cb = ZTOS(b);
45: if ( ca > cb ) return sa;
46: else if ( ca < cb ) return -sa;
47: else return 0;
48: } else
49: return -sa;
50: else if ( IS_IMM(b) )
51: return sa;
52: else {
53: da = SL(a)*sa; db = SL(b)*sa;
54: if ( da > db ) return sa;
55: else if ( da < db ) return -sa;
56: else {
57: for ( i = da-1, ma = BD(a)+i, mb = BD(b)+i; i >= 0; i-- )
58: if ( *ma > *mb ) return sa;
59: else if ( *ma < *mb ) return -sa;
60: return 0;
61: }
62: }
63: }
64: }
1.1 noro 65:
1.5 noro 66: Z stoz(int c)
1.3 noro 67: {
68: Z z;
69:
70: z = ZALLOC(1);
71: if ( c < 0 ) {
72: SL(z) = -1; BD(z)[0] = -c;
73: } else {
74: SL(z) = 1; BD(z)[0] = c;
75: }
76: return z;
77: }
78:
1.5 noro 79: Z utoz(unsigned int c)
80: {
81: Z z;
82:
83: z = ZALLOC(1);
84: SL(z) = 1; BD(z)[0] = c;
85: return z;
86: }
87:
88: Z simpz(Z n)
89: {
90: Z n1;
91:
92: if ( !n ) return 0;
93: else if ( IS_IMM(n) ) return n;
94: else if ( IS_SZ(n) ) {
95: SZTOZ(n,n1); return n1;
96: } else
97: return n;
98: }
99:
1.3 noro 100: int sgnz(Z n)
1.2 noro 101: {
102: if ( !n ) return 0;
1.5 noro 103: else if ( IS_IMM(n) ) return ZTOS(n)>0?1:-1;
1.2 noro 104: else if ( SL(n) < 0 ) return -1;
105: else return 1;
106: }
107:
1.1 noro 108: z_mag(Z n)
109: {
1.4 noro 110: int c,i;
1.3 noro 111:
112: if ( !n ) return 0;
113: else if ( IS_IMM(n) ) {
1.5 noro 114: c = ZTOS(n);
1.3 noro 115: if ( c < 0 ) c = -c;
116: for ( i = 0; c; c >>= 1, i++ );
117: return i;
118: }
119: else return n_bits((N)n);
1.1 noro 120: }
121:
122: Z qtoz(Q n)
123: {
1.3 noro 124: Z r,t;
1.4 noro 125: int c;
1.1 noro 126:
127: if ( !n ) return 0;
128: else if ( !INT(n) )
129: error("qtoz : invalid input");
130: else {
1.3 noro 131: t = (Z)NM(n);
1.5 noro 132: if ( IS_SZ(t) ) {
1.4 noro 133: c = SGN(n) < 0 ? -BD(t)[0] : BD(t)[0];
1.5 noro 134: STOZ(c,r);
1.3 noro 135: } else {
136: r = dupz((Z)t);
1.4 noro 137: if ( SGN(n) < 0 ) SL(r) = -SL(r);
1.3 noro 138: }
1.1 noro 139: return r;
140: }
141: }
142:
143: Q ztoq(Z n)
144: {
145: Q r;
146: Z nm;
1.3 noro 147: int sgn,c;
1.1 noro 148:
149: if ( !n ) return 0;
1.3 noro 150: else if ( IS_IMM(n) ) {
1.5 noro 151: c = ZTOS(n);
1.3 noro 152: STOQ(c,r);
153: return r;
154: } else {
1.1 noro 155: nm = dupz(n);
156: if ( SL(nm) < 0 ) {
157: sgn = -1;
158: SL(nm) = -SL(nm);
159: } else
160: sgn = 1;
161: NTOQ((N)nm,sgn,r);
162: return r;
163: }
164: }
165:
166: Z dupz(Z n)
167: {
1.5 noro 168: Z r;
1.1 noro 169: int sd,i;
170:
171: if ( !n ) return 0;
1.3 noro 172: else if ( IS_IMM(n) ) return n;
173: else {
174: if ( (sd = SL(n)) < 0 ) sd = -sd;
1.5 noro 175: r = ZALLOC(sd);
176: SL(r) = SL(n);
177: for ( i = 0; i < sd; i++ ) BD(r)[i] = BD(n)[i];
178: return r;
1.3 noro 179: }
1.1 noro 180: }
181:
182: Z chsgnz(Z n)
183: {
184: Z r;
1.3 noro 185: int c;
1.1 noro 186:
187: if ( !n ) return 0;
1.3 noro 188: else if ( IS_IMM(n) ) {
1.5 noro 189: c = -ZTOS(n);
190: STOZ(c,r);
1.3 noro 191: return r;
192: } else {
1.1 noro 193: r = dupz(n);
194: SL(r) = -SL(r);
195: return r;
196: }
197: }
198:
1.5 noro 199: Z absz(Z n)
200: {
201: Z r;
202: int c;
203:
204: if ( !n ) return 0;
205: else if ( sgnz(n) > 0 )
206: return n;
207: else
208: return chsgnz(n);
209: }
1.3 noro 210:
1.1 noro 211: Z addz(Z n1,Z n2)
212: {
1.3 noro 213: int d1,d2,d,c;
214: Z r,r1;
215: struct oZ t;
1.1 noro 216:
217: if ( !n1 ) return dupz(n2);
218: else if ( !n2 ) return dupz(n1);
1.3 noro 219: else if ( IS_IMM(n1) ) {
220: if ( IS_IMM(n2) ) {
1.5 noro 221: c = ZTOS(n1)+ZTOS(n2);
222: STOZ(c,r);
1.3 noro 223: } else {
1.5 noro 224: c = ZTOS(n1);
1.3 noro 225: if ( c < 0 ) {
226: t.p = -1; t.b[0] = -c;
1.1 noro 227: } else {
1.3 noro 228: t.p = 1; t.b[0] = c;
1.1 noro 229: }
1.3 noro 230: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
231: r = ZALLOC(d2+1);
232: _addz(&t,n2,r);
1.5 noro 233: if ( !SL(r) ) r = 0;
1.3 noro 234: }
235: } else if ( IS_IMM(n2) ) {
1.5 noro 236: c = ZTOS(n2);
1.3 noro 237: if ( c < 0 ) {
238: t.p = -1; t.b[0] = -c;
1.1 noro 239: } else {
1.3 noro 240: t.p = 1; t.b[0] = c;
241: }
242: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
243: r = ZALLOC(d1+1);
244: _addz(n1,&t,r);
1.5 noro 245: if ( !SL(r) ) r = 0;
1.3 noro 246: } else {
247: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
248: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
249: d = MAX(d1,d2)+1;
250: r = ZALLOC(d);
251: _addz(n1,n2,r);
252: if ( !SL(r) ) r = 0;
1.1 noro 253: }
1.5 noro 254: if ( r && !((int)r&1) && IS_SZ(r) ) {
255: SZTOZ(r,r1); r = r1;
256: }
257: return r;
1.1 noro 258: }
259:
260: Z subz(Z n1,Z n2)
261: {
1.3 noro 262: int d1,d2,d,c;
263: Z r,r1;
264: struct oZ t;
1.1 noro 265:
1.5 noro 266: if ( !n1 )
267: return chsgnz(n2);
268: else if ( !n2 ) return n1;
1.3 noro 269: else if ( IS_IMM(n1) ) {
270: if ( IS_IMM(n2) ) {
1.5 noro 271: c = ZTOS(n1)-ZTOS(n2);
272: STOZ(c,r);
1.3 noro 273: } else {
1.5 noro 274: c = ZTOS(n1);
1.3 noro 275: if ( c < 0 ) {
276: t.p = -1; t.b[0] = -c;
277: } else {
278: t.p = 1; t.b[0] = c;
279: }
280: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
281: r = ZALLOC(d2+1);
282: _subz(&t,n2,r);
1.5 noro 283: if ( !SL(r) ) r = 0;
1.3 noro 284: }
285: } else if ( IS_IMM(n2) ) {
1.5 noro 286: c = ZTOS(n2);
1.3 noro 287: if ( c < 0 ) {
288: t.p = -1; t.b[0] = -c;
289: } else {
290: t.p = 1; t.b[0] = c;
291: }
292: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
293: r = ZALLOC(d1+1);
294: _subz(n1,&t,r);
1.5 noro 295: if ( !SL(r) ) r = 0;
1.3 noro 296: } else {
297: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
298: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
299: d = MAX(d1,d2)+1;
300: r = ZALLOC(d);
301: _subz(n1,n2,r);
302: if ( !SL(r) ) r = 0;
1.1 noro 303: }
1.5 noro 304: if ( r && !((int)r&1) && IS_SZ(r) ) {
305: SZTOZ(r,r1); r = r1;
306: }
307: return r;
1.1 noro 308: }
309:
310: Z mulz(Z n1,Z n2)
311: {
1.3 noro 312: int d1,d2,sgn,i;
1.5 noro 313: int c1,c2;
1.1 noro 314: unsigned int u1,u2,u,l;
1.3 noro 315: Z r;
1.1 noro 316:
317: if ( !n1 || !n2 ) return 0;
1.3 noro 318:
319: if ( IS_IMM(n1) ) {
1.5 noro 320: c1 = ZTOS(n1);
1.3 noro 321: if ( IS_IMM(n2) ) {
1.5 noro 322: c2 = ZTOS(n2);
323: if ( c1 == 1 )
324: return n2;
325: else if ( c1 == -1 )
326: return chsgnz(n2);
327: else if ( c2 == 1 )
328: return n1;
329: else if ( c2 == -1 )
330: return chsgnz(n1);
331: else {
332: sgn = 1;
333: if ( c1 < 0 ) { c1 = -c1; sgn = -sgn; }
334: if ( c2 < 0 ) { c2 = -c2; sgn = -sgn; }
335: u1 = (unsigned int)c1; u2 = (unsigned int)c2;
336: DM(u1,u2,u,l);
337: if ( !u ) {
338: UTOZ(l,r);
339: } else {
340: r = ZALLOC(2); SL(r) = 2; BD(r)[1] = u; BD(r)[0] = l;
341: }
1.1 noro 342: }
343: } else {
1.5 noro 344: sgn = 1;
345: if ( c1 < 0 ) { c1 = -c1; sgn = -sgn; }
346: u1 = (unsigned int)c1;
1.3 noro 347: if ( (d2 = SL(n2)) < 0 ) { sgn = -sgn; d2 = -d2; }
348: r = ZALLOC(d2+1);
349: for ( i = d2; i >= 0; i-- ) BD(r)[i] = 0;
350: muln_1(BD(n2),d2,u1,BD(r));
351: SL(r) = BD(r)[d2]?d2+1:d2;
1.1 noro 352: }
1.3 noro 353: } else if ( IS_IMM(n2) ) {
1.5 noro 354: c2 = ZTOS(n2);
355: if ( c2 == 1 )
356: return n1;
357: else if ( c2 == -1 )
358: return chsgnz(n1);
359:
1.3 noro 360: sgn = 1;
1.5 noro 361: if ( c2 < 0 ) { sgn = -sgn; c2 = -c2; }
362: u2 = (unsigned int)c2;
1.3 noro 363: if ( (d1 = SL(n1)) < 0 ) { sgn = -sgn; d1 = -d1; }
364: r = ZALLOC(d1+1);
365: for ( i = d1; i >= 0; i-- ) BD(r)[i] = 0;
366: muln_1(BD(n1),d1,u2,BD(r));
367: SL(r) = BD(r)[d1]?d1+1:d1;
368: } else {
369: sgn = 1;
1.5 noro 370: if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
371: if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
1.3 noro 372: r = ZALLOC(d1+d2);
373: _mulz(n1,n2,r);
1.1 noro 374: }
1.5 noro 375: if ( sgn < 0 ) r = chsgnz(r);
1.3 noro 376: return r;
1.1 noro 377: }
378:
1.3 noro 379: /* kokokara */
1.4 noro 380: #if 0
1.1 noro 381: Z divsz(Z n1,Z n2)
382: {
1.3 noro 383: int sgn,d1,d2;
384: Z q;
1.1 noro 385:
386: if ( !n2 ) error("divsz : division by 0");
1.3 noro 387: if ( !n1 ) return 0;
388:
389: if ( IS_IMM(n1) ) {
390: if ( !IS_IMM(n2) )
391: error("divsz : cannot happen");
1.5 noro 392: c = ZTOS(n1)/ZTOS(n2);
393: STOZ(c,q);
1.3 noro 394: return q;
1.1 noro 395: }
1.3 noro 396: if ( IS_IMM(n2) ) {
397: sgn = 1;
1.5 noro 398: u2 = ZTOS(n2); if ( u2 < 0 ) { sgn = -sgn; u2 = -u2; }
1.3 noro 399: diviz(n1,u2,&q);
400: if ( sgn < 0 ) SL(q) = -SL(q);
401: return q;
402: }
403:
404: sgn = 1;
405: if ( (d2 = SL(n2)) < 0 ) { sgn = -sgn; d2 = -d2; }
406: if ( d2 == 1 ) {
407: diviz(n1,BD(u2)[0],&q);
408: if ( sgn < 0 ) SL(q) = -SL(q);
409: return q;
410: }
411: if ( (d1 = SL(n1)) < 0 ) { sgn = -sgn; d1 = -d1; }
412: if ( d1 < d2 ) error("divsz : cannot happen");
413: return q;
1.1 noro 414: }
1.5 noro 415: #endif
416:
417: /* XXX */
418: Z divz(Z n1,Z n2,Z *r)
419: {
420: int s1,s2;
421: Q t1,t2,qq,rq;
422: N qn,rn;
423:
424: if ( !n1 ) {
425: *r = 0; return 0;
426: }
427: if ( !n2 )
428: error("divz : division by 0");
429: t1 = ztoq(n1); t2 = ztoq(n2);
430: s1 = sgnz(n1); s2 = sgnz(n2);
431: /* n1 = qn*SGN(n1)*SGN(n2)*n2+SGN(n1)*rn */
432: divn(NM(t1),NM(t2),&qn,&rn);
433: NTOQ(qn,s1*s2,qq);
434: NTOQ(rn,s1,rq);
435: *r = qtoz(rq);
436: return qtoz(qq);
437: }
438:
439: Z divsz(Z n1,Z n2)
440: {
441: int s1,s2;
442: Q t1,t2,qq;
443: N qn;
444:
445: if ( !n1 )
446: return 0;
447: if ( !n2 )
448: error("divsz : division by 0");
449: t1 = ztoq(n1); t2 = ztoq(n2);
450: s1 = sgnz(n1); s2 = sgnz(n2);
451: /* n1 = qn*SGN(n1)*SGN(n2)*n2 */
452: divsn(NM(t1),NM(t2),&qn);
453: NTOQ(qn,s1*s2,qq);
454: return qtoz(qq);
455: }
456:
457: int gcdimm(int c1,int c2)
458: {
459: int r;
460:
461: if ( !c1 ) return c2;
462: else if ( !c2 ) return c1;
463: while ( 1 ) {
464: r = c1%c2;
465: if ( !r ) return c2;
466: c1 = c2; c2 = r;
467: }
468: }
1.1 noro 469:
470: Z gcdz(Z n1,Z n2)
471: {
1.5 noro 472: int c1,c2,g;
473: Z gcd,r;
474: N gn;
1.1 noro 475:
476: if ( !n1 ) return n2;
477: else if ( !n2 ) return n1;
478:
1.5 noro 479: if ( IS_IMM(n1) ) {
480: c1 = ZTOS(n1);
481: if ( c1 < 0 ) c1 = -c1;
482: if ( IS_IMM(n2) )
483: c2 = ZTOS(n2);
484: else
485: c2 = remzi(n2,c1>0?c1:-c1);
486: if ( c2 < 0 ) c2 = -c2;
487: g = gcdimm(c1,c2);
488: STOZ(g,gcd);
489: return gcd;
490: } else if ( IS_IMM(n2) ) {
491: c2 = ZTOS(n2);
492: if ( c2 < 0 ) c2 = -c2;
493: c1 = remzi(n1,c2>0?c2:-c2);
494: if ( c1 < 0 ) c1 = -c1;
495: g = gcdimm(c1,c2);
496: STOZ(g,gcd);
497: return gcd;
498: } else {
499: n1 = dupz(n1);
500: if ( SL(n1) < 0 ) SL(n1) = -SL(n1);
501: n2 = dupz(n2);
502: if ( SL(n2) < 0 ) SL(n2) = -SL(n2);
503: gcdn((N)n1,(N)n2,&gn);
504: gcd = (Z)gn;
505: if ( IS_SZ(gcd) ) {
506: SZTOZ(gcd,r); gcd = r;
507: }
508: return gcd;
509: }
1.1 noro 510: }
511:
512: int remzi(Z n,int m)
513: {
514: unsigned int *x;
515: unsigned int t,r;
1.5 noro 516: int i,c;
1.1 noro 517:
518: if ( !n ) return 0;
1.5 noro 519: if ( IS_IMM(n) ) {
520: c = ZTOS(n)%m;
521: if ( c < 0 ) c += m;
522: return c;
523: }
524:
1.1 noro 525: i = SL(n);
526: if ( i < 0 ) i = -i;
527: for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
528: #if defined(sparc)
529: r = dsa(m,r,*x);
530: #else
531: DSAB(m,r,*x,t,r);
532: #endif
533: }
1.5 noro 534: if ( r && SL(n) < 0 )
535: r = m-r;
1.1 noro 536: return r;
537: }
538:
539: Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)
540: {
541: Z gcd;
542:
543: gcd = gcdz(n1,n2);
544: *c1 = divsz(n1,gcd);
545: *c2 = divsz(n2,gcd);
546: return gcd;
547: }
548:
1.5 noro 549: #if 0
1.1 noro 550: Z estimate_array_gcdz(Z *b,int n)
551: {
552: int m,i,j,sd;
553: Z *a;
554: Z s,t;
555:
556: a = (Z *)ALLOCA(n*sizeof(Z));
557: for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
558: n = j;
559: if ( !n ) return 0;
560: if ( n == 1 ) return a[0];
561:
562: m = n/2;
563: for ( i = 0, s = 0; i < m; i++ ) {
564: if ( !a[i] ) continue;
565: else s = (SL(a[i])<0)?subz(s,a[i]):addz(s,a[i]);
566: }
567: for ( t = 0; i < n; i++ ) {
568: if ( !a[i] ) continue;
569: else t = (SL(a[i])<0)?subz(t,a[i]):addz(t,a[i]);
570: }
571: return gcdz(s,t);
572: }
573:
574: Z array_gcdz(Z *b,int n)
575: {
576: int m,i,j,sd;
577: Z *a;
578: Z gcd;
579:
580: a = (Z *)ALLOCA(n*sizeof(Z));
581: for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
582: n = j;
583: if ( !n ) return 0;
584: if ( n == 1 ) return a[0];
585: gcd = a[0];
586: for ( i = 1; i < n; i++ )
587: gcd = gcdz(gcd,a[i]);
588: return gcd;
589: }
1.4 noro 590: #endif
1.1 noro 591:
592: void _copyz(Z n1,Z n2)
593: {
594: int n,i;
595:
596: if ( !n1 || !SL(n1) ) SL(n2) = 0;
597: else {
598: n = SL(n2) = SL(n1);
599: if ( n < 0 ) n = -n;
600: for ( i = 0; i < n; i++ ) BD(n2)[i] = BD(n1)[i];
601: }
602: }
603:
604: void _addz(Z n1,Z n2,Z nr)
605: {
1.3 noro 606: int d1,d2;
1.1 noro 607:
608: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
609: else if ( !n2 || !SL(n2) ) _copyz(n1,nr);
1.3 noro 610: else if ( (d1=SL(n1)) > 0 )
611: if ( (d2=SL(n2)) > 0 )
612: SL(nr) = _addz_main(BD(n1),d1,BD(n2),d2,BD(nr));
1.1 noro 613: else
1.3 noro 614: SL(nr) = _subz_main(BD(n1),d1,BD(n2),-d2,BD(nr));
615: else if ( (d2=SL(n2)) > 0 )
616: SL(nr) = _subz_main(BD(n2),d2,BD(n1),-d1,BD(nr));
1.1 noro 617: else
1.3 noro 618: SL(nr) = -_addz_main(BD(n1),-d1,BD(n2),-d2,BD(nr));
1.1 noro 619: }
620:
621: void _subz(Z n1,Z n2,Z nr)
622: {
1.3 noro 623: int d1,d2;
1.1 noro 624:
625: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
626: else if ( !n2 || !SL(n2) ) {
627: _copyz(n1,nr);
628: SL(nr) = -SL(nr);
1.3 noro 629: } else if ( (d1=SL(n1)) > 0 )
630: if ( (d2=SL(n2)) > 0 )
631: SL(nr) = _subz_main(BD(n1),d1,BD(n2),d2,BD(nr));
1.1 noro 632: else
1.3 noro 633: SL(nr) = _addz_main(BD(n1),d1,BD(n2),-d2,BD(nr));
634: else if ( (d2=SL(n2)) > 0 )
635: SL(nr) = -_addz_main(BD(n1),-d1,BD(n2),d2,BD(nr));
1.1 noro 636: else
1.3 noro 637: SL(nr) = -_subz_main(BD(n1),-d1,BD(n2),-d2,BD(nr));
1.1 noro 638: }
639:
640: void _mulz(Z n1,Z n2,Z nr)
641: {
1.3 noro 642: int d1,d2;
1.1 noro 643: int i,d,sgn;
644: unsigned int mul;
645: unsigned int *m1,*m2;
646:
647: if ( !n1 || !SL(n1) || !n2 || !SL(n2) )
648: SL(nr) = 0;
649: else {
1.3 noro 650: d1 = SL(n1); d2 = SL(n2);
1.1 noro 651: sgn = 1;
1.3 noro 652: if ( d1 < 0 ) { sgn = -sgn; d1 = -d1; }
653: if ( d2 < 0 ) { sgn = -sgn; d2 = -d2; }
654: d = d1+d2;
1.1 noro 655: for ( i = d-1, m1 = BD(nr); i >= 0; i-- ) *m1++ = 0;
1.3 noro 656: for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
657: if ( mul = *m2 ) muln_1(m1,d1,mul,BD(nr)+i);
1.1 noro 658: SL(nr) = sgn*(BD(nr)[d-1]?d:d-1);
659: }
660: }
661:
662: int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
663: {
664: int d,i;
665: unsigned int tmp,c;
666: unsigned int *t;
667:
668: if ( d2 > d1 ) {
669: t = m1; m1 = m2; m2 = t;
670: d = d1; d1 = d2; d2 = d;
671: }
672: #if defined(VISUAL)
673: __asm {
674: push esi
675: push edi
676: mov esi,m1
677: mov edi,m2
678: mov ebx,mr
679: mov ecx,d2
680: xor eax,eax
681: Lstart__addz:
682: mov eax,DWORD PTR [esi]
683: mov edx,DWORD PTR [edi]
684: adc eax,edx
685: mov DWORD PTR [ebx],eax
686: lea esi,DWORD PTR [esi+4]
687: lea edi,DWORD PTR [edi+4]
688: lea ebx,DWORD PTR [ebx+4]
689: dec ecx
690: jnz Lstart__addz
691: pop edi
692: pop esi
693: mov eax,0
694: adc eax,eax
695: mov c,eax
696: }
697: #elif defined(i386)
698: asm volatile("\
699: movl %1,%%esi;\
700: movl %2,%%edi;\
701: movl %3,%%ebx;\
702: movl %4,%%ecx;\
703: testl %%eax,%%eax;\
704: Lstart__addz:;\
705: movl (%%esi),%%eax;\
706: movl (%%edi),%%edx;\
707: adcl %%edx,%%eax;\
708: movl %%eax,(%%ebx);\
709: leal 4(%%esi),%%esi;\
710: leal 4(%%edi),%%edi;\
711: leal 4(%%ebx),%%ebx;\
712: decl %%ecx;\
713: jnz Lstart__addz;\
714: movl $0,%%eax;\
715: adcl %%eax,%%eax;\
716: movl %%eax,%0"\
717: :"=m"(c)\
718: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
719: :"eax","ebx","ecx","edx","esi","edi");
720: #else
1.6 ! noro 721: for ( i = 0, c = 0; i < d2; i++, m1++, m2++, mr++ ) {
1.1 noro 722: tmp = *m1 + *m2;
723: if ( tmp < *m1 ) {
724: tmp += c;
725: c = 1;
726: } else {
727: tmp += c;
728: c = tmp < c ? 1 : 0;
729: }
730: *mr = tmp;
731: }
732: #endif
733: for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
734: tmp = *m1++ + c;
735: c = tmp < c ? 1 : 0;
736: *mr++ = tmp;
737: }
738: for ( ; i < d1; i++ )
739: *mr++ = *m1++;
740: *mr = c;
741: return (c?d1+1:d1);
742: }
743:
744: int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
745: {
746: int d,i,sgn;
747: unsigned int t,tmp,br;
748: unsigned int *m;
749:
750: if ( d1 > d2 ) sgn = 1;
751: else if ( d1 < d2 ) sgn = -1;
752: else {
753: for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
754: if ( i < 0 ) return 0;
755: if ( m1[i] > m2[i] ) sgn = 1;
756: else if ( m1[i] < m2[i] ) sgn = -1;
757: }
758: if ( sgn < 0 ) {
759: m = m1; m1 = m2; m2 = m;
760: d = d1; d1 = d2; d2 = d;
761: }
762: #if defined(VISUAL)
763: __asm {
764: push esi
765: push edi
766: mov esi,m1
767: mov edi,m2
768: mov ebx,mr
769: mov ecx,d2
770: xor eax,eax
771: Lstart__subz:
772: mov eax,DWORD PTR [esi]
773: mov edx,DWORD PTR [edi]
774: sbb eax,edx
775: mov DWORD PTR [ebx],eax
776: lea esi,DWORD PTR [esi+4]
777: lea edi,DWORD PTR [edi+4]
778: lea ebx,DWORD PTR [ebx+4]
779: dec ecx
780: jnz Lstart__subz
781: pop edi
782: pop esi
783: mov eax,0
784: adc eax,eax
785: mov br,eax
786: }
787: #elif defined(i386)
788: asm volatile("\
789: movl %1,%%esi;\
790: movl %2,%%edi;\
791: movl %3,%%ebx;\
792: movl %4,%%ecx;\
793: testl %%eax,%%eax;\
794: Lstart__subz:;\
795: movl (%%esi),%%eax;\
796: movl (%%edi),%%edx;\
797: sbbl %%edx,%%eax;\
798: movl %%eax,(%%ebx);\
799: leal 4(%%esi),%%esi;\
800: leal 4(%%edi),%%edi;\
801: leal 4(%%ebx),%%ebx;\
802: decl %%ecx;\
803: jnz Lstart__subz;\
804: movl $0,%%eax;\
805: adcl %%eax,%%eax;\
806: movl %%eax,%0"\
807: :"=m"(br)\
808: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
809: :"eax","ebx","ecx","edx","esi","edi");
810: #else
1.6 ! noro 811: for ( i = 0, br = 0; i < d2; i++, mr++ ) {
1.1 noro 812: t = *m1++;
813: tmp = *m2++ + br;
814: if ( br > 0 && !tmp ) {
815: /* tmp = 2^32 => br = 1 */
816: }else {
817: tmp = t-tmp;
818: br = tmp > t ? 1 : 0;
819: *mr = tmp;
820: }
821: }
822: #endif
823: for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
824: t = *m1++;
825: tmp = t - br;
826: br = tmp > t ? 1 : 0;
827: *mr++ = tmp;
828: }
829: for ( ; i < d1; i++ )
830: *mr++ = *m1++;
831: for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
832: return sgn*(i+1);
833: }
834:
835: /* XXX */
836:
837: void printz(Z n)
838: {
1.4 noro 839: int sd,u;
1.1 noro 840:
841: if ( !n )
842: fprintf(asir_out,"0");
1.4 noro 843: else if ( IS_IMM(n) ) {
1.5 noro 844: u = ZTOS(n);
1.4 noro 845: fprintf(asir_out,"%d",u);
846: } else {
1.1 noro 847: if ( (sd = SL(n)) < 0 ) { SL(n) = -SL(n); fprintf(asir_out,"-"); }
848: printn((N)n);
849: if ( sd < 0 ) SL(n) = -SL(n);
1.2 noro 850: }
851: }
852:
853: /*
854: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
855: *
856: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
857: * where W(k,l,i) = i! * kCi * lCi
858: */
859:
860: void mkwcz(int k,int l,Z *t)
861: {
862: int i,n,up,low;
863: N nm,d,c;
864:
865: n = MIN(k,l);
866: for ( t[0] = (Z)ONEN, i = 1; i <= n; i++ ) {
867: DM(k-i+1,l-i+1,up,low);
868: if ( up ) {
869: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
870: } else {
871: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
872: }
873: kmuln((N)t[i-1],nm,&d); divin(d,i,&c); t[i] = (Z)c;
1.1 noro 874: }
875: }