;; Model 100 Pie, TRS-80 Model 100 8085-based Pi/e
;;
;; Jan 21, 2013
;;
;; (c) 2012, All Right Researved, Egan Ford (egan@sense.net)
;;
;; THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY
;; KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
;; IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
;; PARTICULAR PURPOSE.
;;
;; Model 100 Pie will compute 1000 decimal digits of Pi or e.
;;
;; All 8085 Pie computation is base 256 and then converted to base 10 for
;; display (technically base 100, but they both display the same).  Arrays
;; of bytes are used to represent multiprecision numbers.  In the comments
;; below "array" refers to a big endian multiprecision number.
;;
;; C snippets:  J.W. Stumpel, (c) May 1991
;;
;; Assembler: Macro Assembler AS V1.42
;;
;; To compile/link:
;;    asl -L -gnuerrors -u -x -x -U m100pi.s
;;    p2bin -r \$-\$ m100pi.p pi.co
;;
;; To Install: ($D000-$DFFF used, can be changed by editing start =)
;;    From BASIC:
;;       clear256,53248
;;       loadm"pi
;;       menu
;;

				; Assembler: Macro Assembler AS V1.42
	cpu	8080		; Intel 8080

;; constants (R/O)

dec_len	=	1000		; 1000 decimal digits
bin_len	=	418		; ceil(1001 / log(256)) + 1 (+ 1 if odd) = 418
				;   Duff's Device 4x unrolled loops assume
				;   bin_len % 4 = 2
bin_end =	bin_len-1
				; counting up is easier
bcd_len	=	65536-(dec_len/2)

home	=	0422Dh		; home the cursor
cls	=	04231h		; clear screen
crout	=	04222h		; CR/LF out
cout	=	04B44h		; char out
getkey	=	012CBh		; get key
menu	=	05797h		; menu
gettime	=	0190Fh		; read time
start	=	0D000h		; start here


;; start of global macros/functions

; Description:  Return hi 8-bits of a label
;               e.g.: hi(label)

hi	function x,(x>>8)&255


; Description:  Return lo 8-bits of a label
;               e.g.: lo(label)

lo	function x,x&255


; Description:  Set initial value left of decimal of arg1 to arg2
; Input:        arg1 (16-bit), arg2 (8-bit)
;               arg1 = immediate address at start of array
;               arg2 = immediate integer to be set left of decimal point
; Calls:        set_mp
; Registers:    A, B, C, D, E, H, L

mset	macro	arg1, arg2	; arg1 = arg2
	lxi	d,arg1		; DE = address of a
	mvi	a,arg2		; load A with initial value
	call	set_mp
	endm


; Description:  Convert and print a base 256 number (ptr_a) as base 10/100.
; Input:        arg1 (16-bit) immediate address at start of an array
; Calls:        print_mp
; Registers:    A, B, C, D, E, H, L

mprint	macro	arg1		; print (arg1) base 100
	lxi	h,arg1		; HL <- ptr
	call	print_mp
	endm


; Description:  Divide (arg1) by an 8-bit quantity (arg1 = arg1 / arg2)
; Input:        arg1 (16-bit), arg2 (8-bit)
;               arg1 = immediate address at start of array
;               arg2 = immediate address of divisor
; Calls:        div_mp
; Return:       If arg1 = 0 (array of zeros), then carry is set,
;               otherwise carry is clear.
; Registers:    A, B, C, D, E, H, L

mdiv	macro	arg1, arg2	; arg1 = arg1 / arg2
	lxi	d,arg1		; load DE ptr to array
	lxi	b,arg2		; load BC ptr to divisor
	call	div_mp
	endm


; Description:  Divide (arg1) by a 16-bit quantity (arg1 = arg1 / arg2)
; Input:        arg1 (16-bit), arg2 (16-bit)
;               arg1 = immediate address at start of array
;               arg2 = immediate address of divisor
; Calls:        div16_mp
; Return:       If a = 0 (array of zeros), then carry is set,
;               otherwise carry is clear.
; Registers:    A, B, C, D, E, H, L

mdiv16	macro	arg1, arg2	; arg1 = arg1 / arg2
	lxi	d,arg1		; load DE ptr to array
	lxi	b,arg2		; load BC ptr to divisor
	call	div16_mp
	endm


; Description:  Transfer array (arg2)[0,sizeof(arg2)] to array (arg1)
; Input:        arg1 (16-bit), arg2 (16-bit)
;               arg1 = immediate address at start of array
;               arg2 = immediate address at start of array
; Calls:        copy_mp
; Registers:    A, B, C, D, E, H, L

mcopy	macro	arg1, arg2	; arg1 = arg2
	lxi	d,arg1		; DE = arg1
	lxi	b,arg2		; BC = arg2
	call	copy_mp
	endm


; Description:  Transfer array (arg2)[0,sizeof(arg2)] to array (arg1)
; Input:        arg1 (16-bit), arg2 (16-bit)
;               arg1 = direct address to start of array
;               arg2 = immediate address at start of array
; Calls:        copy_mp
; Registers:    A, B, C, D, E, H, L

mcopym	macro	arg1, arg2	; arg1 = arg2
	lhld	arg1		; HL = (arg1)
	xchg			; DE <-> HL
	lxi	b,arg2		; BC = arg2
	call	copy_mp
	endm


; Description:  Left shift array (arg1)
; Input:        arg1 (16-bit)
;               arg1 = immediate address at start of array
; Calls:        asl_mp
; Registers:    A, D, E, H, L

