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