projects / unix-history / blob
? search:
summary | tags | clone url | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Start development on 386BSD 0.0
[unix-history] / .ref-BSD-4_3_Net_2 / usr / src / usr.bin / lisp / franz / divbig.c
#ifndef lint
static char *rcsid =
"$Header: divbig.c,v 1.4 83/11/26 12:10:16 sklower Exp $";
#endif
/* -[Sat Jan 29 12:22:36 1983 by jkf]-
* divbig.c $Locker: $
* bignum division
*
* (c) copyright 1982, Regents of the University of California
*/
#include "global.h"
#define b 0x40000000
#define toint(p) ((int) (p))
divbig(dividend, divisor, quotient, remainder)
lispval dividend, divisor, *quotient, *remainder;
{
register *ujp, *vip;
int *alloca(), d, negflag = 0, m, n, carry, rem, qhat, j;
int borrow, negrem = 0;
long *utop = sp(), *ubot, *vbot, *qbot;
register lispval work; lispval export();
Keepxs();
/* copy dividend */
for(work = dividend; work; work = work ->s.CDR)
stack(work->s.I);
ubot = sp();
if(*ubot < 0) { /* knuth's division alg works only for pos
bignums */
negflag ^= 1;
negrem = 1;
dsmult(utop-1,ubot,-1);
}
stack(0);
ubot = sp();
/*copy divisor */
for(work = divisor; work; work = work->s.CDR)
stack(work->s.I);
vbot = sp();
stack(0);
if(*vbot < 0) {
negflag ^= 1;
dsmult(ubot-1,vbot,-1);
}
/* check validity of data */
n = ubot - vbot;
m = utop - ubot - n - 1;
if (n == 1) {
/* do destructive division by a single. */
rem = dsdiv(utop-1,ubot,*vbot);
if(negrem)
rem = -rem;
if(negflag)
dsmult(utop-1,ubot,-1);
if(remainder)
*remainder = inewint(rem);
if(quotient)
*quotient = export(utop,ubot);
Freexs();
return;
}
if (m < 0) {
if (remainder)
*remainder = dividend;
if(quotient)
*quotient = inewint(0);
Freexs();
return;
}
qbot = alloca(toint(utop) + toint(vbot) - 2 * toint(ubot));
d1:
d = b /(*vbot +1);
dsmult(utop-1,ubot,d);
dsmult(ubot-1,vbot,d);
d2: for(j=0,ujp=ubot; j <= m; j++,ujp++) {
d3:
qhat = calqhat(ujp,vbot);
d4:
if((borrow = mlsb(ujp + n, ujp, ubot, -qhat)) < 0) {
adback(ujp + n, ujp, ubot);
qhat--;
}
qbot[j] = qhat;
}
d8: if(remainder) {
dsdiv(utop-1, utop - n, d);
if(negrem) dsmult(utop-1,utop-n,-1);
*remainder = export(utop,utop-n);
}
if(quotient) {
if(negflag)
dsmult(qbot+m,qbot,-1);
*quotient = export(qbot + m + 1, qbot);
}
Freexs();
}
/*
* asm code commented out due to optimizer bug
* also, this file is now shared with the 68k version!
calqhat(ujp,v1p)
register int *ujp, *v1p;
{
asm(" cmpl (r10),(r11) # v[1] == u[j] ??");
asm(" beql 2f ");
asm(" # calculate qhat and rhat simultaneously,");
asm(" # qhat in r0");
asm(" # rhat in r1");
asm(" emul (r11),$0x40000000,4(r11),r4 # u[j]b+u[j+1] into r4,r5");
asm(" ediv (r10),r4,r0,r1 # qhat = ((u[j]b+u[j+1])/v[1]) into r0");
asm(" # (u[j]b+u[j+1] -qhat*v[1]) into r1");
asm(" # called rhat");
asm("1:");
asm(" # check if v[2]*qhat > rhat*b+u[j+2]");
asm(" emul r0,4(r10),$0,r2 # qhat*v[2] into r3,r2");
asm(" emul r1,$0x40000000,8(r11),r8 #rhat*b + u[j+2] into r9,r8");
asm(" # give up if r3,r2 <= r9,r8, otherwise iterate");
asm(" subl2 r8,r2 # perform r3,r2 - r9,r8");
asm(" sbwc r9,r3");
asm(" bleq 3f # give up if negative or equal");
asm(" decl r0 # otherwise, qhat = qhat - 1");
asm(" addl2 (r10),r1 # since dec'ed qhat, inc rhat by v[1]");
asm(" jbr 1b");
asm("2: ");
asm(" # get here if v[1]==u[j]");
asm(" # set qhat to b-1");
asm(" # rhat is easily calculated since if we substitute b-1 for qhat in");
asm(" # the formula, then it simplifies to (u[j+1] + v[1])");
asm(" # ");
asm(" addl3 4(r11),(r10),r1 # rhat = u[j+1] + v[1]");
asm(" movl $0x3fffffff,r0 # qhat = b-1");
asm(" jbr 1b");
asm("3:");
}
mlsb(utop,ubot,vtop,nqhat)
register int *utop, *ubot, *vtop;
register int nqhat;
{
asm(" clrl r0");
asm("loop2: addl2 (r11),r0");
asm(" emul r8,-(r9),r0,r2");
asm(" extzv $0,$30,r2,(r11)");
asm(" extv $30,$32,r2,r0");
asm(" acbl r10,$-4,r11,loop2");
}
adback(utop,ubot,vtop)
register int *utop, *ubot, *vtop;
{
asm(" clrl r0");
asm("loop3: addl2 -(r9),r0");
asm(" addl2 (r11),r0");
asm(" extzv $0,$30,r0,(r11)");
asm(" extv $30,$2,r0,r0");
asm(" acbl r10,$-4,r11,loop3");
}
dsdiv(top,bot,div)
register int* bot;
{
asm(" clrl r0");
asm("loop4: emul r0,$0x40000000,(r11),r1");
asm(" ediv 12(ap),r1,(r11),r0");
asm(" acbl 4(ap),$4,r11,loop4");
}
dsmult(top,bot,mult)
register int* top;
{
asm(" clrl r0");
asm("loop5: emul 12(ap),(r11),r0,r1");
asm(" extzv $0,$30,r1,(r11)");
asm(" extv $30,$32,r1,r0");
asm(" acbl 8(ap),$-4,r11,loop5");
asm(" movl r1,4(r11)");
}
*/
lispval
export(top,bot)
register long *top, *bot;
{
register lispval p;
lispval result;
top--; /* screwey convention matches original
vax assembler convenience */
while(bot < top)
{
if(*bot==0)
bot++;
else if(*bot==-1)
*++bot |= 0xc0000000;
else break;
}
if(bot==top) return(inewint(*bot));
result = p = newsdot();
protect(p);
p->s.I = *top--;
while(top >= bot) {
p = p->s.CDR = newdot();
p->s.I = *top--;
}
p->s.CDR = 0;
np--;
return(result);
}
#define MAXINT 0x80000000L
Ihau(fix)
register int fix;
{
register count;
if(fix==MAXINT)
return(32);
if(fix < 0)
fix = -fix;
for(count = 0; fix; count++)
fix /= 2;
return(count);
}
lispval
Lhau()
{
register count;
register lispval handy;
register dum1,dum2;
lispval Labsval();
handy = lbot->val;
top:
switch(TYPE(handy)) {
case INT:
count = Ihau(handy->i);
break;
case SDOT:
handy = Labsval();
for(count = 0; handy->s.CDR!=((lispval) 0); handy = handy->s.CDR)
count += 30;
count += Ihau(handy->s.I);
break;
default:
handy = errorh1(Vermisc,"Haulong: bad argument",nil,
TRUE,997,handy);
goto top;
}
return(inewint(count));
}
lispval
Lhaipar()
{
register lispval work;
register n;
register int *top = sp() - 1;
register int *bot;
int mylen;
/*chkarg(2);*/
work = lbot->val;
/* copy data onto stack */
on1:
switch(TYPE(work)) {
case INT:
stack(work->i);
break;
case SDOT:
for(; work!=((lispval) 0); work = work->s.CDR)
stack(work->s.I);
break;
default:
work = errorh1(Vermisc,"Haipart: bad first argument",nil,
TRUE,996,work);
goto on1;
}
bot = sp();
if(*bot < 0) {
stack(0);
dsmult(top,bot,-1);
bot--;
}
for(; *bot==0 && bot < top; bot++);
/* recalculate haulong internally */
mylen = (top - bot) * 30 + Ihau(*bot);
/* get second argument */
work = lbot[1].val;
while(TYPE(work)!=INT)
work = errorh1(Vermisc,"Haipart: 2nd arg not int",nil,
TRUE,995,work);
n = work->i;
if(n >= mylen || -n >= mylen)
goto done;
if(n==0) return(inewint(0));
if(n > 0) {
/* Here we want n most significant bits
so chop off mylen - n bits */
stack(0);
n = mylen - n;
for(n; n >= 30; n -= 30)
top--;
if(top < bot)
error("Internal error in haipart #1",FALSE);
dsdiv(top,bot,1<<n);
} else {
/* here we want abs(n) low order bits */
stack(0);
bot = top + 1;
for(; n <= 0; n += 30)
bot--;
n = 30 - n;
*bot &= ~ (-1<<n);
}
done:
return(export(top + 1,bot));
}
#define STICKY 1
#define TOEVEN 2
lispval
Ibiglsh(bignum,count,mode)
lispval bignum, count;
{
register lispval work;
register n;
register int *top = sp() - 1;
register int *bot;
int mylen, guard = 0, sticky = 0, round = 0;
lispval export();
/* get second argument */
work = count;
while(TYPE(work)!=INT)
work = errorh1(Vermisc,"Bignum-shift: 2nd arg not int",nil,
TRUE,995,work);
n = work->i;
if(n==0) return(bignum);
for(; n >= 30; n -= 30) {/* Here we want to multiply by 2^n
so start by copying n/30 zeroes
onto stack */
stack(0);
}
work = bignum; /* copy data onto stack */
on1:
switch(TYPE(work)) {
case INT:
stack(work->i);
break;
case SDOT:
for(; work!=((lispval) 0); work = work->s.CDR)
stack(work->s.I);
break;
default:
work = errorh1(Vermisc,"Bignum-shift: bad bignum argument",nil,
TRUE,996,work);
goto on1;
}
bot = sp();
if(n >= 0) {
stack(0);
bot--;
dsmult(top,bot,1<<n);
} else {
/* Trimming will only work without leading
zeroes without my having to think
a lot harder about it, if the inputs
are canonical */
for(n = -n; n > 30; n -= 30) {
if(guard) sticky |= 1;
guard = round;
if(top > bot) {
round = *top;
top --;
} else {
round = *top;
*top >>= 30;
}
}
if(n > 0) {
if(guard) sticky |= 1;
guard = round;
round = dsrsh(top,bot,-n,-1<<n);
}
stack(0); /*so that dsadd1 will work;*/
if (mode==STICKY) {
if(((*top&1)==0) && (round | guard | sticky))
dsadd1(top,bot);
} else if (mode==TOEVEN) {
int mask;
if(n==0) n = 30;
mask = (1<<(n-1));
if(! (round & mask) ) goto chop;
mask -= 1;
if( ((round&mask)==0)
&& guard==0
&& sticky==0
&& (*top&1)==0 ) goto chop;
dsadd1(top,bot);
}
chop:;
}
work = export(top + 1,bot);
return(work);
}
/*From drb Mon Jul 27 01:25:56 1981
To: sklower
The idea is that the answer/2
is equal to the exact answer/2 rounded towards - infinity. The final bit
of the answer is the "or" of the true final bit, together with all true
bits after the binary point. In other words, the 1's bit of the answer
is almost always 1. THE FINAL BIT OF THE ANSWER IS 0 IFF n*2^i = THE
ANSWER RETURNED EXACTLY, WITH A 0 FINAL BIT.
To try again, more succintly: the answer is correct to within 1, and
the 1's bit of the answer will be 0 only if the answer is exactly
correct. */
lispval
Lsbiglsh()
{
register struct argent *mylbot = lbot;
chkarg(2,"sticky-bignum-leftshift");
return(Ibiglsh(lbot->val,lbot[1].val,STICKY));
}
lispval
Lbiglsh()
{
register struct argent *mylbot = lbot;
chkarg(2,"bignum-leftshift");
return(Ibiglsh(lbot->val,lbot[1].val,TOEVEN));
}
lispval
HackHex() /* this is a one minute function so drb and kls can debug biglsh */
/* (HackHex i) returns a string which is the result of printing i in hex */
{
register struct argent *mylbot = lbot;
char buf[32];
sprintf(buf,"%lx",lbot->val->i);
return((lispval)inewstr(buf));
}
Continuous source code history from original PDP-7 UNIX to FreeBSD 2, the first fully unencumbered FreeBSD release.