masl	macro	arg1		; arg1 = arg1 * 2
	lxi	h,arg1		; HL = arg1
	call	asl_mp
	endm


; Description:  Add array arg2 to array arg1 (arg1 = arg1 + arg2)
; Input:        arg1 (16-bit), arg2 (16-bit)
;               arg1 = immediate address at start of array
;               arg2 = immediate address at start of array
; Calls:        add_mp
; Registers:    A, B, C, D, E, H, L

madd	macro	arg1, arg2	; arg1 = arg1 + arg2
	lxi	b,arg1		; BC = arg1
	lxi	h,arg2		; HL = arg2
	call	add_mp
	endm


; Description:  Add array arg2 to array arg1 (arg1 = arg1 + arg2)
; Input:        arg1 (16-bit), arg2 (16-bit)
;               arg1 = direct address to start of array
;               arg2 = immediate address at start of array
; Calls:        add_mp
; Registers:    A, B, C, D, E, H, L

maddm	macro	arg1, arg2	; arg1 = arg1 + arg2
	lhld	arg1		; HL = (arg1)
	mov	b,h		; BC <- HL
	mov	c,l
	lxi	h,arg2		; HL = arg2
	call	add_mp
	endm


; Description:  Subtract array b from array a (a = a - b)
; Input:        arg1 (16-bit), arg2 (16-bit)
;               arg1 = immediate address at start of array
;               arg2 = immediate address at start of array
; Calls:        sub_mp
; Registers:    A, B, C, D, E, H, L

msub	macro	arg1, arg2	; arg1 = arg1 - arg2
	lxi	b,arg1		; BC = arg1
	lxi	h,arg2		; HL = arg2
	call	sub_mp
	endm


; Description:  Subtract array b from array a (a = a - b)
; Input:        arg1 (16-bit), arg2 (16-bit)
;               arg1 = direct address to start of array
;               arg2 = immediate address at start of array
; Calls:        sub_mp
; Registers:    A, B, C, D, E, H, L

msubm	macro	arg1, arg2	; arg1 = arg1 - arg2
	lhld	arg1		; HL = (arg1)
	mov	b,h		; BC <- HL
	mov	c,l
	lxi	h,arg2		; HL = arg2
	call	sub_mp
	endm


; Description:  Compute arctan(1/arg2) and store at (arg1).
; Input:        arg1 (16-bit), arg2 (8-bit)
;               arg1 = immediate address at start of array
;               arg2 = immediate integer to be set left of decimal point
; Calls:        atan_mp
; Registers:    A, B, C, D, E, H, L

matan	macro	arg1, arg2	; arg1 = atan(1/arg2)
	lxi	h,arg1		; HL = arg1
	mvi	a,arg2		; A = arg2
	call	atan_mp
	endm


;; header m100

	org	start-6

	dw	start	
	dw	end-main
	dw	start

;; start of main code

	org	start
main:
	call	cls
	call	home

	lxi	h,stime
	call	gettime

	lxi	h,pitext
	call	print

	call	pi
	;call	e
	mprint	mp_a
	call	crout
	call	crout

	lxi	h,etime
	call	gettime

	lxi	h,stimel
	call	print
	lxi	h,etimel
	call	print
	call	crout

	lxi	h,msg
	call	print

	call	getkey
	jmp	menu


;; start of variables

bcd_cnt:			; counter used for print_mp
	db	0, 0		; 0-1
stk_ptr:			; 16-bits used to backup stack ptr
beg_ptr:			; 16-bits used to store pointers to arrays
	db	0, 0		; 2-3
end_ptr:			; 16-bits used to store pointers to arrays
	db	0, 0		; 4-5
n:	db	0, 0		; 6-7, used by atan_mp
regx:	db	0		; B, the x in atan(1/x)
x2:	db	0		; C, cached x*x
ptr_a:	db	0, 0		; D-E, pointer to array used by atan_mp
				; table of squares
sqrtbl:	db	0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100
	db	121, 144, 169, 196, 225
				; Change PI to E for, well, e
pitext:	db	"1000 DIGITS OF PI = \0"
stimel:	db	"Start Time: "
stime:	db	0,0,0,0,0,0,0,0
	db	"\r\n\0"
etimel:	db	"  End Time: "
etime:	db	0,0,0,0,0,0,0,0
	db	"\r\n\0"
msg:	db	"Send results to egan@sense.net,\r\npress any key to exit. \0"


;; start of subs

; Description:  Compute e using the infinite series:
;
;               oo
;               =====
;               \      1
;       e   =    >    --  = 1 + 1/1! + 1/2! + 1/3! + 1/4! + ...
;               /     k!
;               =====
;               k = 0
;
;               and save in array mp_a.
;
; Input:        None
;
; Output:       mp_a = e = 1 + 1/1! + 1/2! + 1/3! + 1/4! + ...
;
; Registers:    A, B, C, D, E, H, L
;
; Calls:        mset/set_mp, mdiv/div_mp, mdiv16/div16_mp, madd/add
;
; Globals:      n, mp_a, mp_x
;
; C Algorithm:
;
; int n = 1;
; setbig(A, 1, 0);
; setbig(X, 1, 0);
; while(divbig(X, n++))         // while(dividend != 0)
;     addbig(A, X);

e:	lxi	h,0		; n = 0 little endian
	shld	n		; n = HL
	mset	mp_a,1		; a = 1
	mset	mp_x,1		; x = 1

