Found some Scheme code that I wrote when I was reading Concrete Mathematics
The Stern–Brocot tree was discovered independently by Moritz Stern (1858) and Achille Brocot (1861). Stern was a German number theorist; Brocot was a French clockmaker who used the Stern–Brocot tree to design systems of gears with a gear ratio close to some desired value by finding a ratio of smooth numbers near that value.
(text below from Concrete Mathematics)
The idea is to start with the two fractions (0/1, 1/0) and then to repeat the following opeation as many times as desired:
Insert (m + m')/(n + n') between two adjacent fractions m/n and m'/n'
The new fraction (m + m')/(n + n') is called the _mediant_ of m/n and m'/n'. For example, the first step gives us one new entry between 0/1 and 1/0:
0/1, 1/1, 1/0
and the next step gives two more
0/1, 1/2, 1/1, 2/1, 1/0
The next gives four more,
0/1, 1/3, 1/2, 2/3, 1/1, 3/2, 2/1, 3/1, 1/0
The entire array can be regarded as an infinite binary tree structure whose top levels look like this:
1 1 0
— - - - - — - - - - —
0 1 1
/ \
/ \
/ \
/ \
1 / \ 2
— —
2 1
/ \ / \
/ \ / \
/ \ / \
1 2 3 3
— — — —
3 3 2 1
/ \ / \ / \ / \
1 2 3 3 4 5 5 4
— — — — — — — —
4 5 5 4 3 3 2 1We can, in fact, regard the Stern-Brocot tree as a number system for representing rational numbers, because each positive, reduced fraction occurs exactly once. Let's use the letters L and R to stand for going down to the left or right branch as we proceed from the root of the tree to a particular fraction; then a string of L's and R's uniquely identifies a place in the tree. For example, LRRL means that we go left from 1/1 down to 1/2, then right to 2/3, then right to 3/4, then left to 5/7. The fraction 1/1 corresponds to the empty string, let's agree to call it I.
Graham, Knuth, Patashnik, Concrete Mathematics, Addison Wesley, 2nd ed. pp 115-123
We define
S = a string of L's and R's
f(S) = fraction corresponding to S
M(S) = ((n m) (n' m')) = a 2x2 matrix that holds the four quantities involved in the ancestral fractions m/n and m'/n'
M(I) = ((1 0) (0 1))
M(SL) = M(S) . ((1 0) (1 1))
M(SR) = M(S) . ((1 1) (0 1))
Therefore, if we define L and R as 2x2 matrices,
L = ((1 0) (1 1))
R = ((1 1) (0 1))
In Scheme:
(define I '((1 0) (0 1)))
(define R '((1 0) (1 1)))
(define L '((1 1) (0 1)))While we're at it let's define a few useful operations.
Since we're representing matrices as lists, we'll need matrix multiplication
(define (mat-mul m1 m2)
(map
(lambda (row)
(apply map (lambda column (apply + (map * row column))) m2))
m1))and conversion from the 2x2 matrix representation to rational number
(define (mat->rat mat)
`(/ ,(+ (caadr mat) (cadadr mat)) ,(+ (caar mat) (cadar mat))))Irrational numbers don't appear in the Stern-Brocot tree, but all the rational numbers that are ``close'' to them do. We can find the infinite representation of an irrational number α with this procedure:
if α < 1 then (output(L); α := α / (1 - α))
else (output(R); α := α - 1)In Scheme
(define (brocot-string num steps)
(letrec ((helper (lambda (num acc steps)
(cond ((equal? steps 0) (reverse acc))
((< num 1) (helper (/ num (- 1 num)) (cons 'L acc) (- steps 1)))
(else (helper (- num 1) (cons 'R acc) (- steps 1)))))))
(helper num '() steps)))Now we can look for a regular pattern
> (brocot-string 2.7182818284590452353599104282424311632484 32)
'(R R L R R L R L L L L R L R R R R R R L R L L L L L L L L R L R)In the Stern-Brocot system e = R²LR²LRL⁴RLR⁶LRL⁸RLR¹⁰LRL¹²....
We'll use infinite streams to represent this string. We start with the sequence of bases RLRLRL...
(require srfi/41)
(define bases (stream-cons R (stream-cons L bases)))The sequence of exponents is (2,1,1,2,1,1,4,1,1,6,1,1,8,1,1,10...) You may recognize Euler's continued fraction https://oeis.org/A003417
(define exponents
(stream-map
(lambda (x)
(cond ((equal? x 1) 2)
((equal? (modulo x 3) 0) (- x (/ x 3)))
(else 1)))
(stream-from 1)))Now we just have to zip the two streams
(define exp-pairs (stream-zip bases exponents))and we build the infinite sequence of rational approximations of e
(define-stream (exponentiate a b)
(if (equal? b 0)
empty-stream
(stream-cons a (exponentiate a (- b 1)))))
(define-stream (make-e-seq pairs)
(let ((a (cadr (stream-car pairs)))
(b (car (stream-car pairs))))
(if (equal? a 1)
(stream-cons b (make-e-seq (stream-cdr pairs)))
(stream-append (exponentiate b a) (make-e-seq (stream-cdr pairs))))))
(define e-seq (make-e-seq exp-pairs))We define a function to calculate the nth rational approximation of e in our Stern-Brocot representation like this
(define (e ref)
(mat->rat (stream-ref (stream-scan mat-mul I e-seq) ref)))now let's test our code
> (require math/bigfloat)
> (bf-precision 100000)
> (time (bf (eval (e 1000000))))
cpu time: 3258 real time: 3260 gc time: 1711
(time (bf (eval (e 1000000))))
(bf #e2.7182818284590452353...)Yields a result correct to 6331 digits. While far from being the most efficient way to approximate e, that was fun!
What's the fraction 1/0 supposed to be?
>>4
According to Donald Knuth
I guess 1/0 is infinity "in lowest terms"
0/1 and 1/0 do not belong to the Stern-Brocot language, they're only the ancestral fractions of the empty element 1/1. Scaffolding.
This is not your personal blog.
>>6
By far more /prog/ than the last 20 posts.
In a similar spirit but not as advanced, including a small challenge at the end:
https://textboard.org/prog/34#t34p165
>>7
Who cares, it is a wall of text and should be posted elsewhere.
>>6,8
There were discussions here about Concrete Mathematics in the past, I thought someone might be interested in a Scheme implementation of cool clockmaker technologies.
But ok, I should have turned this into some kind of challenge like "let's post original ways to approximate e".
>>8
Missed that. Thanks.
>>1
erratum
0 1 1
— - - - - — - - - - —
1 1 0
/ \
/ \
/ \
/ \
1 / \ 2
— —
2 1
/ \ / \
/ \ / \
/ \ / \
1 2 3 3
— — — —
3 3 2 1>>9
You mean code walls of text that belong on /prog/ or refuting >>1-kun for having shit code walls of text that could be simpler and cleaner elegant code five liners at most?
Some other nice patterns in the Stern-Brocot tree:
> (brocot-string (sqrt 2) 32)
'(R L L R R L L R R L L R R L L R R L L R R L L R R L L R R L L R)
> (brocot-string (sqrt 3) 32)
'(R L R R L R R L R R L R R L R R L R R L R R L R R L R R L R R L)
> (define phi (/ (+ 1 (sqrt 5)) 2))
> (brocot-string phi 32)
'(R L R L R L R L R L R L R L R L R L R L R L R L R L R L R L R L)>>13
What is the characterization of the subset of real numbers whose SB string expansion can be recognized by a FSM?
class SBtree:
@ staticmethod
def string (x):
while x > 0:
if x < 1:
yield 'L'
x = x / (1 - x)
else:
yield 'R'
x = x - 1
@ staticmethod
def nstring (x, n):
return ''.join (itertools.islice (SBtree.string (x), n))
>>> sb.SBtree.nstring (math.sqrt (2), 32)
'RLLRRLLRRLLRRLLRRLLRRLLRRLLRRLLR'
>>> sb.SBtree.nstring (math.sqrt (3), 32)
'RLRRLRRLRRLRRLRRLRRLRRLRRLRRLRRL'
>>> sb.SBtree.nstring (math.sqrt (5), 32)
'RRLLLLRRRRLLLLRRRRLLLLRRRRLLLLRR'
>>> sb.SBtree.nstring (math.sqrt (7), 32)
'RRLRLRRRRLRLRRRRLRLRRRRLRLRRRRLR'
>>> sb.SBtree.nstring ((1 + math.sqrt (5)) / 2, 32)
'RLRLRLRLRLRLRLRLRLRLRLRLRLRLRLRL'
>>> sb.SBtree.nstring (math.e, 32)
'RRLRRLRLLLLRLRRRRRRLRLLLLLLLLRLR'
>>> sb.SBtree.nstring (math.pi, 32)
'RRRLLLLLLLRRRRRRRRRRRRRRRLRRRRRR'>>14
I don't know if they have a name. The pattern of irrational number approximations in the Stern-Brocot tree is similar to the continued fraction, a more common notation than the Stern-Brocot tree sequences.
sqrt(2) = [1;2,2,2,2,2,...] (algebraic)
phi = [1;1,1,1,1,1,1,...] (algebraic)
e = [2;1,2,1,1,4,1,1,6,1,1,8,...] (transcendental)
cube root of 2 = [1;3,1,5,1,1,4,1,1,8,1,14,1,10,2,1,4,12,2,3,2...] (algebraic, non periodic CF)
π = [3;7,15,1,292,1,1,1,2,1,3,1,...] (transcendental, no known pattern)
My guess is that there's an infinite but countable number of regular patterns you can use to generate an irrational number. Since the set of irrational numbers is uncountable almost all continued fractions of irrational numbers have no pattern.
I do appreciate this and would like to see more of these. But I would change the format a bit.
You could write the problem in the OP as a motivator, and a link to a literary .scm file with the solution.
>>16
The guy doesn't care, he just writes down his stream of consciousness, and dumps it on the board.
maybe verbose but I would rather see this than recycled memes complaining about rust
rewrite Stern-Brocot tree in rust(lacks example)
http://rosettacode.org/wiki/Stern-Brocot_sequence
>>19
And then try to do that once you have the sequence. This rational is an approximation of e.
> (time (e 1000000))
cpu time: 2995 real time: 2997 gc time: 1174
'(/
340017086135953955237517024368071023319667336443292804743952295811218300447359859240275563254009942677431489707597147675747035681413226512627011095930560268265456381178976392711767087225557220527288583530246999554026111396756681847411034069691700458364305110849779785985479469458249422285136155670021333443718590111827401706503305002386802304743268641986242675018801066689637090293781703589884119797742595395970458469615457316464652161835058300625085365074839652205963909281728057007645884413972470559421095590641084923213749096077911891464157447244982848623998799902248659508491205347662407542487598477506176854976430107224759805723469500480207518560065060801139517286177001208163592054763542622442833088837014707870098680498790101406422989090225950831187269091905076078387996734933063315154563648844011105458348826524339618765836090847570234445324485311773505775960915934136456454684723065101560864875759072426638458977878476314281969987155314153269630993616187041722772926774199317277275408014196669434574127395283760934088304504454441826341472487028029542307844361121860169948845000863967607247499379959564922815018909507873427673764317550105676071503852510957626027913732521804089379035997601445017705811769908093393745107753664844766235361183695170527428489809363428463651419304155765399729238548589316792609372558848237780044857470539820865022763311896953342264917235482949907211105305961474199254134498023181661330808198705378200836895070871333229848624381644393862053739146946714031748697068394959781227898824577572717091477283125470687771132502660646193012227555089150091316446335434009362342321577385694092293349561432035238800308129504247688449971272567391760416476691624159595214779684041072400751170742699454671514147869612822126793025761288002185791555487625331998621292869836022457378328419035361513580174231044804345605326811465916238059290186883564724761265694886919748527987916817620628526399573566506494861425004396642032056780484911383385189989513516563131778797299221345966468968956122130938820067631844304729073653295278994294270111849584185735076742525508006069816526335483941312168466769941517563643711961854946980989175852514399703705680462519538562629390325716019370164673384882170485495116424556766549082251149470708965255450905049612453587859276129878129571158069679466773167708728881129438071391603459651281058662550128959761153549617728975030404264330537814833057633963701181617935774238987394896893837492212211374342977762904418298881285931234958814978103651305047589748736431989381870422496158486260942264340691466697262368432164893697234466153250046703893250452429139844480273241829106099238215033421990835160228005567092855709150858705097715014615303931593755514322351810007680897794492648609028769501447879692294277237362011716820403857042445580674191735659531498184966721518718160689953505802281847373678940279750467265390416599246074195127826056496476091242971985446723655542587173824287311438381591741110560862978217691099053628393417793105301826047059921235216516021504425461929120276377753800102305173326425003579123379528380303069690154701626211729549743364833186844412203484432054917804678030588228660410069011019359882579002
125085295636436909683809909842556856885178429583466250503633118421625758576495571982549603221529371824249449947226736257956572986278847452327489918555893035992800158136526570747790366607976473974086490509833266140241624411002520378548860256844094890877496352409672412181476991240644660422957359119247812657697715895808101749155909929315560672524600833028834500581928675977242236292980202330265409027931557849262071183680232807134194663977226617526758164092581075236564819647024212853553521938166754805064971997583680157504280085693771587742124565102644454392184568192455131506363793673061697667809133470803909462294541130705110877585153688820760830173021946921156853394209265699621359388107424247279451921507178896866252151258182423463805654515203717737130171743215112991658603794371604988588191145881380551692440558754609581883970182436373351721827090765483767025558740593461170481647907645229692750706772617083777913648346892911751008260597754177442242884897528540552656654142051828032526143451856483844303282993772253401076047130518814970561671090962371290205074655269040331486455093397311641724716513914711853050343863410602433828186055109883845169168186241477852812819511630717226031782266261200943940615962914287809900772211847247623544281170759686802721987002900335098524747135307554944852406985163571494649355659613567838732869723194876669393386105031678437801228285473666870189915779831674530812838317813362749429334047728206554077502413410488307816529862238758762434688416460199490950423539362723027760398951127569670882323069657999986406666638536904433413611732504237085435613428187295322449920616488897190353761498046868750584302985283042516291072517661563543124917633435349040950215732545008269904452842072489224708167231629873428753864275250495686428820211229340636698055009463068010066247920722428973426571945857005028864658276369838400927087527667218274517604615390540969915685158805164446004960391637638952325069774576683392245223799253721839340049299692557960655127171880922366931873834678505507980929392029365039716768900191201134404082041954658457421412098953281935182415736165076549914653156936288250673567627314590360080813173455844097539997651557001995149793696142371449923446980427725864245921082517781640485031251339553931144656900391589822516702794312100705930445806007992360069962039035567107884172990813990752000736096406721872831171253639698564151167421533401665295128817139455342798715270038691942997621086685910211094868754454516817862989819051640051770874497216852913390188488061385723926005176343183394443750361968944989563225128221788985578992361222409045215774195773008647206524738774180554917621147514549299392328337681997169072439473578702176284273236998017086945818255227347039675511648057103771074013479256786923260298195015744892867688675269801389398308311777208591853723170164583690632297703927756753727378438541760265853399635262151837406060792526386670515259891008695811004879840637123479657107736250722757524021127684175490377871823327703940937023628221667804570593378320396524365465248186987504652311179460359354823215116827906074424119159403212171240501026466951623914101612248957295893464404775833237162234624918433172999)Also note that the if the stream exponents is Euler's continued fraction, the stream bases (RLRLRLRLRLRLR...) is nothing else than the golden ratio in our Stern-Brocot tree number system. So, if you want a rational approximation of Φ instead, just do this
> (define (Φ ref)
(mat->rat (stream-ref (stream-scan mat-mul I bases) ref)))
> (time (Φ 10000))
cpu time: 33 real time: 34 gc time: 15
'(/
88083137989997064605355872998857923445691333015376030932812485815888664307789011385238647061572694566755888008658862476758094375234981509702215595106015601812940878487465890539696395631360292400123725490667987980947195761919733084221263262792135552511961663188744083262743015393903228035182529922900769207624088879893951554938584166812233127685528968882435827903110743620870056104022290494963321073406865860606579792362403866826411642270661211435590340090149458419810817251120025713501918959350654895682804718752319215892119222907223279849851227166387954139546662644064653804466345416102543306712688251378793506564112970620367672131344559199027717813404940431009754143637417645359401155245658646088296578097547699141284451819782703782878668237441026255023475279003880007450550868002409533068098127495095667313120369142331519140185017719214501847645741030739351025342932514280625453085775191996236343792432215700850773568257988920265539647922172315902209901079830195949058505943508013044450503826167880993094540503572266189964694973263576375908606977788395730196227274629745722872833622300472769312273603346624292690875697438264265712313123637644491367875538847442013130532147345613099333195400845560466085176375175045485046787815133225349388996334014329318304865656815129208586686515835880811316065788759195646547703631454040090435955879604123186007481842117640574158367996845627012099571008761776991075470991386301988104753915798231741447012236434261594666985397841758348337030914623617101746431922708522824868155612811426016775968762121429282582582088871795463467796927317452368633552346819405423359738696980252707545944266042764236577381721803749442538053900196250284054406347238606575093877669323501452512412179883698552204038865069179867773579705703841178650618818357366165649529547898801198617541432893443650952033983923542592952070864044249738338089778163986683069566736505126466886304227253105034231716761535350441178724210841830855527586882822093246545813120624113290391593897765219320931179697869997243770533719319530526369830529543842405655495229382251039116426750156771132964376
54438373113565281338734260993750380135389184554695967026247715841208582865622349017083051547938960541173822675978026317384359584751116241439174702642959169925586334117906063048089793531476108466259072759367899150677960088306597966641965824937721800381441158841042480997984696487375337180028163763317781927941101369262750979509800713596718023814710669912644214775254478587674568963808002962265133111359929762726679441400101575800043510777465935805362502461707918059226414679005690752321895868142367849593880756423483754386342639635970733756260098962462668746112041739819404875062443709868654315626847186195620146126642232711815040367018825205314845875817193533529827837800351902529239517836689467661917953884712441028463935449484614450778762529520961887597272889220768537396475869543159172434537193611263743926337313005896167248051737986306368115003088396749587102619524631352447499505204198305187168321623283859794627245919771454628218399695789223798912199431775469705216131081096559950638297261253848242007897109054754028438149611930465061866170122983288964352733750792786069444761853525144421077928045979904561298129423809156055033032338919609162236698759922782923191896688017718575555520994653320128446502371153715141749290913104897203455577507196645425232862022019506091483585223882711016708433051169942115775151255510251655931888164048344129557038825477521111577395780115868397072602565614824956460538700280331311861485399805397031555727529693399586079850381581446276433858828529535803424850845426446471681531001533180479567436396815653326152509571127480411928196022148849148284389124178520174507305538928717857923509417743383331506898239354421988805429332440371194867215543576548565499134519271098919802665184564927827827212957649240235507595558205647569365394873317659000206373126570643509709482649710038733517477713403319028105575667931789470024118803094604034362953471997461392274791549730356412633074230824051999996101549784667340458326852960388301120765629245998136251652347093963049734046445106365304163630823669242257761468288461791843224793434406079917883360676846711185597501)This has very little to do with the code on the rosettacode's page, it's something fun from Concrete Mathematics
cpu time: 33 real time: 34 gc time: 15
target 1.618033988749895
result 1.618033988749895
88083137989997064605355872998857923445691333015376030932812485815888664307789011385238647061572694566755888008658862476758094375234981509702215595106015601812940878487465890539696395631360292400123725490667987980947195761919733084221263262792135552511961663188744083262743015393903228035182529922900769207624088879893951554938584166812233127685528968882435827903110743620870056104022290494963321073406865860606579792362403866826411642270661211435590340090149458419810817251120025713501918959350654895682804718752319215892119222907223279849851227166387954139546662644064653804466345416102543306712688251378793506564112970620367672131344559199027717813404940431009754143637417645359401155245658646088296578097547699141284451819782703782878668237441026255023475279003880007450550868002409533068098127495095667313120369142331519140185017719214501847645741030739351025342932514280625453085775191996236343792432215700850773568257988920265539647922172315902209901079830195949058505943508013044450503826167880993094540503572266189964694973263576375908606977788395730196227274629745722872833622300472769312273603346624292690875697438264265712313123637644491367875538847442013130532147345613099333195400845560466085176375175045485046787815133225349388996334014329318304865656815129208586686515835880811316065788759195646547703631454040090435955879604123186007481842117640574158367996845627012099571008761776991075470991386301988104753915798231741447012236434261594666985397841758348337030914623617101746431922708522824868155612811426016775968762121429282582582088871795463467796927317452368633552346819405423359738696980252707545944266042764236577381721803749442538053900196250284054406347238606575093877669323501452512412179883698552204038865069179867773579705703841178650618818357366165649529547898801198617541432893443650952033983923542592952070864044249738338089778163986683069566736505126466886304227253105034231716761535350441178724210841830855527586882822093246545813120624113290391593897765219320931179697869997243770533719319530526369830529543842405655495229382251039116426750156771132964376
54438373113565281338734260993750380135389184554695967026247715841208582865622349017083051547938960541173822675978026317384359584751116241439174702642959169925586334117906063048089793531476108466259072759367899150677960088306597966641965824937721800381441158841042480997984696487375337180028163763317781927941101369262750979509800713596718023814710669912644214775254478587674568963808002962265133111359929762726679441400101575800043510777465935805362502461707918059226414679005690752321895868142367849593880756423483754386342639635970733756260098962462668746112041739819404875062443709868654315626847186195620146126642232711815040367018825205314845875817193533529827837800351902529239517836689467661917953884712441028463935449484614450778762529520961887597272889220768537396475869543159172434537193611263743926337313005896167248051737986306368115003088396749587102619524631352447499505204198305187168321623283859794627245919771454628218399695789223798912199431775469705216131081096559950638297261253848242007897109054754028438149611930465061866170122983288964352733750792786069444761853525144421077928045979904561298129423809156055033032338919609162236698759922782923191896688017718575555520994653320128446502371153715141749290913104897203455577507196645425232862022019506091483585223882711016708433051169942115775151255510251655931888164048344129557038825477521111577395780115868397072602565614824956460538700280331311861485399805397031555727529693399586079850381581446276433858828529535803424850845426446471681531001533180479567436396815653326152509571127480411928196022148849148284389124178520174507305538928717857923509417743383331506898239354421988805429332440371194867215543576548565499134519271098919802665184564927827827212957649240235507595558205647569365394873317659000206373126570643509709482649710038733517477713403319028105575667931789470024118803094604034362953471997461392274791549730356412633074230824051999996101549784667340458326852960388301120765629245998136251652347093963049734046445106365304163630823669242257761468288461791843224793434406079917883360676846711185597501
stack 6
$ python3 -m timeit -s 'import flake.sbtree as mod' 'mod.work_speclimit (mod.Specs.Phi (), 10000)'
10000 loops, best of 3: 137 usec per looptarget 1.414213562373095
result 1.414213562373095
912372904077792820333022931514683905196095335332461243563566371630623725027756174840958798586146594892914095401037642377610054228482576727704110500231020662439134133236781589366415205524878535314722341540748126631282687462031097448656828645302597037940779004487605981247667669583729123395867583899548179259053671636594464946726564337037930617138893433544440429360350179319034433043179386253653890778056959564090099584131042792609171485551446220287193670064157175032113243945323799379335206189982084006721103924837454561272670000817253624301004165698525457461249517756759114499889703222149412160665301694111168795976530493998573585466679141412035628294794135220268925566327084096717726075716622819757810639450486297919988578046158833828525325469062759563222180775691512072690399043415592162621696366117965993334661822419744848270408206137218896579046774810149448397093592171548356413795629132503305148903487947415107183019346627904110216664835865582002041482574392065817955357317802858225690925207376993605542582914690980260418563086534219739082870274044629514790108045021964648588957147877553335242189622725230861151908928597780357076966509346517444159661044195373532127149501448465623900205359156173053520045633506497748488413010211676218187024245944636731897289083327716265252592579848721079948056891537057074054078404059256984291844030570801386300484362165471294405751570533709092619045671164240921912086888637832641806419666332796920984444963799864385789395922731763421593489233820681623305468082385959617126639107953866786715082681199375295298398342516439245368301352147472945010786588515500962019666031356507887402263555290987070386890083958938515638981902578172203283552384977890173501578909020338496992971143893003367573490861994803868623180409853340292013727017820739449406566091157326121807628137523647618653406547187947774090879311484132159187408769424091732434374817993103851978756820565200612073505393
645145067444270761356716481857445357641649553114979126977678963754439951082038984632678061418130657858763297544529821984444895464626862603691739802477466631503873767432694953865148667058442664481918025121653166100343954869841332731542457362834113318740074258352418751185500785071084474652030341872751114517691124799298946476356463130421256001538326419686685024658804199090849748412873447604535750572396245530424783088241899918727047006913682721495224443247600409827124799470185578822892780248099201063388167346497317474442071198019623551456189630759811865831617485719078716492727748680645636681156261142523348830419155615693078563602052074308536791767398023563062654117446029432097791614030894600205918745588000027749241499556416620236383481141542139618267239995908984466111486647077115808307902368887069463321555768022390997259037097147923161550796376448633639191058160377727679526847368225829716455421516498274839007505526603632830905156575689466421113905804893041617407988995957582361699110853657191118269861280913215596304247181727470222880786012786174601277208667732522083157085633502053869442681988756433014137552047292645249667623080960158029307041797974528005411231083666859195730472054478916752735859020924743323704976823318050126745409118196223466719223339522073412802366352813857415809629646773071321021322915809378632028416028413299468319362402530515502588589399067652778615784607050131186014530406266363654361196873709332572288106710280508411833502919921220462451339963101921144001361401160661249966470252485720782595383560914505654141717325381005598727773686382562943798678618507388289518297731421477986248887295505359290084507833070865662747168546322649930506385522948488952975319520205239729675192003571242819758609748037543641240137577371696785255726670121334640733501179715052678269583483097831741247328687771258281063372501649165518987809812008399308777781407031699396066250798983206988348120665
stack 6
$ python3 -m timeit -s 'import flake.sbtree as mod' 'mod.work_speclimit (mod.Specs.Sqrt2 (), 10000)'
10000 loops, best of 3: 121 usec per looptarget 1.732050807568877
result 1.732050807568877
37333050739547620284741346177443232768419828695576483192119227709011874720867558339004144048908608643954187575967122760016462436597149422351883158769394191023679913703570251600882068481387414359397619517424061973492935942203954445483788494431222578988497811330117180887243303058687612115480389195262103976386982801547226241736184765087358487101684839240435787719763032119054686805648308753341087388646696533367579871775911044407929199778666970896621939895063461817937810347316850879892599866435706231219754394830144519751453447850442566201884566425006242081130213811073324017618229313206662953896117018758439730158342409695682744685435055507681964742791636723618042926447421832733978027037356893003642400795083277917572826017744015949356157886964966927490914775661127518930358507018785066186644947170327972057946083883044850852277818698118365103260660825088738036827240142991301488107222873963750819769771904654194225287688794059905383024205725563175121133535944182786005199341008982313788920025582392443629294875181601118908847120206222270984968236498764086338625250980712437368193355490705121246779227656912224111729593933533272780650483740863985497980606931536674447177718482803533891496589890376318229724885453329711793624837676433843135324083030155160267058536032981571743325390625805043153062883213991856782364670780605229622859726138556007967627528794502411945065099156358124051734768863065914792552486457279065577226417382681617126063892480294139719789670858396610321459720941260501653332523601524041423566189024845550000941398132271344064174395659569621527186657227534096104268698881663551635501366925200893982635208764415188639236599527755755220572368830848307065718490436530312669439969677316707180302233189875768562448229382203076191884103142375500528208253188467810619697091413712152142298355067830147049981234594952906086026836206347827204014637887719424320921238973858066084698096953687711047
21554246894147775797650832578303319263028404773092017926791774965040179802055156425720667289665504525915863815507020606442822591122943147707530277952196331686389340468052522281951900764255429081274790186257286268230621724824017922276471739177768334695384953697871877837881155031465359055028538862532846056226236947786925864165196534881747615168645767076159981751248460569675652894545741239891833300198605315822337729504784372003441206061759695669624704595618478359792979354122596033333248015656636731973572897856628025650797008826972131622163470243837577196348969655965538229305259367543808617363175028855903945171232070605035695614607015184689952951586312421207546531341892706569271601923463595043791163292377687138580802118407026421147932224182625632590924856698714323013638768197141356547467582293140034400356391974440346504311412439829638421716251542120579618728973231768592667750073773590267419218843931830221497884407792589017184003540992191405403034918783995177085651898934538952701703076197715837649791991627458931196137767975686601141239200515090942654152115666835346260740768607950637521722151011656654374992897516330417865006506372510681289630018666381838777552937120944730239512806638205231951609646164304538015782410981802569959188173997031855775431753188170602372552800771539717805865719257089560105167973807987552658427273583942730449293553775055987686609702926083577493656009021484804414735012057096501093436794050801716970153395976595414469165021583793220009073307292114697925615393562085411284136700452395376891838361636625895532086290413315475775102695473088715260812363265582034603291944111937827403421961367979542615557894681375864831794298679205586154133397902722523733839564466051081942627041120700837220749704695868172599837356090872581186433799610358187908547841059944683830596399048437391260574434794830712097077693113271446484457383044766723608216668485987749656447980544044121356
stack 7
$ python3 -m timeit -s 'import flake.sbtree as mod' 'mod.work_speclimit (mod.Specs.Sqrt3 (), 10000)'
10000 loops, best of 3: 127 usec per looptarget 2.236067977499790
result 2.236067977499790
41332294493938545531741898189124475002595083290489742271982042410001179649590638352971327178662782342658715702623030919737457907830190691919815009012258506376641135023825853565815239172591664165125213495590102205194513046640807952860130344653843839739939014372203482062252430269007112899956699379555747989512246507615679658302303597205403785255823841691991082530008399436256706957923656044401564042746883350613913909297196677594960296197551329608870980667422294986691904386738935339267820859141352189405773882440494245555502768587027398841266564236644931168995594325256768307938601042996121095218241154575833662343043702732992390353103692853602670343958004075610587622455232284763903928028648731969351298672275652254240329774371782479484528772043935076738418647432348313216320003488034702089990883852569243503909749609181481316530786070722847681935414366840677535822771205760207869055158774687324451786213766720466147055533425266225306535307310450569765855050917117310080021623332198083686457831950357034237316279546694295104672547709923461351337599635696874489909962163561569998895704141643514246636570145413107594734738158346892744164182479359093299630818413585724011399825669588321391668100430142680865137048656972350481497356791937266235756440335503077080610088030645733619795947172505340600403178545449010069030779945613043875596021032181052174542895568049342840041819628616283091446792264365775194267951629158048818995367725670881481269134411330911655910252452531346900379867176605513744986701481901397038137684255542940091075132568064399281896930450391998725001
18484364030897371437300920331810666605620710023854620676063981218524652362414027825042685017628799482783457930163920239444169022341244320998258403142770892372144915632343436000992920792602361123265031091519899666493465827806254868440765822317542750159888888783785299054848253120329811978818108194558415103032630528570493513177986315310466471480532082979588479340324654453261594178812021489384865926085887335446285369050537228888493984139295507672130875197799909665706437635949989954884651880162170248896892053994937480676824100235778078258089007324143972798519797707032130647969884610321425997189568119602840375894038155972908248004382498468916557706979429998716578160611084323023757946616513714721274482895333064156415741467388984618053241491703589401980162709229995190315146112057114148528371268073240697867276018185731510023421667460054167043946189996319239451804850583222207623709979312120485235589607032199701387734728901493075402111739944241663097291682065128415467408794012133045576298823323231644319444771309101979838714607268362885013353236236279028547075826693271711076508425978753112672873471483112345162587327323284559170894801540946256076581492385877327920830212228117773961574323020309037154581962604279919770924565083410929733620139585525373839607694142835671049697791915617732496134164500439645102443987515076670635346496272691715990210824938942259376967327736103437367741646251904790549748690944569347164625686510650942902612493691144502089014530540918554198933184653476942446080685583723088707560004626657585412385543184301378502185399284958449745001
stack 6
$ python3 -m timeit -s 'import flake.sbtree as mod' 'mod.work_speclimit (mod.Specs.Sqrt5 (), 10000)'
10000 loops, best of 3: 93.6 usec per looptarget 2.645751311064591
result 2.645751311064591
88685109179915130341887104907585686521847537035156964074461380628433449732421486996699064062444540064158389056209363824310211054734553092226299945457021199112700497855553114289504071487912489870285448166897774268426045437736838881291135766978714310577197146720994928482690407035258346866833661597831052596381967292664054207176771296004452023561901916859428680387946648537327818721735456944347866760599816478742091731402036923914749011584466135017369187496638702604704904099567360882401992490615065358974108247116798594573520556021596510768434830377703395683058096577446210575258033599847116433595329505939001912158627483541548413908052794044326067277738229833351839833732621914148391648757352594454100888285735908886235371717556569350093976323499838384402503010849693022056133132459889455774624180654855740421068215212274465740819478763364431816220680747532316564598964601291404834137936045821906603748156817494328651241872320246841849963968760441611551680171915171870107132895649133384251789520357115932060244542818093510776757680312455843828454609567341874249691543056825373273269204071745454267556076788350083854242116696619839481745027404336920399635197690200093935910164689111311839586735442346668112486814384577793261763215032118149485337833479163988065715333603367509804484339306442650722582405639197507565528586654028311449710723867519836823306113282478980440642161464033492837646551961460316452507922119182266791487068039633701371545426504396858443073658611354269445043000296653262607719089598658043789113954673633318703927713435275149460809489273849788240351191345839326411839428019773563791946938637627130882929019426915426450059256851311863911361172986681994529441158570079902956177659085523863706292659368
33519820554952402077344952617115996430062602096298132910527390831423668803224786792409482114736441464414434466600804779984731735712635377260497095812209531262037812938763480625651274971130960829211892715611645468602169701448987251532767391210039108389525077224714582291644331715154178763831397834298494044533585708820980314756995071855712127833971373932571326174459953709149532727724636735804897739228824429719381518007647530374783592205329481298090989741894174116420193651705325603383474062569368008883036209486064571228798270612570245995726989027981078917804184187499121621704880192862070899839335763600538337612503028675222083143790259516383967477391105199604260172432581074516295386295207443781246839223963989741293797783635987581484624039801831036560454434936317427599590411496224537788272692769212532372138459731579486311965568423854036095144587150585857725218696351547762091779862472592476271252396689972340657268628812154795429091581895847878967677044235215574631131931485165870936334656348188917546216466160109951993701418977615343196637062373754515878578286200786667708184435930028647016106141546838211875075480787744084715947501964293329801736675201320126556780597972877755791737287370240036747650626088972123685671173383351988672241516353767982686940196092077710683695881056812441315674315369699222038356904379301779532128607451321996820067317087887933988110334316185513431274490816808887993078291114719390594741039617560489171670668108767799640798842133852400136686687954287968401830668881914079204844843027519987932569505672423814843276175818680085566684747600593430391192915923373560623197517985154594912997008324835735525214421317644885339320934402944417757163610788379320383776665001229589830530087683
stack 7
$ python3 -m timeit -s 'import flake.sbtree as mod' 'mod.work_speclimit (mod.Specs.Sqrt7 (), 10000)'
10000 loops, best of 3: 120 usec per loop[1/2]
import math
# 0 1 a c a a+c 1 1 a+c c 1 0
# 1 0 b d b b+d 0 1 b+d d 1 1
class Matrix:
I = [1, 0, 0, 1]
S = [0, 1, 1, 0]
L = [1, 1, 0, 1]
R = [1, 0, 1, 1]
@ staticmethod
def new ():
return [None] * 4
@ staticmethod
def copy (m):
return m [:]
@ staticmethod
def copyinto (msrc, mdst):
mdst [:] = msrc
@ staticmethod
def pair (m):
return (m [0] + m [1], m [2] + m [3])
@ staticmethod
def mul (m1, m2, mout):
mout [0] = m1 [0] * m2 [0] + m1 [1] * m2 [2]
mout [1] = m1 [0] * m2 [1] + m1 [1] * m2 [3]
mout [2] = m1 [2] * m2 [0] + m1 [3] * m2 [2]
mout [3] = m1 [2] * m2 [1] + m1 [3] * m2 [3]
class Stack:
@ staticmethod
def alloc (size, matrixclass):
new = matrixclass.new
return [new () for k in range (size)]
@ staticmethod
def ensure (stack, need, matrixclass):
have = len (stack)
if have < need:
stack.extend (Stack.alloc (need - have, matrixclass))
class Expo:
@ staticmethod
def bits (m, exp, mstack, stackpos, matrixclass):
MC = matrixclass
if exp < 1:
MC.copyinto (MC.I, mstack [stackpos])
return
if exp == 1:
MC.copyinto (m, mstack [stackpos])
return
outnow = stackpos
outother = stackpos + 1
powernow = stackpos + 2
powerother = stackpos + 3
div, mod = divmod (exp, 2)
left = div
mul = MC.mul
Stack.ensure (mstack, stackpos + 4, MC)
MC.copyinto (m, mstack [powernow])
MC.copyinto (m if mod else MC.I, mstack [outnow])
while left:
div, mod = divmod (left, 2)
left = div
mul (mstack [powernow], mstack [powernow], mstack [powerother])
powernow, powerother = powerother, powernow
if mod:
mul (mstack [outnow], mstack [powernow], mstack [outother])
outnow, outother = outother, outnow
if outnow != stackpos:
MC.copyinto (mstack [outnow], mstack [stackpos])
# context
# stack
# matrixclass
class Group:
@ classmethod
def kind (cls):
pass
def eval (self, ctx, stackpos):
pass
class EmptyGroup (Group):
@ classmethod
def kind (cls):
return "empty"
def eval (self, ctx, stackpos):
MC = ctx ["matrixclass"]
stack = ctx ["stack" ]
MC.copyinto (MC.I, stack [stackpos])
class FixedGroup (Group):
@ classmethod
def kind (cls):
return "fixed"
def __init__ (self, fixedstr):
self.fixedstr = fixedstr
def eval (self, ctx, stackpos):
MC = ctx ["matrixclass"]
stack = ctx ["stack" ]
fs = self.fixedstr
count = len (fs)
if count < 1:
MC.copyinto (MC.I, stack [stackpos])
return
sym = {'L': MC.L, 'R': MC.R}
if count == 1:
MC.copyinto (sym [fs], stack [stackpos])
return
div, mod = divmod (count - 1, 2)
if mod:
now = stackpos + 1
other = stackpos
else:
now = stackpos
other = stackpos + 1
Stack.ensure (stack, stackpos + 2, MC)
mnow = stack [now ]
mother = stack [other]
MC.copyinto (sym [fs [0]], mnow)
mul = MC.mul
fspos = 1
for k in range (div):
mul (mnow, sym [fs [fspos ]], mother)
mul (mother, sym [fs [fspos + 1]], mnow)
fspos += 2
if mod:
mul (mnow, sym [fs [fspos]], mother)
[2/2]
class ListGroup (Group):
@ classmethod
def kind (cls):
return "list"
def __init__ (self, groups):
self.groups = groups
def eval (self, ctx, stackpos):
MC = ctx ["matrixclass"]
stack = ctx ["stack" ]
gs = self.groups
count = len (gs)
if count < 1:
MC.copyinto (MC.I, stack [stackpos])
return
if count == 1:
gs [0].eval (ctx, stackpos)
return
div, mod = divmod (count - 1, 2)
if mod:
now = stackpos + 1
other = stackpos
else:
now = stackpos
other = stackpos + 1
Stack.ensure (stack, stackpos + 3, MC)
factor = stackpos + 2
mfactor = stack [factor]
mnow = stack [now ]
mother = stack [other ]
gs [0].eval (ctx, now)
mul = MC.mul
gspos = 1
for k in range (div):
gs [gspos ].eval (ctx, factor)
mul (mnow, mfactor, mother)
gs [gspos + 1].eval (ctx, factor)
mul (mother, mfactor, mnow)
gspos += 2
if mod:
gs [gspos].eval (ctx, factor)
mul (mnow, mfactor, mother)
class ExpoGroup (Group):
@ classmethod
def kind (cls):
return "expo"
def __init__ (self, base, expo, *, alwaysevalbase = False):
self.base = base
self.expo = expo
self.alwaysevalbase = alwaysevalbase
def eval (self, ctx, stackpos):
MC = ctx ["matrixclass"]
stack = ctx ["stack" ]
base = self.base
expo = self.expo
always = self.alwaysevalbase
if (expo < 1) and not always:
MC.copyinto (MC.I, stack [stackpos])
return
Stack.ensure (stack, stackpos + 2, MC)
base.eval (ctx, stackpos)
Expo.bits (stack [stackpos], expo, stack, stackpos + 1, MC)
MC.copyinto (stack [stackpos + 1], stack [stackpos])
class Spec:
@ classmethod
def kind (cls):
pass
@ classmethod
def targetfloat (cls):
pass
class StarSpec (Spec):
@ classmethod
def kind (cls):
return "star"
def __init__ (self, group):
self.group = group
class Compiler:
@ staticmethod
def fixedlimit (spec, limit):
what = spec.kind ()
ret = getattr (Compiler, "fixedlimit_" + what) (spec, limit)
return ret
@ staticmethod
def fixedlimit_star (spec, limit):
if limit < 1:
return {
"group": EmptyGroup (),
"count": 0,
}
g = spec.group
what = g.kind ()
if what != "fixed":
raise ValueError (what)
fs = g.fixedstr
fslen = len (fs)
if fslen < 1:
raise ValueError (fslen)
div, mod = divmod (limit, fslen)
if mod:
h = ListGroup ([
ExpoGroup (g, div),
FixedGroup (fs [:mod])])
else:
h = ExpoGroup (g, div)
return {
"group": h,
"count": limit,
}
class Specs:
class Phi (StarSpec):
@ classmethod
def targetfloat (cls):
return (1 + math.sqrt (5)) / 2
def __init__ (self):
StarSpec.__init__ (self, FixedGroup ("RL"))
class Sqrt2 (StarSpec):
@ classmethod
def targetfloat (cls):
return math.sqrt (2)
def __init__ (self):
StarSpec.__init__ (self, FixedGroup ("RLLR"))
class Sqrt3 (StarSpec):
@ classmethod
def targetfloat (cls):
return math.sqrt (3)
def __init__ (self):
StarSpec.__init__ (self, FixedGroup ("RLR"))
class Sqrt5 (StarSpec):
@ classmethod
def targetfloat (cls):
return math.sqrt (5)
def __init__ (self):
StarSpec.__init__ (self, FixedGroup ("RRLLLLRR"))
class Sqrt7 (StarSpec):
@ classmethod
def targetfloat (cls):
return math.sqrt (7)
def __init__ (self):
StarSpec.__init__ (self, FixedGroup ("RRLRLRR"))
def work_speclimit (spec, limit):
log = print
M = Matrix
C = Compiler
slots = 2
mstack = Stack.alloc (slots, M)
target = spec.targetfloat ()
ctx = {
"matrixclass": M,
"stack": mstack,
}
log ("target {:.15f}".format (target))
ret = C.fixedlimit (spec, limit)
group = ret ["group"]
group.eval (ctx, 1)
M.mul (M.S, mstack [1], mstack [0])
a, b = M.pair (mstack [0])
log ("result {:.15f}".format (a / b))
log ("{:d}".format (a))
log ("{:d}".format (b))
log ("stack {:d}".format (len (mstack)))Is there anything more ugly than inserting whitespace before the opening parenthesis?
>>29
not using lisp scheme and instead algol
>>30
True
cpu time: 2995 real time: 2997 gc time: 1174
group List (Fixed RRLRR, List (Loop (498, List (Store (a, List (Load (a, Fixed LR), Load (l4, Fixed LLLL))), Store (b, List (Load (b, Fixed RLRR), Load (r4, Fixed RRRR))))), List (Fixed LR, Expo (Fixed L, 1996), Fixed RL, Expo (Fixed R, 999))))
target 2.718281828459045
result 2.718281828459045
340017086135953955237517024368071023319667336443292804743952295811218300447359859240275563254009942677431489707597147675747035681413226512627011095930560268265456381178976392711767087225557220527288583530246999554026111396756681847411034069691700458364305110849779785985479469458249422285136155670021333443718590111827401706503305002386802304743268641986242675018801066689637090293781703589884119797742595395970458469615457316464652161835058300625085365074839652205963909281728057007645884413972470559421095590641084923213749096077911891464157447244982848623998799902248659508491205347662407542487598477506176854976430107224759805723469500480207518560065060801139517286177001208163592054763542622442833088837014707870098680498790101406422989090225950831187269091905076078387996734933063315154563648844011105458348826524339618765836090847570234445324485311773505775960915934136456454684723065101560864875759072426638458977878476314281969987155314153269630993616187041722772926774199317277275408014196669434574127395283760934088304504454441826341472487028029542307844361121860169948845000863967607247499379959564922815018909507873427673764317550105676071503852510957626027913732521804089379035997601445017705811769908093393745107753664844766235361183695170527428489809363428463651419304155765399729238548589316792609372558848237780044857470539820865022763311896953342264917235482949907211105305961474199254134498023181661330808198705378200836895070871333229848624381644393862053739146946714031748697068394959781227898824577572717091477283125470687771132502660646193012227555089150091316446335434009362342321577385694092293349561432035238800308129504247688449971272567391760416476691624159595214779684041072400751170742699454671514147869612822126793025761288002185791555487625331998621292869836022457378328419035361513580174231044804345605326811465916238059290186883564724761265694886919748527987916817620628526399573566506494861425004396642032056780484911383385189989513516563131778797299221345966468968956122130938820067631844304729073653295278994294270111849584185735076742525508006069816526335483941312168466769941517563643711961854946980989175852514399703705680462519538562629390325716019370164673384882170485495116424556766549082251149470708965255450905049612453587859276129878129571158069679466773167708728881129438071391603459651281058662550128959761153549617728975030404264330537814833057633963701181617935774238987394896893837492212211374342977762904418298881285931234958814978103651305047589748736431989381870422496158486260942264340691466697262368432164893697234466153250046703893250452429139844480273241829106099238215033421990835160228005567092855709150858705097715014615303931593755514322351810007680897794492648609028769501447879692294277237362011716820403857042445580674191735659531498184966721518718160689953505802281847373678940279750467265390416599246074195127826056496476091242971985446723655542587173824287311438381591741110560862978217691099053628393417793105301826047059921235216516021504425461929120276377753800102305173326425003579123379528380303069690154701626211729549743364833186844412203484432054917804678030588228660410069011019359882579002
125085295636436909683809909842556856885178429583466250503633118421625758576495571982549603221529371824249449947226736257956572986278847452327489918555893035992800158136526570747790366607976473974086490509833266140241624411002520378548860256844094890877496352409672412181476991240644660422957359119247812657697715895808101749155909929315560672524600833028834500581928675977242236292980202330265409027931557849262071183680232807134194663977226617526758164092581075236564819647024212853553521938166754805064971997583680157504280085693771587742124565102644454392184568192455131506363793673061697667809133470803909462294541130705110877585153688820760830173021946921156853394209265699621359388107424247279451921507178896866252151258182423463805654515203717737130171743215112991658603794371604988588191145881380551692440558754609581883970182436373351721827090765483767025558740593461170481647907645229692750706772617083777913648346892911751008260597754177442242884897528540552656654142051828032526143451856483844303282993772253401076047130518814970561671090962371290205074655269040331486455093397311641724716513914711853050343863410602433828186055109883845169168186241477852812819511630717226031782266261200943940615962914287809900772211847247623544281170759686802721987002900335098524747135307554944852406985163571494649355659613567838732869723194876669393386105031678437801228285473666870189915779831674530812838317813362749429334047728206554077502413410488307816529862238758762434688416460199490950423539362723027760398951127569670882323069657999986406666638536904433413611732504237085435613428187295322449920616488897190353761498046868750584302985283042516291072517661563543124917633435349040950215732545008269904452842072489224708167231629873428753864275250495686428820211229340636698055009463068010066247920722428973426571945857005028864658276369838400927087527667218274517604615390540969915685158805164446004960391637638952325069774576683392245223799253721839340049299692557960655127171880922366931873834678505507980929392029365039716768900191201134404082041954658457421412098953281935182415736165076549914653156936288250673567627314590360080813173455844097539997651557001995149793696142371449923446980427725864245921082517781640485031251339553931144656900391589822516702794312100705930445806007992360069962039035567107884172990813990752000736096406721872831171253639698564151167421533401665295128817139455342798715270038691942997621086685910211094868754454516817862989819051640051770874497216852913390188488061385723926005176343183394443750361968944989563225128221788985578992361222409045215774195773008647206524738774180554917621147514549299392328337681997169072439473578702176284273236998017086945818255227347039675511648057103771074013479256786923260298195015744892867688675269801389398308311777208591853723170164583690632297703927756753727378438541760265853399635262151837406060792526386670515259891008695811004879840637123479657107736250722757524021127684175490377871823327703940937023628221667804570593378320396524365465248186987504652311179460359354823215116827906074424119159403212171240501026466951623914101612248957295893464404775833237162234624918433172999
stack 12
$ python3 -m timeit -s 'import flake.sbtree as mod' 'mod.work_speclimit (mod.Specs.E (), 1000000)'
100 loops, best of 3: 4.5 msec per loopSpec for the e stream:
class StamploopSpec (Spec):
@ classmethod
def kind (cls):
return "stamploop"
def __init__ (self, stamp):
self.stamp = stamp
def stamplength (self, index):
pass
def tailgroup (self, index):
pass
def work_loadstore_repeat ():
return ListGroup ([
StoreGroup ("a", ListGroup ([
LoadGroup ("a", FixedGroup ("LR")),
LoadGroup ("l4", FixedGroup ("LLLL"))])),
StoreGroup ("b", ListGroup ([
LoadGroup ("b", FixedGroup ("RLRR")),
LoadGroup ("r4", FixedGroup ("RRRR"))]))])
class EStamp (StamploopSpec):
def __init__ (self):
StamploopSpec.__init__ (self,
work_loadstore_repeat ())
def stamplength (self, index):
return 14 + 8 * index
def tailgroup (self, index):
return ListGroup ([
FixedGroup ("LR"),
ExpoGroup (FixedGroup ("L"), 4 + 4 * index),
FixedGroup ("RL"),
ExpoGroup (FixedGroup ("R"), 6 + 4 * index)])
class Specs:
class E (ListGroup):
@ classmethod
def targetfloat (cls):
return math.e
def __init__ (self):
ListGroup.__init__ (self, [
FixedGroup ("RRLRR"), EStamp ()])The daily code-dump is much appreciated.
The additional building blocks used by the e spec:
# context
# stack
# matrixclass
# storage
class Group:
@ classmethod
def kind (cls):
pass
def eval (self, ctx, stackpos):
pass
def show (self):
pass
def isfixed (self):
pass
def fixedlength (self):
pass
class LoopGroup (ListGroup):
class Loop:
def __init__ (self, count, item):
self.count = count
self.item = item
def __len__ (self):
return self.count
def __getitem__ (self, key):
return self.item
# no slicing
@ classmethod
def kind (cls):
return "loop"
def __init__ (self, count, group):
ListGroup.__init__ (self, LoopGroup.Loop (count, group))
self.count = count
self.group = group
def show (self):
return "Loop ({}, {})".format (self.count, self.group.show ())
def isfixed (self):
return self.group.isfixed ()
def fixedlength (self):
return self.count * self.group.fixedlength ()
class LoadGroup (Group):
@ classmethod
def kind (cls):
return "load"
def __init__ (self, name, fallback):
self.name = name
self.fallback = fallback
def eval (self, ctx, stackpos):
MC = ctx ["matrixclass"]
stack = ctx ["stack" ]
storage = ctx ["storage" ]
name = self.name
value = storage.get (name)
if value is None:
self.fallback.eval (ctx, stackpos)
storage [name] = MC.copy (stack [stackpos])
else:
MC.copyinto (value, stack [stackpos])
def show (self):
return "Load ({}, {})".format (self.name, self.fallback.show ())
def isfixed (self):
return False
def fixedlength (self):
raise Wrong
class StoreGroup (Group):
@ classmethod
def kind (cls):
return "store"
def __init__ (self, name, group):
self.name = name
self.group = group
def eval (self, ctx, stackpos):
MC = ctx ["matrixclass"]
stack = ctx ["stack" ]
storage = ctx ["storage" ]
name = self.name
group = self.group
group.eval (ctx, stackpos)
# after
m = storage.get (name)
if m is None:
storage [name] = MC.copy (stack [stackpos])
else:
MC.copyinto (stack [stackpos], m)
def show (self):
return "Store ({}, {})".format (self.name, self.group.show ())
def isfixed (self):
return self.group.isfixed ()
def fixedlength (self):
return self.group.fixedlength ()The evaluation context updated with storage:
def work_fixedgroup (target, group, *, debug = True, matrixclass = None):
log = print
M = Matrix if matrixclass is None else matrixclass
slots = 2
mstack = Stack.alloc (slots, M)
ctx = {
"matrixclass": M,
"stack": mstack,
"storage": {},
}
if debug:
log (" group {}".format (group.show ()))
log ("target {:.15f}".format (target))
group.eval (ctx, 1)
M.mul (M.S, mstack [1], mstack [0])
a, b = M.pair (mstack [0])
if debug:
log ("result {:.15f}".format (a / b))
log ("{:d}".format (a))
log ("{:d}".format (b))
log ("stack {:d}".format (len (mstack)))
def work_speclimit (spec, limit, *, debug = True, matrixclass = None, compiler = None):
C = Compiler if compiler is None else compiler
target = spec.targetfloat ()
ret = C.fixedlimit (spec, limit)
group = ret ["group"]
work_fixedgroup (target, group, debug = debug, matrixclass = matrixclass)The updated compiler:
class Compiler:
@ staticmethod
def fixedlimit (spec, limit):
what = spec.kind ()
fun = getattr (Compiler, "fixedlimit_" + what)
# after
if limit < 1:
return {
"group": EmptyGroup (),
"count": 0,
}
ret = fun (spec, limit)
return ret
@ staticmethod
def fixedlimit_fixed (spec, limit):
fs = spec.fixedstr
fslen = len (fs)
if fslen < 1:
return {
"group": EmptyGroup (),
"count": 0,
}
if fslen <= limit:
return {
"group": spec,
"count": fslen,
}
return {
"group": FixedGroup (fs [:limit]),
"count": limit,
}
@ staticmethod
def fixedlimit_expo (spec, limit):
base = spec.base
baselen = base.fixedlength ()
expo = spec.expo
have = expo * baselen
if have < 1:
return {
"group": EmptyGroup (),
"count": 0,
}
if have <= limit:
return {
"group": spec,
"count": have,
}
div, mod = divmod (limit, baselen)
if mod:
ret = Compiler.fixedlimit (base, mod)
gmod = ret ["group"]
if ret ["count"] != mod:
raise Wrong (mod, ret ["count"])
if div:
if div > 1:
h = ListGroup ([ExpoGroup (base, div), gmod])
else:
h = ListGroup ([base, gmod])
else:
h = gmod
else:
if div > 1:
h = ExpoGroup (base, div)
else:
h = base
return {
"group": h,
"count": limit,
}
@ staticmethod
def fixedlimit_list (spec, limit):
gs = spec.groups
gslen = len (gs)
if gslen < 1:
return {
"group": EmptyGroup (),
"count": 0,
}
take = []
left = limit
used = 0
flim = Compiler.fixedlimit
while (left > 0) and (used < gslen):
ret = flim (gs [used], left)
used += 1
now = ret ["count"]
if now > 0:
take.append (ret ["group"])
left -= now
taken = len (take)
if taken < 1:
return {
"group": EmptyGroup (),
"count": 0,
}
return {
"group": take [0] if taken == 1 else ListGroup (take),
"count": limit - left,
}
@ staticmethod
def fixedlimit_star (spec, limit):
g = spec.group
if not g.isfixed ():
raise Wrong (g.kind ())
have = g.fixedlength ()
if have < 1:
raise Wrong
div, mod = divmod (limit, have)
if mod:
ret = Compiler.fixedlimit (g, mod)
gmod = ret ["group"]
if ret ["count"] != mod:
raise Wrong (mod, ret ["count"])
if div:
if div > 1:
h = ListGroup ([ExpoGroup (g, div), gmod])
else:
h = ListGroup ([g, gmod])
else:
h = gmod
else:
if div > 1:
h = ExpoGroup (g, div)
else:
h = g
return {
"group": h,
"count": limit,
}
@ staticmethod
def fixedlimit_stamploop (spec, limit):
more = True
index = 0
left = limit
stal = spec.stamplength
while more:
now = stal (index)
if now <= left:
index += 1
left -= now
more = (left > 0)
else:
more = False
out = []
if index > 1:
out.append (LoopGroup (index, spec.stamp))
elif index == 1:
out.append (spec.stamp)
if left > 0:
g = spec.tailgroup (index)
ret = Compiler.fixedlimit (g, left)
gmod = ret ["group"]
if ret ["count"] != left:
raise Wrong (left, ret ["count"])
out.append (gmod)
outlen = len (out)
if outlen < 1:
raise Wrong
return {
"group": out [0] if outlen == 1 else ListGroup (out),
"count": limit,
}Ratios from OP's timings, using vanilla python without gmpy2. For (phi, 10^4) >>21 >>22:
>>> 33000 / 137
240.87591240875912>>> 2995 / 4.5
665.5555555555555Both ratios will increase when the stream index increases.
Mhhh, undocumented spam-code... love it when it takes up more than half the space on the frontpage.
Moving the integer arithmetic to gmpy2 yields a modest speedup for (e, 10^6), shaving off less than 8% of the run time.
$ python3 -m timeit -s 'import numsci.sbtree as mod' 'mod.work_speclimit (mod.fsb.Specs.E (), 1000000)'
100 loops, best of 3: 4.17 msec per loopThe gains should increase for higher stream positions as multiplications take over the bulk of the run time.
class Matrix:
z0 = gmpy2.mpz (0)
z1 = gmpy2.mpz (1)
I = [z1, z0, z0, z1]
S = [z0, z1, z1, z0]
L = [z1, z1, z0, z1]
R = [z1, z0, z1, z1]>>40
We all know that that is not how code documentation works.
Here is stream index 10^12 in the e stream, to check against other implementations.
http://0x0.st/-K_o.gz
$ gawk '{ if (length ($0) < 1000) { print NR ": " $0 } else { print NR "! " length ($0) } }' e12.txt
1: group List (Fixed RRLRR, List (Loop (499998, List (Store (a, List (Load (a, Fixed LR), Load (l4, Fixed LLLL))), Store (b, List (Load (b, Fixed RLRR), Load (r4, Fixed RRRR))))), List (Fixed LR, Expo (Fixed L, 1999996), Fixed RL, Expo (Fixed R, 999999))))
2: target 2.718281828459045
3: result 2.718281828459045
4! 6167766
5! 6167765
6: stack 12The download is nearly 6 megabytes, the uncompressed file just over 12 megabytes and the numerator has 6167766 digits, so text editors might not like it.
>>40
Do not feed.
The phi stream grows much faster in numerator and denominator digit counts, since they are essentially the fibs, and by stream index 10^9 we are already past 200 million.
$ gawk '{ if (length ($0) < 1000) { print NR ": " $0 } else { print NR "! " length ($0) } }' phi9.txt
1: group Expo (Fixed RL, 500000000)
2: target 1.618033988749895
3: result 1.618033988749895
4! 208987641
5! 208987641
6: stack 6The e stream was still just a little over 6 million digits >>43 at a thousand times the stream index. For verification here is the hash of the phi 10^9 numerator:
$ sha1sum <(gawk 'NR == 4 { printf "%s", $0 }' phi9.txt)
4a74c3993775189c754cfcc6ecf5de166bb718c7 /dev/fd/63And the hash of the denominator:
$ sha1sum <(gawk 'NR == 5 { printf "%s", $0 }' phi9.txt)
24e5e3350b839eea494e9daec46eec323f0f32dd /dev/fd/63>>43
We all know who the real troll here is.
If there is no substantial engagement in the thread, we all know who the real spammer is. Just get your own site.
The (e, 10^12) link >>43 will expire at some point so here is the hash of the numerator:
$ sha1sum <(gawk 'NR == 4 { printf "%s", $0 }' e12.txt)
2a4d66da41f5564d04a5c4356159dbc7ccdd9827 /dev/fd/63And the hash of the denominator:
$ sha1sum <(gawk 'NR == 5 { printf "%s", $0 }' e12.txt)
e405d307f115213561485d57d5e4d9b8f658b047 /dev/fd/63Stream index 10^9 for the Sqrt2 >>28 stream:
$ gawk '{ if (length ($0) < 1000) { print NR ": " $0 } else { print NR "! " length ($0) } }' sqrt29.txt
1: group Expo (Fixed RLLR, 250000000)
2: target 1.414213562373095
3: result 1.414213562373095
4! 191387843
5! 191387843
6: stack 6The hash of the numerator:
$ sha1sum <(gawk 'NR == 4 { printf "%s", $0 }' sqrt29.txt)
e726c20ee05cebf894df98462d5c793a08318c66 /dev/fd/63The hash of the denominator:
$ sha1sum <(gawk 'NR == 5 { printf "%s", $0 }' sqrt29.txt)
86e33050005cc854e7ba9b734b9f2a2f1333e6a0 /dev/fd/63A hypothesis >>43 >>44 presents itself. It would seem that the number of L -> R and R -> L switches is an indicator of the rate at which the numerator and denominator will grow.
spec index digits group
e 10^12 6167766 RRLRR(LRL[4+4*k]RLR[8+4*k])*
phi 10^9 208987641 (RL)*
sqrt2 10^9 191387843 (RLLR)*The phi stream has the maximum number of switches a stream can have, and it grows the fastest. The sqrt2 stream has half the number of switches and it grows slower than phi. The e stream has contiguous runs of increasing length between switches, and it grows very slowly. This should predict the rate of growth of the other >>28 streams based on their generating pattern.
* cricket sounds *
The [8+4*k] (eight) in >>47 is a typo. The correct length of the R run is [6+4*k] (six) as in >>33.
An equivalent group for e is:
R(RLR[2+4*k]LRL[4+4*k])*(Sqrt3, 10^9) >>28:
$ gawk '{ if (length ($0) < 1000) { print NR ": " $0 } else { print NR "! " length ($0) } }' sqrt39.txt
1: group List (Expo (Fixed RLR, 333333333), Fixed R)
2: target 1.732050807568877
3: result 1.732050807568877
4! 190649183
5! 190649183
6: stack 7Numerator/denominator hashes:
$ sha1sum <(gawk 'NR == 4 { printf "%s", $0 }' sqrt39.txt)
d464586d1eccfb68056109d70d486a142f9d8992 /dev/fd/63
$ sha1sum <(gawk 'NR == 5 { printf "%s", $0 }' sqrt39.txt)
98fcfa9aada30daa813b423f30f79a1d6390b902 /dev/fd/63(Sqrt5, 10^9) >>28:
$ gawk '{ if (length ($0) < 1000) { print NR ": " $0 } else { print NR "! " length ($0) } }' sqrt59.txt
1: group Expo (Fixed RRLLLLRR, 125000000)
2: target 2.236067977499790
3: result 2.236067977499790
4! 156740731
5! 156740731
6: stack 6Numerator/denominator hashes:
$ sha1sum <(gawk 'NR == 4 { printf "%s", $0 }' sqrt59.txt)
c6e680622d444ae51479e5f4f8e2f2902436674f /dev/fd/63
$ sha1sum <(gawk 'NR == 5 { printf "%s", $0 }' sqrt59.txt)
7396e2d4afd5060c98e11267bbb4c80389325a35 /dev/fd/63(Sqrt7, 10^9) >>28:
$ gawk '{ if (length ($0) < 1000) { print NR ": " $0 } else { print NR "! " length ($0) } }' sqrt79.txt
1: group List (Expo (Fixed RRLRLRR, 142857142), Fixed RRLRLR)
2: target 2.645751311064591
3: result 2.645751311064591
4! 171773357
5! 171773356
6: stack 7Numerator/denominator hashes:
$ sha1sum <(gawk 'NR == 4 { printf "%s", $0 }' sqrt79.txt)
06253fc6faf29a64a2cc1bd2a0b159e2968865ab /dev/fd/63
$ sha1sum <(gawk 'NR == 5 { printf "%s", $0 }' sqrt79.txt)
e91a7792c7c4c91e0c57b432036306c4c72b83c1 /dev/fd/63The hypothesis of >>47 does not hold.
spec index digits group
e 10^12 6167766 RRLRR(LRL[4+4*k]RLR[6+4*k])*
sqrt5 10^9 156740731 (RRLLLLRR)*
sqrt7 10^9 171773357 (RRLRLRR)*
sqrt3 10^9 190649183 (RLR)*
sqrt2 10^9 191387843 (RLLR)*
phi 10^9 208987641 (RL)*Sqrt3 has more switches (2N/3) than Sqrt2 (N/2) but grows at a slower rate. Sqrt7 also has more switches (4N/7) than Sqrt2 and slower growth.
Just realized the Stern-Brocot tree being the list of all coprime numbers could be used to generate the pythagorean triples (uint a,b,c such that a²+b²=c²). We need to take only the right subtree where every m > n.
Euclide's formula to generate a pythagorean triple, with m and n coprimes:
a = m²-n², b = 2mn, c = m²+n²If we need only primitive pythagorean triples then we just need to filter out the fractions where both m and n are odd.