From 983d66f0f2e48ee500ec7eee2da9244c001c5064 Mon Sep 17 00:00:00 2001 From: Charles Bethin Date: Tue, 13 Oct 2026 01:16:09 -0400 Subject: [PATCH] Prove sharp inverse-square rounding bounds for cyclic scans --- .../cyclic-scan-sharp-bound/EXPECTED.json | 390 ++++++++++++++++++ research/cyclic-scan-sharp-bound/PROOF.md | 177 ++++++++ research/cyclic-scan-sharp-bound/README.md | 97 +++++ research/cyclic-scan-sharp-bound/SOURCE.json | 29 ++ .../cyclic-scan-sharp-bound/VALIDATION.json | 99 +++++ .../cyclic-scan-sharp-bound/fixed_scan.py | 59 +++ .../cyclic-scan-sharp-bound/vendor/NOTICE.md | 24 ++ .../vendor/check_budgets.py | 59 +++ .../cyclic-scan-sharp-bound/vendor/scan.py | 123 ++++++ .../vendor/weighted.py | 118 ++++++ research/cyclic-scan-sharp-bound/verify.py | 257 ++++++++++++ 11 files changed, 1432 insertions(+) create mode 100644 research/cyclic-scan-sharp-bound/EXPECTED.json create mode 100644 research/cyclic-scan-sharp-bound/PROOF.md create mode 100644 research/cyclic-scan-sharp-bound/README.md create mode 100644 research/cyclic-scan-sharp-bound/SOURCE.json create mode 100644 research/cyclic-scan-sharp-bound/VALIDATION.json create mode 100644 research/cyclic-scan-sharp-bound/fixed_scan.py create mode 100644 research/cyclic-scan-sharp-bound/vendor/NOTICE.md create mode 100644 research/cyclic-scan-sharp-bound/vendor/check_budgets.py create mode 100644 research/cyclic-scan-sharp-bound/vendor/scan.py create mode 100644 research/cyclic-scan-sharp-bound/vendor/weighted.py create mode 100644 research/cyclic-scan-sharp-bound/verify.py diff --git a/research/cyclic-scan-sharp-bound/EXPECTED.json b/research/cyclic-scan-sharp-bound/EXPECTED.json new file mode 100644 index 0000000000..ea3687ca67 --- /dev/null +++ b/research/cyclic-scan-sharp-bound/EXPECTED.json @@ -0,0 +1,390 @@ +{ + "status": "PASS_SHARP_CYCLIC_SCAN_BOUND", + "source_head": "0768f97d06c623b8211e651543c66135796b62bf", + "seed": 68, + "counts": { + "closed_form_cases": 5024, + "residual_cases": 24, + "original_weighted_matches": 3362, + "original_constant_matches": 482, + "expanded_admission_cases": 1686 + }, + "scan_transcript_sha256": "5880549eac0a0b8deace1c2a7c15c543d4f51c9a8c37f7cbce5578e36f3059b1", + "sharp_family": [ + { + "k": 2, + "beta": "33/128", + "error_over_u_beta_squared": "3135/8192" + }, + { + "k": 3, + "beta": "129/1024", + "error_over_u_beta_squared": "115455/262144" + }, + { + "k": 4, + "beta": "513/8192", + "error_over_u_beta_squared": "3939327/8388608" + }, + { + "k": 5, + "beta": "2049/65536", + "error_over_u_beta_squared": "130084863/268435456" + }, + { + "k": 6, + "beta": "8193/524288", + "error_over_u_beta_squared": "4228366335/8589934592" + }, + { + "k": 7, + "beta": "32769/4194304", + "error_over_u_beta_squared": "136369340415/274877906944" + }, + { + "k": 8, + "beta": "131073/33554432", + "error_over_u_beta_squared": "4380899934207/8796093022208" + }, + { + "k": 9, + "beta": "524289/268435456", + "error_over_u_beta_squared": "140462877835263/281474976710656" + }, + { + "k": 10, + "beta": "2097153/2147483648", + "error_over_u_beta_squared": "4499203724148735/9007199254740992" + }, + { + "k": 11, + "beta": "8388609/17179869184", + "error_over_u_beta_squared": "144044836494770175/288230376151711744" + }, + { + "k": 12, + "beta": "33554433/137438953472", + "error_over_u_beta_squared": "4610560255892389887/9223372036854775808" + }, + { + "k": 13, + "beta": "134217729/1099511627776", + "error_over_u_beta_squared": "147555939290410123263/295147905179352825856" + }, + { + "k": 14, + "beta": "536870913/8796093022208", + "error_over_u_beta_squared": "4722078261288512782335/9444732965739290427392" + }, + { + "k": 15, + "beta": "2147483649/70368744177664", + "error_over_u_beta_squared": "151111115836174668660735/302231454903657293676544" + }, + { + "k": 16, + "beta": "8589934593/562949953421312", + "error_over_u_beta_squared": "4835629492045154634170367/9671406556917033397649408" + }, + { + "k": 17, + "beta": "34359738369/4503599627370496", + "error_over_u_beta_squared": "154741324323555347858980863/309485009821345068724781056" + }, + { + "k": 18, + "beta": "137438953473/36028797018963968", + "error_over_u_beta_squared": "4951741267711618143156699135/9903520314283042199192993792" + }, + { + "k": 19, + "beta": "549755813889/288230376151711744", + "error_over_u_beta_squared": "158456022797362000806434308095/316912650057057350374175801344" + }, + { + "k": 20, + "beta": "2199023255553/2305843009213693952", + "error_over_u_beta_squared": "5070597565211944986081281179647/10141204801825835211973625643008" + }, + { + "k": 21, + "beta": "8796093022209/18446744073709551616", + "error_over_u_beta_squared": "162259199457979354781792352600063/324518553658426726783156020576256" + }, + { + "k": 22, + "beta": "35184372088833/147573952589676412928", + "error_over_u_beta_squared": "5192295620594935917032442362331135/10384593717069655257060992658440192" + }, + { + "k": 23, + "beta": "140737488355329/1180591620717411303424", + "error_over_u_beta_squared": "166153479666075036138230726583648255/332306998946228968225951765070086144" + }, + { + "k": 24, + "beta": "562949953421313/9444732965739290427392", + "error_over_u_beta_squared": "5316911666227022879289717706329161727/10633823966279326983230456482242756608" + }, + { + "k": 25, + "beta": "2251799813685249/75557863725914323419136", + "error_over_u_beta_squared": "170141178389866906376628920043767332863/340282366920938463463374607431768211456" + }, + { + "k": 26, + "beta": "9007199254740993/604462909807314587353088", + "error_over_u_beta_squared": "5444517789605377605270203816035364110335/10889035741470030830827987437816582766592" + }, + { + "k": 27, + "beta": "36028797018963969/4835703278458516698824704", + "error_over_u_beta_squared": "174224570565446283495244098273363902857215/348449143727040986586495598010130648530944" + }, + { + "k": 28, + "beta": "144115188075855873/38685626227668133590597632", + "error_over_u_beta_squared": "5575186278863468389930244993477862498500607/11150372599265311570767859136324180752990208" + }, + { + "k": 29, + "beta": "576460752303423489/309485009821345068724781056", + "error_over_u_beta_squared": "178405961255937986495541786623658691095691263/356811923176489970264571492362373784095686656" + }, + { + "k": 30, + "beta": "2305843009213693953/2475880078570760549798248448", + "error_over_u_beta_squared": "5708990765506927543569360460141826835780468735/11417981541647679048466287755595961091061972992" + }, + { + "k": 31, + "beta": "9223372036854775809/19807040628566084398385987584", + "error_over_u_beta_squared": "182687704581292273065033028833811065923725950975/365375409332725729550921208179070754913983135744" + }, + { + "k": 32, + "beta": "36893488147419103233/158456325028528675187087900672", + "error_over_u_beta_squared": "5846006547962482205289441801966375347788907020287/11692013098647223345629478661730264157247460343808" + }, + { + "k": 33, + "beta": "147573952589676412929/1267650600228229401496703205376", + "error_over_u_beta_squared": "187072209556577502048399247525961333136643550347263/374144419156711147060143317175368453031918731001856" + }, + { + "k": 34, + "beta": "590295810358705651713/10141204801825835211973625643008", + "error_over_u_beta_squared": "5986310706158929209245393293020044894091825261838335/11972621413014756705924586149611790497021399392059392" + }, + { + "k": 35, + "beta": "2361183241434822606849/81129638414606681695789005144064", + "error_over_u_beta_squared": "191561942602660920995241852246814602699393219256385535/383123885216472214589586756787577295904684780545900544" + }, + { + "k": 36, + "beta": "9444732965739290427393/649037107316853453566312041152512", + "error_over_u_beta_squared": "6129982163374352452639914653142391825371997876170784767/12259964326927110866866776217202473468949912977468817408" + }, + { + "k": 37, + "beta": "37778931862957161709569/5192296858534827628530496329220096", + "error_over_u_beta_squared": "196159429229406526177167651891039748803513917085124132863/392318858461667547739736838950479151006397215279002157056" + }, + { + "k": 38, + "beta": "151115727451828646838273/41538374868278621028243970633760768", + "error_over_u_beta_squared": "6277101735363844800752535603485601816980736862595250651135/12554203470773361527671578846415332832204710888928069025792" + }, + { + "k": 39, + "beta": "604462909807314587353089/332306998946228968225951765070086144", + "error_over_u_beta_squared": "200867255532008409033412868120093349127138709805071016656895/401734511064747568885490523085290650630550748445698208825344" + }, + { + "k": 40, + "beta": "2417851639229258349412353/2658455991569831745807614120560689152", + "error_over_u_beta_squared": "6427752177030115095618527416147827235753989372208151472898047/12855504354071922204335696738729300820177623950262342682411008" + }, + { + "k": 41, + "beta": "9671406556917033397649409/21267647932558653966460912964485513216", + "error_over_u_beta_squared": "205688069665057219164581991300551709832859294438018400926040063/411376139330301510538742295639337626245683966408394965837152256" + }, + { + "k": 42, + "beta": "38685626227668133590597633/170141183460469231731687303715884105728", + "error_over_u_beta_squared": "6582018229283327590943250055782344907054102467954159227173339135/13164036458569648337239753460458804039861886925068638906788872192" + }, + { + "k": 43, + "beta": "154742504910672534362390529/1361129467683753853853498429727072845824", + "error_over_u_beta_squared": "210624583337090428153010027215058483148935336063781264368185901055/421249166674228746791672110734681729275580381602196445017243910144" + }, + { + "k": 44, + "beta": "618970019642690137449562113/10889035741470030830827987437816582766592", + "error_over_u_beta_squared": "6739986666787276824781537310429353823121291419127940493037040304127/13479973333575319897333507543509815336818572211270286240551805124608" + }, + { + "k": 45, + "beta": "2475880078570760549798248449/87112285931760246646623899502532662132736", + "error_over_u_beta_squared": "215679573337198988375172657227835897932743849030056862873773213220863/431359146674410236714672241392314090778194310760649159697657763987456" + }, + { + "k": 46, + "beta": "9903520314283042199192993793/696898287454081973172991196020261297061888", + "error_over_u_beta_squared": "6901746346790465707720140446086988805695433518514997585794567384334335/13803492693581127574869511724554050904902217944340773110325048447598592" + }, + { + "k": 47, + "beta": "39614081257132168796771975169/5575186299632655785383929568162090376495104", + "error_over_u_beta_squared": "220855883097296471922478340928249041830712261750069759480164155985494015/441711766194596082395824375185729628956870974218904739530401550323154944" + }, + { + "k": 48, + "beta": "158456325028528675187087900673/44601490397061246283071436545296723011960832", + "error_over_u_beta_squared": "7067388259113512209926248456293220210549303686207232802552593297868587007/14134776518227074636666380005943348126619871175004951664972849610340958208" + }, + { + "k": 49, + "beta": "633825300114114700748351602689/356811923176489970264571492362373784095686656", + "error_over_u_beta_squared": "226156424291632792452151015347881496458571342212042567352950120354638987263/452312848583266388373324160190187140051835877600158453279131187530910662656" + }, + { + "k": 50, + "beta": "2535301200456458802993406410753/2854495385411919762116571938898990272765493248", + "error_over_u_beta_squared": "7237005577332255786221009527084746568366416591643639334680105253609182068735/14474011154664524427946373126085988481658748083205070504932198000989141204992" + }, + { + "k": 51, + "beta": "10141204801825835211973625643009/22835963083295358096932575511191922182123945984", + "error_over_u_beta_squared": "231584178474632288003107137442021016984049354834689089629783087415319672651775/463168356949264781694283940034751631413079938662562256157830336031652518559744" + }, + { + "k": 52, + "beta": "40564819207303340847894502572033/182687704666362864775460604089535377456991567872", + "error_over_u_beta_squared": "7410693711188234861603985719350166635344762824034136937978902596608244934770687/14821387422376473014217086081112052205218558037201992197050570753012880593911808" + }, + { + "k": 53, + "beta": "162259276829213363391578010288129/1461501637330902918203684832716283019655932542976", + "error_over_u_beta_squared": "237142198758023541899400460158497622305627338580217480560209571467164775836811263/474284397516047136454946754595585670566993857190463750305618264096412179005177856" + }, + { + "k": 54, + "beta": "649037107316853453566312041152513/11692013098647223345629478661730264157247460343808", + "error_over_u_beta_squared": "7588550360256753762029981399300635629412889627587738284573538646703174470216974335/15177100720513508366558296147058741458143803430094840009779784451085189728165691392" + }, + { + "k": 55, + "beta": "2596148429267413814265248164610049/93536104789177786765035829293842113257979682750464", + "error_over_u_beta_squared": "242833611528216127124946071565280008199651872304391344486161204596704642858102030335/485667223056432267729865476705879726660601709763034880312953102434726071301302124544" + }, + { + "k": 56, + "beta": "10384593717069655257060992658440193/748288838313422294120286634350736906063837462003712", + "error_over_u_beta_squared": "7770675568902916175838060958691517196190405321552308741557872160219719158258538119167/15541351137805832567355695254588151253139254712417116170014499277911234281641667985408" + }, + { + "k": 57, + "beta": "41538374868278621028243970633760769/5986310706507378352962293074805895248510699696029696", + "error_over_u_beta_squared": "248661618204893319352254537375769479177849816336795765492998083652636301738394060324863/497323236409786642155382248146820840100456150797347717440463976893159497012533375533056" + }, + { + "k": 58, + "beta": "166153499473114484112975882535043073/47890485652059026823698344598447161988085597568237568", + "error_over_u_beta_squared": "7957171782556586246879130583186878339758760615708488160634333393827199493684269611483135/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792" + }, + { + "k": 59, + "beta": "664613997892457936451903530140172289/383123885216472214589586756787577295904684780545900544", + "error_over_u_beta_squared": "254629497041810760341843944856576188118733059238984615633005342953054345940893258215325695/509258994083621521567111422102344540262867098416484062659035112338595324940834176545849344" + }, + { + "k": 60, + "beta": "2658455991569831745807614120560689153/3064991081731777716716694054300618367237478244367204352", + "error_over_u_beta_squared": "8148143905337944338006394494523975328937674653423848650634408281076004242189172022808936447/16296287810675888690147565507275025288411747149327490005089123594835050398106693649467179008" + }, + { + "k": 61, + "beta": "10633823966279326983230456482242756609/24519928653854221733733552434404946937899825954937634816", + "error_over_u_beta_squared": "260740604970814218929282835970583807545776842995547276780932912619167486021363214499199320063/521481209941628438084722096232800809229175908778479680162851955034721612739414196782949728256" + }, + { + "k": 62, + "beta": "42535295865117307932921825928971026433/196159429230833773869868419475239575503198607639501078528", + "error_over_u_beta_squared": "8343699359066055007546302145391747394369677328925760096183085597899916950351457354185112027135/16687398718132110018711107079449625895333629080911349765211262561111091607661254297054391304192" + }, + { + "k": 63, + "beta": "170141183460469231731687303715884105729/1569275433846670190958947355801916604025588861116008628224", + "error_over_u_beta_squared": "266998379490113760270429690961864965471014594476256289470681369800248302853277887004576218873855/533996758980227520598755426542388028650676130589163192486760401955554931445160137505740521734144" + }, + { + "k": 64, + "beta": "680564733841876926926749214863536422913/12554203470773361527671578846415332832204710888928069025792", + "error_over_u_beta_squared": "8543948143683640329116918467728943676729088352862632809901396942005048225943083076519396825366527/17087896287367280659160173649356416916821636178853222159576332862577757806245124400183696695492608" + }, + { + "k": 65, + "beta": "2722258935367507707706996859454145691649/100433627766186892221372630771322662657637687111424552206336", + "error_over_u_beta_squared": "273406340597876490539152084678514434162138069448861715342827974895543716448638928719825223091748863/546812681195752981093125556779405341338292357723303109106442651602488249799843980805878294255763456" + } + ], + "cascade_recurrence_cases": 384, + "budget_cases": 276, + "budget_examples": [ + { + "radix": 3, + "p": 128, + "g": 27, + "old": 244, + "new": 217, + "saved": 27 + }, + { + "radix": 4, + "p": 128, + "g": 41, + "old": 342, + "new": 301, + "saved": 41 + }, + { + "radix": 9, + "p": 128, + "g": 465, + "old": 4790, + "new": 4324, + "saved": 466 + }, + { + "radix": 16, + "p": 128, + "g": 1967, + "old": 33580, + "new": 31613, + "saved": 1967 + }, + { + "radix": 256, + "p": 128, + "g": 8356379, + "old": 2147589552, + "new": 2139233173, + "saved": 8356379 + } + ], + "budget_transcript_sha256": "68ff59cc4a71b22745c7c00726b0b09dd75abe0943f0f4f116b9703127b8b02b", + "controls": [ + "zero_seed", + "plus_denominator", + "constant_coefficient", + "wrong_orientation", + "unstable_coefficient", + "below_closure_guard", + "false_inverse_linear_margin_bound" + ], + "scope": "Exact finite support for the written all-length proof; no factorizer, full kernel, tape runtime, or multiplication theorem validation." +} diff --git a/research/cyclic-scan-sharp-bound/PROOF.md b/research/cyclic-scan-sharp-bound/PROOF.md new file mode 100644 index 0000000000..6fa0440ed5 --- /dev/null +++ b/research/cyclic-scan-sharp-bound/PROOF.md @@ -0,0 +1,177 @@ +# A sharp stability-margin power for periodic scan rounding + +The two-pass periodic scan in [PR628](https://github.com/CrocSwap/integer-mult-bounds/pull/628) +is sound under its stated hypotheses. Its proof gives the sufficient bound +`6(B+1)u/beta^3` and requires `u <= beta^2/2`. +Keeping a common geometric factor improves both conclusions. +This result concerns that scan, including its varying-coefficient extension. +It does not change the selected multiplication construction. + +## Statement + +Let `L >= 1`, `0 < beta <= 1`, and `B >= 0`. +For every cyclic sequence of exact complex coefficients and inputs, assume + +``` +|b_i| <= 1-beta, |x_i| <= B. +``` + +Choose either cyclic orientation. Compute the zero-seed pass, the running +product, the closure division, and the second pass as in PR628. +Each completed complex operation has absolute rounding error at most `u`. +Assume `0 < u <= beta/2`. Then the rounded closure denominator is nonzero and + +``` +max_i |yhat_i - y_i| <= u (2B/beta^2 + 3/beta + 1) + <= 4(B+1)u/beta^2, +``` + +where `y_i-b_i y_(i-d)=x_i` is the unique periodic solution. +The bound is uniform in `L`, the orientation, and the coefficient pattern. +Nearest rounding of both components on a grid of spacing `u` meets the +operation-error hypothesis, since its complex error is at most `u/sqrt(2)`. +Coefficients are exact; coefficient approximation requires a separate budget. + +## Proof + +Relabel the chosen traversal as `1,...,L`. +Put `q_i=|b_i|`, `Q=product_i q_i`, and + +``` +G = sum_(j=1)^L product_(k=j+1)^L q_k. +``` + +An empty product is one. The exact affine return map is `s -> z+c*s`, +where `c=product_i b_i` and `z` is the last zero-seed value. +Thus the exact closure denominator is `d=1-c`. +The following identity is the key estimate: + +``` +1-Q = sum_(j=1)^L (1-q_j) product_(k=j+1)^L q_k >= beta*G. +``` + +Also `|d| >= 1-Q > 0`. Both rounded first-pass recurrences multiply +their previous errors by `b_i` and add one error of modulus at most `u`. +Consequently + +``` +|zhat-z| <= u*G, |chat-c| <= u*G <= (u/beta)|d|. +``` + +With `dhat=1-chat`, the hypothesis `u <= beta/2` gives +`|dhat| >= |d|/2`. Computing `1-chat` is exact on the dyadic grid. +This is the step that avoids bounding the numerator and denominator with +two unrelated worst-case geometric sums. + +The cyclic weighted shift has norm at most `1-beta` in the maximum norm. +Its Neumann inverse therefore has norm at most `1/beta`. +In particular, the exact periodic seed `s=z/d` satisfies `|s| <= B/beta`. +Before rounding the closure division, + +``` +|zhat/dhat-s| + <= (|zhat-z| + |s| |dhat-d|)/|dhat| + <= 2u/beta + 2B*u/beta^2. +``` + +Division rounding adds at most `u`. +At each output, the second pass attenuates the seed error by a factor at +most one and contributes at most `u/beta` from its own rounding errors. +This proves the first bound. The second follows from `beta <= 1`: + +``` +2B/beta^2 + 3/beta + 1 <= (2B+4)/beta^2 <= 4(B+1)/beta^2. +``` + +All statements include `L=1` and zero coefficients. No diagonalization, +probabilistic rounding assumption, or limit in `L` enters the proof. + +## The inverse-square dependence is necessary for this algorithm + +For every integer `k >= 2`, take a real one-entry input with + +``` +L=1, x=1, u=2^(-3k), a=2^(-k), beta=a+u/2, b=1-beta. +``` + +Use the implemented nearest rule, with ties toward positive infinity. +All inputs and coefficients are exact dyadics. The working grid is `u`. +The first-pass value is exactly one. Rounding the product `b*1` gives +`chat=1-a`, because `b` is exactly halfway between adjacent grid points. +The closure seed is exactly `1/a=2^k`. +The second-pass value is also exactly on the grid: + +``` +yhat = 1+b/a = 2^k-2^(k-1)u, y = 1/beta. +yhat-y = b*u/(2*a*beta) >= u/(4*beta^2). +``` + +The inequality uses `a <= beta` and `b >= 1/2`. +Moreover `u <= beta^2/2`, so every example satisfies even PR628's older +precision admission. As `k` increases, `beta` tends to zero. +For any `q < 2`, no constant independent of `beta` can replace the error +bound by `C*u/beta^q` for this algorithm, even with `B=1` and `L=1`. +This is not a lower bound for other periodic solvers or a claim about an +optimally rounded exact closure product. + +## Consequence for the existing factor cascade + +Retain PR628's hypotheses and its exact certified factors. +An outer scalar or inner unit phase contributes at most one extra rounding +error `u` and has modulus at most one. +Thus each factor has error at most `5(B+1)u/beta^2`. +Put `n=r-1`, as in PR628. The unchanged exact prefix norm is bounded by +`r*beta^(-j)` after `j` factors. If `E_j` is the accumulated error, then + +``` +E_j <= E_(j-1)/beta + + 5(r*beta^(-(j-1)) + E_(j-1) + 1)u/beta^2. +T_j = beta^j E_j + <= (1+5u/beta)T_(j-1) + 5(r+1)u/beta. +``` + +For `u <= beta/(10n)`, the growth factor over `n` stages is at most +`(1+1/(2n))^n < 2`. +One elementary proof expands the binomial and bounds the degree-`j` +coefficient by `1/(2^j j!)`; the finite sum is strictly below +`sum_(j>=0) 2^(-j)=2`. +It follows that + +``` +E_n <= 10n(r+1)u*beta^(-(n+1)). +``` + +For `beta=2^-g`, the following sufficient fractional precision replaces +the execution budget in PR628's PROOF.md: + +``` +w_new = max(g + ceil(log2(10n)), + p + (n+1)g + ceil(log2(20n(r+1)))). +``` + +It ensures `E_n <= 2^(-p-1)`. All coefficient, root, phase-oracle and +inverse-perturbation budgets stay as stated in PR628. +The original execution budget is + +``` +w_old = max(2g + ceil(log2(14n)), + p + (n+2)g + ceil(log2(28n(r+1)))). +``` + +Both new terms are at most their old counterparts minus `g`. +Thus `w_new <= w_old-g`. This is a proved saving in a sufficient +execution-grid budget. It is not a measured speedup or a change to the +`O(p+r^4)` asymptotic precision and preparation bound. +The factor-cascade consequence inherits PR628's exact-factor and phase +assumptions; this package does not independently validate its factorizer. + +## Attribution and status + +PR628 by Chafik Boukhalfa supplies the algorithm, the original rounding +argument, its varying-phase extension, and the surrounding factor cascade. +Periodic affine closure and geometric error bounds are standard methods. +This note adds the shared-factor bound, the sharp family for this algorithm, +and the consequent execution guard. No literature-priority claim is made. +Prepared for Charles Bethin with substantial OpenAI Codex assistance. +No human peer review or proof-assistant verification is claimed. +New material uses Apache-2.0. Vendored notices remain intact. diff --git a/research/cyclic-scan-sharp-bound/README.md b/research/cyclic-scan-sharp-bound/README.md new file mode 100644 index 0000000000..94c006751f --- /dev/null +++ b/research/cyclic-scan-sharp-bound/README.md @@ -0,0 +1,97 @@ +# Sharp rounding bounds for cyclic scans + +**Mathematical refinement, no new kappa.** The clean two-pass periodic scan +from [PR628](https://github.com/CrocSwap/integer-mult-bounds/pull/628) has +length-independent error at most `4(B+1)u/beta^2`, under `u <= beta/2`. +The existing sufficient bound is `6(B+1)u/beta^3`, under `u <= beta^2/2`. +Here `u` is the working grid spacing, `B` bounds the input magnitude, +and all exact complex coefficients have modulus at most `1-beta`. + +[PROOF.md](PROOF.md) proves the result for every cycle length, either +orientation, and arbitrary varying coefficients within the margin. +A one-entry dyadic family has error at least `u/(4*beta^2)` for the same +algorithm and rounding rule. Thus the inverse-square margin dependence is +necessary for a uniform bound on this algorithm. + +The proof also lowers PR628's sufficient factor-cascade execution precision +by at least `g` bits when `beta=2^-g`. It preserves the coefficient and +phase-preparation assumptions. The asymptotic `O(p+r^4)` precision bound +stays the same. No runtime improvement, factorization implementation, +supplier admission, selected frontier change, or full multiplication proof +follows. This package proposes a guard refinement; it does not activate it +inside PR628 or modify that contribution's files. + +## Reproduce + +From any working directory, run with Python 3.11 or later: + +```sh +python3 -B /path/to/repository/research/cyclic-scan-sharp-bound/verify.py +``` + +Expected status: `PASS_SHARP_CYCLIC_SCAN_BOUND`. +The command checks source hashes and compares its full result with +[EXPECTED.json](EXPECTED.json). It uses only the standard library. +To preserve a fresh receipt, use `--write /absolute/new/result.json`. +The output path must not exist. The default command never overwrites evidence. + +## Evidence and independence + +- 5,024 cases compare the rounded output with an independently evaluated + finite geometric inverse. This oracle evaluates each matrix row directly. +- 24 longer cases, through length 2,048, use exact residual certificates. + Dividing the residual norm by `beta` bounds the forward error. +- 3,362 cases match PR628's unchanged varying-coefficient scan exactly. + Of these, 482 also match its unchanged constant-coefficient scan. +- 1,686 cases meet the new guard but lie outside the original admission. + The new integer-grid implementation executes these cases. +- 64 exact members of the sharp family match the symbolic output and error + formulas in the original scan and the integer-grid implementation. +- 384 exact cascade recurrences and 276 inherited-budget comparisons check + the precision formulas. The written inequalities establish their scope. +- Seven controls reject missing closure, wrong denominator sign, substituted + coefficients, wrong orientation, unstable coefficients, insufficient + admitted precision, and a false inverse-linear margin bound. + +The original-scan comparisons establish implementation agreement, not +independent validation of a shared traversal. The direct finite inverse, +residual certificates, and symbolic one-entry family provide separate +mathematical checks. None of these finite checks proves the all-length +theorem by enumeration. There is no independent human or formal review. + +The exhaustive domain uses lengths 1–3, all real coefficient sequences over +`{-1/2,0,1/2}`, all inputs over `{-1,0,1}`, both orientations, and three grids. +Additional complex cases use the recorded seed 68 and exact rational unit +directions. Long cases use two margins, two grids, and both orientations. +These are fixed validation domains, not an optimization search. + +## Sources and scope + +[SOURCE.json](SOURCE.json) pins PR628 at +`0768f97d06c623b8211e651543c66135796b62bf`, with construction base +`66849e1475c930fd7732b1a0067ffee7e1e95f60`. +The vendored files match that Git source byte for byte. +`verify.py` executes the unchanged `scan` function AST from `weighted.py` +with its original arithmetic module. It does not execute the unrelated +kernel factorization imports or claim their validation. +The vendored weighted file is an archival source, not a standalone command +inside this reduced package. + +PR663's nonlinear restoration audit and PR673's dirty-echo construction and +alias repair address separate questions. Their results are not dependencies +of this clean-array rounding theorem. Existing dirty restoration, arbitrary +stride routing, factor preparation, and full-transform obligations remain. + +Joshua Snider's [PR692](https://github.com/CrocSwap/integer-mult-bounds/pull/692), +at `c5b53a5e80372a1fae9830d8dcf21aaabce3f729`, already uses the standard finite +geometric inverse to audit weighted stride cycles and their bilinear dual. +This package does not claim that oracle identity or method as new. +Its distinct result is the uniform sharp rounding bound and execution guard. +The oracle here was written independently; no PR692 code is consumed. + +Chafik Boukhalfa's PR628 supplies the scan and its prior proof. +Its predecessor attribution is preserved in [vendor/NOTICE.md](vendor/NOTICE.md). +The shared-factor estimate, matching family, guard derivation, and new +validation code were prepared for Charles Bethin with substantial OpenAI +Codex assistance. No literature-priority claim is made. +New material is Apache-2.0, under the repository license. diff --git a/research/cyclic-scan-sharp-bound/SOURCE.json b/research/cyclic-scan-sharp-bound/SOURCE.json new file mode 100644 index 0000000000..57a408f5c5 --- /dev/null +++ b/research/cyclic-scan-sharp-bound/SOURCE.json @@ -0,0 +1,29 @@ +{ + "base_commit": "66849e1475c930fd7732b1a0067ffee7e1e95f60", + "pr": 628, + "head": "0768f97d06c623b8211e651543c66135796b62bf", + "vendor": { + "vendor/NOTICE.md": { + "source_path": "research/prime-power-cyclic-scan/NOTICE.md", + "sha256": "abeacf1ad1491d5797cfcf1f3832799bbf3b24b99597229d7bc884059a51ec68" + }, + "vendor/check_budgets.py": { + "source_path": "research/prime-power-cyclic-scan/check_budgets.py", + "sha256": "581ba58a30e05cb09dadb35be6e25bd92232b2137a0cc7eb5b83e2ee20c0bf02" + }, + "vendor/scan.py": { + "source_path": "research/prime-power-cyclic-scan/scan.py", + "sha256": "1a9c01ca1a29a9be1e756858161202621db29d5d332a4266b4017a9502a50d7e" + }, + "vendor/weighted.py": { + "source_path": "research/prime-power-cyclic-scan/weighted.py", + "sha256": "8b0f93d131a30e64eb0ccd893779daffea79c0c00f1d2a30a91d23eb2bf0bff5" + } + }, + "proof_sources": { + "research/prime-power-cyclic-scan/PROOF.md": "595893609764c64a90d9b21b9cbd6d9e61980703f3deb43361105cf8c3d1d093", + "research/prime-power-cyclic-scan/SCAN-PROOF.md": "6be99f0dce58d82ecb07967664e4a23f202ae3858f124a1bdcdcac3e660548db", + "research/prime-power-cyclic-scan/WEIGHTED-PROOF.md": "b0643ad6ae448ca2fe11ee9a9597b5656e5e2feaaa89b16939e1d1bc7ba7437d" + }, + "seed": 68 +} diff --git a/research/cyclic-scan-sharp-bound/VALIDATION.json b/research/cyclic-scan-sharp-bound/VALIDATION.json new file mode 100644 index 0000000000..5171f0792c --- /dev/null +++ b/research/cyclic-scan-sharp-bound/VALIDATION.json @@ -0,0 +1,99 @@ +{ + "status": "local_checks_passed", + "recorded_at": "2026-10-13T05:15:32.964616+00:00", + "platform": "macOS-27.2-arm64-arm-64bit", + "interpreters": { + "python3.13": "Python 3.13.11", + "python3.14": "Python 3.14.6" + }, + "checks": [ + { + "command": "python3.13 -B research/cyclic-scan-sharp-bound/verify.py --write /absolute/new/result.json", + "status": "passed; receipt preserved as EXPECTED.json" + }, + { + "command": "python3.14 -I -O -B research/cyclic-scan-sharp-bound/verify.py", + "status": "passed; full expected receipt comparison" + }, + { + "command": "make selected-record-check", + "status": "passed; selected result unchanged" + }, + { + "command": "git diff --check", + "status": "passed" + } + ], + "packaging": { + "status": "PASS_PORTABLE_AND_PACKAGING_CONTROLS", + "records": [ + { + "name": "minimal_package", + "command": [ + "python3.13", + "-I", + "-B", + "/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py" + ], + "exit": 0, + "output": "{\"status\": \"PASS_SHARP_CYCLIC_SCAN_BOUND\", \"counts\": {\"closed_form_cases\": 5024, \"residual_cases\": 24, \"original_weighted_matches\": 3362, \"original_constant_matches\": 482, \"expanded_admission_cases\": 1686}, \"cascade_recurrence_cases\": 384, \"budget_cases\": 276, \"controls\": [\"zero_seed\", \"plus_denominator\", \"constant_coefficient\", \"wrong_orientation\", \"unstable_coefficient\", \"below_closure_guard\", \"false_inverse_linear_margin_bound\"]}\n" + }, + { + "name": "source_drift", + "command": [ + "python3.13", + "-I", + "-B", + "/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py" + ], + "exit": 1, + "output": "Traceback (most recent call last):\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 257, in \n main()\n ~~~~^^\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 104, in main\n source, original, weighted, old_budgets = source_check()\n ~~~~~~~~~~~~^^\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 83, in source_check\n require(hashlib.sha256((ROOT/path).read_bytes()).hexdigest() == pin['sha256'],\n ~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n 'source drift: '+path)\n ^^^^^^^^^^^^^^^^^^^^^^\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 26, in require\n raise ValueError(message)\nValueError: source drift: vendor/scan.py\n" + }, + { + "name": "evidence_drift", + "command": [ + "python3.13", + "-I", + "-B", + "/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py" + ], + "exit": 1, + "output": "Traceback (most recent call last):\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 257, in \n main()\n ~~~~^^\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 252, in main\n require(json.loads((ROOT/'EXPECTED.json').read_text()) == result, 'expected evidence mismatch')\n ~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 26, in require\n raise ValueError(message)\nValueError: expected evidence mismatch\n" + }, + { + "name": "refuse_existing_output", + "command": [ + "python3.13", + "-I", + "-B", + "/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py", + "--write", + "/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/existing.json" + ], + "exit": 1, + "output": "Traceback (most recent call last):\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 257, in \n main()\n ~~~~^^\n File \"/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/package/verify.py\", line 249, in main\n with args.write.open('x') as f:\n ~~~~~~~~~~~~~~~^^^^^\n File \"/opt/homebrew/Cellar/python@3.13/3.13.11_1/Frameworks/Python.framework/Versions/3.13/lib/python3.13/pathlib/_local.py\", line 537, in open\n return io.open(self, mode, buffering, encoding, errors, newline)\n ~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\nFileExistsError: [Errno 17] File exists: '/var/folders/g0/3_pt5b7911vc69wz960x6gq40000gn/T/hourly68 minimal path ifgfxn1h/existing.json'\n" + } + ] + }, + "proof_review": { + "precise_scope": "All L>=1, B>=0, 0 0 and len(coefficients) == n, "coefficient length") + require(direction in (-1, 1), "direction") + require(type(bits) is int and bits >= 0, "precision") + require(0 < beta <= 1, "margin") + require(all(a*a+b*b <= (1-beta)**2 for a, b in coefficients), "stability") + scale = 1 << bits + require(F(1, scale) <= beta / 2, "closure precision") + order = range(n) if direction == 1 else range(n-1, -1, -1) + + def update(value, coefficient, offset=(0, 0)): + a, b = coefficient + v, w = value + return (nearest(a*v-b*w+scale*offset[0]), + nearest(a*w+b*v+scale*offset[1])) + + z, power = (0, 0), (scale, 0) + for i in order: + z = update(z, coefficients[i], x[i]) + power = update(power, coefficients[i]) + sign = -1 if mutation != "plus_denominator" else 1 + den = (scale+sign*power[0], sign*power[1]) + den_norm = den[0]**2 + den[1]**2 + require(den_norm > 0, "closure denominator") + seed = (nearest(F(scale*(z[0]*den[0]+z[1]*den[1]), den_norm)), + nearest(F(scale*(z[1]*den[0]-z[0]*den[1]), den_norm))) + if mutation == "zero_seed": + seed = (0, 0) + output, previous = [(0, 0)]*n, seed + for i in order: + coefficient = coefficients[i if mutation != "constant_coefficient" else 0] + previous = update(previous, coefficient, x[i]) + output[i] = previous + return [(F(a, scale), F(b, scale)) for a, b in output] + + +def error_bound(B, beta, bits): + """The non-simplified universal bound from PROOF.md.""" + return F(1, 1 << bits) * (2*B/beta**2 + 3/beta + 1) diff --git a/research/cyclic-scan-sharp-bound/vendor/NOTICE.md b/research/cyclic-scan-sharp-bound/vendor/NOTICE.md new file mode 100644 index 0000000000..22fbeb8ef0 --- /dev/null +++ b/research/cyclic-scan-sharp-bound/vendor/NOTICE.md @@ -0,0 +1,24 @@ +# Attribution and scope + +Ali Romman's PR621, exact head 520774462b19b29af41554b74def31f9d01e997e, +supplies the prime-power kernel circle bound. Its kernel attribution is +OpenAI's exact-DFT preprint at adc7f1241b42e322a6451854ab7e4b4c146bf78a +(`loc:newton`, `loc:fourier-factorization`) and classical q-Pochhammer/Newton +factorization, including Alexey Kuznetsov's cited q-calculus work. This +package does not claim to originate those constructions. + +Pan's approximate-factorization theorem is used in the form stated by +Mehlhorn, Sagraloff and Wang; Brent's work supplies fast elementary-function +evaluation. Causal/anticausal recursive filtering is standard, as illustrated +by the Briand-Davy source listed in SOURCE.json. Sources are cited for these +specific tools, not as proofs of this package's complete claim. + +The scan rounding analysis, explicit budget composition, interval and +factor certificates, routed reference execution, and packaging were +prepared with substantial OpenAI Codex assistance for chafreaky. The dense +matrix solver and algebraic coefficient crosschecks provide independent +computational paths, not human peer review or formal Lean verification. + +No new multiplication exponent, literature priority, maintainer acceptance, +physical tape execution, or production performance result is claimed. +The code and documents use the repository's Apache-2.0 license. diff --git a/research/cyclic-scan-sharp-bound/vendor/check_budgets.py b/research/cyclic-scan-sharp-bound/vendor/check_budgets.py new file mode 100644 index 0000000000..ff21b6e190 --- /dev/null +++ b/research/cyclic-scan-sharp-bound/vendor/check_budgets.py @@ -0,0 +1,59 @@ +"""Exact integer checks of the explicit uniform precision budgets.""" +from pathlib import Path +import argparse,json,math + +def need(ok,msg): + if not ok:raise ValueError(msg) +def ceil_log2(n): + need(type(n)==int and n>0,'positive integer log argument') + return (n-1).bit_length() +def prime_power_data(r): + need(type(r)==int and r>=3,'radix') + prime=next((d for d in range(2,math.isqrt(r)+1) if r%d==0),r) + q=r + while q%prime==0:q//=prime + need(q==1,'prime power') + return r-r//prime +def budgets(r,p): + D=prime_power_data(r);need(type(p)==int and p>=1,'target bits') + n=r-1;ell=ceil_log2(r);d=1+2*ell+r*n*D;g=d+ell+2*r+2;a=ell+n*g + e=max(d+1,a+1,p+2*a+3);b=max(n*ell,e+r+3) + w=max(2*g+ceil_log2(14*n),p+(n+2)*g+ceil_log2(28*n*(r+1))) + return dict(radix=r,target_bits=p,n=n,D=D,ell=ell,d=d,g=g,a=a,e=e,b=b,w=w) +def check(x): + r,p,n,D,ell,d,g,a,e,b,w=(x[k] for k in ('radix','target_bits','n','D','ell','d','g','a','e','b','w')) + need(D==prime_power_data(r) and n==r-1 and ell==ceil_log2(r),'kernel parameters') + need(d>=1+2*ell+r*(n*D-1),'delta bound exponent') + need(g>=d+ell+2*r and g>=2,'root gap exponent') + need(a>=ell+n*g,'inverse prefix bound') + need(e>=d+1 and e>=a+1 and e>=p+2*a+3,'factor perturbation budget') + need(b>=n*ell and b>=e+r+3,'approximate factorization precision') + need(w-2*g>=ceil_log2(14*n),'rounded recurrence growth budget') + gap=w-p-(n+2)*g-1 + need(gap>=0 and 14*n*(r+1)<=1< 0, 'empty cycle') + need(direction in (-1, 1), 'invalid direction') + need(0 < beta <= 1, 'invalid margin') + need(norm2(b) <= (1-beta)**2, 'uncertified stable factor') + u = F(1, 1 << bits) + need(u <= beta*beta/2, 'precision too low for closure') + order = list(range(len(x))) if direction == 1 else list(reversed(range(len(x)))) + z, power = ZERO, ONE + for i in order: + z = rounded(add(x[i], mul(b, z)), bits) + power = rounded(mul(b, power), bits) + den = sub(ONE, power) + need(norm2(den) >= (beta/2)**2, 'closure lower bound failed') + seed = rounded(div(z, den), bits) + y, previous = [ZERO]*len(x), seed + for i in order: + previous = rounded(add(x[i], mul(b, previous)), bits) + y[i] = previous + return y + +def dense_solve(x, b, direction): + """Independent exact Gaussian elimination on the cyclic matrix.""" + n = len(x) + a = [[ZERO for _ in range(n)] + [x[i]] for i in range(n)] + for i in range(n): + a[i][i] = add(a[i][i], ONE) + j = (i-direction) % n + a[i][j] = sub(a[i][j], b) + for k in range(n): + pivot = next(i for i in range(k, n) if a[i][k] != ZERO) + a[k], a[pivot] = a[pivot], a[k] + p = a[k][k] + a[k] = [div(v, p) for v in a[k]] + for i in range(n): + if i == k: continue + p = a[i][k] + a[i] = [sub(v, mul(p, w)) for v, w in zip(a[i], a[k])] + return [row[-1] for row in a] + +def residual_bound(x, y, b, beta, direction, bound): + # Neumann inverse norm <= 1/beta converts residual to forward error. + for i in range(len(x)): + residual = sub(sub(y[i], mul(b, y[(i-direction) % len(x)])), x[i]) + need(norm2(residual) <= (beta*bound)**2, 'residual error bound failed') + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument('--out', type=Path, required=True) + args = parser.parse_args() + cases = [] + for beta in (F(1, 2), F(1, 16), F(1, 4096)): + for unit in (ONE, (F(-1), F(0)), (F(3, 5), F(4, 5))): + b = mul((1-beta, F(0)), unit) + for n in (1, 2, 3, 8, 17, 64): + x = [(F((i*7)%13-6, 8), F((i*11)%17-8, 16)) for i in range(n)] + B = max(abs(v[0])+abs(v[1]) for v in x) + for direction in (-1, 1): + bits = 64 + bound = 6*(B+1)*F(1, 1 << bits)/beta**3 + y = scan(x, b, beta, bits, direction) + residual_bound(x, y, b, beta, direction, bound) + if n <= 8: + exact = dense_solve(x, b, direction) + need(all(norm2(sub(v,w)) <= bound**2 for v,w in zip(y,exact)), 'dense comparison failed') + cases.append(dict(beta=str(beta), b=[str(v) for v in b], length=n, + direction=direction, bits=bits, forward_error_bound=str(bound), + dense_crosscheck=n<=8, residual_check=True)) + controls = [] + for name, x, b, beta, bits, direction in ( + ('unstable', [ONE], ONE, F(1,2), 64, 1), + ('underprecision', [ONE], (F(15,16),F(0)), F(1,16), 3, 1), + ('empty', [], ZERO, F(1,2), 64, 1), + ('bad_direction', [ONE], ZERO, F(1,2), 64, 0), + ): + try: scan(x,b,beta,bits,direction) + except ValueError as error: controls.append(dict(name=name,rejected=str(error))) + else: raise ValueError('negative control accepted: '+name) + x = [ONE, ZERO, ONE] + b, beta = (F(1,2),F(0)), F(1,2) + y = scan(x,b,beta,64,1) + y[0] = add(y[0],ONE) + try: residual_bound(x,y,b,beta,1,F(1,1<<40)) + except ValueError as error: controls.append(dict(name='corrupt_output',rejected=str(error))) + else: raise ValueError('corrupt output accepted') + args.out.write_text(json.dumps(dict(status='passed',cases=cases,controls=controls),indent=2)+'\n') + print(json.dumps(dict(status='passed',cases=len(cases),dense_cases=sum(c['dense_crosscheck'] for c in cases),controls=len(controls)))) + +if __name__ == '__main__': main() diff --git a/research/cyclic-scan-sharp-bound/vendor/weighted.py b/research/cyclic-scan-sharp-bound/vendor/weighted.py new file mode 100644 index 0000000000..da90ee6e83 --- /dev/null +++ b/research/cyclic-scan-sharp-bound/vendor/weighted.py @@ -0,0 +1,118 @@ +"""Exact certificates for varying unit phases on a cyclic neighbor shift.""" +from pathlib import Path +from fractions import Fraction as F +import argparse, importlib.util, json +P=Path(__file__).resolve().parent +s=importlib.util.spec_from_file_location('kernel',P/'check.py');k=importlib.util.module_from_spec(s);s.loader.exec_module(k) +c=k.c + +def scan(x,bs,beta,bits,direction): + c.need(len(x)>0 and len(bs)==len(x),'coefficient length') + c.need(direction in (-1,1),'direction') + c.need(0=(beta/2)**2,'closure denominator') + previous=c.rounded(c.div(z,den),bits);y=[c.ZERO]*len(x) + for i in order: + previous=c.rounded(c.add(x[i],c.mul(bs[i],previous)),bits);y[i]=previous + return y + +def apply(x,state,phases,direction,wrong_phase=False): + n=len(x);c.need(n>0 and len(phases)==n,'phase length') + c.need(all(c.norm2(t)==1 for t in phases),'unit phases') + y=[(v[0]*state['radix'],v[1]*state['radix'])for v in x] + for inside,b in state['factors']: + if inside: + inv=[(phases[(i+direction)%n][0],-phases[(i+direction)%n][1])for i in range(n)] + if wrong_phase:inv=[(v[0],-v[1])for v in phases] + bs=[c.mul(b,t)for t in inv] + y=scan(y,bs,state['beta'],state['bits'],-direction) + y=[c.rounded(c.mul(inv[i],y[(i+direction)%n]),state['bits'])for i in range(n)] + else: + y=scan(y,[c.mul(b,t)for t in phases],state['beta'],state['bits'],direction) + y=[c.rounded(c.neg(c.mul(b,v)),state['bits'])for v in y] + return y + +def matrix(q,phases,direction): + # Independently expand W^j entry-by-entry, including wrap-around phases. + n=len(phases);a=[[c.ZERO for _ in range(n)]for _ in range(n)] + for i in range(n): + weight=c.ONE + for j,v in enumerate(q): + col=(i-direction*j)%n;a[i][col]=c.add(a[i][col],c.mul(v,weight)) + weight=c.mul(weight,phases[col]) + return a + +def residual(x,y,a,bound): + for i,row in enumerate(a): + v=c.ZERO + for z,t in zip(row,y):v=c.add(v,c.mul(z,t)) + c.need(c.norm2(c.sub(v,x[i]))<=bound**2,'expanded weighted residual') + +def solve(a,x): + a=[row[:]+[v]for row,v in zip(a,x)];n=len(a) + for j in range(n): + i=next(i for i in range(j,n)if a[i][j]!=c.ZERO);a[i],a[j]=a[j],a[i] + z=a[j][j];a[j]=[c.div(v,z)for v in a[j]] + for i in range(n): + if i!=j: + z=a[i][j];a[i]=[c.sub(v,c.mul(z,t))for v,t in zip(a[i],a[j])] + return [row[-1]for row in a] + +def main(): + ap=argparse.ArgumentParser();ap.add_argument('--out',type=Path,required=True);a=ap.parse_args() + units=[c.ONE,(F(0),F(1)),(F(3,5),F(4,5)),(F(-5,13),F(12,13))] + cases=[];certs=[];generic=[];controls=[] + for beta in (F(1,2),F(1,16),F(1,4096)): + for n in (1,2,3,8,17,64): + bs=[c.mul((1-beta,F(0)),units[(i*i+3*i+1)%4])for i in range(n)] + # Alternate magnitudes as well as phases. + bs=[(z[0]/(1+i%2),z[1]/(1+i%2))for i,z in enumerate(bs)] + x=[(F(i%5-2,4),F(i%7-3,8))for i in range(n)] + B=max(abs(z[0])+abs(z[1])for z in x);bound=6*(B+1)*F(1,1<<64)/beta**3 + for d in (-1,1): + y=scan(x,bs,beta,64,d) + m=[[c.ZERO for _ in range(n)]for _ in range(n)] + for i in range(n):m[i][i]=c.ONE;m[i][(i-d)%n]=c.sub(m[i][(i-d)%n],bs[i]) + residual(x,y,m,beta*bound) + if n<=8:c.need(all(c.norm2(c.sub(v,t))<=bound**2 for v,t in zip(y,solve(m,x))),'generic dense') + generic.append(dict(beta=str(beta),length=n,direction=d,dense=n<=8)) + sigma,nu,phase_cert=k.irrational_gauge() + for row in json.loads((P/'FACTORS.json').read_text()): + if row['radix'] not in (4,9,16):continue + st=k.certify(row);A=F(st['certificate']['inverse_norm_bound']) + eps=F(st['certificate']['coefficient_error_bound'])+sum(j*(abs(v[0])+abs(v[1]))for j,v in enumerate(st['q']))*nu + c.need(A*eps<1 and A*A*eps/(1-A*eps)+st['budgets'][-1]<=F(1,1<= u/(4*beta**2), 'sharp lower bound') + require(y == fixed.run([ONE], [(b, F(0))], beta, bits)[0], 'sharp integer replay') + sharp.append(dict(k=k, beta=str(beta), error_over_u_beta_squared=str(error*beta**2/u))) + + # General cascade inequality, checked with exact recurrence on moderate sizes. + recurrence_cases = 0 + for n in range(1, 33): + r = n+1 + for g in (1, 2, 4, 8): + beta = F(1, 1 << g) + for p in (1, 32, 128): + w = max(g+ceil_log2(10*n), p+(n+1)*g+ceil_log2(20*n*(r+1))) + u = F(1, 1 << w) + E, amplification = F(0), F(r) + for _ in range(n): + E = E/beta+5*(amplification+E+1)*u/beta**2 + amplification /= beta + require(E <= 10*n*(r+1)*u/beta**(n+1) <= F(1, 1 << (p+1)), 'cascade bound') + recurrence_cases += 1 + budgets = [] + for r in range(3, 257): + try: + old_budgets.prime_power_data(r) + except ValueError: + continue + for p in (1, 128, 1024, 1048576): + old = old_budgets.budgets(r, p) + old_budgets.check(old) + n, g = r-1, old['g'] + new = max(g+ceil_log2(10*n), p+(n+1)*g+ceil_log2(20*n*(r+1))) + require(new <= old['w']-g, 'execution guard saving') + require(new-g >= ceil_log2(10*n), 'new growth guard') + require(new-p-(n+1)*g >= ceil_log2(20*n*(r+1)), 'new target guard') + budgets.append(dict(radix=r, p=p, g=g, old=old['w'], new=new, saved=old['w']-new)) + + controls = [] + def rejected(name, fn): + try: + fn() + except (ValueError, ZeroDivisionError): + controls.append(name) + else: + raise ValueError('control accepted: '+name) + x, bs, beta, bits = [ONE, ZERO, (F(1, 2), F(1, 4))], [(F(1, 2), F(0)), (F(0), F(1, 2)), (F(-1, 4), F(0))], F(1, 2), 24 + exact = exact_inverse(x, bs, 1) + bound = fixed.error_bound(F(1), beta, bits) + for mutation in ('zero_seed', 'plus_denominator', 'constant_coefficient'): + def test(mutation=mutation): + y = fixed.run(x, bs, beta, bits, 1, mutation) + require(all(norm2(minus(a, b)) <= bound**2 for a, b in zip(y, exact)), mutation) + rejected(mutation, test) + rejected('wrong_orientation', lambda: require(all(norm2(minus(a,b)) <= bound**2 for a,b in zip(fixed.run(x,bs,beta,bits,-1),exact)), 'orientation')) + rejected('unstable_coefficient', lambda: fixed.run([ONE], [ONE], beta, bits)) + rejected('below_closure_guard', lambda: fixed.run([ONE], [ZERO], F(1, 16), 4)) + k=8; u=F(1,1<<(3*k)); a=F(1,1<