/* $Header: bigmath.c 1.4 83/06/09 00:50:06 sklower Exp $ */

#include "config.h"

	.globl	_dmlad
/*
	routine for destructive multiplication and addition to a bignum by
	two fixnums.

	from C, the invocation is dmlad(sdot,mul,add);
	where sdot is the address of the first special cell of the bignum
	mul is the multiplier, add is the fixnum to be added (The latter
	being passed by value, as is the usual case.


	Register assignments:

	r11 = current sdot
	r10 = carry
	r9  = previous sdot, for relinking.
*/
_dmlad:	.word	0x0e00
	movl	4(ap),r11		#initialize cell pointer
	movl	12(ap),r10		#initialize carry
loop:	emul	8(ap),(r11),r10,r0	#r0 gets cell->car times mul + carry
/*	ediv	$0x40000000,r0,r10,(r11)#cell->car gets prod % 2**30
					#carry gets quotient
*/
	extzv	$0,$30,r0,(r11)
	extv	$30,$32,r0,r10
	movl	r11,r9			#save last cell for fixup at end.
	movl	4(r11),r11		#move to next cell
	bneq	loop			#done indicated by 0 for next sdot
	tstl	r10			#if carry zero no need to allocate
	beql	done			#new bigit
	mcoml	r10,r3			#test to see if neg 1.
	bneq	alloc			#if not must allocate new cell.
	tstl	(r9)			#make sure product isn't -2**30
	beql	alloc
	movl	r0,(r9)			#save old lower half of product.
	brb	done
alloc:	jsb	_qnewdot			#otherwise allocate new bigit
	movl	r10,(r0)		#store carry
	movl	r0,4(r9)		#save new link cell
done:	movl	4(ap),r0
	ret
	.globl _dodiv
/*
	routine to destructively divide array representation of a bignum by 
	1000000000

	invocation:
		remainder = dodiv(top,bottom)
		int *top, *bottom;
	where *bottom is the address of the biggning of the array, *top is
	the top of the array.

	register assignments:
	r0 = carry
	r1 & r2 = 64bit temporary
	r3 = pointer
*/
_dodiv:	.word	0
	clrl		r0		#no carry to begin.
	movl		8(ap),r3	#get pointer to array.
loop2:	emul		$0x40000000,r0,(r3),r1
	ediv		$1000000000,r1,(r3),r0
	acbl		4(ap),$4,r3,loop2
	ret
	.globl	_dsneg
/*
	dsneg(top, bot);
	int *top, *bot;

	routine to destructively negate a bignum stored in array format
	lower order stuff at higher addresses. It is assume that the
	result will be positive.
*/
_dsneg:	.word	0
	movl	4(ap),r1	#load up address.
	clrl	r2		#set carry
loop3:	mnegl	(r1),r0		#negate and take carry into account.
	addl2	r2,r0
	extzv	$0,$30,r0,(r1)
	extv	$30,$2,r0,r2
	acbl	8(ap),$-4,r1,loop3
				#decrease r1, and branch back if appropriate.
	ret

/*
	bignum add routine
	basic data representation is each bigit is a positive number
	less than 2^30, except for the leading bigit, which is in
	the range -2^30 < x < 2^30.
*/

	.globl	_adbig
	.globl	Bexport
	.globl	backfr
/*
	Initialization section
*/
_adbig:	.word	0x0fc0		#save registers 6-11
	movl	4(ap),r1	#arg1 = addr of 1st bignum
	movl	8(ap),r2	#arg2 = addr of 2nd bignum
	clrl	r5		#r5   = carry
	movl	$0xC0000000,r4	#r4   = clear constant.
	movl	sp,r10		#save start address of bignum on stack.
				#note well that this is 4 above the actual
				#low order word.
/*
	first loop is to waltz through both bignums adding
	bigits, pushing them onto stack. 
*/
loop4:	addl3	(r1),(r2),r0	#add bigits
	addl2	r5,r0		#add carry
	bicl3	r4,r0,-(sp)	#save sum, no overflow possible
	extv	$30,$2,r0,r5	#sign extend two high order bits
				#to be next carry.
	movl	4(r1),r1	#get cdr
	bleq	out1		#negative indicates end of list.
	movl	4(r2),r2	#get cdr of second bignum
	bgtr	loop4		#if neither list at end, do it again
