Commit | Line | Data |
---|---|---|
bc51592a C |
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 | } |