e_loop:
				; n = n + 1
	lhld	n		; HL = n
	inx	h		; n++
	shld	n		; n = HL
	lda	n+1		; if n < 256 then div else div16
	ana	a		; update zero flag
	jnz	+		; >= 256 use 16 bit div
	mdiv	mp_x,n		; < 256 use 8 bit div
	rc			; dividend = 0, done
	jmp	++		; still working on it
/	mdiv16	mp_x,n		; div16
	rc			; dividend = 0, done
/	madd	mp_a,mp_x	; a = a + x
	jmp	e_loop


; Description:  Compute pi using the Gregory expansion of Machin's arctan
;               formula and save in array mp_a.
;
;   pi = 4 * (4 *      atan(1/5)           -           atan(1/239)          )
;
;
;   __      /   / 1     1      1       \       / 1       1        1        \\
;   || = 4 | 4 |  - - ---- + ---- - ... |  -  | --- - ------ + ------ - ... ||
;          |   |  5      3      5       |     | 239        3        5       ||
;           \   \     3x5    5x5       /       \      3x239    5x239       //
;
;
; Input:        None
;
; Output:       mp_a = pi = 4 * (4 * atan(1/5) - atan(1/239))
;
; Registers:    A, B, C, D, E, H, L
;
; Calls:        matan/atan_mp, masl/asl_mp, msub/sub_mp

pi:	matan	mp_a,5		; a = atan(1/5)
	masl	mp_a		; a = a * 4
	masl	mp_a
	matan	mp_b,239	; b = atan(1/239)
	msub	mp_a,mp_b	; a = a - b
	masl	mp_a		; a = a * 4
	masl	mp_a
	ret


; Description:  Compute arctan(1/N) using the Gregory expansion of Machin's
;               arctan formula and save in array (ptr_a).
;
;
;                               / 1     1      1       \
;               arctan(1/N) =  |  - - ---- + ---- - ... |
;                              |  N      3      5       |
;                               \     3xN    5xN       /
;
;
; Input:        H/L (hi/lo) pointer to array
;               B = N (8 bit)
;
; Output:       H/L (hi/lo) pointer to array = arctan(1/N)
;
; Registers:    A, B, C, D, E, H, L
;
; Calls:        mset/set_mp, mdiv/div_mp, maddm/add_mp, msubm/sub_mp,
;               mdiv16/div16_mp, mcopy(m)/copy_mp
;
; Globals:      mp_x (16-bit), mp_y (16-bit)
;
; C Algorithm:
;
; void atanbig(bignum A, unsigned short x)
; {
;     bignum X, Y;
;     unsigned short n = 1;
;
;     setbig(X, 1, 0);
;     divbig(X, x);
;     copybig(A, X);
;     x *= x;
;     while (1) {
;         n += 2;
;         divbig(X, x);
;         copybig(Y, X);
;         if (!divbig(Y, n))      // dividend = 0
;             break;
;         if (n & 2)
;             subbig(A, Y);
;         else
;             addbig(A, Y);
;     }
; }

atan_mp:
	shld	ptr_a		; ptr_a to point to array
	sta	regx		; save A

	cpi	16		; if A < 16 then get square
	jnc	+
	lxi	h,sqrtbl	; HL points to sqrtbl
	add	l		; A += L
	mov	l,a		; L <- A
	mov	a,m		; load A with square
/	sta	x2		; x2 = A

	lxi	h,1		; n = 1 little endian
	shld	n

	mset	mp_x,1		; x = 1
	mdiv	mp_x,regx	; x /= regx
	mcopym	ptr_a,mp_x	; a = x

atan_mp_loop:			; main loop
				; n = n + 2
	lhld	n		; HL = n
	inx	h		; n++
	inx	h		; n++
	shld	n		; n = HL

	lda	regx		; A <- regx
	cpi	16		; if > 16 /x/x othersize /x2
	jc	+		; x already x*x, one div required
	mdiv	mp_x,x2		; x >= 16, then do two div
/	mdiv	mp_x,x2		;   (faster than 16 bit div)

	mcopy	mp_y,mp_x	; y = x

	lda	n+1		; if n < 256 then div else div16
	ana	a		; update zero flag
	jnz	+		; >= 256 use 16 bit div
	mdiv	mp_y,n		; < 256 use 8 bit div
	rc			; dividend = 0, done
	jmp	++		; still working on it
/	mdiv16	mp_y,n		; div16
	rc			; dividend = 0, done

/	lda	n		; add or sub?
	ani	2		; odd/even check on 2nd bit
	jz	+		; add it
	msubm	ptr_a,mp_y	; a = a - y
	jmp	atan_mp_loop	; back to top
/	maddm	ptr_a,mp_y	; a = a + y
	jmp	atan_mp_loop	; back to top


; Description:  Multiprecision subtraction:  a = a - b
;
; Input:        H/L (hi/lo) address to array (a)
;               B/C (hi/lo) address to array (b)
;
; Output:       a = a - b
;
; Registers:    A, B, C, D, E, H, L

sub_mp:
				; BC = BC - HL
	lxi	d,bin_len	; DE = bin_len
	mov	a,c		; BC += DE (end address)
	add	e
	mov	c,a
	mov	a,b
	adc	d
	mov	b,a
	dad	d		; HL += DE (end address)
	lxi	d,65536-bin_len	; set DE to roll to 0 when done
	ora	a		; clear carry
	jmp	++		; Duff's Device hardcoded for fixed array length