/*
	second loop propagates carries through higher order words.
	It assumes remaining list in r1.
*/
loop5:	addl3	(r1),r5,r0	#add bigits and carry
	bicl3	r4,r0,-(sp)	#save sum, no overflow possible
	extv	$30,$2,r0,r5	#sign extend two high order bits
				#to be next carry.
	movl	4(r1),r1	#get cdr
out2:	bgtr	loop5		#negative indicates end of list.
out2a:	pushl	r5
/*
	suppress unnecessary leading zeroes and -1's

	WARNING:  this code is duplicated in C in divbig.c
*/
iexport:movl	sp,r11		#more set up for output routine
ckloop:	
Bexport:tstl	(r11)		#look at leading bigit
	bgtr	copyit		#if positive, can allocate storage etc.
	blss	negchk		#if neg, still a chance we can get by
	cmpl	r11,r10		#check to see that
	bgeq	copyit		#we don't pop everything off of stack
	tstl	(r11)+		#incr r11
	brb	ckloop		#examine next
negchk:
	mcoml	(r11),r3		#r3 is junk register
	bneq	copyit		#short test for -1
	tstl	4(r11)		#examine next bigit
	beql	copyit		#if zero must must leave as is.
	cmpl	r11,r10		#check to see that
	bgeq	copyit		#we don't pop everything off of stack
	tstl	(r11)+		#incr r11
	bisl2	r4,(r11)	#set high order two bits
	brb	negchk		#try to supress more leading -1's
/*
	The following code is an error exit from the first loop
 	and is out of place to avoid a jump around a jump.
*/
out1:	movl	4(r2),r1	#get next addr of list to continue.
	brb	out2		#if second list simult. exhausted, do
				#right thing.
/*
	loop6 is a faster version of loop5 when carries are no
	longer necessary.
*/
loop6a: pushl	(r1)		#get datum
loop6:	movl	4(r1),r1	#get cdr
	bgtr	loop6a		#if not at end get next cell
	brb	out2a

/*
	create linked list representation of bignum
*/
copyit:	subl3	r11,r10,r2	#see if we can get away with allocating an int
	bneq	on1		#test for having popped everything
	subl3	$4,r10,r11	#if so, fix up pointer to bottom
	brb	intout		#and allocate int.
on1:	cmpl	r2,$4		#if = 4, then can do
	beql	intout
	calls	$0,_newsdot	#get new cell for new bignum
backfr:
#ifdef PORTABLE
	movl	r0,*_np
	addl2	$4,_np
#else
	movl	r0,(r6)+	#push address of cell on
				#arg stack to save from garbage collection.
				#There is guaranteed to be slop for a least 1
				#push without checking.
#endif
loop7:	movl	-(r10),(r0)	#save bigit
	movl	r0,r9		#r9 = old cell, to link
	cmpl	r10,r11		#have we copy'ed all the bigits?
	bleq	Edone
	jsb	_qnewdot	#get new cell for new bignum
	movl	r0,4(r9)	#link new cell to old
	brb	loop7
Edone:	
	clrl	4(r9)		#indicate end of list with 0
#ifdef PORTABLE
	subl2	$4,_np
	movl	*_np,r0
#else
	movl	-(r6),r0	#give resultant address.
#endif
	ret
/*
	export integer
*/
intout: pushl	(r11)
	calls	$1,_inewint
	ret
	.globl	_mulbig
/*
	bignum multiplication routine

	Initialization section
*/
_mulbig:.word	0x0fc0		#save regs 6-11
	movl	4(ap),r1	#get address of first bignum
	movl	sp,r11		#save top of 1st bignum
mloop1:	pushl	(r1)		#get bigit
	movl	4(r1),r1	#get cdr
	bgtr	mloop1		#repeat if not done
	movl	sp,r10		#save bottom of 1st bignum, top of 2nd bignum
	
	movl	8(ap),r1	#get address of 2nd bignum
