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