```   1: /* \$Header: bigmath.c 1.4 83/06/09 00:50:06 sklower Exp \$ */
2:
3: #include "config.h"
4:
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
81:     clrl    r2      #set carry
82: loop3:  mnegl   (r1),r0     #negate and take carry into account.
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:
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: */
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
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
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
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
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

pushl defined in line 5; used 1 times
tstl defined in line 5; used 2 times
 Last modified: 1985-08-14 Generated: 2016-12-26 Generated by src2html V0.67 page hit count: 741