mloop2:	pushl	(r1)		#get bigit
	movl	4(r1),r1	#get cdr
	bgtr	mloop2		#repeat if not done
	movl	sp,r9		#save bottom of 2nd bignum
	subl3	r9,r11,r6	#r6 contains sum of lengths of bignums
	subl2	r6,sp
	movl	sp,r8		#save bottom of product bignum
/*
	Actual multiplication
*/
m1:	movc5	$0,(r8),$0,r6,(r8)#zap out stack space
	movl	r9,r7		#r7 = &w[j +n] (+4 for a.d.) through calculation
	subl3	$4,r10,r4	#r4 = &v[j]

m3:	movl	r7,r5		#r7 = &w[j+n]
	subl3	$4,r11,r3	#r3 = &u[i]
	clrl	r2		#clear carry.

m4:	addl2	-(r5),r2	#add w[i + j] to carry (no ofl poss)
	emul	(r3),(r4),r2,r0 #r0 = u[i] * v[j] + sext(carry)
	extzv	$0,$30,r0,(r5)	#get new bigit
	extv	$30,$32,r0,r2	#get new carry

m5:	acbl	r10,$-4,r3,m4	#r3 =- 4; if(r3 >= r10) goto m4; r10 = &[u1];
	movl	r2,-(r5)	#save w[j] = carry

m6:	subl2	$4,r7		#add just &w[j+n] (+4 for autodec)
	acbl	r9,$-4,r4,m3	#r4 =- 4; if(r4>=r9) goto m5; r9 = &v[1]

	movl	r9,r10		#set up for output routine
	movl	$0xC0000000,r4	#r4   = clear constant.
	movq	20(fp),r6	#restor _np and _lbot !
	brw	iexport		#do it!
/*
 The remainder of this file are routines used in bignum division.
 Interested parties should consult Knuth, Vol 2, and divbig.c.
 These files are here only due to an optimizer bug.
*/
	.align	1
	.globl	_calqhat
_calqhat:
	.word	0xf00
	movl	4(ap),r11		# &u[j] into r11
	movl	8(ap),r10		# &v[1] into r10
	cmpl	(r10),(r11)		# v[1] == u[j] ??
	beql	L102			
	# calculate qhat and rhat simultaneously,
	#  qhat in r0
	#  rhat in r1
	emul	(r11),$0x40000000,4(r11),r4 # u[j]b+u[j+1] into r4,r5
	ediv	(r10),r4,r0,r1		# qhat = ((u[j]b+u[j+1])/v[1]) into r0
					# (u[j]b+u[j+1] -qhat*v[1]) into r1
					# called rhat
L101:
	# check if v[2]*qhat > rhat*b+u[j+2]
	emul	r0,4(r10),$0,r2		# qhat*v[2] into r3,r2
	emul	r1,$0x40000000,8(r11),r8 #rhat*b + u[j+2] into r9,r8
	# give up if r3,r2 <= r9,r8, otherwise iterate
	subl2	r8,r2			# perform r3,r2 - r9,r8
	sbwc	r9,r3
	bleq	L103			# give up if negative or equal
	decl	r0			# otherwise, qhat = qhat - 1
	addl2	(r10),r1		# since dec'ed qhat, inc rhat by v[1]
	jbr	L101
L102:	
	# get here if v[1]==u[j]
	# set qhat to b-1
	# rhat is easily calculated since if we substitute b-1 for qhat in
	# the formula, then it simplifies to (u[j+1] + v[1])
	# 
	addl3	4(r11),(r10),r1		# rhat = u[j+1] + v[1]
	movl	$0x3fffffff,r0		# qhat = b-1
	jbr	L101
	
L103:
	ret

	.align	1
	.globl	_mlsb
_mlsb:
	.word	.R2
	movl	4(ap),r11
	movl	8(ap),r10
	movl	12(ap),r9
	movl	16(ap),r8
	clrl	r0
loop8:	addl2	(r11),r0
	emul	r8,-(r9),r0,r2
	extzv	$0,$30,r2,(r11)
	extv	$30,$32,r2,r0
	acbl	r10,$-4,r11,loop8
	ret
	.set	.R2,0xf00
	.align	1
	.globl	_adback
_adback:
	.word	.R3
	movl	4(ap),r11
	movl	8(ap),r10
	movl	12(ap),r9
	clrl	r0