/	dcx	b		; BC--
	ldax	b		; A <- (BC)
	dcx	h		; HL--
	sbb	m		; A - (HL) - carry
	stax	b		; (BC) <- A
	inr	e		; E += 1
	dcr	c		; C--
	ldax	b		; A <- (BC)
	dcr	l		; L--
	sbb	m		; A - (HL) - carry
	stax	b		; (BC) <- A
	inr	e		; E += 1
/	dcr	c		; C--
	ldax	b		; A <- (BC)
	dcr	l		; L--
	sbb	m		; A - (HL) - carry
	stax	b		; (BC) <- A
	inr	e		; E += 1
	dcr	c		; C--
	ldax	b		; A <- (BC)
	dcr	l		; L--
	sbb	m		; A - (HL) - carry
	stax	b		; (BC) <- A
	inr	e		; E += 1
	jnz	--		; E != 0 jmp back to //
	inr	d		; D += 1
	jnz	--		; D != 0 jmp back to //
	ret


; Description:  Multiprecision addition:  a = a + b
;
; Input:        B/C (hi/lo) address to array (a)
;               H/L (hi/lo) address to array (b)
;
; Output:       a = a + b
;
; Registers:    A, B, C, D, E, H, L

add_mp:
				; BC = BC + HL
	lxi	d,bin_len	; DE = bin_len
	mov	a,c		; BC += DE (end address)
	add	e
	mov	c,a
	mov	a,b
	adc	d
	mov	b,a
	dad	d		; HL += DE (end address)
	lxi	d,65536-bin_len	; set DE to roll to 0 when done
	ora	a		; clear carry
	jmp	++		; Duff's Device hardcoded for fixed array length
/	dcx	b		; BC--
	ldax	b		; A <- (BC)
	dcx	h		; HL--
	adc	m		; A + (HL) + carry
	stax	b		; (BC) <- A
	inr	e		; E += 1
	dcr	c		; C--
	ldax	b		; A <- (BC)
	dcr	l		; L--
	adc	m		; A + (HL) + carry
	stax	b		; (BC) <- A
	inr	e		; E += 1
/	dcr	c		; C--
	ldax	b		; A <- (BC)
	dcr	l		; L--
	adc	m		; A + (HL) + carry
	stax	b		; (BC) <- A
	inr	e		; E += 1
	dcr	c		; C--
	ldax	b		; A <- (BC)
	dcr	l		; L--
	adc	m		; A + (HL) + carry
	stax	b		; (BC) <- A
	inr	e		; E += 1
	jnz	--		; E != 0 jmp back to //
	inr	d		; D += 1
	jnz	--		; D != 0 jmp back to //
	ret


; Description:  Multiprecision left shift:  a = a * 2
;
; Input:        H/L (hi/lo) address to array (a)
;
; Output:       a = a * 2
;
; Registers:    A, D, E, H, L

asl_mp:
	lxi	d,bin_len	; DE = bin_len
	dad	d		; HL += DE (end address)
	lxi	d,65536-bin_len	; set DE to roll to 0 when done
	ora	a		; clear carry
	jmp	++		; Duff's Device hardcoded for fixed array length
/	dcx	h		; HL--
	mov	a,m		; A <- (HL)
	ral			; rol
	mov	m,a		; (HL) <- A
	inr	e		; E += 1
	dcr	l		; L--
	mov	a,m		; A <- (HL)
	ral			; rol
	mov	m,a		; (HL) <- A
	inr	e		; E += 1
/	dcr	l		; L--
	mov	a,m		; A <- (HL)
	ral			; rol
	mov	m,a		; (HL) <- A
	inr	e		; E += 1
	dcr	l		; L--
	mov	a,m		; A <- (HL)
	ral			; rol
	mov	m,a		; (HL) <- A
	inr	e		; E += 1
	jnz	--		; E != 0 jmp back to //
	inr	d		; D += 1
	jnz	--		; D != 0 jmp back to //
	ret


; Description:  Multiprecision transfer:  a = b
;
; Input:        B/C (hi/lo) address to array (a)
;               H/L (hi/lo) address to array (b)
;
; Output:       a = b
;
; Registers:    A, B, C, D, E, H, L

copy_mp:
	lxi	h,bin_len	; HL = bin_len
	dad	d		; HL += DE (end address)
	xchg			; DE <-> HL, HL points to array a
	jmp	++		; Duff's Device hardcoded for fixed array length
				;   bin_len % 4 = 2
/	ldax	b		; A <- (BC)
	mov	m,a		; (HL) <- A
	inr	c		; C++
	inr	l		; L++
	ldax	b		; A <- (BC)
	mov	m,a		; (HL) <- A
	inx	b		; BC++
	inx	h		; HL++
/	ldax	b		; A <- (BC)
	mov	m,a		; (HL) <- A
	inr	c		; C++
	inr	l		; L++
	ldax	b		; A <- (BC)
	mov	m,a		; (HL) <- A
	inr	c		; C++
	inr	l		; L++
	mov	a,e		; A <- E
	cmp	l		; E = L?
	jnz	--		; E != L jmp back to //
	mov	a,d		; A <- D
	cmp	h		; D = H?
	jnz	--		; D != H jmp back to //
	ret


; Description:  Skip leading zeros (used by div_mp and div16_mp)
;
; Input:        H/L (hi/lo) address to array (a)
;
; Output:       None
;
; Registers:    A, D, E, H, L
;
; Globals:      end_ptr

