![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ddN.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / asm|
Return to ddN.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / asm |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/asm/ddN.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */
2: #ifndef FBASE
3: #define FBASE
4: #endif
5:
6: #include "ca.h"
7: #include "base.h"
8: #include "inline.h"
9:
10: void bxprintn(N),bprintn(N);
11: void divnmain_special(int,int,unsigned int *,unsigned int *,unsigned int *);
12: void dupn(N,N);
13:
14: void divn(n1,n2,nq,nr)
15: N n1,n2,*nq,*nr;
16: {
17: int tmp,b;
18: int i,j,d1,d2,dd;
19: unsigned int *m1,*m2;
20: unsigned int uq,ur,msw;
21: N q,r,t1,t2;
22:
23: if ( !n2 )
24: error("divn: divsion by 0");
25: else if ( !n1 ) {
26: *nq = 0; *nr = 0;
27: } else if ( UNIN(n2) ) {
28: COPY(n1,*nq); *nr = 0;
29: } else if ( (d1 = PL(n1)) < (d2 = PL(n2)) ) {
30: COPY(n1,*nr); *nq = 0;
31: } else if ( d2 == 1 ) {
32: if ( d1 == 1 ) {
33: DQR(BD(n1)[0],BD(n2)[0],uq,ur)
34: STON(uq,*nq); STON(ur,*nr);
35: } else {
36: ur = divin(n1,BD(n2)[0],nq); STON(ur,*nr);
37: }
38: return;
39: } else if ( (d1 == d2) && ((tmp = cmpn(n1,n2)) <= 0) )
40: if ( !tmp ) {
41: *nr = 0; COPY(ONEN,*nq);
42: } else {
43: COPY(n1,*nr); *nq = 0;
44: }
45: else {
46: msw = BD(n2)[d2-1];
47: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
48: b = BSH-i;
49: if ( b ) {
50: bshiftn(n1,-b,&t1); bshiftn(n2,-b,&t2);
51: } else {
52: COPY(n1,t1); COPY(n2,t2);
53: }
54: m1 = (unsigned int *)ALLOCA((PL(t1)+1)*sizeof(unsigned int));
55: bcopy(BD(t1),m1,PL(t1)*sizeof(unsigned int)); m1[PL(t1)] = 0;
56: m2 = BD(t2); dd = d1-d2;
57: *nq = q = NALLOC(dd+1); INITRC(q);
58: divnmain(d1,d2,m1,m2,BD(q)); FREEN(t1); FREEN(t2);
59: PL(q) = (BD(q)[dd]?dd+1:dd);
60: if ( !PL(q) ) {
61: FREEN(q);
62: }
63: for ( j = d2-1; (j >= 0) && (m1[j] == 0); j--);
64: if ( j == -1 )
65: *nr = 0;
66: else {
67: r = NALLOC(j+1); INITRC(r); PL(r) = j+1;
68: bcopy(m1,BD(r),(j+1)*sizeof(unsigned int));
69: if ( b ) {
70: bshiftn(r,b,nr); FREEN(r);
71: } else
72: *nr = r;
73: }
74: }
75: }
76:
77: void divsn(n1,n2,nq)
78: N n1,n2,*nq;
79: {
80: int d1,d2,dd,i,b;
81: unsigned int *m1,*m2;
82: unsigned int uq,msw;
83: N q,t1,t2;
84:
85: if ( !n1 )
86: *nq = 0;
87: else if ( UNIN(n2) ) {
88: COPY(n1,*nq);
89: } else if ( (d1 = PL(n1)) < (d2 = PL(n2)) )
90: error("divsn: cannot happen");
91: else if ( d2 == 1 ) {
92: if ( d1 == 1 ) {
93: uq = BD(n1)[0] / BD(n2)[0]; STON(uq,*nq);
94: } else
95: divin(n1,BD(n2)[0],nq);
96: } else {
97: msw = BD(n2)[d2-1];
98: for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
99: b = BSH-i;
100: if ( b ) {
101: bshiftn(n1,-b,&t1); bshiftn(n2,-b,&t2);
102: } else {
103: COPY(n1,t1); COPY(n2,t2);
104: }
105: m1 = (unsigned int *)ALLOCA((PL(t1)+1)*sizeof(unsigned int));
106: bcopy(BD(t1),m1,PL(t1)*sizeof(unsigned int)); m1[PL(t1)] = 0;
107: m2 = BD(t2); dd = d1-d2;
108: *nq = q = NALLOC(dd+1); INITRC(q);
109: divnmain(d1,d2,m1,m2,BD(q));
110: FREEN(t1); FREEN(t2);
111: PL(q) = (BD(q)[dd]?dd+1:dd);
112: if ( !PL(q) )
113: FREEN(q);
114: }
115: }
116:
117: void remn(n1,n2,nr)
118: N n1,n2,*nr;
119: {
120: int d1,d2,tmp;
121: unsigned int uq,ur;
122: N q;
123:
124: if ( !n2 )
125: error("remn: divsion by 0");
126: else if ( !n1 || UNIN(n2) )
127: *nr = 0;
128: else if ( (d1 = PL(n1)) < (d2 = PL(n2)) ) {
129: COPY(n1,*nr);
130: } else if ( d2 == 1 ) {
131: if ( d1 == 1 ) {
132: DQR(BD(n1)[0],BD(n2)[0],uq,ur)
133: STON(ur,*nr);
134: } else {
135: ur = rem(n1,BD(n2)[0]); UTON(ur,*nr);
136: }
137: return;
138: } else if ( (d1 == d2) && ((tmp = cmpn(n1,n2)) <= 0) )
139: if ( !tmp ) {
140: *nr = 0;
141: } else {
142: COPY(n1,*nr);
143: }
144: else
145: divn(n1,n2,&q,nr);
146: }
147:
148: /* d = 2^(32*words)-lower */
149:
150: void remn_special(a,d,bits,lower,b)
151: N a,d;
152: int bits;
153: unsigned int lower;
154: N *b;
155: {
156: int words;
157: N r;
158: int rl,i;
159: unsigned int c,tmp;
160: unsigned int *rb,*src,*dst;
161:
162: if ( cmpn(a,d) < 0 ) {
163: *b = a;
164: return;
165: }
166: words = bits>>5;
167: r = NALLOC(PL(a)); dupn(a,r);
168: rl = PL(r); rb = BD(r);
169: while ( rl > words ) {
170: src = rb+words;
171: dst = rb;
172: i = rl-words;
173: for ( c = 0; i > 0; i--, src++, dst++ ) {
174: DMA2(*src,lower,c,*dst,c,*dst)
175: *src = 0;
176: }
177: for ( ; c; dst++ ) {
178: tmp = *dst + c;
179: c = tmp < c ? 1 : 0;
180: *dst = tmp;
181: }
182: rl = dst-rb;
183: }
184: for ( i = words-1; i >= 0 && !rb[i]; i-- );
185: PL(r) = i + 1;
186: if ( cmpn(r,d) >= 0 ) {
187: tmp = rb[0]-BD(d)[0];
188: UTON(tmp,*b);
189: } else
190: *b = r;
191: }
192: void mulin(n,d,p)
193: N n;
194: unsigned int d;
195: unsigned int *p;
196: {
197: unsigned int carry;
198: unsigned int *m,*me;
199:
200: bzero((char *)p,(int)((PL(n)+1)*sizeof(int)));
201: for ( carry = 0, m = BD(n), me = m+PL(n); m < me; p++, m++ ) {
202: DMA2(*m,d,*p,carry,carry,*p)
203: }
204: *p = carry;
205: }
206:
207: unsigned int divin(n,dvr,q)
208: N n;
209: unsigned int dvr;
210: N *q;
211: {
212: int d;
213: unsigned int up;
214: unsigned int *m,*mq,*md;
215: N nq;
216:
217: d = PL(n); m = BD(n);
218: if ( ( d == 1 ) && ( dvr > *m ) ) {
219: *q = 0;
220: return *m;
221: }
222: *q = nq = NALLOC(d); INITRC(nq);
223: for ( md = m+d-1, mq = BD(nq)+d-1, up = 0; md >= m; md--, mq-- ) {
224: DSAB(dvr,up,*md,*mq,up)
225: }
226: PL(nq) = (BD(nq)[d-1]?d:d-1);
227: if ( !PL(nq) )
228: FREEN(nq);
229: return ( up );
230: }
231:
232: void bprintn(n)
233: N n;
234: {
235: int l,i;
236: unsigned int *b;
237:
238: if ( !n )
239: printf("0");
240: else {
241: l = PL(n); b = BD(n);
242: for ( i = l-1; i >= 0; i-- )
243: printf("%010u|",b[i]);
244: }
245: }
246:
247: void bxprintn(n)
248: N n;
249: {
250: int l,i;
251: unsigned int *b;
252:
253: if ( !n )
254: printf("0");
255: else {
256: l = PL(n); b = BD(n);
257: for ( i = l-1; i >= 0; i-- )
258: printf("%08x|",b[i]);
259: }
260: }
261:
262: #if defined(VISUAL) || defined(i386)
263: void muln(n1,n2,nr)
264: N n1,n2,*nr;
265: {
266: unsigned int tmp,carry,mul;
267: unsigned int *p1,*m1,*m2;
268: int i,d1,d2,d;
269: N r;
270:
271: if ( !n1 || !n2 )
272: *nr = 0;
273: else if ( UNIN(n1) )
274: COPY(n2,*nr);
275: else if ( UNIN(n2) )
276: COPY(n1,*nr);
277: else if ( (PL(n1) == 1) && (PL(n2) == 1) ) {
278: DM(*BD(n1),*BD(n2),carry,tmp)
279: if ( carry ) {
280: *nr = r = NALLOC(2); INITRC(r);
281: PL(r) = 2; p1 = BD(r); *p1++ = tmp; *p1 = carry;
282: } else
283: STON(tmp,*nr);
284: } else {
285: d1 = PL(n1); d2 = PL(n2);
286: d = d1+d2;
287: *nr = r = NALLOC(d); INITRC(r);
288: bzero((char *)BD(r),(int)((d1+d2)*sizeof(int)));
289: for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
290: if ( mul = *m2 )
291: muln_1(m1,d1,mul,BD(r)+i);
292: PL(r) = (BD(r)[d-1]?d:d-1);
293: }
294: }
295:
296: void _muln(n1,n2,nr)
297: N n1,n2,nr;
298: {
299: unsigned int mul;
300: unsigned int *m1,*m2;
301: int i,d1,d2,d;
302:
303: if ( !n1 || !PL(n1) || !n2 || !PL(n2) )
304: PL(nr) = 0;
305: else if ( UNIN(n1) )
306: dupn(n2,nr);
307: else if ( UNIN(n2) )
308: dupn(n1,nr);
309: else {
310: d1 = PL(n1); d2 = PL(n2);
311: d = d1+d2;
312: bzero((char *)BD(nr),(int)((d1+d2)*sizeof(int)));
313: for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
314: if ( mul = *m2 )
315: muln_1(m1,d1,mul,BD(nr)+i);
316: PL(nr) = (BD(nr)[d-1]?d:d-1);
317: }
318: }
319:
320: void muln_1(p,s,d,r)
321: unsigned int *p;
322: int s;
323: unsigned int d;
324: unsigned int *r;
325: {
326: /* esi : p, edi : r, carry : ebx, s : ecx */
327: #if defined(VISUAL)
328: __asm {
329: push esi
330: push edi
331: mov esi,p
332: mov edi,r
333: mov ecx,s
334: xor ebx,ebx
335: Lstart_muln:
336: mov eax,DWORD PTR [esi]
337: mul d
338: add eax,DWORD PTR [edi]
339: adc edx,0
340: add eax,ebx
341: adc edx,0
342: mov DWORD PTR[edi],eax
343: mov ebx,edx
344: lea esi,DWORD PTR [esi+4]
345: lea edi,DWORD PTR [edi+4]
346: dec ecx
347: jnz Lstart_muln
348: mov DWORD PTR [edi],ebx
349: pop edi
350: pop esi
351: }
352: #else
353: asm volatile("\
354: movl %0,%%esi;\
355: movl %1,%%edi;\
356: movl $0,%%ebx;\
357: Lstart_muln:;\
358: movl (%%esi),%%eax;\
359: mull %2;\
360: addl (%%edi),%%eax;\
361: adcl $0,%%edx;\
362: addl %%ebx,%%eax;\
363: adcl $0,%%edx;\
364: movl %%eax,(%%edi);\
365: movl %%edx,%%ebx;\
366: leal 4(%%esi),%%esi;\
367: leal 4(%%edi),%%edi;\
368: decl %3;\
369: jnz Lstart_muln;\
370: movl %%ebx,(%%edi)"\
371: :\
372: :"m"(p),"m"(r),"m"(d),"m"(s)\
373: :"eax","ebx","edx","esi","edi");
374: #endif
375: }
376:
377: void divnmain(d1,d2,m1,m2,q)
378: int d1,d2;
379: unsigned int *m1,*m2,*q;
380: {
381: int i,j;
382: UL r,ltmp;
383: unsigned int l,ur;
384: unsigned int *n1,*n2;
385: unsigned int u,qhat;
386: unsigned int v1,v2;
387:
388: v1 = m2[d2-1]; v2 = m2[d2-2];
389: #if 0
390: if ( v1 == (1<<31) && !v2 ) {
391: divnmain_special(d1,d2,m1,m2,q);
392: return;
393: }
394: #endif
395: for ( j = d1-d2, m1 += j, q += j; j >= 0; j--, q--, m1-- ) {
396: n1 = m1+d2; n2 = m1+d2-1;
397: if ( *n1 == v1 ) {
398: qhat = ULBASE - 1;
399: r = (UL)*n1 + (UL)*n2;
400: } else {
401: DSAB(v1,*n1,*n2,qhat,ur)
402: r = (UL)ur;
403: }
404: DM(v2,qhat,u,l)
405: while ( 1 ) {
406: if ((r > (UL)u) || ((r == (UL)u) && (*(n1-2) >= l)))
407: break;
408: if ( l >= v2 )
409: l -= v2;
410: else {
411: l = (unsigned int)((UL)l+(ULBASE-(UL)v2)); u--;
412: }
413: r += v1; qhat--;
414: }
415: if ( qhat ) {
416: u = divn_1(m2,d2,qhat,m1);
417: if ( *(m1+d2) < u ) {
418: for ( i = d2, qhat--, ur = 0, n1 = m1, n2 = m2;
419: i > 0; i--, n1++, n2++ ) {
420: ltmp = (UL)*n1 + (UL)*n2 + (UL)ur;
421: *n1 = (unsigned int)(ltmp & BMASK);
422: ur = (unsigned int)(ltmp >> BSH);
423: }
424: }
425: *n1 = 0;
426: }
427: *q = qhat;
428: }
429: }
430:
431: void divnmain_special(d1,d2,m1,m2,q)
432: int d1,d2;
433: unsigned int *m1,*m2,*q;
434: {
435: int i,j;
436: UL ltmp;
437: unsigned int ur;
438: unsigned int *n1,*n2;
439: unsigned int u,qhat;
440: unsigned int v1,v2;
441:
442: v1 = m2[d2-1]; v2 = 0;
443: for ( j = d1-d2, m1 += j, q += j; j >= 0; j--, q--, m1-- ) {
444: n1 = m1+d2; n2 = m1+d2-1;
445: qhat = ((*n1)<<1)|((*n2)>>31);
446: if ( qhat ) {
447: u = divn_1(m2,d2,qhat,m1);
448: if ( *(m1+d2) < u ) {
449: for ( i = d2, qhat--, ur = 0, n1 = m1, n2 = m2;
450: i > 0; i--, n1++, n2++ ) {
451: ltmp = (UL)*n1 + (UL)*n2 + (UL)ur;
452: *n1 = (unsigned int)(ltmp & BMASK);
453: ur = (unsigned int)(ltmp >> BSH);
454: }
455: }
456: *n1 = 0;
457: }
458: *q = qhat;
459: }
460: }
461:
462: unsigned int divn_1(p,s,d,r)
463: unsigned int *p;
464: int s;
465: unsigned int d;
466: unsigned int *r;
467: {
468: /*
469: unsigned int borrow,l;
470:
471: for ( borrow = 0; s--; p++, r++ ) {
472: DMA(*p,d,borrow,borrow,l)
473: if ( *r >= l )
474: *r -= l;
475: else {
476: *r = (unsigned int)((UL)*r+(ULBASE-(UL)l)); borrow++;
477: }
478: }
479: return borrow;
480: */
481: /* esi : p, edi : r, borrow : ebx, s : ecx */
482: #if defined(VISUAL)
483: __asm {
484: push esi
485: push edi
486: mov esi,p
487: mov edi,r
488: mov ecx,s
489: xor ebx,ebx
490: Lstart_divn:
491: mov eax,DWORD PTR [esi]
492: mul d
493: add eax,ebx
494: adc edx,0
495: sub DWORD PTR [edi],eax
496: adc edx,0
497: mov ebx,edx
498: lea esi,DWORD PTR [esi+4]
499: lea edi,DWORD PTR [edi+4]
500: dec ecx
501: jnz Lstart_divn
502: mov eax,ebx
503: pop edi
504: pop esi
505: }
506: /* return value is in eax. */
507: #else
508: unsigned int borrow;
509:
510: asm volatile("\
511: movl %1,%%esi;\
512: movl %2,%%edi;\
513: movl $0,%%ebx;\
514: Lstart_divn:;\
515: movl (%%esi),%%eax;\
516: mull %3;\
517: addl %%ebx,%%eax;\
518: adcl $0,%%edx;\
519: subl %%eax,(%%edi);\
520: adcl $0,%%edx;\
521: movl %%edx,%%ebx;\
522: leal 4(%%esi),%%esi;\
523: leal 4(%%edi),%%edi;\
524: decl %4;\
525: jnz Lstart_divn;\
526: movl %%ebx,%0"\
527: :"=m"(borrow)\
528: :"m"(p),"m"(r),"m"(d),"m"(s)\
529: :"eax","ebx","edx","esi","edi");
530:
531: return borrow;
532: #endif
533: }
534:
535: #else
536:
537: void muln(n1,n2,nr)
538: N n1,n2,*nr;
539: {
540: unsigned int tmp,carry,mul;
541: unsigned int *p1,*pp,*m1,*m2;
542: int i,j,d1,d2;
543: N r;
544:
545: if ( !n1 || !n2 )
546: *nr = 0;
547: else if ( UNIN(n1) )
548: COPY(n2,*nr);
549: else if ( UNIN(n2) )
550: COPY(n1,*nr);
551: else if ( (PL(n1) == 1) && (PL(n2) == 1) ) {
552: DM(*BD(n1),*BD(n2),carry,tmp)
553: if ( carry ) {
554: *nr = r = NALLOC(2); INITRC(r);
555: PL(r) = 2; p1 = BD(r); *p1++ = tmp; *p1 = carry;
556: } else
557: STON(tmp,*nr);
558: } else {
559: d1 = PL(n1); d2 = PL(n2);
560: *nr = r = NALLOC(d1+d2); INITRC(r);
561: bzero((char *)BD(r),(int)((d1+d2)*sizeof(int)));
562: for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
563: if ( mul = *m2 ) {
564: for ( j = d1, carry = 0, p1 = m1, pp = BD(r)+i;
565: j--; pp++ ) {
566: DMA2(*p1++,mul,*pp,carry,carry,*pp)
567: }
568: *pp = carry;
569: }
570: PL(r) = (carry?d1+d2:d1+d2-1);
571: }
572: }
573:
574: void _muln(n1,n2,nr)
575: N n1,n2,nr;
576: {
577: unsigned int tmp,carry,mul;
578: unsigned int *p1,*pp,*m1,*m2;
579: int i,j,d1,d2;
580:
581: if ( !n1 || !PL(n1) || !n2 || !PL(n2) )
582: PL(nr) = 0;
583: else if ( UNIN(n1) )
584: dupn(n2,nr);
585: else if ( UNIN(n2) )
586: dupn(n1,nr);
587: else {
588: d1 = PL(n1); d2 = PL(n2);
589: bzero((char *)BD(nr),(int)((d1+d2)*sizeof(int)));
590: for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
591: if ( mul = *m2 ) {
592: for ( j = d1, carry = 0, p1 = m1, pp = BD(nr)+i;
593: j--; pp++ ) {
594: DMA2(*p1++,mul,*pp,carry,carry,*pp)
595: }
596: *pp = carry;
597: }
598: PL(nr) = (carry?d1+d2:d1+d2-1);
599: }
600: }
601:
602: /* r[0...s] = p[0...s-1]*d */
603:
604: void muln_1(p,s,d,r)
605: unsigned int *p;
606: int s;
607: unsigned int d;
608: unsigned int *r;
609: {
610: unsigned int carry;
611:
612: for ( carry = 0; s--; p++, r++ ) {
613: DMA2(*p,d,carry,*r,carry,*r)
614: }
615: *r = carry;
616: }
617:
618: void divnmain(d1,d2,m1,m2,q)
619: int d1,d2;
620: unsigned int *m1,*m2,*q;
621: {
622: int i,j;
623: UL r,ltmp;
624: unsigned int l,ur,tmp;
625: unsigned int *n1,*n2;
626: unsigned int u,qhat;
627: unsigned int v1,v2;
628:
629: v1 = m2[d2-1]; v2 = m2[d2-2];
630: for ( j = d1-d2, m1 += j, q += j; j >= 0; j--, q--, m1-- ) {
631: n1 = m1+d2; n2 = m1+d2-1;
632: if ( *n1 == v1 ) {
633: qhat = ULBASE - 1;
634: r = (UL)*n1 + (UL)*n2;
635: } else {
636: DSAB(v1,*n1,*n2,qhat,ur)
637: r = (UL)ur;
638: }
639: DM(v2,qhat,u,l)
640: while ( 1 ) {
641: if ((r > (UL)u) || ((r == (UL)u) && (*(n1-2) >= l)))
642: break;
643: if ( l >= v2 )
644: l -= v2;
645: else {
646: l = (unsigned int)((UL)l+(ULBASE-(UL)v2)); u--;
647: }
648: r += v1; qhat--;
649: }
650: if ( qhat ) {
651: u = divn_1(m2,d2,qhat,m1);
652: if ( *(m1+d2) < u ) {
653: for ( i = d2, qhat--, ur = 0, n1 = m1, n2 = m2;
654: i > 0; i--, n1++, n2++ ) {
655: ltmp = (UL)*n1 + (UL)*n2 + (UL)ur;
656: *n1 = (unsigned int)(ltmp & BMASK);
657: ur = (unsigned int)(ltmp >> BSH);
658: }
659: }
660: *n1 = 0;
661: }
662: *q = qhat;
663: }
664: }
665:
666: unsigned int divn_1(p,s,d,r)
667: unsigned int *p;
668: int s;
669: unsigned int d;
670: unsigned int *r;
671: {
672: unsigned int borrow,l;
673:
674: for ( borrow = 0; s--; p++, r++ ) {
675: DMA(*p,d,borrow,borrow,l)
676: if ( *r >= l )
677: *r -= l;
678: else {
679: *r = (unsigned int)((UL)*r+(ULBASE-(UL)l)); borrow++;
680: }
681: }
682: return borrow;
683: }
684: #endif