float_finite_product_quotient.py (10577B)
1 #!/usr/bin/env python3 2 """Finite products and nonzero-divisor quotients, with exact rational oracles.""" 3 import argparse 4 import hashlib 5 import json 6 from pathlib import Path 7 import random 8 import statistics 9 10 from float_aligned_quantize import nearest 11 from float_quantize import ROOT, value, dyadic 12 from float_multiply import Runner, sha 13 from base_workloads import encode as nat 14 from http_headers import LEAF, TRUE, pair, tree_list, expect 15 16 from float_finite import SIGN, LIMIT, MASK, finite_value, encoded, special_words, normalized 17 18 19 def oracle(a,b,divide=False): 20 if a>=LIMIT or b>=LIMIT:return False,2 21 if (a & MASK)>>52==2047 or (b & MASK)>>52==2047:return False,1 22 if divide and not (b & MASK):return False,3 23 x=finite_value(a)/finite_value(b) if divide else finite_value(a)*finite_value(b) 24 bits,inexact,underflow=nearest(abs(x)) 25 return True,(bits | ((a^b)&SIGN),inexact,underflow) 26 27 28 def cases(): 29 words=special_words();rows=[(a,b) for a in words for b in words] 30 # Around half-minimum-subnormal, normal/subnormal and overflow boundaries. 31 for a in [1,2,3,2**52-1,2**52,2**52+1,0x7fefffffffffffff]: 32 for b in [0x3fdfffffffffffff,0x3fe0000000000000,0x3fe0000000000001, 33 0x3fefffffffffffff,0x3ff0000000000000,0x3ff0000000000001, 34 0x3fffffffffffffff,0x4000000000000000,0x4000000000000001]: 35 for sa in [0,SIGN]: 36 for sb in [0,SIGN]:rows.extend([(a|sa,b|sb),(b|sb,a|sa)]) 37 for a,b in double_round_pairs()+overflow_pairs(): 38 for sa in [0,SIGN]: 39 for sb in [0,SIGN]:rows.extend([(a|sa,b|sb),(b|sb,a|sa)]) 40 rng=random.Random(64054) 41 def word():return (rng.randrange(2047)<<52)|rng.getrandbits(52)|(rng.randrange(2)<<63) 42 rows += [(word(),word()) for _ in range(384)] 43 # Correlated exponents exercise finite normal and subnormal answers. 44 for _ in range(192): 45 ea=rng.randrange(1,2047) 46 eb=max(1,min(2046,rng.choice([1023-ea,2046-ea,ea-1023,ea])+rng.randrange(-3,4))) 47 rows.append(((ea<<52)|rng.getrandbits(52)|(rng.randrange(2)<<63), 48 (eb<<52)|rng.getrandbits(52)|(rng.randrange(2)<<63))) 49 return list(dict.fromkeys(rows)) 50 51 52 def raw_cases(): 53 lo,hi=2**52,2**53-1 54 values=[lo,lo+1,lo+2,3*2**51,hi-1,hi] 55 rng=random.Random(5653) 56 return [(a,b) for a in values for b in values]+[(rng.randrange(lo,hi+1),rng.randrange(lo,hi+1)) for _ in range(96)] 57 58 59 def overflow_pairs(): 60 # Exact overflow midpoint: (2^27-1)*(2^27+1) = 2^54-1. 61 a=(1534<<52)|(((2**27-1)<<26)-2**52) 62 b=(1535<<52)|(((2**27+1)<<25)-2**52) 63 return [(a+i,b) for i in [-1,0,1]] 64 65 66 def double_round_pairs(): 67 return [((512<<52)|1,(511<<52)|4),((512<<52)|1,(511<<52)|(2**52-3))] 68 69 70 def double_round_wrong(a,b): 71 _,sa,ma=normalized(a);_,sb,mb=normalized(b) 72 n=ma*mb;s=sa+sb-2304;k=max(n.bit_length()-53,0) 73 q,r=divmod(n,2**k) 74 q+=2*r>2**k or (2*r==2**k and q%2==1) 75 bits,ix,uf=nearest(dyadic(q,s+k-2304));ix=ix or bool(r) 76 return True,(bits|((a^b)&SIGN),ix,ix and bits<2**52) 77 78 79 def negative_controls(): 80 result=[('doubleProducts',a,b,False,double_round_wrong(a,b)) for a,b in double_round_pairs()] 81 for a,b in [(0x3ff0000000000000,0x3fffffffffffffff),(0x3ff0000000000000,0x3ff0000000000001)]: 82 _,sa,ma=normalized(a);_,sb,mb=normalized(b) 83 wrong=nearest(dyadic((ma<<56)//mb,sa-sb-56)) 84 result.append(('truncateDivides',a,b,True,(True,wrong))) 85 return result 86 87 88 def main(): 89 p=argparse.ArgumentParser(description=__doc__) 90 p.add_argument('tricu',type=Path);p.add_argument('--directory',required=True,type=Path) 91 args=p.parse_args();root=args.directory.resolve();root.mkdir(parents=True,exist_ok=False) 92 work=root/'workspace';work.mkdir() 93 sources={'driver.tri':'bench/float-finite-product-quotient.tri','finite.tri':'bench/float-finite.tri', 94 'adders.tri':'bench/float-adders.tri','multiply.tri':'bench/float-multiply.tri','near.tri':'bench/float-near-jam.tri', 95 'aligned.tri':'bench/float-aligned-quantize.tri','fast.tri':'bench/float-quantize-fast.tri', 96 'quantize.tri':'bench/float-quantize.tri','pipelines.tri':'bench/float-pipelines.tri', 97 'kernels.tri':'bench/float-significands.tri','layouts.tri':'bench/float-layouts.tri', 98 'native_math.tri':'lib/math/native.tri'} 99 for name,path in sources.items():(work/name).write_bytes((ROOT/path).read_bytes()) 100 report=dict(complete=False,scope='finite signed products and nonzero-divisor quotients; RNE; rejects nonfinite inputs', 101 runtime_sha256=sha(args.tricu),compiler_sha256=sha(ROOT/'bootstrap/build-stable.arboricx'), 102 harness_sha256=sha(Path(__file__)),sources={n:sha(ROOT/p) for n,p in sources.items()}, 103 validation=[],demand=[],samples=[],summaries=[]) 104 runner=Runner(args.tricu,root,2400,memory_mib=2048) 105 try: 106 tool=runner.imported(ROOT/'bootstrap/build-stable.arboricx','compiler');bundle=root/'finite.arboricx' 107 runner.run('compile',['compile','--file',work/'driver.tri','--toolchain',tool,'--tool-store',runner.store, 108 '--store',runner.store,'--output',bundle,'--compile-timeout','240'],270) 109 pin=runner.imported(bundle,'finite');report['bundle_sha256']=sha(bundle) 110 def run(label,name,raw,wanted,count,measure=False): 111 file=root/(label+'.input.tree');file.write_bytes(raw) 112 actual=runner.invoke(label,pin,name,['--arg-tree-file',file]+(['--profile-full-result'] if measure else [])) 113 expect(actual,wanted,label) 114 if not measure: 115 report['validation'].append(dict(label=label,cases=count,input_sha256=sha(file), 116 output_sha256=hashlib.sha256(actual).hexdigest()));return 117 points={v['profile']:v for line in (root/(label+'.log')).read_text().splitlines() 118 if line.startswith('{') for v in [json.loads(line)]} 119 a,b=points['invoke-begin-1'],points['force-end-1'] 120 return dict(wall_ns=b['wall_ns']-a['wall_ns'],allocated_bytes=b['allocated_bytes']-a['allocated_bytes']) 121 def material(rows,sub): 122 return tree_list([pair(nat(a),nat(b)) for a,b in rows]),tree_list([encoded(oracle(a,b,sub)) for a,b in rows]) 123 rows=cases();report['finite_pairs']=len(rows) 124 for name,sub in [('nativeMultiplies',False),('orderedMultiplies',False),('nativeDivides',True),('boundedDivides',True)]: 125 for offset in range(0,len(rows),12): 126 batch=rows[offset:offset+12];raw,wanted=material(batch,sub) 127 run(f'check-{name}-{offset}',name,raw,wanted,len(batch)) 128 raw_pairs=raw_cases();report['raw_quotient_pairs']=len(raw_pairs) 129 for offset in range(0,len(raw_pairs),12): 130 batch=raw_pairs[offset:offset+12] 131 raw=tree_list([pair(nat(a),nat(b)) for a,b in batch]) 132 wanted=tree_list([pair(nat((a<<56)//b),nat((a<<56)%b)) for a,b in batch]) 133 for name in ['nativeRaw','boundedRaw']: 134 run(f'raw-{name}-{offset}',name,raw,wanted,len(batch)) 135 bad=[0x7ff0000000000000,0xfff0000000000000,0x7ff8000000000001, 136 0xfff8000000000001,0x7ff0000000000001,0xfff0000000000001,LIMIT,LIMIT+1,2**128] 137 rejects=[(a,b) for a in bad for b in [0,SIGN,0x3ff0000000000000]] 138 rejects += [(b,a) for a,b in rejects]+[(bad[0],LIMIT),(LIMIT,bad[0])] 139 report['rejection_pairs']=len(rejects) 140 for name,sub in [('nativeMultiplies',False),('orderedMultiplies',False),('nativeDivides',True),('boundedDivides',True)]: 141 raw,wanted=material(rejects,sub);run('reject-'+name,name,raw,wanted,len(rejects)) 142 for name,wanted in [ 143 ('rangeDemand',encoded((False,2))), 144 ('productZeroDemand',pair(nat(0),pair(LEAF,LEAF))), 145 ('quotientZeroDemand',pair(nat(0),pair(LEAF,LEAF))), 146 ('jamOddDemand',nat(7)), 147 ('walkExactDemand',pair(nat(2**56),nat(0))), 148 ('divisionRejectDemand',encoded((False,3)))]: 149 actual=runner.invoke('demand-'+name,pin,name,['--arg','0']) 150 expect(actual,wanted,name);report['demand'].append(dict(name=name,outcome='expected-result')) 151 for i,(name,a,b,divide,wrong) in enumerate(negative_controls()): 152 correct=oracle(a,b,divide);assert wrong!=correct 153 run('negative-'+str(i),name,tree_list([pair(nat(a),nat(b))]),tree_list([encoded(wrong)]),1) 154 report.setdefault('negative_controls',[]).append(dict(name=name,words=[a,b],correct=correct,wrong=wrong)) 155 rng=random.Random(64054) 156 groups={ 157 'dense':[(0x3fffffffffffffff-i,0x3fefffffffffffff-i) for i in range(16)], 158 'random':[(0x3ff0000000000000|rng.getrandbits(52),0x3ff0000000000000|rng.getrandbits(52)) for _ in range(16)], 159 'exact':[(0xbff0000000000000+i,0xbff0000000000000+i) for i in range(16)], 160 'subnormal':[(SIGN|(1+i),0x3fe0000000000000+(i%3)-1) for i in range(16)], 161 'tiny':[(SIGN|(1+i),0x0010000000000000+i) for i in range(16)], 162 'zero':[(SIGN,0xbff0000000000000+i) for i in range(16)], 163 'overflow':[(0x7fefffffffffffff-i,0x4000000000000000+i) for i in range(16)]} 164 jobs=[] 165 for group,batch in groups.items(): 166 for name,sub in [('nativeMultiplies',False),('orderedMultiplies',False),('nativeDivides',True),('boundedDivides',True)]: 167 raw,wanted=material(batch,sub);jobs.append((group,name,raw,wanted,len(batch))) 168 rng=random.Random(64053) 169 for rep in range(5): 170 schedule=list(jobs);rng.shuffle(schedule) 171 for group,name,raw,wanted,size in schedule: 172 sample=run(f'measure-{group}-{name}-{rep}',name,raw,wanted,size,True) 173 report['samples'].append(dict(group=group,name=name,size=size,repeat=rep,warmup=rep==0,**sample)) 174 for group,name,_,_,size in jobs: 175 samples=[s for s in report['samples'] if s['group']==group and s['name']==name and not s['warmup']] 176 row=dict(group=group,name=name,size=size) 177 for key in ['wall_ns','allocated_bytes']: 178 row[key]=statistics.median(s[key] for s in samples) 179 row[key+'_min']=min(s[key] for s in samples);row[key+'_max']=max(s[key] for s in samples) 180 report['summaries'].append(row);print(json.dumps(row),flush=True) 181 report['complete']=True 182 finally:(root/'report.json').write_text(json.dumps(report,indent=2)+'\n') 183 184 185 if __name__=='__main__':main()