skipzeros:
				; DE points to array
	lxi	h,bin_len	; HL = bin_len
	dad	d		; HL += DE
	shld	end_ptr		; store HL to end_ptr
	xchg			; DE <-> HL
				; HL points to array
				; DE = end of array
	jmp	++		; Duff's Device hardcoded for fixed array length
				;   bin_len % 4 = 2
/	mov	a,m		; A <- M
	ana	a		; update zero flag
	jnz	even
	inr	l		; L += 1
	mov	a,m		; A <- M
	ana	a		; update zero flag
	jnz	odd
	inx	h		; HL += 1
/	mov	a,m		; A <- M
	ana	a		; update zero flag
	jnz	even
	inr	l		; L += 1
	mov	a,m		; A <- M
	ana	a		; update zero flag
	jnz	odd
	inr	l		; L += 1
	mov	a,e		; A <- E
	cmp	l		; E = L?
	jnz	--		; E != L jmp back to //
	mov	a,d		; A <- D
	cmp	h		; D = H?
	jnz	--		; D != H jmp back to //
	stc			; set carry
	ret			;   all zeros, return
odd:	dcx	h		; HL--
even:	ora	a		; clear carry
	ret


; Description:  Multiprecision 16-bit division:  a = a / b
;
; Input:        H/L (hi/lo) address to array (a)
;               B/C (hi/lo) address to divisor (b)
;
; Output:       a = a / b
;
; Registers:    A, B, C, D, E, H, L, SP
;
; Globals:      end_ptr, stk_ptr
;
; Calls:        skipzeros
;
; C Algorithm:  see div_mp

div16_mp:
	call	skipzeros	; skip leading zeros for speed
	rc			; return if carry set (all zeros)

				; HL pointing to array MS 16-bit Digit
				; A, D, E disposable
				; BC has pointer to divisor

	xchg			; backup HL to DE
	di			; disable interrupts
				; backup stack ptr
	lxi	h,0		; HL = 0
	dad	sp		; HL += SP
	shld	stk_ptr		; store SP to stk_ptr
	xchg			; restore HL from DE
	sphl			; store HL to stack pointer
				; load divisor in DE
	ldax	b		; A <- (BC)
	mov	e,a		; E <- A lo(BC)
	inx	b		; BC++
	ldax	b		; A <- (BC)
	mov	d,a		; D <- A hi(BC)

	mvi	b,0		; set carry/remainder to 0, dividend+0
	mvi	c,0		; set carry/remainder to 0, dividend+1

				; SP pointing to array MS 16-bit digit
				; BC to be 16-bit carry/remainder/dividend
				; DE to be divisor
				; BCDE off limits, need for chain
				; BCLH = dividend+0-3

/	pop	h		; load 16-bit pair into HL
				; 32/16 division
	rept	16		; do it, 16 times
	mov	a,h		; A <- H dividend+3
	add	a		; asl (A = A + A)
	mov	h,a		; H <- A dividend+3
	mov	a,l		; A <- L dividend+2
	ral			; rol
	mov	l,a		; L <- A dividend+2
	mov	a,c		; A <- C dividend+1
	ral			; rol
	mov	c,a		; C <- A dividend+1
	mov	a,b		; A <- B dividend+0
	ral			; rol
	mov	b,a		; B <- A dividend+0
				; trial division, nothing changed
	mov	a,c		; A <- C
	sub	e		; C = C - E (dividend+1 - divisor(lo))
	mov	a,b		; A <- B
	sbb	d		; B = B - D (dividend+0 - divisor(hi)) - carry
	jc	m0		; too small
				; do it again, but this time save it
	mov	b,a		; B <- A, (hi) done, just save it
	mov	a,c		; A <- C, (lo)
	sub	e		; C = C - E (dividend+1 - divisor(lo))
	mov	c,a		; C <- A
	inr	h		; H++, dividend+3
m0:
	endm			; end 32/16 division

	push	h		; store 16-bit LSB pair HL
	pop	h		; SP += 2
				; check SP = end_ptr, load SP in HL
	lxi	h,0		; HL = 0
	dad	sp		; HL += SP
				; compare HL to end_ptr
	lda	end_ptr		; A <- lo(end_ptr)
	cmp	l		; A = L?
	jnz	-
	lda	end_ptr+1	; A <- hi(end_ptr)
	cmp	h		; A = H?
	jnz	-
				; restore stack ptr
	lhld	stk_ptr		; HL = (stk_ptr)
	sphl			; SP = HL
	ei			; enable interrupts
	ora	a		; clear carry
	ret


; Description:  Multiprecision 8-bit division:  a = a / b
;
; Input:        H/L (hi/lo) address to array (a)
;               B/C (hi/lo) address to divisor (b)
;
; Output:       a = a / b
;
; Registers:    A, B, C, D, E, H, L
;
; Globals:      end_ptr
;
; Calls:        skipzeros, div8(macro)
;
; C Algorithm:
;
; short divbig(number, x)
; bignum number;
; unsigned short x;
; {
;     dword result;
;     short j = 0;
;     unsigned short rest = 0;
;
;     while (number[j] == 0 && j < MAXSIZE)
;         j++;
;     if (j == MAXSIZE)
;         return (0);
;     while (j < MAXSIZE) {
;         result.w.lo = number[j];
;         result.w.hi = rest;
;         number[j] = result.L / x;
;         rest = result.L % x;
;         j++;
;     }
;     return (1);
; }

