4.3BSD - /usr/src/ucb/lisp/franz/divbig.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: }

Defined functions

Defined variables

Defined macros