![[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.2 ! noro 1: #if 0
1.1 noro 2: #include "ca.h"
3: #include "inline.h"
1.2 ! noro 4: #endif
1.1 noro 5:
6: typedef struct oZ {
7: int p;
8: unsigned int b[1];
9: } *Z;
10:
11: #define ZALLOC(d) ((Z)MALLOC_ATOMIC(TRUESIZE(oZ,(d)-1,int)))
12: #define SL(n) ((n)->p)
13:
14: Z qtoz(Q n);
15: Q ztoq(Z n);
16: Z chsgnz(Z n);
17: Z dupz(Z n);
18: Z addz(Z n1,Z n2);
19: Z subz(Z n1,Z n2);
20: Z mulz(Z n1,Z n2);
21: Z divsz(Z n1,Z n2);
22: Z divz(Z n1,Z n2,Z *rem);
23: Z gcdz(Z n1,Z n2);
24: Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2);
25: Z estimate_array_gcdz(Z *a,int n);
26: Z array_gcdz(Z *a,int n);
1.2 ! noro 27: void mkwcz(int k,int l,Z *t);
1.1 noro 28: int remzi(Z n,int m);
29: inline void _addz(Z n1,Z n2,Z nr);
30: inline void _subz(Z n1,Z n2,Z nr);
31: inline void _mulz(Z n1,Z n2,Z nr);
32: inline int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
33: inline int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
34:
1.2 ! noro 35: sgnz(Z n)
! 36: {
! 37: if ( !n ) return 0;
! 38: else if ( SL(n) < 0 ) return -1;
! 39: else return 1;
! 40: }
! 41:
1.1 noro 42: z_mag(Z n)
43: {
44: return n_bits((N)n);
45: }
46:
47: Z qtoz(Q n)
48: {
49: Z r;
50:
51: if ( !n ) return 0;
52: else if ( !INT(n) )
53: error("qtoz : invalid input");
54: else {
55: r = dupz((Z)NM(n));
56: if ( SGN(n) < 0 ) SL(r) = -SL(r);
57: return r;
58: }
59: }
60:
61: Q ztoq(Z n)
62: {
63: Q r;
64: Z nm;
65: int sgn;
66:
67: if ( !n ) return 0;
68: else {
69: nm = dupz(n);
70: if ( SL(nm) < 0 ) {
71: sgn = -1;
72: SL(nm) = -SL(nm);
73: } else
74: sgn = 1;
75: NTOQ((N)nm,sgn,r);
76: return r;
77: }
78: }
79:
80: Z dupz(Z n)
81: {
82: Z nr;
83: int sd,i;
84:
85: if ( !n ) return 0;
86: if ( (sd = SL(n)) < 0 ) sd = -sd;
87: nr = ZALLOC(sd);
88: SL(nr) = SL(n);
89: for ( i = 0; i < sd; i++ ) BD(nr)[i] = BD(n)[i];
90: return nr;
91: }
92:
93: Z chsgnz(Z n)
94: {
95: Z r;
96:
97: if ( !n ) return 0;
98: else {
99: r = dupz(n);
100: SL(r) = -SL(r);
101: return r;
102: }
103: }
104:
105: Z addz(Z n1,Z n2)
106: {
107: unsigned int u1,u2,t;
108: int sd1,sd2,d,d1,d2;
109: Z nr;
110:
111: if ( !n1 ) return dupz(n2);
112: else if ( !n2 ) return dupz(n1);
113: else {
114: d1 = ((sd1 = SL(n1))< 0)?-sd1:sd1;
115: d2 = ((sd2 = SL(n2))< 0)?-sd2:sd2;
116: if ( d1 == 1 && d2 == 1 ) {
117: u1 = BD(n1)[0]; u2 = BD(n2)[0];
118: if ( sd1 > 0 ) {
119: if ( sd2 > 0 ) {
120: t = u1+u2;
121: if ( t < u1 ) {
122: nr=ZALLOC(2); SL(nr) = 2; BD(nr)[1] = 1;
123: } else {
124: nr=ZALLOC(1); SL(nr) = 1;
125: }
126: BD(nr)[0] = t;
127: } else {
128: if ( u1 == u2 ) nr = 0;
129: else if ( u1 > u2 ) {
130: nr=ZALLOC(1); SL(nr) = 1; BD(nr)[0] = u1-u2;
131: } else {
132: nr=ZALLOC(1); SL(nr) = -1; BD(nr)[0] = u2-u1;
133: }
134: }
135: } else {
136: if ( sd2 > 0 ) {
137: if ( u2 == u1 ) nr = 0;
138: else if ( u2 > u1 ) {
139: nr=ZALLOC(1); SL(nr) = 1; BD(nr)[0] = u2-u1;
140: } else {
141: nr=ZALLOC(1); SL(nr) = -1; BD(nr)[0] = u1-u2;
142: }
143: } else {
144: t = u1+u2;
145: if ( t < u1 ) {
146: nr=ZALLOC(2); SL(nr) = -2; BD(nr)[1] = 1;
147: } else {
148: nr=ZALLOC(1); SL(nr) = -1;
149: }
150: BD(nr)[0] = t;
151: }
152: }
153: } else {
154: d = MAX(d1,d2)+1;
155: nr = ZALLOC(d);
156: _addz(n1,n2,nr);
157: if ( !SL(nr) ) nr = 0;
158: }
159: return nr;
160: }
161: }
162:
163: Z subz(Z n1,Z n2)
164: {
165: int sd1,sd2,d;
166: Z nr;
167:
168: if ( !n1 )
169: if ( !n2 ) return 0;
170: else {
171: nr = dupz(n2);
172: SL(nr) = -SL(nr);
173: }
174: else if ( !n2 ) return dupz(n1);
175: else {
176: if ( (sd1 = SL(n1)) < 0 ) sd1 = -sd1;
177: if ( (sd2 = SL(n2)) < 0 ) sd2 = -sd2;
178: d = MAX(sd1,sd2)+1;
179: nr = ZALLOC(d);
180: _subz(n1,n2,nr);
181: if ( !SL(nr) ) nr = 0;
182: return nr;
183: }
184: }
185:
186: Z mulz(Z n1,Z n2)
187: {
188: int sd1,sd2,d1,d2;
189: unsigned int u1,u2,u,l;
190: Z nr;
191:
192: if ( !n1 || !n2 ) return 0;
193: else {
194: d1 = ((sd1 = SL(n1))< 0)?-sd1:sd1;
195: d2 = ((sd2 = SL(n2))< 0)?-sd2:sd2;
196: if ( d1 == 1 && d2 == 1 ) {
197: u1 = BD(n1)[0]; u2 = BD(n2)[0];
198: DM(u1,u2,u,l);
199: if ( u ) {
200: nr = ZALLOC(2); SL(nr) = sd1*sd2*2; BD(nr)[1] = u;
201: } else {
202: nr = ZALLOC(1); SL(nr) = sd1*sd2;
203: }
204: BD(nr)[0] = l;
205: } else {
206: nr = ZALLOC(d1+d2);
207: _mulz(n1,n2,nr);
208: }
209: return nr;
210: }
211: }
212:
213: /* XXX */
214: Z divz(Z n1,Z n2,Z *rem)
215: {
216: int sgn,sd1,sd2;
217: N q,r;
218:
219: if ( !n2 ) error("divz : division by 0");
220: else if ( !n1 ) {
221: *rem = 0; return 0;
222: } else {
223: sd2 = SL(n2);
224: if ( sd2 == 1 && BD(n2)[0] == 1 ) {
225: *rem = 0; return n1;
226: } else if ( sd2 == -1 && BD(n2)[0] == 1 ) {
227: *rem = 0; return chsgnz(n1);
228: }
229:
230: sd1 = SL(n1);
231: n1 = dupz(n1);
1.2 ! noro 232: if ( sd1 < 0 ) SL(n1) = -sd1;
! 233: if ( sd2 < 0 ) SL(n2) = -sd2;
1.1 noro 234: divn((N)n1,(N)n2,&q,&r);
235: if ( q && ((sd1>0&&sd2<0)||(sd1<0&&sd2>0)) ) SL((Z)q) = -SL((Z)q);
236: if ( r && sd1 < 0 ) SL((Z)r) = -SL((Z)r);
1.2 ! noro 237: SL(n2) = sd2;
1.1 noro 238: *rem = (Z)r;
239: return (Z)q;
240: }
241: }
242:
243: Z divsz(Z n1,Z n2)
244: {
245: int sgn,sd1,sd2;
246: N q;
247:
248: if ( !n2 ) error("divsz : division by 0");
249: else if ( !n1 ) return 0;
250: else {
251: sd2 = SL(n2);
252: if ( sd2 == 1 && BD(n2)[0] == 1 ) return n1;
253: else if ( sd2 == -1 && BD(n2)[0] == 1 ) return chsgnz(n1);
254:
255: sd1 = SL(n1);
256: n1 = dupz(n1);
1.2 ! noro 257: if ( sd1 < 0 ) SL(n1) = -sd1;
! 258: if ( sd2 < 0 ) SL(n2) = -sd2;
1.1 noro 259: divsn((N)n1,(N)n2,&q);
260: if ( q && ((sd1>0&&sd2<0)||(sd1<0&&sd2>0)) ) SL((Z)q) = -SL((Z)q);
1.2 ! noro 261: SL(n2) = sd2;
1.1 noro 262: return (Z)q;
263: }
264: }
265:
266: Z gcdz(Z n1,Z n2)
267: {
268: int sd1,sd2;
269: N gcd;
270:
271: if ( !n1 ) return n2;
272: else if ( !n2 ) return n1;
273:
274: n1 = dupz(n1);
275: if ( SL(n1) < 0 ) SL(n1) = -SL(n1);
276: n2 = dupz(n2);
277: if ( SL(n2) < 0 ) SL(n2) = -SL(n2);
278: gcdn((N)n1,(N)n2,&gcd);
279: return (Z)gcd;
280: }
281:
282: int remzi(Z n,int m)
283: {
284: unsigned int *x;
285: unsigned int t,r;
286: int i;
287:
288: if ( !n ) return 0;
289: i = SL(n);
290: if ( i < 0 ) i = -i;
291: for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
292: #if defined(sparc)
293: r = dsa(m,r,*x);
294: #else
295: DSAB(m,r,*x,t,r);
296: #endif
297: }
298: return r;
299: }
300:
301: Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)
302: {
303: Z gcd;
304:
305: gcd = gcdz(n1,n2);
306: *c1 = divsz(n1,gcd);
307: *c2 = divsz(n2,gcd);
308: return gcd;
309: }
310:
311: Z estimate_array_gcdz(Z *b,int n)
312: {
313: int m,i,j,sd;
314: Z *a;
315: Z s,t;
316:
317: a = (Z *)ALLOCA(n*sizeof(Z));
318: for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
319: n = j;
320: if ( !n ) return 0;
321: if ( n == 1 ) return a[0];
322:
323: m = n/2;
324: for ( i = 0, s = 0; i < m; i++ ) {
325: if ( !a[i] ) continue;
326: else s = (SL(a[i])<0)?subz(s,a[i]):addz(s,a[i]);
327: }
328: for ( t = 0; i < n; i++ ) {
329: if ( !a[i] ) continue;
330: else t = (SL(a[i])<0)?subz(t,a[i]):addz(t,a[i]);
331: }
332: return gcdz(s,t);
333: }
334:
335: Z array_gcdz(Z *b,int n)
336: {
337: int m,i,j,sd;
338: Z *a;
339: Z gcd;
340:
341: a = (Z *)ALLOCA(n*sizeof(Z));
342: for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
343: n = j;
344: if ( !n ) return 0;
345: if ( n == 1 ) return a[0];
346: gcd = a[0];
347: for ( i = 1; i < n; i++ )
348: gcd = gcdz(gcd,a[i]);
349: return gcd;
350: }
351:
352: void _copyz(Z n1,Z n2)
353: {
354: int n,i;
355:
356: if ( !n1 || !SL(n1) ) SL(n2) = 0;
357: else {
358: n = SL(n2) = SL(n1);
359: if ( n < 0 ) n = -n;
360: for ( i = 0; i < n; i++ ) BD(n2)[i] = BD(n1)[i];
361: }
362: }
363:
364: void _addz(Z n1,Z n2,Z nr)
365: {
366: int sd1,sd2;
367:
368: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
369: else if ( !n2 || !SL(n2) ) _copyz(n1,nr);
370: else if ( (sd1=SL(n1)) > 0 )
371: if ( (sd2=SL(n2)) > 0 )
372: SL(nr) = _addz_main(BD(n1),sd1,BD(n2),sd2,BD(nr));
373: else
374: SL(nr) = _subz_main(BD(n1),sd1,BD(n2),-sd2,BD(nr));
375: else if ( (sd2=SL(n2)) > 0 )
376: SL(nr) = _subz_main(BD(n2),sd2,BD(n1),-sd1,BD(nr));
377: else
378: SL(nr) = -_addz_main(BD(n1),-sd1,BD(n2),-sd2,BD(nr));
379: }
380:
381: void _subz(Z n1,Z n2,Z nr)
382: {
383: int sd1,sd2;
384:
385: if ( !n1 || !SL(n1) ) _copyz(n2,nr);
386: else if ( !n2 || !SL(n2) ) {
387: _copyz(n1,nr);
388: SL(nr) = -SL(nr);
389: } else if ( (sd1=SL(n1)) > 0 )
390: if ( (sd2=SL(n2)) > 0 )
391: SL(nr) = _subz_main(BD(n1),sd1,BD(n2),sd2,BD(nr));
392: else
393: SL(nr) = _addz_main(BD(n1),sd1,BD(n2),-sd2,BD(nr));
394: else if ( (sd2=SL(n2)) > 0 )
395: SL(nr) = -_addz_main(BD(n1),-sd1,BD(n2),sd2,BD(nr));
396: else
397: SL(nr) = -_subz_main(BD(n1),-sd1,BD(n2),-sd2,BD(nr));
398: }
399:
400: void _mulz(Z n1,Z n2,Z nr)
401: {
402: int sd1,sd2;
403: int i,d,sgn;
404: unsigned int mul;
405: unsigned int *m1,*m2;
406:
407: if ( !n1 || !SL(n1) || !n2 || !SL(n2) )
408: SL(nr) = 0;
409: else {
410: sd1 = SL(n1); sd2 = SL(n2);
411: sgn = 1;
412: if ( sd1 < 0 ) { sgn = -sgn; sd1 = -sd1; }
413: if ( sd2 < 0 ) { sgn = -sgn; sd2 = -sd2; }
414: d = sd1+sd2;
415: for ( i = d-1, m1 = BD(nr); i >= 0; i-- ) *m1++ = 0;
416: for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < sd2; i++, m2++ )
417: if ( mul = *m2 ) muln_1(m1,sd1,mul,BD(nr)+i);
418: SL(nr) = sgn*(BD(nr)[d-1]?d:d-1);
419: }
420: }
421:
422: int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
423: {
424: int d,i;
425: unsigned int tmp,c;
426: unsigned int *t;
427:
428: if ( d2 > d1 ) {
429: t = m1; m1 = m2; m2 = t;
430: d = d1; d1 = d2; d2 = d;
431: }
432: #if defined(VISUAL)
433: __asm {
434: push esi
435: push edi
436: mov esi,m1
437: mov edi,m2
438: mov ebx,mr
439: mov ecx,d2
440: xor eax,eax
441: Lstart__addz:
442: mov eax,DWORD PTR [esi]
443: mov edx,DWORD PTR [edi]
444: adc eax,edx
445: mov DWORD PTR [ebx],eax
446: lea esi,DWORD PTR [esi+4]
447: lea edi,DWORD PTR [edi+4]
448: lea ebx,DWORD PTR [ebx+4]
449: dec ecx
450: jnz Lstart__addz
451: pop edi
452: pop esi
453: mov eax,0
454: adc eax,eax
455: mov c,eax
456: }
457: #elif defined(i386)
458: asm volatile("\
459: movl %1,%%esi;\
460: movl %2,%%edi;\
461: movl %3,%%ebx;\
462: movl %4,%%ecx;\
463: testl %%eax,%%eax;\
464: Lstart__addz:;\
465: movl (%%esi),%%eax;\
466: movl (%%edi),%%edx;\
467: adcl %%edx,%%eax;\
468: movl %%eax,(%%ebx);\
469: leal 4(%%esi),%%esi;\
470: leal 4(%%edi),%%edi;\
471: leal 4(%%ebx),%%ebx;\
472: decl %%ecx;\
473: jnz Lstart__addz;\
474: movl $0,%%eax;\
475: adcl %%eax,%%eax;\
476: movl %%eax,%0"\
477: :"=m"(c)\
478: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
479: :"eax","ebx","ecx","edx","esi","edi");
480: #else
481: for ( i = 0, c = 0, mr = BD(nr); i < d2; i++, m1++, m2++, mr++ ) {
482: tmp = *m1 + *m2;
483: if ( tmp < *m1 ) {
484: tmp += c;
485: c = 1;
486: } else {
487: tmp += c;
488: c = tmp < c ? 1 : 0;
489: }
490: *mr = tmp;
491: }
492: #endif
493: for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
494: tmp = *m1++ + c;
495: c = tmp < c ? 1 : 0;
496: *mr++ = tmp;
497: }
498: for ( ; i < d1; i++ )
499: *mr++ = *m1++;
500: *mr = c;
501: return (c?d1+1:d1);
502: }
503:
504: int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
505: {
506: int d,i,sgn;
507: unsigned int t,tmp,br;
508: unsigned int *m;
509:
510: if ( d1 > d2 ) sgn = 1;
511: else if ( d1 < d2 ) sgn = -1;
512: else {
513: for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
514: if ( i < 0 ) return 0;
515: if ( m1[i] > m2[i] ) sgn = 1;
516: else if ( m1[i] < m2[i] ) sgn = -1;
517: }
518: if ( sgn < 0 ) {
519: m = m1; m1 = m2; m2 = m;
520: d = d1; d1 = d2; d2 = d;
521: }
522: #if defined(VISUAL)
523: __asm {
524: push esi
525: push edi
526: mov esi,m1
527: mov edi,m2
528: mov ebx,mr
529: mov ecx,d2
530: xor eax,eax
531: Lstart__subz:
532: mov eax,DWORD PTR [esi]
533: mov edx,DWORD PTR [edi]
534: sbb eax,edx
535: mov DWORD PTR [ebx],eax
536: lea esi,DWORD PTR [esi+4]
537: lea edi,DWORD PTR [edi+4]
538: lea ebx,DWORD PTR [ebx+4]
539: dec ecx
540: jnz Lstart__subz
541: pop edi
542: pop esi
543: mov eax,0
544: adc eax,eax
545: mov br,eax
546: }
547: #elif defined(i386)
548: asm volatile("\
549: movl %1,%%esi;\
550: movl %2,%%edi;\
551: movl %3,%%ebx;\
552: movl %4,%%ecx;\
553: testl %%eax,%%eax;\
554: Lstart__subz:;\
555: movl (%%esi),%%eax;\
556: movl (%%edi),%%edx;\
557: sbbl %%edx,%%eax;\
558: movl %%eax,(%%ebx);\
559: leal 4(%%esi),%%esi;\
560: leal 4(%%edi),%%edi;\
561: leal 4(%%ebx),%%ebx;\
562: decl %%ecx;\
563: jnz Lstart__subz;\
564: movl $0,%%eax;\
565: adcl %%eax,%%eax;\
566: movl %%eax,%0"\
567: :"=m"(br)\
568: :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
569: :"eax","ebx","ecx","edx","esi","edi");
570: #else
571: for ( i = 0, br = 0, mr = BD(nr); i < d2; i++, mr++ ) {
572: t = *m1++;
573: tmp = *m2++ + br;
574: if ( br > 0 && !tmp ) {
575: /* tmp = 2^32 => br = 1 */
576: }else {
577: tmp = t-tmp;
578: br = tmp > t ? 1 : 0;
579: *mr = tmp;
580: }
581: }
582: #endif
583: for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
584: t = *m1++;
585: tmp = t - br;
586: br = tmp > t ? 1 : 0;
587: *mr++ = tmp;
588: }
589: for ( ; i < d1; i++ )
590: *mr++ = *m1++;
591: for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
592: return sgn*(i+1);
593: }
594:
595: /* XXX */
596:
597: void printz(Z n)
598: {
599: int sd;
600:
601: if ( !n )
602: fprintf(asir_out,"0");
603: else {
604: if ( (sd = SL(n)) < 0 ) { SL(n) = -SL(n); fprintf(asir_out,"-"); }
605: printn((N)n);
606: if ( sd < 0 ) SL(n) = -SL(n);
1.2 ! noro 607: }
! 608: }
! 609:
! 610: /*
! 611: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
! 612: *
! 613: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
! 614: * where W(k,l,i) = i! * kCi * lCi
! 615: */
! 616:
! 617: void mkwcz(int k,int l,Z *t)
! 618: {
! 619: int i,n,up,low;
! 620: N nm,d,c;
! 621:
! 622: n = MIN(k,l);
! 623: for ( t[0] = (Z)ONEN, i = 1; i <= n; i++ ) {
! 624: DM(k-i+1,l-i+1,up,low);
! 625: if ( up ) {
! 626: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
! 627: } else {
! 628: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
! 629: }
! 630: kmuln((N)t[i-1],nm,&d); divin(d,i,&c); t[i] = (Z)c;
1.1 noro 631: }
632: }