nx_linalg_affine_domain_candidate_t342.nx source
↩ module page · 140 lines · 8780 B
1// Private additive linear-algebra extension for scale-preserving affine skinning.
2// Uses existing Q30 arithmetic; no new persistent asset format.
3import "nx_linalg.nx"
4const LA_AFFINE_RANGE:i64=0-60
5const LA_AFFINE_SINGULAR:i64=0-61
6const LA_AFFINE_RESIDUAL:i64=0-62
7const LA_AFFINE_MAX:i64=70368744177664
8const LA_AFFINE_I64_MAX:i64=9223372036854775807
9func la_affine_domain(v:i64)->i64{if v<=0-LA_AFFINE_MAX||v>=LA_AFFINE_MAX{return LA_AFFINE_RANGE};return 0}
10// Row-major 3x3 Q30. Caller output stays unchanged on domain/singularity/residual refusal.
11func la_inverse3_q30(a:*Mat,out:*Mat,residual_budget:i64,report:*i64)->i64{
12 if a.rows!=3||a.cols!=3||out.rows!=3||out.cols!=3||residual_budget<0{return LA_AFFINE_RANGE}
13 let cofactor:*i64=sys_mmap(72) as *i64;let inverse:*i64=sys_mmap(72) as *i64
14 var rc:i64=0;var r:i64=0
15 while r<3&&rc==0{var c:i64=0;while c<3&&rc==0{if la_affine_domain(nx_mat_get(a,r,c))!=0{rc=LA_AFFINE_RANGE};c=c+1};r=r+1}
16 r=0;while r<3&&rc==0{var c:i64=0;while c<3&&rc==0{
17 let x:i64=fq_mul(nx_mat_get(a,(r+1)%3,(c+1)%3),nx_mat_get(a,(r+2)%3,(c+2)%3))
18 let y:i64=fq_mul(nx_mat_get(a,(r+1)%3,(c+2)%3),nx_mat_get(a,(r+2)%3,(c+1)%3))
19 // Each product is below 2^62 by the shared primitive's admitted domain.
20 let v:i64=x-y;if la_affine_domain(v)!=0{rc=LA_AFFINE_RANGE}else{cofactor[r*3+c]=v};c=c+1
21 };r=r+1}
22 var determinant:i64=0;var c:i64=0
23 while c<3&&rc==0{let term:i64=fq_mul(nx_mat_get(a,0,c),cofactor[c]);if term>0&&determinant>LA_AFFINE_I64_MAX-term{rc=LA_AFFINE_RANGE}else{if term<0&&determinant<(0-LA_AFFINE_I64_MAX)-term{rc=LA_AFFINE_RANGE}else{determinant=determinant+term}};c=c+1}
24 if determinant==0&&rc==0{rc=LA_AFFINE_SINGULAR}
25 r=0;while r<3&&rc==0{c=0;while c<3&&rc==0{
26 let numerator:i64=cofactor[c*3+r]
27 if fq_abs(numerator)/fq_abs(determinant)>=LA_AFFINE_MAX/FQ_ONE{rc=LA_AFFINE_RANGE}else{let v:i64=fq_div(numerator,determinant);if la_affine_domain(v)!=0{rc=LA_AFFINE_RANGE}else{inverse[r*3+c]=v}};c=c+1
28 };r=r+1}
29 var maximum:i64=0;r=0
30 while r<3&&rc==0{c=0;while c<3&&rc==0{var value:i64=0;var k:i64=0
31 while k<3{let term:i64=fq_mul(nx_mat_get(a,r,k),inverse[k*3+c]);if term>0&&value>LA_AFFINE_I64_MAX-term{rc=LA_AFFINE_RANGE}else{if term<0&&value<(0-LA_AFFINE_I64_MAX)-term{rc=LA_AFFINE_RANGE}else{value=value+term}};k=k+1}
32 var target:i64=0;if r==c{target=FQ_ONE};if value<0-LA_AFFINE_MAX||value>LA_AFFINE_MAX{rc=LA_AFFINE_RANGE}else{let difference:i64=fq_abs(value-target);if difference>maximum{maximum=difference}};c=c+1
33 };r=r+1}
34 report[0]=determinant;report[1]=maximum
35 if rc==0&&maximum>residual_budget{rc=LA_AFFINE_RESIDUAL}
36 if rc==0{r=0;while r<3{c=0;while c<3{nx_mat_set(out,r,c,inverse[r*3+c]);c=c+1};r=r+1}}
37 sys_munmap(cofactor as *u8,72);sys_munmap(inverse as *u8,72);return rc
38}
39
40// Column-major 4x4 Q30; domain follows the existing exact fixed-point multiply owner.
41func la_affine_compose_q30(a:*i64,b:*i64,out:*i64)->i64{
42 var i:i64=0;while i<16{if la_affine_domain(a[i])!=0||la_affine_domain(b[i])!=0{return LA_AFFINE_RANGE};i=i+1}
43 let temporary:*i64=sys_mmap(128) as *i64;var rc:i64=0;var c:i64=0
44 while c<4&&rc==0{var r:i64=0;while r<4&&rc==0{var sum:i64=0;var k:i64=0
45 while k<4&&rc==0{let term:i64=fq_mul(a[k*4+r],b[c*4+k]);if term>0&&sum>LA_AFFINE_I64_MAX-term{rc=LA_AFFINE_RANGE}else{if term<0&&sum<(0-LA_AFFINE_I64_MAX)-term{rc=LA_AFFINE_RANGE}else{sum=sum+term}};k=k+1}
46 if la_affine_domain(sum)!=0{rc=LA_AFFINE_RANGE};temporary[c*4+r]=sum;r=r+1
47 };c=c+1}
48 if rc==0{i=0;while i<16{out[i]=temporary[i];i=i+1}}
49 sys_munmap(temporary as *u8,128);return rc
50}
51
52// Static yaw draw enclosure, binary32 input ABI: bounds loXYZ/hiXYZ, sin, cos,
53// base scale, instance scale, translation XYZ. Deformation is a caller admission
54// prerequisite. The result is a world-axis Q8 AABB, not anatomical geometry.
55// Q30 interval endpoints are outward rounded; source / 2^10 and base * 2^10
56// preserve the exact real expression while staying in fq_mul's proved domain.
57const LA_STATIC_INPUT_WORDS:i64=13
58const LA_STATIC_OUTPUT_WORDS:i64=6
59const LA_STATIC_SCRATCH_WORDS:i64=48
60const LA_BINARY32_FRACTION_BITS:i64=23
61const LA_BINARY32_ULP_DEN:i64=8388608
62const LA_STATIC_SOURCE_SHIFT:i64=10
63const LA_STATIC_PATH_ROUNDS:i64=5
64const LA_STATIC_GAMMA_FACTOR:i64=2*LA_STATIC_PATH_ROUNDS
65const LA_STATIC_Q8_SHIFT:i64=FQ_SHIFT-8
66func la_iv(p:*i64,index:i64)->*i64{return ((p as i64)+index*16) as *i64}
67func la_iv_abs(p:*i64)->i64{let a:i64=fq_abs(p[0]);let b:i64=fq_abs(p[1]);if a>b{return a};return b}
68func la_f32_ordered(a:i64,b:i64)->i64{
69 if (a&2147483647)==0&&(b&2147483647)==0{return 1}
70 var x:i64=2147483648+a;var y:i64=2147483648+b
71 if (a&2147483648)!=0{x=4294967295-a};if (b&2147483648)!=0{y=4294967295-b}
72 return (x<=y) as i64
73}
74func la_iv_f32(bits:i64,power:i64,out:*i64)->i64{
75 if bits<0||bits>4294967295{return LA_AFFINE_RANGE}
76 let exp:i64=(bits>>23)&255;var mant:i64=bits&8388607
77 if exp==255{return LA_AFFINE_RANGE}
78 var shift:i64=exp-120+power
79 if exp==0{shift=0-119+power}else{mant=mant+8388608}
80 var lo:i64=0;var hi:i64=0
81 if mant!=0{
82 if shift>=0{if shift>=63{return LA_AFFINE_RANGE};if mant>((LA_AFFINE_MAX-1)>>shift){return LA_AFFINE_RANGE};lo=mant<<shift;hi=lo}
83 else{let n:i64=0-shift;if n>=63{hi=1}else{lo=mant>>n;hi=lo;if (lo<<n)!=mant{hi=hi+1}}}
84 }
85 if hi>=LA_AFFINE_MAX{return LA_AFFINE_RANGE}
86 if (bits&2147483648)!=0{out[0]=0-hi;out[1]=0-lo}else{out[0]=lo;out[1]=hi}
87 return 0
88}
89func la_iv_mul(a:*i64,b:*i64,out:*i64)->i64{
90 var low:i64=LA_AFFINE_I64_MAX;var high:i64=0-LA_AFFINE_I64_MAX
91 var i:i64=0;while i<2{var j:i64=0;while j<2{
92 if la_affine_domain(a[i])!=0||la_affine_domain(b[j])!=0{return LA_AFFINE_RANGE}
93 let value:i64=fq_mul(a[i],b[j]);if value<low{low=value};if value>high{high=value};j=j+1
94 };i=i+1}
95 // fq_mul truncates toward zero. One Q30 unit on either side encloses its remainder.
96 if low<=0-LA_AFFINE_MAX+1||high>=LA_AFFINE_MAX-1{return LA_AFFINE_RANGE}
97 out[0]=low-1;out[1]=high+1;return 0
98}
99func la_static_storage_disjoint(a:*i64,aw:i64,b:*i64,bw:i64)->i64{
100 let x:i64=a as i64;let y:i64=b as i64
101 if x<=0||y<=0{return 0}
102 if x<=y{if y-x<aw*8{return 0}}else{if x-y<bw*8{return 0}}
103 return 1
104}
105func la_static_yaw_bounds_f32(input:*i64,input_words:i64,out:*i64,out_words:i64,scratch:*i64,scratch_words:i64)->i64{
106 if input_words!=LA_STATIC_INPUT_WORDS||out_words!=LA_STATIC_OUTPUT_WORDS||scratch_words<LA_STATIC_SCRATCH_WORDS{return LA_AFFINE_RANGE}
107 if la_static_storage_disjoint(input,input_words,out,out_words)==0||la_static_storage_disjoint(input,input_words,scratch,LA_STATIC_SCRATCH_WORDS)==0||la_static_storage_disjoint(out,out_words,scratch,LA_STATIC_SCRATCH_WORDS)==0{return LA_AFFINE_RANGE}
108 var i:i64=0
109 while i<LA_STATIC_INPUT_WORDS{var power:i64=0;if i<6{power=0-LA_STATIC_SOURCE_SHIFT};if i==8{power=LA_STATIC_SOURCE_SHIFT};let rc:i64=la_iv_f32(input[i],power,la_iv(scratch,i));if rc!=0{return rc};i=i+1}
110 i=0;while i<3{if la_f32_ordered(input[i],input[i+3])==0{return LA_AFFINE_RANGE};let lo:*i64=la_iv(scratch,i);let hi:*i64=la_iv(scratch,i+3);if lo[0]>hi[1]{return LA_AFFINE_RANGE};lo[1]=hi[1];i=i+1}
111 let base:*i64=la_iv(scratch,8);let scale:*i64=la_iv(scratch,9)
112 if base[0]<=0||scale[0]<=0{return LA_AFFINE_RANGE}
113 let combined:*i64=la_iv(scratch,13);if la_iv_mul(base,scale,combined)!=0{return LA_AFFINE_RANGE}
114 var axis:i64=0
115 while axis<3{
116 let first:*i64=la_iv(scratch,14);let second:*i64=la_iv(scratch,15)
117 let term1:*i64=la_iv(scratch,16);let term2:*i64=la_iv(scratch,17)
118 if axis==1{if la_iv_mul(la_iv(scratch,2),combined,term1)!=0{return LA_AFFINE_RANGE};term2[0]=0;term2[1]=0}
119 else{
120 var a:i64=7;var b:i64=6;if axis==2{a=6;b=7}
121 if la_iv_mul(la_iv(scratch,0),la_iv(scratch,a),first)!=0||la_iv_mul(la_iv(scratch,1),la_iv(scratch,b),second)!=0{return LA_AFFINE_RANGE}
122 if la_iv_mul(first,combined,term1)!=0||la_iv_mul(second,combined,term2)!=0{return LA_AFFINE_RANGE}
123 if axis==2{let old:i64=term2[0];term2[0]=0-term2[1];term2[1]=0-old}
124 }
125 let translation:*i64=la_iv(scratch,10+axis)
126 let mass:i64=la_iv_abs(term1)+la_iv_abs(term2)+la_iv_abs(translation)
127 // Each monomial reaches the output through at most 3 multiplies and 2 adds.
128 // With u=2^-23 (either adjacent f32), gamma5=5u/(1-5u)<10u.
129 // Absolute monomial mass also covers cancellation and permitted reassociation.
130 let error:i64=(mass*LA_STATIC_GAMMA_FACTOR+LA_BINARY32_ULP_DEN-1)/LA_BINARY32_ULP_DEN+1
131 let low:i64=term1[0]+term2[0]+translation[0]-error
132 let high:i64=term1[1]+term2[1]+translation[1]+error
133 if la_affine_domain(low)!=0||la_affine_domain(high)!=0{return LA_AFFINE_RANGE}
134 // Arithmetic shifts give floor; negated floor of negative gives ceil.
135 scratch[40+axis]=low>>LA_STATIC_Q8_SHIFT
136 scratch[43+axis]=0-((0-high)>>LA_STATIC_Q8_SHIFT)
137 axis=axis+1
138 }
139 i=0;while i<LA_STATIC_OUTPUT_WORDS{out[i]=scratch[40+i];i=i+1};return 0
140}