div8	macro {noexpand}	; Input:  B = lo, C = hi, D = divisor
				; Output: B = lo, C = hi
	mov	b,m		; B <- M, dividend lo
	rept	8		; do it 8 times
	mov	a,b		; A <- B
	add	a		; asl (A = A + A)
	mov	b,a		; B <- A, dividend lo
	mov	a,c		; A <- C, dividend hi
	ral			; rol A (hi)
	jc	m0		; 9th bit hi?  Sub it.
	cmp	d		; A - D, if A < D set carry
	jc	m1		; too small
m0:	sub	d		; A = A - D
	inr	b		; B++
m1:	mov	c,a		; C <- A
	endm
	mov	m,b		; M <- B save it
	endm

div_mp:
	call	skipzeros	; skip leading zeros for speed
	rc			; return if carry set (all zeros)
				; BC points to divisor
	ldax	b		; A = lo(divisor)
	mov	d,a		; D = divisor
	mvi	c,0		; set carry to 0
				; we know HL is even and bin_len is also even
				;   so unroll twice
/	div8			; div M / D
	inr	l		; L++
	div8			; div M / D
	inx	h		; HL++
	mov	a,e		; A <- E
	cmp	l		; E = L? (lo)
	jnz	-		;   if not back to /
				; check if H = end_ptr hi
	lda	end_ptr+1	; A <- end_ptr+1 (hi)
	cmp	h		; does A = end_ptr(hi)
	jnz	-		; if not match back //
	ora	a		; clear carry
	ret


; Description:  Print an mp array base 10/100
;
; Input:        H/L (hi/lo) address to array (a)
;
; Output:       (a) base 10/100 out to screen
;
; Registers:    A, B, C, D, E, H, L
;
; Calls:        bintobcd, prbcd, mult100(macro)
;
; Globals:      beg_ptr, end_ptr, bcd_cnt, mt100_lo, mt100_hi

aslb	macro {noexpand}
	mov	a,b		; A <- B
	add	a		; asl (A = A + A)
	mov	b,a		; B <- A
	endm

rolc	macro {noexpand}
	mov	a,c		; A <- C
	ral			; rol
	mov	c,a		; C <- A
	endm

mult100_slow	macro {noexpand}	; 100 = 1100100B, so
				;   a = ((x*2+x)*8+x)*4
				; IOW, shift and add,
				; 1 first MSB, is the initial load
				; 1 shift (2x) and add +x
				; 0 shift (2x)
				; 0 shift (2x)
				; 1 shift (2x) and add +x
				; 0 shift (2x)
				; 0 shift (2x)
				; end of digits

				; 1 initial load
	mvi	c,0		; C = MSB
	mov	a,m		; B/A = LSB
				; 1 shift (*2) ...
	add	a		; *2, asl (A = A + A)
	mov	b,a		; B <- A
	rolc 			; A <- C, * 2,  C <- A
				; ... and add +x
	mov	a,b		; A <- B
	add	m		; A = A + M
	jnc	m1		; no carry
	inr	c		; inc MSB
m1:				; 0 shift (*2)
	add	a		; *2, asl (A = A + A)
	mov	b,a		; B <- A
	rolc 			; A <- C, * 2,  C <- A
				; 0 shift (*2)
	aslb			; A <- B, A & A (clear carry), * 2, B <- A
	rolc 			; A <- C, * 2,  C <- A
				; 1 shift (*2) ...
	aslb			; A <- B, A & A (clear carry), * 2, B <- A
	rolc 			; A <- C, * 2,  C <- A
				; ... and add +x	
	mov	a,b		; A <- B
	add	m		; A = A + M
	jnc	m2		; no carry
	inr	c		; inc MSB
m2:				; 0 shift (*2)
	add	a		; *2, asl (A = A + A)
	mov	b,a		; B <- A
	rolc 			; A <- C, * 2,  C <- A
				; 0 shift (*2)
	aslb			; A <- B, A & A (clear carry), * 2, B <- A
	rolc 			; A <- C, * 2,  C <- A
				; end of x 100
				; add carry for chain computation
	mov	a,b		; A <- B
	add	e		; A = A + E (carry_mp)
	jnc	m3		; no carry
	inr	c		; inc MSB
m3:	mov	m,a		; M <- A, store LSB
	mov	e,c		; E <- C, tranfer MSB
	endm

mult100	macro
	mov	a,m		; get A (page aligned table offset)
	mvi	b,hi(mt100_lo)	; BC to point to low table
	mov	c,a		; C = A (position of A*100)
	ldax	b		; A <- (BC), LSB A = A*100
	add	e		; LSB A = A + E, add mp carry
	mov	m,a		; M <- A, store LSB
	inr	b		; B++ to get MSB from table
	ldax	b		; A <- (BC), MSB A = A*100
	mov	e,a		; E <- A, store mp carry
	jnc	m0		; if carry set from add, then
	inr	e		;   E++
m0:
	endm

print_mp:
				; HL now set to MSB
	mov	a,m		; load A with MSB
	cpi	10		; is < 10, then use ASCII method
	jc	+		;   otherwise, use 2 digit BCD method
	call	bintobcd	; convert to BCD
	call	prbcd		; print it
	jmp	++		; skip over ASCII method
/	adi	030h		; make ASCII
	call	echo		; print it
