// ================================
// Nihonium
//		SFPU: sqrt function test
//		(c) 2023	1YEN Toru
// ================================
asm		"mcoc_irom.mem"
incl	"mcoc115.incl"
# ================================
equ		chr_lf,0x0a
equ		chr_spc,0x20
equ		reso,127
# ================================


// initialize
sesr	sreg_b_nh
ldc		sp,ramtop+0x0100
ldbiu	r0,simc_noto
ldwi	r7,simctrl
stw		[r7],r0

// check SFPU
ldbiu	r0,0
fneg	r0,r0
ldli	r1,sngl_zer_n
cmp		r0,r1
beq		chk_ok
// error: no sfpu
ldwi	r1,lab_chk_str_no_sfpu
jalwr7	print
bra		tstfail

chk_str_no_sfpu:
dats	"ERR: no SFPU found\r\n\0"
chk_ok:


// ================================
// sqrt test
// ================================
ldwi	r1,lab_func_name
jalwr7	print
ldwi	r7,princtl
ldbiu	r0,prin_flt_3|prin_nld
stw		[r7],r0							// not print NL

// sqrt (x)
ldwi	r6,lab_test_tab
ldbil	r7,prinflt&0xff
#
# print hexadecimal, for debug
#@ldbil	r7,prinhxl&0xff
#
mov		r5,r7
ldbil	r5,princhr&0xff
test_loop:
ldl		r1,[r6]
addi	r6,4
// break condition
cmpi	r1,1
beq		tstpass
cmpi	r1,2
beq		test_break
// test
pushl	r6
stl		[r7],r1
ldbiu	r0,chr_spc
stw		[r5],r0
jalwr6	math_sqrt
stl		[r7],r0
ldbiu	r1,chr_spc
stw		[r5],r1
popl	r6
ldl		r2,[r6]
addi	r6,4
stl		[r7],r2
ldbiu	r1,chr_lf
stw		[r5],r1
// check
cmp		r0,r2
bne		tstfail
bra		test_loop
test_break:


// ================================
// calculate plot data
// ================================
// header
ldwi	r1,lab_test_head
jalwr7	print
// print long hexadecimal number format
ldbiu	r0,prin_nld
ldwi	r1,princtl
stw		[r1],r0
ldwi	r7,prinhxl

// 0<=x<reso
// plot range
fldi	r4,10
ldwi	r5,reso
fsngl	r1,r5
fdiv	r4,r1
ldbiu	r5,0
// loop
test_pl_lp:
fsngl	r1,r5
fmul	r1,r4
// print r1+" "
ldbil	r7,prinhxl&0xff
stl		[r7],r1
ldbil	r7,princhr&0xff
ldbiu	r0,chr_spc
stw		[r7],r0
// sqrt (x)
pushlm	r4,r5
jalwr6	math_sqrt
// print r0+LF
ldbil	r7,prinhxl&0xff
stl		[r7],r0
ldbil	r7,princhr&0xff
ldbiu	r0,chr_lf
stw		[r7],r0
// end of loop
poplm	r5,r4
addi	r5,1
ldwi	r6,reso
cmp		r6,r5
bhs		test_pl_lp
// end
ldwi	r1,lab_test_end
jalwr7	print

// end of test
bra		tstpass


// ================================
// print r1 ; disturbed r0,r1,r7
// ================================
print:
ldwi	r7,princhr
prn_loop:
ldb		r0,[r1]
mov		r0,r0
beq		prn_break
stw		[r7],r0
addi	r1,1
bra		prn_loop
prn_break:
rtnw

// strings
func_name:
dats	"sqrt\n\0"

test_head:
dats	"PLOT\n\0"

test_end:
dats	"END\n\0"


// test table
test_tab:
#datl	0x2,0x0							// break mark
#
datl	sngl_zer,sngl_zer
datl	sngl_zer_n,sngl_zer
datl	sngl_nan,sngl_nan
datl	sngl_inf,sngl_inf
datl	sngl_inf_n,sngl_nan
#
datf	-2
datl	sngl_nan
#
datl	sngl_min
datf	sqrt(1.1754944e-38)
#
datl	sngl_max
datf	sqrt(3.4028235e+38)
#
datf	0.25,sqrt(0.25)
datf	0.5,sqrt(0.5)
datf	0.75,sqrt(0.75)
datf	1,sqrt(1)
datf	2,sqrt(2)
datf	3,sqrt(3)
datf	4,sqrt(4)
datf	5,sqrt(5)
datf	7,sqrt(7)
datf	9,sqrt(9)
datf	11,sqrt(11)
datf	16,sqrt(16)
datf	17,sqrt(17)
datf	25,sqrt(25)
datf	36,sqrt(36)
datf	49,sqrt(49)
datf	64,sqrt(64)
datf	81,sqrt(81)
datf	90,sqrt(90)
datf	10,sqrt(10)
datf	200,sqrt(200)
datf	3_000,sqrt(3_000)
datf	40_000,sqrt(40_000)
datf	500_000,sqrt(500_000)
datf	6_000_000,sqrt(6_000_000)
datf	70_000_000,sqrt(70_000_000)
datf	800_000_000,sqrt(800_000_000)
datf	9_000_000_000,sqrt(8.9999995e+09)
#datl	0x1,0x0							// end mark
datl	0x2,0x0							// break mark


// ================================
// float r0=sqrt (float r1), disturbed r1~r4
// ================================
math_sqrt:
// zero or non-finite value
ldbiu	r0,0
fcmp	r1,r0
beq		math_sqrt_return
ldli	r0,sngl_nan
bmi		math_sqrt_return
ldli	r0,sngl_inf
fcmp	r1,r0
bcs		math_sqrt_return
// sqrt exponent
mov		r2,r1
lsfti	r2,-23
subi	r2,127
ldli	r0,~0x7f80_0000
and		r1,r0
// if ex&0x01==0
mov		r0,r2
andi	r0,0x01
bne		math_sqrt_ex_odd
asfti	r2,-1
ldli	r0,0x3f80_0000
or		r1,r0
bra		math_sqrt_ex_endif
align	4
math_sqrt_ex_odd:
// if ex&0x01==1
subi	r2,1
ldli	r0,0x4000_0000
asfti	r2,-1
or		r1,r0
math_sqrt_ex_endif:
asfti	r2,23
// r0=y<n>
// r1=x
// r2=exponent
// Newton-Raphson method
//	sqrt (x)=y3, y<n+1>=(y<n> + x/y<n>)/2, y0=0.3259*x + 0.7407
ldli	r4,0x0080_0000					// lsb of exponent
// y0
fldi	r0,0.3259
fmul	r0,r1
fldi	r3,0.7407
fadd	r0,r3
// iteration 3 times: y1~y3
rept	3
// y<$(repti)+1>
mov		r3,r1
fdiv	r3,r0
fadd	r0,r3
sub		r0,r4
endr
// calculate exponent
add		r0,r2
math_sqrt_return:
rtnw


// ================================
tstpass:
ldwi	r0,simpass						// pass: LED=blue
bra		tstfnsh
tstfail:
ldwi	r0,simfail						// fail: LED=red
tstfnsh:
ldwi	r7,simfnsh
stw		[r7],r0
bra		pcnt-2
nop
// ================================

