code wiki / _hdl_build / nx_bd_calc.nx
nx_bd_calc.nx source
↩ module page · 130 lines · 6554 B
1// nx_bd_calc.nx -- SOVEREIGN Bjontegaard BD-rate (2026-07-29; retires the bd_*.py lane instruments
2// per the sovereign-no-python law -- python remains ONLY as the cubic ORACLE inside this organ's gate).
3// ESTIMATOR (declared): BD-PWL -- piecewise-LINEAR interpolation of log10(rate) over the overlapping
4// PSNR span, trapezoid-integrated. The classic estimator fits a cubic; on the 4-5 point curves this
5// lane sweeps, PWL is monotone, has no ill-conditioned Vandermonde solve, and its delta vs the cubic
6// oracle is MEASURED and bounded in nx_bd_calc_gate (tolerance derived from measurement, never taste).
7// All math is bits-up sovereign binary64 (nx_f64 add/sub/mul/div/log2/exp/cvt).
8// usage: nx_bd_calc <anchor.pts> <test.pts>
9// each line: "<kbps_milli> <psnr_milli>" (integers; kbps_milli = bytes*5 for 48f@30fps CIF sweeps;
10// the milli scaling cancels: log10(1000x)=log10(x)+3 shifts BOTH curves equally, psnr-milli scales
11// numerator and denominator of the average equally)
12// output: "BD_PWL_PERMILLE=<n>" (test-vs-anchor bits at equal quality; negative = test wins)
13// license: ORIGINAL
14import "nx_syscalls.nx"
15import "nx_tier.nx"
16import "nx_f64.nx"
17import "nx_f64_div.nx"
18import "nx_f64_cvt.nx"
19import "nx_f64_exp.nx"
20import "nx_f64_log2.nx"
21
22const BD_MAXPTS: i64 = 64
23const BD_LOG10_2: i64 = 0x3FD34413509F79FF // log10(2) in binary64
24const BD_LN10: i64 = 0x40026BB1BBB55516 // ln(10) in binary64
25const BD_ONE: i64 = 0x3FF0000000000000
26const BD_HALF: i64 = 0x3FE0000000000000
27const BD_KILO: i64 = 0x408F400000000000 // 1000.0
28
29func gw(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 }
30func gn(v: i64) -> i64 {
31 let b: *u8=sys_mmap(28); var m: i64=v; if m<0{sys_write(1,"-" as *u8,1);m=0-m}
32 let t: *u8=sys_mmap(28); var k: i64=0; if m==0{t[0]=48 as u8;k=1} while m>0{t[k]=(48+(m%10)) as u8;m=m/10;k=k+1}
33 var i: i64=0; while i<k{b[i]=t[k-1-i];i=i+1} sys_write(1,b,k); return 0 }
34
35// parse "<int> <int>\n" rows; returns count. km/pm caller-allocated BD_MAXPTS i64.
36func bd_parse(path: *u8, km: *i64, pm: *i64) -> i64 {
37 let box: *i64 = sys_mmap(16) as *i64
38 let d: *u8 = sys_read_file(path, box)
39 if (d as i64) == 0 { return 0 - 1 }
40 let n: i64 = box[0]
41 var cnt: i64 = 0
42 var i: i64 = 0
43 while i < n {
44 // skip blank/comment lines
45 if d[i] == (10 as u8) { i = i + 1; continue }
46 if d[i] == (13 as u8) { i = i + 1; continue }
47 if d[i] == (35 as u8) { while i < n { if d[i] == (10 as u8) { break } i = i + 1 } continue }
48 var v1: i64 = 0
49 while i < n { let c: i64 = d[i] & 0xff
50 if c < 48 { break } if c > 57 { break }
51 v1 = v1*10 + (c-48); i = i + 1 }
52 while i < n { if d[i] != (32 as u8) { break } i = i + 1 }
53 var v2: i64 = 0
54 while i < n { let c2: i64 = d[i] & 0xff
55 if c2 < 48 { break } if c2 > 57 { break }
56 v2 = v2*10 + (c2-48); i = i + 1 }
57 while i < n { if d[i] == (10 as u8) { i = i + 1; break } i = i + 1 }
58 if v1 > 0 { if v2 > 0 { if cnt < BD_MAXPTS { km[cnt] = v1; pm[cnt] = v2; cnt = cnt + 1 } } }
59 }
60 return cnt }
61
62// insertion sort both arrays by ascending pm
63func bd_sort(km: *i64, pm: *i64, n: i64) -> i64 {
64 var i: i64 = 1
65 while i < n {
66 let kp: i64 = km[i]; let pp: i64 = pm[i]
67 var j: i64 = i - 1
68 while j >= 0 { if pm[j] > pp { km[j+1] = km[j]; pm[j+1] = pm[j]; j = j - 1 } else { break } }
69 km[j+1] = kp; pm[j+1] = pp
70 i = i + 1 }
71 return 0 }
72
73// v[i] = log10(km[i]) as f64 raw (milli offset cancels across curves)
74func bd_logs(km: *i64, v: *i64, n: i64) -> i64 {
75 var i: i64 = 0
76 while i < n { v[i] = nx_f64_mul(nx_f64_log2(nx_i64_to_f64(km[i])), BD_LOG10_2); i = i + 1 }
77 return 0 }
78
79// trapezoid-integrate the piecewise-linear v(p) over [lo,hi] (i64 psnr-milli bounds).
80// returns f64 raw integral (units: log10 * psnr_milli).
81func bd_integ(pm: *i64, v: *i64, n: i64, lo: i64, hi: i64) -> i64 {
82 var acc: i64 = 0 // f64 +0.0
83 var i: i64 = 0
84 while i < n - 1 {
85 let p0: i64 = pm[i]; let p1: i64 = pm[i+1]
86 if p1 > p0 {
87 var a: i64 = p0; if a < lo { a = lo }
88 var b: i64 = p1; if b > hi { b = hi }
89 if b > a {
90 let seg: i64 = nx_i64_to_f64(p1 - p0)
91 let dv: i64 = nx_f64_sub(v[i+1], v[i])
92 let va: i64 = nx_f64_add(v[i], nx_f64_mul(dv, nx_f64_div(nx_i64_to_f64(a - p0), seg)))
93 let vb: i64 = nx_f64_add(v[i], nx_f64_mul(dv, nx_f64_div(nx_i64_to_f64(b - p0), seg)))
94 acc = nx_f64_add(acc, nx_f64_mul(nx_f64_mul(nx_f64_add(va, vb), BD_HALF), nx_i64_to_f64(b - a)))
95 }
96 }
97 i = i + 1 }
98 return acc }
99
100func main(argc: i64, argv: *i64) -> i64 {
101 if argc < 3 { gw("usage: nx_bd_calc <anchor.pts> <test.pts> (lines: kbps_milli psnr_milli)\n" as *u8); return 2 }
102 let ka: *i64 = sys_mmap(BD_MAXPTS*8) as *i64
103 let pa: *i64 = sys_mmap(BD_MAXPTS*8) as *i64
104 let kt: *i64 = sys_mmap(BD_MAXPTS*8) as *i64
105 let pt: *i64 = sys_mmap(BD_MAXPTS*8) as *i64
106 let na: i64 = bd_parse(argv[1] as *u8, ka, pa)
107 let nt: i64 = bd_parse(argv[2] as *u8, kt, pt)
108 if na < 2 { gw("BD-RED anchor needs >=2 points\n" as *u8); return 3 }
109 if nt < 2 { gw("BD-RED test needs >=2 points\n" as *u8); return 3 }
110 bd_sort(ka, pa, na)
111 bd_sort(kt, pt, nt)
112 var lo: i64 = pa[0]; if pt[0] > lo { lo = pt[0] }
113 var hi: i64 = pa[na-1]; if pt[nt-1] < hi { hi = pt[nt-1] }
114 if hi <= lo { gw("BD-RED no PSNR overlap\n" as *u8); return 4 }
115 let va: *i64 = sys_mmap(BD_MAXPTS*8) as *i64
116 let vt: *i64 = sys_mmap(BD_MAXPTS*8) as *i64
117 bd_logs(ka, va, na)
118 bd_logs(kt, vt, nt)
119 let ia: i64 = bd_integ(pa, va, na, lo, hi)
120 let it: i64 = bd_integ(pt, vt, nt, lo, hi)
121 let avg: i64 = nx_f64_div(nx_f64_sub(it, ia), nx_i64_to_f64(hi - lo))
122 // bd = 10^avg - 1 = exp(avg*ln10) - 1; permille with sign-aware half-up rounding
123 let bd: i64 = nx_f64_sub(nx_f64_exp(nx_f64_mul(avg, BD_LN10)), BD_ONE)
124 var scaled: i64 = nx_f64_mul(bd, BD_KILO)
125 if nx_f64_sign(scaled) == 1 { scaled = nx_f64_sub(scaled, BD_HALF) } else { scaled = nx_f64_add(scaled, BD_HALF) }
126 let pm_out: i64 = nx_f64_to_i64(scaled)
127 gw("BD_PWL_PERMILLE=" as *u8); gn(pm_out)
128 gw(" estimator=piecewise-log-linear overlap_milli=" as *u8); gn(lo); gw(".." as *u8); gn(hi)
129 gw(" pts=" as *u8); gn(na); gw("/" as *u8); gn(nt); gw("\n" as *u8)
130 return 0 }