/	mvi	a,'.'		; print decimal point
	call	echo		; print it
	mvi	m,0		; zero MSB
	shld	beg_ptr		; save HL value to location beg_ptr
				; compute end of array since we have
				;   to go backwards
	mov	d,l		; D = lo(beg_ptr)
	lxi	b,bin_len	; BC = bin_len
	dad	b		; HL += BC
	shld	end_ptr		; save HL value to location end_ptr
				; setup bcd counter for outer loop
	lxi	h,bcd_len	; load HL with bcd_cnt
	shld	bcd_cnt		; save HL value to location bcd_cnt

print_mp1:			; main loop
				; multi array x 100
				; loop from LSB to MSB
				; 16-bit product = array[i] * 100 + carry_mp;
				; array[i] = product lo
				; carry_mp = product hi
				;
				; C Algorithm:
				;
				; while (j >= 0) {
				;    result.L = (long) number[j] * 100 + carry;
				;    number[j] = result.w.lo;
				;    carry = result.w.hi;
				;    j--;
				; }

	lhld	end_ptr		; load HL value from location end_ptr
				; backwards loop (LSB to MSB)
	mvi	e,0		; E = carry = 0
				; E is the carry from mult100 to mult100
				; A,B,C,E (carry) used for mult100
				; D and E must not be touched!
	jmp	++		; mp_array * 100 loop
/	dcx	h		; HL--
	mult100			; M = lo(M*100), E = hi(M*100)
	dcr	l		; L--
	mult100			; M = lo(M*100), E = hi(M*100)
/	dcr	l		; L--
	mult100			; M = lo(M*100), E = hi(M*100)
	dcr	l		; HL--
	mult100			; M = lo(M*100), E = hi(M*100)
	mov	a,d		; A <- D
	cmp	l		; D = L? (lo)
	jnz	--
				; check if H = beg_ptr hi
	lda	beg_ptr+1	; A <- beg_ptr+1 (hi)
	cmp	h		; does A = beg_ptr(hi)
	jnz	--		; if not match back //

				; HL now set to MSB
	mov	a,m		; load A with MSB
	call	bintobcd	; convert to BCD
	call	prbcd		; print base 100 digit
	mvi	m,0		; zero MSB

				; check for all decimal digits printed
	lhld	bcd_cnt		; load HL with bcd_cnt value
				;   works with bcd_len = 65536-(dec_len/2)
				;   counting up is easier
	inx	h		; HL++
	shld	bcd_cnt		; store HL into bcd_cnt
	mov	a,l		; A <- L, check lo
	ora	h		; or with hi
	jnz	print_mp1	; if not zero keep going, otherwise
	ret
	

; Description:  Set mp array initial value:  a = b
;
; Input:        H/L (hi/lo) address to array (a)
;               B = value left of decimal
;
; Output:       a[0] = B
;               a[1] - a[length-1] = 0
;
; Registers:    A, B, D, E, H, L

set_mp:
	stax	d		; set inital value, rest of array will be zero
	lxi	h,bin_end	; HL = bin_end
	dad	d		; HL += DE (end address)
	xra	a		; a = 0
	jmp	++		; Duff's Device hardcoded for fixed array length
				;   bin_len % 4 = 2
/	xra	a		; a = 0
	inr	e		; E += 1
	stax	d		; (DE) = 0
	inr	e		; E += 1
	stax	d		; (DE) = 0
	inx	d		; DE += 1
	stax	d		; (DE) = 0
/	inr	e		; E += 1
	stax	d		; (DE) = 0
	mov	a,e		; A <- E
	cmp	l		; E = L?
	jnz	--		; E != L jmp back to //
	mov	a,d		; A <- D
	cmp	h		; D = H?
	jnz	--		; D != H jmp back to //
	ret


; Description:  Output BCD number to display.
;
; Input:        A = BCD number (00 - 99)
;
; Output:       BCD number to display
;
; Registers:    A, B

prbcd:	mov	b,a		; B <- A, backup A
	rrc			; shift right 4 bits
	rrc
	rrc
	rrc
	call	+		; do MSB
	mov	a,b		; then LSB
/	ani	00Fh		; mask off
	ori	030h		; make a number
	call	echo		; print it
	ret


; Description:  Convert BIN/HEX to BCD.
;
; Input:        A = BIN/HEX (00h - 63h)
;
; Output:       A = BCD (00 - 99)
;
; Registers:    A, B, C

bintobcd:
	mov	b,a		; B <- A, backup A
	mvi	c,-1		; start counter at -1
/	inr	c		; C=C+1
	sui	10		; A=A-10
	jp	-		; loop if positive
	mov	a,b		; A <- B, restore A
	jmp	++		; jump to dcc
/	adi	6		; A=A+6
/	dcr	c		; if a < 10, then c will be -1 and end loop
	jp	--		;  othersize add 6 for every 10
	ret


; Description:  Output a string.
;
; Input:        H/L (hi/lo) address to array (a)
;
; Output:       String to screen
;
; Registers:    A, H, L

print:	mov	a,m		; A <- M (HL)
	ana	a		; update zero flag
	rz			; return if zero
	call	echo		; print char
	inx	h		; HL++
	jmp	print		; back to top


echo:	call	cout
	ret


;; tables for faster math must be page aligned

	rept	256-$&0FFh
	db	0
	endm

