![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to bigmath.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDReno / pgrm / lisp / franz / vax|
Return to bigmath.c CVS log | Up to [CSRG BSD Unix] / 43BSDReno / pgrm / lisp / franz / vax |
1.1 root 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
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.