![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to divbig.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 42BSD / ucb / lisp / franz|
Return to divbig.c CVS log | Up to [CSRG BSD Unix] / 42BSD / ucb / lisp / franz |
1.1 root 1: #ifndef lint
2: static char *rcsid =
3: "$Header: divbig.c,v 1.3 83/09/12 14:17:07 sklower Exp $";
4: #endif
5:
6: /* -[Sat Jan 29 12:22:36 1983 by jkf]-
7: * divbig.c $Locker: $
8: * bignum division
9: *
10: * (c) copyright 1982, Regents of the University of California
11: */
12:
13:
14: #include "global.h"
15:
16: #define b 0x40000000
17: #define toint(p) ((int) (p))
18:
19: divbig(dividend, divisor, quotient, remainder)
20: lispval dividend, divisor, *quotient, *remainder;
21: {
22: register *ujp, *vip;
23: int *alloca(), d, negflag = 0, m, n, carry, rem, qhat, j;
24: int borrow, negrem = 0;
25: long *utop = sp(), *ubot, *vbot, *qbot;
26: register lispval work; lispval export();
27: Keepxs();
28:
29: /* copy dividend */
30: for(work = dividend; work; work = work ->s.CDR)
31: stack(work->s.I);
32: ubot = sp();
33: if(*ubot < 0) { /* knuth's division alg works only for pos
34: bignums */
35: negflag ^= 1;
36: negrem = 1;
37: dsmult(utop-1,ubot,-1);
38: }
39: stack(0);
40: ubot = sp();
41:
42:
43: /*copy divisor */
44: for(work = divisor; work; work = work->s.CDR)
45: stack(work->s.I);
46:
47: vbot = sp();
48: stack(0);
49: if(*vbot < 0) {
50: negflag ^= 1;
51: dsmult(ubot-1,vbot,-1);
52: }
53:
54: /* check validity of data */
55: n = ubot - vbot;
56: m = utop - ubot - n - 1;
57: if (n == 1) {
58: /* do destructive division by a single. */
59: rem = dsdiv(utop-1,ubot,*vbot);
60: if(negrem)
61: rem = -rem;
62: if(negflag)
63: dsmult(utop-1,ubot,-1);
64: if(remainder)
65: *remainder = inewint(rem);
66: if(quotient)
67: *quotient = export(utop,ubot);
68: Freexs();
69: return;
70: }
71: if (m < 0) {
72: if (remainder)
73: *remainder = dividend;
74: if(quotient)
75: *quotient = inewint(0);
76: Freexs();
77: return;
78: }
79: qbot = alloca(toint(utop) + toint(vbot) - 2 * toint(ubot));
80: d1:
81: d = b /(*vbot +1);
82: dsmult(utop-1,ubot,d);
83: dsmult(ubot-1,vbot,d);
84:
85: d2: for(j=0,ujp=ubot; j <= m; j++,ujp++) {
86:
87: d3:
88: qhat = calqhat(ujp,vbot);
89: d4:
90: if((borrow = mlsb(ujp + n, ujp, ubot, -qhat)) < 0) {
91: adback(ujp + n, ujp, ubot);
92: qhat--;
93: }
94: qbot[j] = qhat;
95: }
96: d8: if(remainder) {
97: dsdiv(utop-1, utop - n, d);
98: if(negrem) dsmult(utop-1,utop-n,-1);
99: *remainder = export(utop,utop-n);
100: }
101: if(quotient) {
102: if(negflag)
103: dsmult(qbot+m,qbot,-1);
104: *quotient = export(qbot + m + 1, qbot);
105: }
106: Freexs();
107: }
108: /*
109: * asm code commented out due to optimizer bug
110: * also, this file is now shared with the 68k version!
111: calqhat(ujp,v1p)
112: register int *ujp, *v1p;
113: {
114: asm(" cmpl (r10),(r11) # v[1] == u[j] ??");
115: asm(" beql 2f ");
116: asm(" # calculate qhat and rhat simultaneously,");
117: asm(" # qhat in r0");
118: asm(" # rhat in r1");
119: asm(" emul (r11),$0x40000000,4(r11),r4 # u[j]b+u[j+1] into r4,r5");
120: asm(" ediv (r10),r4,r0,r1 # qhat = ((u[j]b+u[j+1])/v[1]) into r0");
121: asm(" # (u[j]b+u[j+1] -qhat*v[1]) into r1");
122: asm(" # called rhat");
123: asm("1:");
124: asm(" # check if v[2]*qhat > rhat*b+u[j+2]");
125: asm(" emul r0,4(r10),$0,r2 # qhat*v[2] into r3,r2");
126: asm(" emul r1,$0x40000000,8(r11),r8 #rhat*b + u[j+2] into r9,r8");
127: asm(" # give up if r3,r2 <= r9,r8, otherwise iterate");
128: asm(" subl2 r8,r2 # perform r3,r2 - r9,r8");
129: asm(" sbwc r9,r3");
130: asm(" bleq 3f # give up if negative or equal");
131: asm(" decl r0 # otherwise, qhat = qhat - 1");
132: asm(" addl2 (r10),r1 # since dec'ed qhat, inc rhat by v[1]");
133: asm(" jbr 1b");
134: asm("2: ");
135: asm(" # get here if v[1]==u[j]");
136: asm(" # set qhat to b-1");
137: asm(" # rhat is easily calculated since if we substitute b-1 for qhat in");
138: asm(" # the formula, then it simplifies to (u[j+1] + v[1])");
139: asm(" # ");
140: asm(" addl3 4(r11),(r10),r1 # rhat = u[j+1] + v[1]");
141: asm(" movl $0x3fffffff,r0 # qhat = b-1");
142: asm(" jbr 1b");
143: asm("3:");
144: }
145: mlsb(utop,ubot,vtop,nqhat)
146: register int *utop, *ubot, *vtop;
147: register int nqhat;
148: {
149: asm(" clrl r0");
150: asm("loop2: addl2 (r11),r0");
151: asm(" emul r8,-(r9),r0,r2");
152: asm(" extzv $0,$30,r2,(r11)");
153: asm(" extv $30,$32,r2,r0");
154: asm(" acbl r10,$-4,r11,loop2");
155: }
156: adback(utop,ubot,vtop)
157: register int *utop, *ubot, *vtop;
158: {
159: asm(" clrl r0");
160: asm("loop3: addl2 -(r9),r0");
161: asm(" addl2 (r11),r0");
162: asm(" extzv $0,$30,r0,(r11)");
163: asm(" extv $30,$2,r0,r0");
164: asm(" acbl r10,$-4,r11,loop3");
165: }
166: dsdiv(top,bot,div)
167: register int* bot;
168: {
169: asm(" clrl r0");
170: asm("loop4: emul r0,$0x40000000,(r11),r1");
171: asm(" ediv 12(ap),r1,(r11),r0");
172: asm(" acbl 4(ap),$4,r11,loop4");
173: }
174: dsmult(top,bot,mult)
175: register int* top;
176: {
177: asm(" clrl r0");
178: asm("loop5: emul 12(ap),(r11),r0,r1");
179: asm(" extzv $0,$30,r1,(r11)");
180: asm(" extv $30,$32,r1,r0");
181: asm(" acbl 8(ap),$-4,r11,loop5");
182: asm(" movl r1,4(r11)");
183: }
184: */
185: lispval
186: export(top,bot)
187: register long *top, *bot;
188: {
189: register lispval p;
190: lispval result;
191:
192: top--; /* screwey convention matches original
193: vax assembler convenience */
194: while(bot < top)
195: {
196: if(*bot==0)
197: bot++;
198: else if(*bot==-1)
199: *++bot |= 0xc0000000;
200: else break;
201: }
202: if(bot==top) return(inewint(*bot));
203: result = p = newsdot();
204: protect(p);
205: p->s.I = *top--;
206: while(top >= bot) {
207: p = p->s.CDR = newdot();
208: p->s.I = *top--;
209: }
210: p->s.CDR = 0;
211: np--;
212: return(result);
213: }
214:
215: #define MAXINT 0x80000000L
216:
217: Ihau(fix)
218: register int fix;
219: {
220: register count;
221: if(fix==MAXINT)
222: return(32);
223: if(fix < 0)
224: fix = -fix;
225: for(count = 0; fix; count++)
226: fix /= 2;
227: return(count);
228: }
229: lispval
230: Lhau()
231: {
232: register count;
233: register lispval handy;
234: register dum1,dum2;
235: lispval Labsval();
236:
237: handy = lbot->val;
238: top:
239: switch(TYPE(handy)) {
240: case INT:
241: count = Ihau(handy->i);
242: break;
243: case SDOT:
244: handy = Labsval();
245: for(count = 0; handy->s.CDR!=((lispval) 0); handy = handy->s.CDR)
246: count += 30;
247: count += Ihau(handy->s.I);
248: break;
249: default:
250: handy = errorh1(Vermisc,"Haulong: bad argument",nil,
251: TRUE,997,handy);
252: goto top;
253: }
254: return(inewint(count));
255: }
256: lispval
257: Lhaipar()
258: {
259: int *sp();
260: register lispval work;
261: register n;
262: register int *top = sp() - 1;
263: register int *bot;
264: int mylen;
265:
266: /*chkarg(2);*/
267: work = lbot->val;
268: /* copy data onto stack */
269: on1:
270: switch(TYPE(work)) {
271: case INT:
272: stack(work->i);
273: break;
274: case SDOT:
275: for(; work!=((lispval) 0); work = work->s.CDR)
276: stack(work->s.I);
277: break;
278: default:
279: work = errorh1(Vermisc,"Haipart: bad first argument",nil,
280: TRUE,996,work);
281: goto on1;
282: }
283: bot = sp();
284: if(*bot < 0) {
285: stack(0);
286: dsmult(top,bot,-1);
287: bot--;
288: }
289: for(; *bot==0 && bot < top; bot++);
290: /* recalculate haulong internally */
291: mylen = (top - bot) * 30 + Ihau(*bot);
292: /* get second argument */
293: work = lbot[1].val;
294: while(TYPE(work)!=INT)
295: work = errorh1(Vermisc,"Haipart: 2nd arg not int",nil,
296: TRUE,995,work);
297: n = work->i;
298: if(n >= mylen || -n >= mylen)
299: goto done;
300: if(n==0) return(inewint(0));
301: if(n > 0) {
302: /* Here we want n most significant bits
303: so chop off mylen - n bits */
304: stack(0);
305: n = mylen - n;
306: for(n; n >= 30; n -= 30)
307: top--;
308: if(top < bot)
309: error("Internal error in haipart #1",FALSE);
310: dsdiv(top,bot,1<<n);
311:
312: } else {
313: /* here we want abs(n) low order bits */
314: stack(0);
315: bot = top + 1;
316: for(; n <= 0; n += 30)
317: bot--;
318: n = 30 - n;
319: *bot &= ~ (-1<<n);
320: }
321: done:
322: return(export(top + 1,bot));
323: }
324: #define STICKY 1
325: #define TOEVEN 2
326: lispval
327: Ibiglsh(bignum,count,mode)
328: lispval bignum, count;
329: {
330: int *sp();
331: register lispval work;
332: register n;
333: register int *top = sp() - 1;
334: register int *bot;
335: int mylen, guard = 0, sticky = 0, round = 0;
336: lispval export();
337:
338: /* get second argument */
339: work = count;
340: while(TYPE(work)!=INT)
341: work = errorh1(Vermisc,"Bignum-shift: 2nd arg not int",nil,
342: TRUE,995,work);
343: n = work->i;
344: if(n==0) return(bignum);
345: for(; n >= 30; n -= 30) {/* Here we want to multiply by 2^n
346: so start by copying n/30 zeroes
347: onto stack */
348: stack(0);
349: }
350:
351: work = bignum; /* copy data onto stack */
352: on1:
353: switch(TYPE(work)) {
354: case INT:
355: stack(work->i);
356: break;
357: case SDOT:
358: for(; work!=((lispval) 0); work = work->s.CDR)
359: stack(work->s.I);
360: break;
361: default:
362: work = errorh1(Vermisc,"Bignum-shift: bad bignum argument",nil,
363: TRUE,996,work);
364: goto on1;
365: }
366: bot = sp();
367: if(n >= 0) {
368: stack(0);
369: bot--;
370: dsmult(top,bot,1<<n);
371: } else {
372: /* Trimming will only work without leading
373: zeroes without my having to think
374: a lot harder about it, if the inputs
375: are canonical */
376: for(n = -n; n > 30; n -= 30) {
377: if(guard) sticky |= 1;
378: guard = round;
379: if(top > bot) {
380: round = *top;
381: top --;
382: } else {
383: round = *top;
384: *top >>= 30;
385: }
386: }
387: if(n > 0) {
388: if(guard) sticky |= 1;
389: guard = round;
390: round = dsrsh(top,bot,-n,-1<<n);
391: }
392: stack(0); /*so that dsadd1 will work;*/
393: if (mode==STICKY) {
394: if(((*top&1)==0) && (round | guard | sticky))
395: dsadd1(top,bot);
396: } else if (mode==TOEVEN) {
397: int mask;
398:
399: if(n==0) n = 30;
400: mask = (1<<(n-1));
401: if(! (round & mask) ) goto chop;
402: mask -= 1;
403: if( ((round&mask)==0)
404: && guard==0
405: && sticky==0
406: && (*top&1)==0 ) goto chop;
407: dsadd1(top,bot);
408: }
409: chop:;
410: }
411: work = export(top + 1,bot);
412: return(work);
413: }
414:
415: /*From drb Mon Jul 27 01:25:56 1981
416: To: sklower
417:
418: The idea is that the answer/2
419: is equal to the exact answer/2 rounded towards - infinity. The final bit
420: of the answer is the "or" of the true final bit, together with all true
421: bits after the binary point. In other words, the 1's bit of the answer
422: is almost always 1. THE FINAL BIT OF THE ANSWER IS 0 IFF n*2^i = THE
423: ANSWER RETURNED EXACTLY, WITH A 0 FINAL BIT.
424:
425:
426: To try again, more succintly: the answer is correct to within 1, and
427: the 1's bit of the answer will be 0 only if the answer is exactly
428: correct. */
429:
430: lispval
431: Lsbiglsh()
432: {
433: register struct argent *mylbot = lbot;
434: chkarg(2,"sticky-bignum-leftshift");
435: return(Ibiglsh(lbot->val,lbot[1].val,STICKY));
436: }
437: lispval
438: Lbiglsh()
439: {
440: register struct argent *mylbot = lbot;
441: chkarg(2,"bignum-leftshift");
442: return(Ibiglsh(lbot->val,lbot[1].val,TOEVEN));
443: }
444: lispval
445: HackHex() /* this is a one minute function so drb and kls can debug biglsh */
446: /* (HackHex i) returns a string which is the result of printing i in hex */
447: {
448: register struct argent *mylbot = lbot;
449: char buf[32];
450: sprintf(buf,"%lx",lbot->val->i);
451: return((lispval)inewstr(buf));
452: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.