loop9:	addl2	-(r9),r0
	addl2	(r11),r0
	extzv	$0,$30,r0,(r11)
	extv	$30,$2,r0,r0
	acbl	r10,$-4,r11,loop9
	ret
	.set	.R3,0xe00
	.align	1
	.globl	_dsdiv
_dsdiv:
	.word	.R4
	movl	8(ap),r11
	clrl	r0
loopa:	emul	r0,$0x40000000,(r11),r1
	ediv	12(ap),r1,(r11),r0
	acbl	4(ap),$4,r11,loopa
	ret
	.set	.R4,0x800
	.align	1
	.globl	_dsmult
_dsmult:
	.word	.R5
	movl	4(ap),r11
	clrl	r0
loopb:	emul	12(ap),(r11),r0,r1
	extzv	$0,$30,r1,(r11)
	extv	$30,$32,r1,r0
	acbl	8(ap),$-4,r11,loopb
	movl	r1,4(r11)
	ret
	.set	.R5,0x800
	.align	1
	.globl	_dsrsh
_dsrsh:
	.word	.R7
	movl	8(ap),r11	# bot
	movl	16(ap),r5	# mask
	movl	12(ap),r4	# shift count
	clrl	r0
L201:	emul	r0,$0x40000000,(r11),r1
	bicl3	r5,r1,r0
	ashq	r4,r1,r1
	movl	r1,(r11)
	acbl	4(ap),$4,r11,L201
	ret
	.set	.R7,0x800
	.align	1
	.globl	_dsadd1
_dsadd1:
	.word	.R8
	movl	4(ap),r11
	movl	$1,r0
L501:	emul	$1,(r11),r0,r1
	extzv	$0,$30,r1,(r11)
	extv	$30,$32,r1,r0
	acbl	8(ap),$-4,r11,L501
	movl	r1,4(r11)
	ret
	.set	.R8,0x800
/*
	myfrexp (value, exp, hi, lo)
		double value;
		int *exp, *hi, *lo;

	myfrexp returns three values, exp, hi, lo,
	Such that value = 2**exp*tmp, where tmp = (hi*2**-23+lo*2**-53)
	is uniquely determined subect to .5< abs(tmp) <= 1.0
	

	Entry point
*/
	.text
	.globl	_myfrexp
_myfrexp:
	.word	0x0000		# We use r2, but do not save it

	clrl	*12(ap)		# Make for easy exit later
	clrl	*16(ap)		# 
	clrl	*20(ap)		# 
	movd	4(ap),r0	# Fetch "value"
	bneq	L301		# if zero return zero exponent + mantissa
	ret
L301:
	extzv	$7,$8,r0,r2	# r2 := biased exponent
	movab	-129(r2),*12(ap)# subtract bias, store exp
	insv	$154,$7,$8,r0	# lie about exponent to get out
				# high 24 bits easily with emodd.
/*
	This instruction does the following:

		Extend the long floating value in r0 with 0, and
		multiply it by 1.0.  Store the integer part of the result
		in hi, and the fractional part of the result in r0-r1.
*/
	emodd	r0,$0,$0f1.0,*16(ap),r0	# How did you like
					# THAT, sports fans? [jfr's comment]

	tstd	r0		# if zero, exit;
	bneq	L401
	ret
L401:
	insv	$158,$7,$8,r0	# lie about exponent to get out
				# next 30 bits easily with emodd.
				# (2^29 takes 30 bits).
	emodd	r0,$0,$0f1.0,*20(ap),r0	# Get last bits out likewise!
	ret				# (r0 should be zero by now!)
	.globl	_inewint
_inewint:.word	0
	movl	4(ap),r0
	cmpl	r0,$1024
	jgeq	Ialloc
	cmpl	r0,$-1024
	jlss	Ialloc
	moval	_Fixzero[r0],r0
	ret
Ialloc:
	calls	$0,_newint
	movl	4(ap),0(r0)
	ret
	.globl	_blzero
_blzero:				# blzero(where,howmuch)
					# char *where;
					# zeroes a block of length howmuch
					# beginning at where.
	.word	0
	movc5	$0,*4(ap),$0,8(ap),*4(ap)
	ret