mt100_lo:
	db	000h,064h,0C8h,02Ch,090h,0F4h,058h,0BCh,020h,084h,0E8h,04Ch,0B0h,014h,078h,0DCh
	db	040h,0A4h,008h,06Ch,0D0h,034h,098h,0FCh,060h,0C4h,028h,08Ch,0F0h,054h,0B8h,01Ch
	db	080h,0E4h,048h,0ACh,010h,074h,0D8h,03Ch,0A0h,004h,068h,0CCh,030h,094h,0F8h,05Ch
	db	0C0h,024h,088h,0ECh,050h,0B4h,018h,07Ch,0E0h,044h,0A8h,00Ch,070h,0D4h,038h,09Ch
	db	000h,064h,0C8h,02Ch,090h,0F4h,058h,0BCh,020h,084h,0E8h,04Ch,0B0h,014h,078h,0DCh
	db	040h,0A4h,008h,06Ch,0D0h,034h,098h,0FCh,060h,0C4h,028h,08Ch,0F0h,054h,0B8h,01Ch
	db	080h,0E4h,048h,0ACh,010h,074h,0D8h,03Ch,0A0h,004h,068h,0CCh,030h,094h,0F8h,05Ch
	db	0C0h,024h,088h,0ECh,050h,0B4h,018h,07Ch,0E0h,044h,0A8h,00Ch,070h,0D4h,038h,09Ch
	db	000h,064h,0C8h,02Ch,090h,0F4h,058h,0BCh,020h,084h,0E8h,04Ch,0B0h,014h,078h,0DCh
	db	040h,0A4h,008h,06Ch,0D0h,034h,098h,0FCh,060h,0C4h,028h,08Ch,0F0h,054h,0B8h,01Ch
	db	080h,0E4h,048h,0ACh,010h,074h,0D8h,03Ch,0A0h,004h,068h,0CCh,030h,094h,0F8h,05Ch
	db	0C0h,024h,088h,0ECh,050h,0B4h,018h,07Ch,0E0h,044h,0A8h,00Ch,070h,0D4h,038h,09Ch
	db	000h,064h,0C8h,02Ch,090h,0F4h,058h,0BCh,020h,084h,0E8h,04Ch,0B0h,014h,078h,0DCh
	db	040h,0A4h,008h,06Ch,0D0h,034h,098h,0FCh,060h,0C4h,028h,08Ch,0F0h,054h,0B8h,01Ch
	db	080h,0E4h,048h,0ACh,010h,074h,0D8h,03Ch,0A0h,004h,068h,0CCh,030h,094h,0F8h,05Ch
	db	0C0h,024h,088h,0ECh,050h,0B4h,018h,07Ch,0E0h,044h,0A8h,00Ch,070h,0D4h,038h,09Ch

mt100_hi:
	db	000h,000h,000h,001h,001h,001h,002h,002h,003h,003h,003h,004h,004h,005h,005h,005h
	db	006h,006h,007h,007h,007h,008h,008h,008h,009h,009h,00Ah,00Ah,00Ah,00Bh,00Bh,00Ch
	db	00Ch,00Ch,00Dh,00Dh,00Eh,00Eh,00Eh,00Fh,00Fh,010h,010h,010h,011h,011h,011h,012h
	db	012h,013h,013h,013h,014h,014h,015h,015h,015h,016h,016h,017h,017h,017h,018h,018h
	db	019h,019h,019h,01Ah,01Ah,01Ah,01Bh,01Bh,01Ch,01Ch,01Ch,01Dh,01Dh,01Eh,01Eh,01Eh
	db	01Fh,01Fh,020h,020h,020h,021h,021h,021h,022h,022h,023h,023h,023h,024h,024h,025h
	db	025h,025h,026h,026h,027h,027h,027h,028h,028h,029h,029h,029h,02Ah,02Ah,02Ah,02Bh
	db	02Bh,02Ch,02Ch,02Ch,02Dh,02Dh,02Eh,02Eh,02Eh,02Fh,02Fh,030h,030h,030h,031h,031h
	db	032h,032h,032h,033h,033h,033h,034h,034h,035h,035h,035h,036h,036h,037h,037h,037h
	db	038h,038h,039h,039h,039h,03Ah,03Ah,03Ah,03Bh,03Bh,03Ch,03Ch,03Ch,03Dh,03Dh,03Eh
	db	03Eh,03Eh,03Fh,03Fh,040h,040h,040h,041h,041h,042h,042h,042h,043h,043h,043h,044h
	db	044h,045h,045h,045h,046h,046h,047h,047h,047h,048h,048h,049h,049h,049h,04Ah,04Ah
	db	04Bh,04Bh,04Bh,04Ch,04Ch,04Ch,04Dh,04Dh,04Eh,04Eh,04Eh,04Fh,04Fh,050h,050h,050h
	db	051h,051h,052h,052h,052h,053h,053h,053h,054h,054h,055h,055h,055h,056h,056h,057h
	db	057h,057h,058h,058h,059h,059h,059h,05Ah,05Ah,05Bh,05Bh,05Bh,05Ch,05Ch,05Ch,05Dh
	db	05Dh,05Eh,05Eh,05Eh,05Fh,05Fh,060h,060h,060h,061h,061h,062h,062h,062h,063h,063h

end:

;; Memory allocation for mp arrays:

	org	$+(4-$&3)	; force arrays on /4 boundary
mp_a:	org	$+bin_len
	org	$+(4-$&3)	; force arrays on /4 boundary
mp_b:	org	$+bin_len
	org	$+(4-$&3)	; force arrays on /4 boundary
mp_x:	org	$+bin_len
	org	$+(4-$&3)	; force arrays on /4 boundary
mp_y:	org	$+bin_len

