code wiki / _hdl_build / nx_recip64_scalar_test.nx
nx_recip64_scalar_test.nx source
↩ module page · 141 lines · 5068 B
1// nx_recip64_scalar_test.nx -- ALGORITHM-FIRST proof of the true 64-bit unsigned
2// Newton-Raphson reciprocal divider (the scalar form, before gate-net emission).
3// Verified bit-exact vs the PROVEN unsigned oracle nx_rv64im_udiv over a hard
4// battery (dense small + 100k random FULL 64-bit incl. bit-63-set + every power
5// of two + edges). Once this is 1:1, the gate-net mirrors it.
6//
7// Method (normalized reciprocal, mulh-only -- no 128-bit datapath needed):
8// s = 63 - msb(D); Dn = D<<s in [2^63,2^64) (dn = Dn/2^64 in [0.5,1))
9// X = 2^63 (seed x=1); repeat: dx=mulhu(Dn,X); t=-dx; Xm=mulhu(X,t);
10// X = (Xm < 2^63) ? Xm<<1 : ALLONES (cap avoids the x->2 overflow)
11// q = mulhu(N,X) >>logical msb(D); then overflow-aware +/- correction.
12// Seed x=1 has error <= 1-dn <= 0.5; quadratic Newton -> 6 iters reach 2^-64.
13//
14// Known answer: "<ok> <total>" with ok==total, exit 0.
15
16import "rv64im_min_alu.nx" // nx_rv64im_mulh_unsigned / udiv / ltu / lshr
17
18const NR64_ITERS: i64 = 7
19const NR64_KCORR: i64 = 6
20const NR64_HALF: i64 = 0 - 9223372036854775808 // 2^63 (as the bit pattern)
21
22func nr64_msb(D: i64) -> i64 { // highest set bit of D>=1 (0..63)
23 var p: i64 = 63
24 while nx_rv64im_lshr(D, p) == 0 { p = p - 1 }
25 return p
26}
27
28func nr64_div(N: i64, D: i64) -> i64 {
29 let p: i64 = nr64_msb(D)
30 let s: i64 = 63 - p
31 let Dn: i64 = D << s
32 var X: i64 = NR64_HALF // 2^63 (x = 1.0 in Q63)
33 var it: i64 = 0
34 while it < NR64_ITERS {
35 let dx: i64 = nx_rv64im_mulh_unsigned(Dn, X)
36 let t: i64 = 0 - dx // 2^64 - dx (= 2 - dn*x in Q63)
37 let Xm: i64 = nx_rv64im_mulh_unsigned(X, t)
38 if nx_rv64im_ltu(Xm, NR64_HALF) == 1 { X = Xm << 1 } else { X = 0 - 1 }
39 it = it + 1
40 }
41 var q: i64 = nx_rv64im_lshr(nx_rv64im_mulh_unsigned(N, X), p)
42
43 // INTERLEAVED overflow-aware correction: each round, if q*D > N step down,
44 // else if (q+1)*D <= N step up, else stable. (mulhu detects the >2^64 case.)
45 var kc: i64 = 0
46 while kc < NR64_KCORR {
47 let phi: i64 = nx_rv64im_mulh_unsigned(q, D)
48 var toobig: i64 = 0
49 if phi != 0 { toobig = 1 }
50 if phi == 0 { if nx_rv64im_ltu(N, q * D) == 1 { toobig = 1 } }
51 if toobig == 1 {
52 q = q - 1
53 } else {
54 let qp1: i64 = q + 1
55 if qp1 != 0 { // guard: q already max (q+1 wraps to 0)
56 let phu: i64 = nx_rv64im_mulh_unsigned(qp1, D)
57 if phu == 0 { if nx_rv64im_ltu(N, qp1 * D) == 0 { q = qp1 } }
58 }
59 }
60 kc = kc + 1
61 }
62 return q
63}
64
65func _emit_dec(v: i64) -> i64 {
66 let b: *u8 = sys_mmap(28); var n: i64 = v; if n < 0 { n = 0 - n }
67 let t: *u8 = sys_mmap(28); var k: i64 = 0
68 if n == 0 { t[0] = 48; k = 1 }
69 while n > 0 { t[k] = 48 + (n % 10); n = n / 10; k = k + 1 }
70 var i: i64 = 0; while i < k { b[i] = t[k - 1 - i]; i = i + 1 }
71 b[k] = 32; sys_write(1, b, k + 1); return 0
72}
73func _nl() -> i64 { let z: *u8 = sys_mmap(2); z[0] = 10; sys_write(1, z, 1); return 0 }
74
75func main() -> i64 {
76 var total: i64 = 0
77 var ok: i64 = 0
78
79 // dense small: D in [1,256], N in [0,4095]
80 var D: i64 = 1
81 while D <= 256 {
82 var N: i64 = 0
83 while N <= 4095 {
84 total = total + 1
85 if nr64_div(N, D) == nx_rv64im_udiv(N, D) { ok = ok + 1 }
86 N = N + 1
87 }
88 D = D + 1
89 }
90
91 // every power of two as the divisor, with assorted dividends
92 var kbit: i64 = 0
93 while kbit < 64 {
94 let pd: i64 = 1 << kbit
95 let nv: *i64 = sys_mmap(8 * 5) as *i64
96 nv[0]=0; nv[1]=1; nv[2]=pd; nv[3]=0-1; nv[4]=(0-1)-pd
97 var j: i64 = 0
98 while j < 5 {
99 total = total + 1
100 if nr64_div(nv[j], pd) == nx_rv64im_udiv(nv[j], pd) { ok = ok + 1 }
101 j = j + 1
102 }
103 kbit = kbit + 1
104 }
105
106 // random FULL 64-bit (N,D incl. bit 63 set)
107 var seed: i64 = 88172645463325252
108 var nr: i64 = 0
109 while nr < 100000 {
110 seed = seed ^ (seed << 13); seed = seed ^ nx_rv64im_lshr(seed, 7); seed = seed ^ (seed << 17)
111 let N: i64 = seed
112 seed = seed ^ (seed << 13); seed = seed ^ nx_rv64im_lshr(seed, 7); seed = seed ^ (seed << 17)
113 var D: i64 = seed
114 if D == 0 { D = 1 }
115 total = total + 1
116 if nr64_div(N, D) == nx_rv64im_udiv(N, D) { ok = ok + 1 }
117 nr = nr + 1
118 }
119
120 // edges
121 let en: *i64 = sys_mmap(8 * 8) as *i64
122 let ed: *i64 = sys_mmap(8 * 8) as *i64
123 en[0]=0-1; ed[0]=1
124 en[1]=0-1; ed[1]=0-1
125 en[2]=0-1; ed[2]=2
126 en[3]=0-1; ed[3]=NR64_HALF
127 en[4]=NR64_HALF; ed[4]=3
128 en[5]=0-1; ed[5]=(0-1)
129 en[6]=12345678901234; ed[6]=99991
130 en[7]=0-1; ed[7]=NR64_HALF + 1
131 var e: i64 = 0
132 while e < 8 {
133 total = total + 1
134 if nr64_div(en[e], ed[e]) == nx_rv64im_udiv(en[e], ed[e]) { ok = ok + 1 }
135 e = e + 1
136 }
137
138 _emit_dec(ok); _emit_dec(total); _nl()
139 if ok != total { sys_exit(1); return 1 }
140 sys_exit(0); return 0
141}