4.3BSD - /usr/src/ucb/lisp/franz/vax/bigmath.c

   1: /* $Header: bigmath.c 1.4 83/06/09 00:50:06 sklower Exp $ */
   2: 
   3: #include "config.h"
   4: 
   5:     .globl  _dmlad
   6: /*
   7: 	routine for destructive multiplication and addition to a bignum by
   8: 	two fixnums.
   9: 
  10: 	from C, the invocation is dmlad(sdot,mul,add);
  11: 	where sdot is the address of the first special cell of the bignum
  12: 	mul is the multiplier, add is the fixnum to be added (The latter
  13: 	being passed by value, as is the usual case.
  14: 
  15: 
  16: 	Register assignments:
  17: 
  18: 	r11 = current sdot
  19: 	r10 = carry
  20: 	r9  = previous sdot, for relinking.
  21: */
  22: _dmlad: .word   0x0e00
  23:     movl    4(ap),r11       #initialize cell pointer
  24:     movl    12(ap),r10      #initialize carry
  25: loop:   emul    8(ap),(r11),r10,r0  #r0 gets cell->car times mul + carry
  26: /*	ediv	$0x40000000,r0,r10,(r11)#cell->car gets prod % 2**30
  27: 					#carry gets quotient
  28: */
  29:     extzv   $0,$30,r0,(r11)
  30:     extv    $30,$32,r0,r10
  31:     movl    r11,r9          #save last cell for fixup at end.
  32:     movl    4(r11),r11      #move to next cell
  33:     bneq    loop            #done indicated by 0 for next sdot
  34:     tstl    r10         #if carry zero no need to allocate
  35:     beql    done            #new bigit
  36:     mcoml   r10,r3          #test to see if neg 1.
  37:     bneq    alloc           #if not must allocate new cell.
  38:     tstl    (r9)            #make sure product isn't -2**30
  39:     beql    alloc
  40:     movl    r0,(r9)         #save old lower half of product.
  41:     brb done
  42: alloc:  jsb _qnewdot            #otherwise allocate new bigit
  43:     movl    r10,(r0)        #store carry
  44:     movl    r0,4(r9)        #save new link cell
  45: done:   movl    4(ap),r0
  46:     ret
  47:     .globl _dodiv
  48: /*
  49: 	routine to destructively divide array representation of a bignum by
  50: 	1000000000
  51: 
  52: 	invocation:
  53: 		remainder = dodiv(top,bottom)
  54: 		int *top, *bottom;
  55: 	where *bottom is the address of the biggning of the array, *top is
  56: 	the top of the array.
  57: 
  58: 	register assignments:
  59: 	r0 = carry
  60: 	r1 & r2 = 64bit temporary
  61: 	r3 = pointer
  62: */
  63: _dodiv: .word   0
  64:     clrl        r0      #no carry to begin.
  65:     movl        8(ap),r3    #get pointer to array.
  66: loop2:  emul        $0x40000000,r0,(r3),r1
  67:     ediv        $1000000000,r1,(r3),r0
  68:     acbl        4(ap),$4,r3,loop2
  69:     ret
  70:     .globl  _dsneg
  71: /*
  72: 	dsneg(top, bot);
  73: 	int *top, *bot;
  74: 
  75: 	routine to destructively negate a bignum stored in array format
  76: 	lower order stuff at higher addresses. It is assume that the
  77: 	result will be positive.
  78: */
  79: _dsneg: .word   0
  80:     movl    4(ap),r1    #load up address.
  81:     clrl    r2      #set carry
  82: loop3:  mnegl   (r1),r0     #negate and take carry into account.
  83:     addl2   r2,r0
  84:     extzv   $0,$30,r0,(r1)
  85:     extv    $30,$2,r0,r2
  86:     acbl    8(ap),$-4,r1,loop3
  87:                 #decrease r1, and branch back if appropriate.
  88:     ret
  89: 
  90: /*
  91: 	bignum add routine
  92: 	basic data representation is each bigit is a positive number
  93: 	less than 2^30, except for the leading bigit, which is in
  94: 	the range -2^30 < x < 2^30.
  95: */
  96: 
  97:     .globl  _adbig
  98:     .globl  Bexport
  99:     .globl  backfr
 100: /*
 101: 	Initialization section
 102: */
 103: _adbig: .word   0x0fc0      #save registers 6-11
 104:     movl    4(ap),r1    #arg1 = addr of 1st bignum
 105:     movl    8(ap),r2    #arg2 = addr of 2nd bignum
 106:     clrl    r5      #r5   = carry
 107:     movl    $0xC0000000,r4  #r4   = clear constant.
 108:     movl    sp,r10      #save start address of bignum on stack.
 109:                 #note well that this is 4 above the actual
 110:                 #low order word.
 111: /*
 112: 	first loop is to waltz through both bignums adding
 113: 	bigits, pushing them onto stack.
 114: */
 115: loop4:  addl3   (r1),(r2),r0    #add bigits
 116:     addl2   r5,r0       #add carry
 117:     bicl3   r4,r0,-(sp) #save sum, no overflow possible
 118:     extv    $30,$2,r0,r5    #sign extend two high order bits
 119:                 #to be next carry.
 120:     movl    4(r1),r1    #get cdr
 121:     bleq    out1        #negative indicates end of list.
 122:     movl    4(r2),r2    #get cdr of second bignum
 123:     bgtr    loop4       #if neither list at end, do it again
 124: /*
 125: 	second loop propagates carries through higher order words.
 126: 	It assumes remaining list in r1.
 127: */
 128: loop5:  addl3   (r1),r5,r0  #add bigits and carry
 129:     bicl3   r4,r0,-(sp) #save sum, no overflow possible
 130:     extv    $30,$2,r0,r5    #sign extend two high order bits
 131:                 #to be next carry.
 132:     movl    4(r1),r1    #get cdr
 133: out2:   bgtr    loop5       #negative indicates end of list.
 134: out2a:  pushl   r5
 135: /*
 136: 	suppress unnecessary leading zeroes and -1's
 137: 
 138: 	WARNING:  this code is duplicated in C in divbig.c
 139: */
 140: iexport:movl    sp,r11      #more set up for output routine
 141: ckloop:
 142: Bexport:tstl    (r11)       #look at leading bigit
 143:     bgtr    copyit      #if positive, can allocate storage etc.
 144:     blss    negchk      #if neg, still a chance we can get by
 145:     cmpl    r11,r10     #check to see that
 146:     bgeq    copyit      #we don't pop everything off of stack
 147:     tstl    (r11)+      #incr r11
 148:     brb ckloop      #examine next
 149: negchk:
 150:     mcoml   (r11),r3        #r3 is junk register
 151:     bneq    copyit      #short test for -1
 152:     tstl    4(r11)      #examine next bigit
 153:     beql    copyit      #if zero must must leave as is.
 154:     cmpl    r11,r10     #check to see that
 155:     bgeq    copyit      #we don't pop everything off of stack
 156:     tstl    (r11)+      #incr r11
 157:     bisl2   r4,(r11)    #set high order two bits
 158:     brb negchk      #try to supress more leading -1's
 159: /*
 160: 	The following code is an error exit from the first loop
 161:  	and is out of place to avoid a jump around a jump.
 162: */
 163: out1:   movl    4(r2),r1    #get next addr of list to continue.
 164:     brb out2        #if second list simult. exhausted, do
 165:                 #right thing.
 166: /*
 167: 	loop6 is a faster version of loop5 when carries are no
 168: 	longer necessary.
 169: */
 170: loop6a: pushl   (r1)        #get datum
 171: loop6:  movl    4(r1),r1    #get cdr
 172:     bgtr    loop6a      #if not at end get next cell
 173:     brb out2a
 174: 
 175: /*
 176: 	create linked list representation of bignum
 177: */
 178: copyit: subl3   r11,r10,r2  #see if we can get away with allocating an int
 179:     bneq    on1     #test for having popped everything
 180:     subl3   $4,r10,r11  #if so, fix up pointer to bottom
 181:     brb intout      #and allocate int.
 182: on1:    cmpl    r2,$4       #if = 4, then can do
 183:     beql    intout
 184:     calls   $0,_newsdot #get new cell for new bignum
 185: backfr:
 186: #ifdef PORTABLE
 187:     movl    r0,*_np
 188:     addl2   $4,_np
 189: #else
 190:     movl    r0,(r6)+    #push address of cell on
 191:                 #arg stack to save from garbage collection.
 192:                 #There is guaranteed to be slop for a least 1
 193:                 #push without checking.
 194: #endif
 195: loop7:  movl    -(r10),(r0) #save bigit
 196:     movl    r0,r9       #r9 = old cell, to link
 197:     cmpl    r10,r11     #have we copy'ed all the bigits?
 198:     bleq    Edone
 199:     jsb _qnewdot    #get new cell for new bignum
 200:     movl    r0,4(r9)    #link new cell to old
 201:     brb loop7
 202: Edone:
 203:     clrl    4(r9)       #indicate end of list with 0
 204: #ifdef PORTABLE
 205:     subl2   $4,_np
 206:     movl    *_np,r0
 207: #else
 208:     movl    -(r6),r0    #give resultant address.
 209: #endif
 210:     ret
 211: /*
 212: 	export integer
 213: */
 214: intout: pushl   (r11)
 215:     calls   $1,_inewint
 216:     ret
 217:     .globl  _mulbig
 218: /*
 219: 	bignum multiplication routine
 220: 
 221: 	Initialization section
 222: */
 223: _mulbig:.word   0x0fc0      #save regs 6-11
 224:     movl    4(ap),r1    #get address of first bignum
 225:     movl    sp,r11      #save top of 1st bignum
 226: mloop1: pushl   (r1)        #get bigit
 227:     movl    4(r1),r1    #get cdr
 228:     bgtr    mloop1      #repeat if not done
 229:     movl    sp,r10      #save bottom of 1st bignum, top of 2nd bignum
 230: 
 231:     movl    8(ap),r1    #get address of 2nd bignum
 232: mloop2: pushl   (r1)        #get bigit
 233:     movl    4(r1),r1    #get cdr
 234:     bgtr    mloop2      #repeat if not done
 235:     movl    sp,r9       #save bottom of 2nd bignum
 236:     subl3   r9,r11,r6   #r6 contains sum of lengths of bignums
 237:     subl2   r6,sp
 238:     movl    sp,r8       #save bottom of product bignum
 239: /*
 240: 	Actual multiplication
 241: */
 242: m1: movc5   $0,(r8),$0,r6,(r8)#zap out stack space
 243:     movl    r9,r7       #r7 = &w[j +n] (+4 for a.d.) through calculation
 244:     subl3   $4,r10,r4   #r4 = &v[j]
 245: 
 246: m3: movl    r7,r5       #r7 = &w[j+n]
 247:     subl3   $4,r11,r3   #r3 = &u[i]
 248:     clrl    r2      #clear carry.
 249: 
 250: m4: addl2   -(r5),r2    #add w[i + j] to carry (no ofl poss)
 251:     emul    (r3),(r4),r2,r0 #r0 = u[i] * v[j] + sext(carry)
 252:     extzv   $0,$30,r0,(r5)  #get new bigit
 253:     extv    $30,$32,r0,r2   #get new carry
 254: 
 255: m5: acbl    r10,$-4,r3,m4   #r3 =- 4; if(r3 >= r10) goto m4; r10 = &[u1];
 256:     movl    r2,-(r5)    #save w[j] = carry
 257: 
 258: m6: subl2   $4,r7       #add just &w[j+n] (+4 for autodec)
 259:     acbl    r9,$-4,r4,m3    #r4 =- 4; if(r4>=r9) goto m5; r9 = &v[1]
 260: 
 261:     movl    r9,r10      #set up for output routine
 262:     movl    $0xC0000000,r4  #r4   = clear constant.
 263:     movq    20(fp),r6   #restor _np and _lbot !
 264:     brw iexport     #do it!
 265: /*
 266:  The remainder of this file are routines used in bignum division.
 267:  Interested parties should consult Knuth, Vol 2, and divbig.c.
 268:  These files are here only due to an optimizer bug.
 269: */
 270:     .align  1
 271:     .globl  _calqhat
 272: _calqhat:
 273:     .word   0xf00
 274:     movl    4(ap),r11       # &u[j] into r11
 275:     movl    8(ap),r10       # &v[1] into r10
 276:     cmpl    (r10),(r11)     # v[1] == u[j] ??
 277:     beql    L102
 278:     # calculate qhat and rhat simultaneously,
 279:     #  qhat in r0
 280:     #  rhat in r1
 281:     emul    (r11),$0x40000000,4(r11),r4 # u[j]b+u[j+1] into r4,r5
 282:     ediv    (r10),r4,r0,r1      # qhat = ((u[j]b+u[j+1])/v[1]) into r0
 283:                     # (u[j]b+u[j+1] -qhat*v[1]) into r1
 284:                     # called rhat
 285: L101:
 286:     # check if v[2]*qhat > rhat*b+u[j+2]
 287:     emul    r0,4(r10),$0,r2     # qhat*v[2] into r3,r2
 288:     emul    r1,$0x40000000,8(r11),r8 #rhat*b + u[j+2] into r9,r8
 289:     # give up if r3,r2 <= r9,r8, otherwise iterate
 290:     subl2   r8,r2           # perform r3,r2 - r9,r8
 291:     sbwc    r9,r3
 292:     bleq    L103            # give up if negative or equal
 293:     decl    r0          # otherwise, qhat = qhat - 1
 294:     addl2   (r10),r1        # since dec'ed qhat, inc rhat by v[1]
 295:     jbr L101
 296: L102:
 297:     # get here if v[1]==u[j]
 298:     # set qhat to b-1
 299:     # rhat is easily calculated since if we substitute b-1 for qhat in
 300:     # the formula, then it simplifies to (u[j+1] + v[1])
 301:     #
 302:     addl3   4(r11),(r10),r1     # rhat = u[j+1] + v[1]
 303:     movl    $0x3fffffff,r0      # qhat = b-1
 304:     jbr L101
 305: 
 306: L103:
 307:     ret
 308: 
 309:     .align  1
 310:     .globl  _mlsb
 311: _mlsb:
 312:     .word   .R2
 313:     movl    4(ap),r11
 314:     movl    8(ap),r10
 315:     movl    12(ap),r9
 316:     movl    16(ap),r8
 317:     clrl    r0
 318: loop8:  addl2   (r11),r0
 319:     emul    r8,-(r9),r0,r2
 320:     extzv   $0,$30,r2,(r11)
 321:     extv    $30,$32,r2,r0
 322:     acbl    r10,$-4,r11,loop8
 323:     ret
 324:     .set    .R2,0xf00
 325:     .align  1
 326:     .globl  _adback
 327: _adback:
 328:     .word   .R3
 329:     movl    4(ap),r11
 330:     movl    8(ap),r10
 331:     movl    12(ap),r9
 332:     clrl    r0
 333: loop9:  addl2   -(r9),r0
 334:     addl2   (r11),r0
 335:     extzv   $0,$30,r0,(r11)
 336:     extv    $30,$2,r0,r0
 337:     acbl    r10,$-4,r11,loop9
 338:     ret
 339:     .set    .R3,0xe00
 340:     .align  1
 341:     .globl  _dsdiv
 342: _dsdiv:
 343:     .word   .R4
 344:     movl    8(ap),r11
 345:     clrl    r0
 346: loopa:  emul    r0,$0x40000000,(r11),r1
 347:     ediv    12(ap),r1,(r11),r0
 348:     acbl    4(ap),$4,r11,loopa
 349:     ret
 350:     .set    .R4,0x800
 351:     .align  1
 352:     .globl  _dsmult
 353: _dsmult:
 354:     .word   .R5
 355:     movl    4(ap),r11
 356:     clrl    r0
 357: loopb:  emul    12(ap),(r11),r0,r1
 358:     extzv   $0,$30,r1,(r11)
 359:     extv    $30,$32,r1,r0
 360:     acbl    8(ap),$-4,r11,loopb
 361:     movl    r1,4(r11)
 362:     ret
 363:     .set    .R5,0x800
 364:     .align  1
 365:     .globl  _dsrsh
 366: _dsrsh:
 367:     .word   .R7
 368:     movl    8(ap),r11   # bot
 369:     movl    16(ap),r5   # mask
 370:     movl    12(ap),r4   # shift count
 371:     clrl    r0
 372: L201:   emul    r0,$0x40000000,(r11),r1
 373:     bicl3   r5,r1,r0
 374:     ashq    r4,r1,r1
 375:     movl    r1,(r11)
 376:     acbl    4(ap),$4,r11,L201
 377:     ret
 378:     .set    .R7,0x800
 379:     .align  1
 380:     .globl  _dsadd1
 381: _dsadd1:
 382:     .word   .R8
 383:     movl    4(ap),r11
 384:     movl    $1,r0
 385: L501:   emul    $1,(r11),r0,r1
 386:     extzv   $0,$30,r1,(r11)
 387:     extv    $30,$32,r1,r0
 388:     acbl    8(ap),$-4,r11,L501
 389:     movl    r1,4(r11)
 390:     ret
 391:     .set    .R8,0x800
 392: /*
 393: 	myfrexp (value, exp, hi, lo)
 394: 		double value;
 395: 		int *exp, *hi, *lo;
 396: 
 397: 	myfrexp returns three values, exp, hi, lo,
 398: 	Such that value = 2**exp*tmp, where tmp = (hi*2**-23+lo*2**-53)
 399: 	is uniquely determined subect to .5< abs(tmp) <= 1.0
 400: 
 401: 
 402: 	Entry point
 403: */
 404:     .text
 405:     .globl  _myfrexp
 406: _myfrexp:
 407:     .word   0x0000      # We use r2, but do not save it
 408: 
 409:     clrl    *12(ap)     # Make for easy exit later
 410:     clrl    *16(ap)     #
 411:     clrl    *20(ap)     #
 412:     movd    4(ap),r0    # Fetch "value"
 413:     bneq    L301        # if zero return zero exponent + mantissa
 414:     ret
 415: L301:
 416:     extzv   $7,$8,r0,r2 # r2 := biased exponent
 417:     movab   -129(r2),*12(ap)# subtract bias, store exp
 418:     insv    $154,$7,$8,r0   # lie about exponent to get out
 419:                 # high 24 bits easily with emodd.
 420: /*
 421: 	This instruction does the following:
 422: 
 423: 		Extend the long floating value in r0 with 0, and
 424: 		multiply it by 1.0.  Store the integer part of the result
 425: 		in hi, and the fractional part of the result in r0-r1.
 426: */
 427:     emodd   r0,$0,$0f1.0,*16(ap),r0 # How did you like
 428:                     # THAT, sports fans? [jfr's comment]
 429: 
 430:     tstd    r0      # if zero, exit;
 431:     bneq    L401
 432:     ret
 433: L401:
 434:     insv    $158,$7,$8,r0   # lie about exponent to get out
 435:                 # next 30 bits easily with emodd.
 436:                 # (2^29 takes 30 bits).
 437:     emodd   r0,$0,$0f1.0,*20(ap),r0 # Get last bits out likewise!
 438:     ret             # (r0 should be zero by now!)
 439:     .globl  _inewint
 440: _inewint:.word  0
 441:     movl    4(ap),r0
 442:     cmpl    r0,$1024
 443:     jgeq    Ialloc
 444:     cmpl    r0,$-1024
 445:     jlss    Ialloc
 446:     moval   _Fixzero[r0],r0
 447:     ret
 448: Ialloc:
 449:     calls   $0,_newint
 450:     movl    4(ap),0(r0)
 451:     ret
 452:     .globl  _blzero
 453: _blzero:                # blzero(where,howmuch)
 454:                     # char *where;
 455:                     # zeroes a block of length howmuch
 456:                     # beginning at where.
 457:     .word   0
 458:     movc5   $0,*4(ap),$0,8(ap),*4(ap)
 459:     ret

Defined functions