![[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] / 43BSD / ucb / lisp / franz|
Return to divbig.c CVS log | Up to [CSRG BSD Unix] / 43BSD / ucb / lisp / franz |
1.1 root 1: #ifndef lint
2: static char *rcsid =
3: "$Header: divbig.c,v 1.4 83/11/26 12:10:16 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: register lispval work;
260: register n;
261: register int *top = sp() - 1;
262: register int *bot;
263: int mylen;
264:
265: /*chkarg(2);*/
266: work = lbot->val;
267: /* copy data onto stack */
268: on1:
269: switch(TYPE(work)) {
270: case INT:
271: stack(work->i);
272: break;
273: case SDOT:
274: for(; work!=((lispval) 0); work = work->s.CDR)
275: stack(work->s.I);
276: break;
277: default:
278: work = errorh1(Vermisc,"Haipart: bad first argument",nil,
279: TRUE,996,work);
280: goto on1;
281: }
282: bot = sp();
283: if(*bot < 0) {
284: stack(0);
285: dsmult(top,bot,-1);
286: bot--;
287: }
288: for(; *bot==0 && bot < top; bot++);
289: /* recalculate haulong internally */
290: mylen = (top - bot) * 30 + Ihau(*bot);
291: /* get second argument */
292: work = lbot[1].val;
293: while(TYPE(work)!=INT)
294: work = errorh1(Vermisc,"Haipart: 2nd arg not int",nil,
295: TRUE,995,work);
296: n = work->i;
297: if(n >= mylen || -n >= mylen)
298: goto done;
299: if(n==0) return(inewint(0));
300: if(n > 0) {
301: /* Here we want n most significant bits
302: so chop off mylen - n bits */
303: stack(0);
304: n = mylen - n;
305: for(n; n >= 30; n -= 30)
306: top--;
307: if(top < bot)
308: error("Internal error in haipart #1",FALSE);
309: dsdiv(top,bot,1<<n);
310:
311: } else {
312: /* here we want abs(n) low order bits */
313: stack(0);
314: bot = top + 1;
315: for(; n <= 0; n += 30)
316: bot--;
317: n = 30 - n;
318: *bot &= ~ (-1<<n);
319: }
320: done:
321: return(export(top + 1,bot));
322: }
323: #define STICKY 1
324: #define TOEVEN 2
325: lispval
326: Ibiglsh(bignum,count,mode)
327: lispval bignum, count;
328: {
329: register lispval work;
330: register n;
331: register int *top = sp() - 1;
332: register int *bot;
333: int mylen, guard = 0, sticky = 0, round = 0;
334: lispval export();
335:
336: /* get second argument */
337: work = count;
338: while(TYPE(work)!=INT)
339: work = errorh1(Vermisc,"Bignum-shift: 2nd arg not int",nil,
340: TRUE,995,work);
341: n = work->i;
342: if(n==0) return(bignum);
343: for(; n >= 30; n -= 30) {/* Here we want to multiply by 2^n
344: so start by copying n/30 zeroes
345: onto stack */
346: stack(0);
347: }
348:
349: work = bignum; /* copy data onto stack */
350: on1:
351: switch(TYPE(work)) {
352: case INT:
353: stack(work->i);
354: break;
355: case SDOT:
356: for(; work!=((lispval) 0); work = work->s.CDR)
357: stack(work->s.I);
358: break;
359: default:
360: work = errorh1(Vermisc,"Bignum-shift: bad bignum argument",nil,
361: TRUE,996,work);
362: goto on1;
363: }
364: bot = sp();
365: if(n >= 0) {
366: stack(0);
367: bot--;
368: dsmult(top,bot,1<<n);
369: } else {
370: /* Trimming will only work without leading
371: zeroes without my having to think
372: a lot harder about it, if the inputs
373: are canonical */
374: for(n = -n; n > 30; n -= 30) {
375: if(guard) sticky |= 1;
376: guard = round;
377: if(top > bot) {
378: round = *top;
379: top --;
380: } else {
381: round = *top;
382: *top >>= 30;
383: }
384: }
385: if(n > 0) {
386: if(guard) sticky |= 1;
387: guard = round;
388: round = dsrsh(top,bot,-n,-1<<n);
389: }
390: stack(0); /*so that dsadd1 will work;*/
391: if (mode==STICKY) {
392: if(((*top&1)==0) && (round | guard | sticky))
393: dsadd1(top,bot);
394: } else if (mode==TOEVEN) {
395: int mask;
396:
397: if(n==0) n = 30;
398: mask = (1<<(n-1));
399: if(! (round & mask) ) goto chop;
400: mask -= 1;
401: if( ((round&mask)==0)
402: && guard==0
403: && sticky==0
404: && (*top&1)==0 ) goto chop;
405: dsadd1(top,bot);
406: }
407: chop:;
408: }
409: work = export(top + 1,bot);
410: return(work);
411: }
412:
413: /*From drb Mon Jul 27 01:25:56 1981
414: To: sklower
415:
416: The idea is that the answer/2
417: is equal to the exact answer/2 rounded towards - infinity. The final bit
418: of the answer is the "or" of the true final bit, together with all true
419: bits after the binary point. In other words, the 1's bit of the answer
420: is almost always 1. THE FINAL BIT OF THE ANSWER IS 0 IFF n*2^i = THE
421: ANSWER RETURNED EXACTLY, WITH A 0 FINAL BIT.
422:
423:
424: To try again, more succintly: the answer is correct to within 1, and
425: the 1's bit of the answer will be 0 only if the answer is exactly
426: correct. */
427:
428: lispval
429: Lsbiglsh()
430: {
431: register struct argent *mylbot = lbot;
432: chkarg(2,"sticky-bignum-leftshift");
433: return(Ibiglsh(lbot->val,lbot[1].val,STICKY));
434: }
435: lispval
436: Lbiglsh()
437: {
438: register struct argent *mylbot = lbot;
439: chkarg(2,"bignum-leftshift");
440: return(Ibiglsh(lbot->val,lbot[1].val,TOEVEN));
441: }
442: lispval
443: HackHex() /* this is a one minute function so drb and kls can debug biglsh */
444: /* (HackHex i) returns a string which is the result of printing i in hex */
445: {
446: register struct argent *mylbot = lbot;
447: char buf[32];
448: sprintf(buf,"%lx",lbot->val->i);
449: return((lispval)inewstr(buf));
450: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.