Challenge Weekly 2.0 (took a year)

Time for a new Challenge thread where each week a challenge will be set for all to follow. There's different challenge levels with 1 being the easiest and 5 being the most difficult. C3W1: Implement a bytebeat synth playing audio directly with either OpenAL or your platform's audio library (e.g alsa, oss, sndio) C2W1: Implement the fibonacci sequence with the tree recursive solution but avoid exponential time complexity. The program must be well defined for the first 100 fibonacci terms. This is the typical interview question. Haskell, Idris, and Agda may not be used for this challenge, bonus points for figuring out why. C1W1: Implement a fizzbuzz program without using the modulo function (%, mod, rem, &c...) Fun task: Implement a variation of the little man computer in any two languages. Memory: An array of 0xffff bytes, this is an 8bit machine due to ease of implementation. Index must be a short int. Instruction set: HLT - halt the emulator and return register value ADD - add value from memory to register SUB - subtract memory value from register value MUL - multiply memory value with register value REM - remainder DIV - integer division SHL - shift left SHR - shift right LDA - load value from memory to register STA - store value from register to memory BRA - branch BRZ - branch if zero BRP - branch if positive INP - read character from user and store in register OUT - print character from register Next week preview: #define uppercase(x) x&(~0x20) #define between(x,y,z) y<x&&x<z #define guard(test,state,msg) if (test) { state.kind=Left;\ state.left=msg; return state; }

>>2036 Memofibtree is linear, which is boring. Here's the matrix method instead, which is logarithmic. def bitsgen (n): while n > 0: di, mo = divmod (n, 2) yield mo n = di def applyfibjump (a, b, p, q): return ( a * p + b * q, a * q + b * (p + q) ) def fibredstep (acc, bit): a, b, p, q = acc p, q = applyfibjump (p, q, p, q) if bit: a, b = applyfibjump (a, b, p, q) return (a, b, p, q) def fibredbits (n): if n <= 0: return 0 if n == 1: return 1 if n == 2: return 1 di, mo = divmod (n - 1, 2) if mo: state = (1, 1, 0, 1) else: state = (0, 1, 0, 1) state = functools.reduce ( fibredstep, bitsgen (di), state) # print (state) return state [1] Terms 100 and 101 with 0 at index 0 and 1 at index 1: 354224848179261915075 573147844013817084101

Have this solution lying about from going through SICP (define (fib n) (fib-iter 1 0 0 1 n)) (define (fib-iter a b p q count) (cond ((= count 0) b) ((even? count) (fib-iter a b (+ (* p p) (* q q)) (* q (+ q (* 2 p))) (/ count 2))) (else (fib-iter (+ (* b q) (* a q) (* a p)) (+ (* b p) (* a q)) p q (- count 1)))))

#lang racket ;; C1W1 (define (fizzbuzz til) (let loop ([n 1] [fizz 2] [buzz 4]) (when (<= n til) (cond [(and (zero? fizz) (zero? buzz)) (displayln 'fizzbuzz) (loop (+ n 1) 2 4)] [(zero? fizz) (displayln 'fizz) (loop (+ n 1) 2 (- buzz 1))] [(zero? buzz) (displayln 'buzz) (loop (+ n 1) (- fizz 1) 4)] [else (displayln n) (loop (+ n 1) (- fizz 1) (- buzz 1))])))) ;; C2W1 (define max-size 100) (define memo (make-vector (+ max-size 1) #f)) (vector-set! memo 0 0) (vector-set! memo 1 1) (define (fibs n) (if (<= n max-size) (or (vector-ref memo n) (let ([res (+ (fibs (- n 1)) (fibs (- n 2)))]) (vector-set! memo n res) res)) (+ (fibs (- n 1)) (fibs (- n 2))))) Challenge 3 coming later, probably. >>2048 >fib-iter >Implement the fibonacci sequence with the tree recursive solution

>>2052 > >Implement the fibonacci sequence with the tree recursive solution I should have posted it in reply to >>2041 I have the sliding solution written somewhere too... Is your solution not exponential? There are two calls to fibs in each iteration

>>2053 it's so-called "memoized", so if you run (fibs 10) then it will calculate it, and store (10 -> 55) in the vector. But since it's recursive, it also stores all the intermediate ones that it calculated, so if you call (fibs 11) it already knows the results of (fibs 9) and (fibs 10) and it doesn't need to recurse. When it was calculating (fibs 10), it needs to calculate (fibs 9) and (fibs 8), and if it runs them in that order, then (fibs 9) will only recurse twice (to get 7 and 8) instead of 2^9 times or whatever it would otherwise be, and because of that it's not exponential. I think

>>2053 >Is your solution not exponential? No, all the results of the computations from (fibs 1) to (fibs 100) are memoized. The reason why Agda, Haskell etc. are disqualified is because they memoize by default.

>>2057 How can the entire language memoize? I mean I understand not having internal state, but what if it tried to memoize "in.nextLine()" or something? I don't know exactly what it looks like because I don't know any of those languages, but it sounds like IO would be impossible

>>2058 In a purely functional language, the same inputs always map to the same outputs, so it can be lazy and only perform a computation when it needs to. I don't know any of those languages well either, so I couldn't tell you how they handle IO, but I know they make a very big deal out of it.

>>2058 >>2060 I also forgot to mention that all those languages make a huge distinction between functions that use IO and functions that don't in the type system, so obviously it can damn well infer what to memoize and what not to.

fizz = cycle ["", "", "Fizz"] buzz = cycle ["", "", "", "", "Buzz"] fizzbuzz :: [String] -> [String] -> Int -> [String] fizzbuzz (f:fs) (b:bs) n = (if null f && null b then show n else f ++ b) : fizzbuzz fs bs (n + 1) main = mapM_ putStrLn $ take 100 (fizzbuzz fizz buzz 1)

>>2056 >>2057 Well, TIL Thanks

>>2057 Haskell does not memoize by default.

>>2036 oh hey, this thread. I did the LMC in java a while ago, with some modifications and liberties taken to improve the programming. I added some additional instructions: LDR - Loads value stored in memory at the address given by the argument PRT - Prints register value as char LIT - Set register to value given by argument INC - Get input as character As well as some extra syntax stuff for the assembler: : - address labels, the assembler replaces them with the correct address for the corresponding label . - interpret the next "word" (token separated by space) as an integer (so the char 0 becomes value 0) $ - interpret the line as bare ascii values # - comments (ignore line) the . integers and $ lines are just treated as data in the program, they are useful if for example you want to store strings, you can use $ and write your string into the program, and use .0 to null-terminate. The first thing I noticed writing some LMC asm, is that its annoying as fuarrrk, you really miss the niceness of strings in high level languages. For example, the following is Hello World. (with my language modifications) LIT :string STA :loc LDR :loop BRZ :end PRT LDA :loc ADD :one STA :loc LDR BRA :loop :end HLT :one .1 :loc .0 :string $Hello, World! .0 You very quickly end up needing self-rewriting code and such. Fibonacci sequence (I think this works, I cant remember though) :loop LDA :f1 ADD :f2 OUT STA :f3 SUB :max LDA 1 ADD :one STA 1 LDA 3 ADD :one STA 3 LDA 7 ADD :one STA 7 BRA :loop :end HLT :one .1 :f1 .0 :f2 .1 :f3 .0 :max .127

>>2036 Is the bytebeat synth challenge mostly about interfacing with the libraries or what?

>>2094 Just pipe stdout to dev/dsp/. Otherwise too much lib fiddling (in my opinion).

>>2095 But then where's the challenge?

>>2096 Figuring out which bytes to print out to make it sound good

>>2036 if i implement the LMC in VHDL (or at least i try, barely know VHDL) it counts?

>fizzbuzz without modulo Well, here you go

>>2096 It's actually libs as I'm looking to expand on it later. It's listed as a C3 due to a later challenge. Think of it as more of a toy for now that could be implemented in one line of idris if it wasn't for the library requirement. >>2095 This doesn't actually pass as it's not portable to my system where I test solutions. >>2098 This is the fun side of it

>>2103 Yikes. How does one even come up with something like that? Seeing it prompted me to do the challenge. #include <iostream> bool divisible(float dividend, float divisor) { float quotient = dividend / divisor; return quotient == static_cast<int>(quotient); } int main() { for(int i = 1; i <= 100; ++i) { if(divisible(i, 15)) std::cout << "FizzBuzz\n"; else if(divisible(i, 3)) std::cout << "Fizz\n"; else if(divisible(i, 5)) std::cout << "Buzz\n"; else std::cout << i << "\n"; } } My heart is at ease once more. Simple truly is best.

>>2115 I like.

Week 2 is here and parsers are the new toys we'll get to play with. C3W2: Implement a parser and interpreter for the following. The language supports both functional and procedural style with only procedural function having access to IO. The only types in this language will be Nat, IO, String, and Bool. The shorter the solution the better. fun fib:Nat (n:Nat) // function declaration with type and argument | 0 => 0 // pattern match | 1 => 1 | n => fib(n-1) + fib(n-2) // call self if our number is above 1 proc main:IO (args:List String) var number = cast(Nat, (first(args))) // bind a value to the variable 'number' if (number > 2) then print(fib(number)) // display the result of the computation endwith(0) // exit with success C2W2: Implement a brainfuck interpreter with optimization and JIT, x64 is the target here. The standard instructions > , . < + - [ ] are to be implemented with loops terminating on '\0' to allow cat program.bf | ./bf "[,.]" - to terminate. C1W2: Implement a plotting program using ascii/ansi graphics. The program must come with it's own small repl for easy plotting of functions such as y = kx+m. Fun task will be posted later during the week.

Since OP only posted coursework this week, here's a standard interview question as well. Write a program that takes a row and column count as arguments and outputs a matrix of that size of properly sized and aligned cells, filled with consecutive numbers starting with 1 in the top left corner and proceeding in a clockwise inwards spiral. Example: output for (5, 7) 1 2 3 4 5 6 7 20 21 22 23 24 25 8 19 32 33 34 35 26 9 18 31 30 29 28 27 10 17 16 15 14 13 12 11

>>2182 That was harder than I expected, there's got to be an easier way than the garbage code I just wrote. http://0x0.st/0Il.c

>>2174 The rules of C3W2 are soykaf. Either give me BNF or don't ask me to write an interpreter.

>>2174 >The only types in this language will be Nat, IO, String, and Bool. >clearly has a List type

>>2174 I'm already planning on writing a language that is a tad more sophisticated than this (ML-like, rank-n polymorphism, existentials + typecase). Should I link the source here when I make progress?

>>2182 Needed parsers for a future challenge which is why they are a bit closer to the usual coursework style questions. Sould have had a better challenge but trying to keep it within beginner grounds. >>2190 I should correct that, Lists too. >>2186 Fell free to skip it. >>2199 Sure

>>2217 I'm really curious what the fun week 3 challenege was going to be.

>>2182 >>2184 Ugly Haskell here: import Prelude hiding (Left,Right) import Data.List --if using GHCI, to try it out call $ putStr $ finally $m $n finally n m = concat $ splitIntoRows (2*n) $ intersperse " " $ renderLocations n m splitIntoRows :: Int -> [String] -> [String] splitIntoRows n [] = [] splitIntoRows n xs = ((take n xs) ++ ["\n"]) ++ (splitIntoRows n (drop n xs)) renderLocations n m = fmap (pad.show.fst) (sortOn (snd) ([1..(n*m)] `zip` (numeralLocations n m))) where pad = \x -> (if length x == 1 then " "++x else if length x == 2 then " "++x else x) numeralLocations :: Int -> Int -> [(Int,Int)] numeralLocations n m = fmap foldPath (fmap (\i -> pathToNumeral n m i Right) [1..(n*m)]) data Direction = Right | Down | Left | Up deriving (Show, Eq) rotateLeft Right = Down rotateLeft Down = Left rotateLeft Left = Up rotateLeft Up = Right pathToNumeral n m i direction = if i > (n*m) then error "index out of range" else if i <= n then [step] else step:(pathToNumeral (m-1) n (i-n) (rotateLeft direction)) where step = (min i n, direction) foldPath :: [(Int,Direction)] -> (Int,Int) foldPath path = foldr stepPath (1, 0) path stepPath :: (Int, Direction) -> (Int,Int) -> (Int,Int) stepPath (i, d) | d == Up = \(a,b) -> (a-i,b) | d == Down = \(a,b) -> (a+i,b) | d == Right = \(a,b) -> (a,b+i) | d == Left = \(a,b) -> (a,b-i) >>2184

>>2182 Here is some rust. The spiral : ( >>2182 ) https://pastebin.com/4z4M7tJr The fizzbuzz : https://pastebin.com/U77rMBz7

>>2182 https://pastebin.com/UPAMi1fs There is my "bounds-checking" solution, which i admitted stupid. I thinking of better solution without any arrays or buffers.

>>2182 J box=: (([:>:[:i.*/+{.)$:]) : ((({.~{.) ; (}.~{.)$:[:(}.,[:<:{.)] )`([:<({.~{.))@.([:(=&0)[:(1&{)]))&, NB. F A B (-C----------) D E NB. A. first y[1] items of x NB. B. x, after dropping first y[1] items NB. C. (|.y), with the shifted element decremented NB. D. see A NB. E. if next cycle would get 0 for y[1], terminate on current cycle NB. F. monad definition: use integers, starting at 1 NB. box successively decreasing lengths of (,x) NB. lengths cycle through items of (,y) decrementing with each cycle transform=: [:,((,:&.>)`(([:|:,:)&.>)`(([:,:|.)&.>)`(([:|:[:,:|.)&.>))/. NB. A B C D NB. A. adds leading axis (see E) final direction: > NB. B. transpose only final direction: v NB. C. reverse and ,: (see E) final direction: < NB. D. reverse and transpose final direction: ^ NB. E. all boxes must contain items with 2 axes combine=: (((,)&.>)`((,~"1)&.>)`((,~)&.>)`((,"1)&.>))/ NB. top right bottom left NB. adds a list to one side of the growing matrix with each verb from the gerund spiral=: (([:>[:combine[:transform[:box{:,[:<:{.)`(''"_)@.([:0&e.])) : ([:$:[,]) NB. unbox (-A-----) (-B--------------) C NB. A. fixes input up for box: converts (x y) to (y (x-1)) NB. B. box is finnicky with 0 length axis; returns '' always NB. C. input is accepted as (spiral r c) or as (r spiral c)mod a little messy, but algorithmically clean

codility.com has some good challenges. What's nice is that it tests the solution against a range of inputs.

>>5450 I'm just a newbie but wouldn't this be more idiomatic for the "sum" function? Instead of your: fn sum(v: Vec<i32>) -> i32 { let mut s = 0; for i in v { s+= i}; s } This: fn sum(v: Vec<i32>) -> i32 { v.iter().fold(0, |acc, x| acc + x) } In any case, why?

>>2182 Another J solution, this time in one line: spiral =: 4 : '>: (x*y) - (|:@|.)^:2 ((, >./@:{: + >:@i.@#@{:)@|:@|.)^:(((x>y)&+)@((-&2)^:(-.@(2&|)@(x<.y)))@<:@(4&*)@<.@-:@(x<.y)) (>:($~ (,&1)@$) i.(-&(x>y)) 2+ | (x-y))'

>>2182 Here's a C99 solution: #include <stdbool.h> #include <stdio.h> #include <stdlib.h> #include <math.h> #define eprintf(...) fprintf(stderr, __VA_ARGS__) static int * spiral_matrix(int row_cnt, int col_cnt); static int safe_atoi(const char * s); int main(int argc, char ** argv) { if(argc != 3) { eprintf("usage: %s <row-count> <column-count>\n", argv[0]); exit(EXIT_FAILURE); } int row_cnt = safe_atoi(argv[1]); int col_cnt = safe_atoi(argv[2]); int * mat = spiral_matrix(row_cnt, col_cnt); int max_num = row_cnt * col_cnt; int max_digits = floor(log10(max_num)) + 1; for(int iy = 0; iy < row_cnt; iy++) { const char * prefix = ""; for(int ix = 0; ix < col_cnt; ix++) { printf("%s%*d", prefix, max_digits, mat[iy * col_cnt + ix]); prefix = " "; } printf("\n"); } free(mat); return EXIT_SUCCESS; } static int * spiral_matrix(int row_cnt, int col_cnt) { int * mat = calloc(row_cnt * col_cnt, sizeof(*mat)); if(!mat) { eprintf("error: out of memory (dimensions too big)\n"); exit(EXIT_FAILURE); } enum { R, D, L, U, NUM_DIRS }; int vecs[NUM_DIRS][2] = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}}; int arit[NUM_DIRS] = {0, 1, 0, 1}; int lens[2] = {col_cnt, row_cnt-1}; int pos[2] = {-1, 0}; int dir = R; int num = 1; while(true) { int n = lens[arit[dir]]; if(n == 0) { break; } for(int i = 0; i < n; i++) { pos[arit[dir]] += vecs[dir][arit[dir]]; mat[pos[1] * col_cnt + pos[0]] = num++; } lens[arit[dir]]--; dir = (dir + 1) % NUM_DIRS; } return mat; } static int safe_atoi(const char * s) { int res; int read_cnt = sscanf(s, "%d", &res); if(read_cnt != 1) { eprintf("error: badly formatted number: '%s'\n", s); exit(EXIT_FAILURE); } return res; }

C1W1: fzbz = [ "" for _ in range(101)] for a in range(3, 101, 3): fzbz[a] += "Fizz" for b in range(5, 101, 5): fzbz[b] += "Buzz" for c in range(1, 101): if fzbz[n] == "": print(n) else: print(fzbz[n])

>>5771 I guess it would, I just didn't know there was a fold function in rust.

>>2182 dumb haskell solution. spiral :: (Integral a) => a -> a -> [[a]] spiral rows cols = spiralRange 1 rows cols spiralRange :: Integral a => a -> a -> a -> [[a]] spiralRange _ 0 _ = [] spiralRange _ _ 0 = [] spiralRange start 1 cols = [[start..cols+start-1]] spiralRange start rows cols = spiralRange start 1 cols ++ zipWith3 surround [upperLeft,upperLeft-1..lowerLeft+1] innerSpiral [upperRight+1..lowerRight-1] ++ [[lowerLeft,lowerLeft-1..lowerRight]] where upperLeft = 2*cols + 2*rows + start - 5 upperRight = cols+start-1 lowerLeft = 2*cols + rows + start - 3 lowerRight = start+cols+rows-2 innerSpiral = spiralRange (upperLeft+1) (rows-2) (cols-2) surround :: a -> [a] -> a -> [a] surround x1 xs x2 = x1 : xs ++ [x2]

>>2182 https://pastebin.com/YS6iMZTk Took me like 25-30 minutes debugging and all Do I win

>>7015 To be completely honest, I half assed the output formatting, just imagine that I wasn't too lazy to use my homebrew integer to string conversion function to convert the numbers to strings and then pad with spaces

Surprised I haven't seen this one yet. #;2> (define (fizz n) (case (gcd n 15) ((15) 'fizzbuzz) ((3) 'fizz) ((5) 'buzz) (else n))) #;3> (define (fizzbuzz from to) (cond ((< from to) (display (fizz from)) (newline) (fizzbuzz (+ from 1) to))))

>>5829 Improved this and forgot to post it: spiral =: 4 : '(x*y) - |."1|.(|:@|. , {:@{: + >:@i.@#)^:(+:<:y<.x) (,.`,: @. (y>x)) i. >: | x-y' I wish lainchan was more active :(

(honestly not sure if this is the correct thread but whatevs) Are any of you doing the Advent of Code challenge this year?

>>11609 I am. Worked solutions in Awk and Scheme

>>11609 >>11632 And am I, in RetroForth and regular Forth

>>11635 >>11609 >>11632 Was just thinking of starting a new thread around this. See you all there!

here's my non-terminating functional fizzbuzz without modulo: (define (fizzbuzz x y pos) (fizzbuzz (if (= x 3) (begin (display "fizz") 1) (+ x 1)) (if (= y 5) (begin (display "buzz") 1) (+ y 1)) (if (or (= x 3) (= y 5)) (begin (display " ") (+ pos 1)) (begin (display pos) (display " ") (+ pos 1))))) (fizzbuzz 1 1 1)

>>11639 I realised I could unroll the loop in my fizzbuzz implementation to remove all the conditionals, the fizz (x) and buzz (y) counters, and reduce my overall amount of (display) calls. turns out to be pretty close to a 14x improvement in performance. #!r7rs (import (scheme base) (scheme write)) ; srfi-48 is unfortunately poorly optimised and using quasiquote with a for-each ; is 14% slower (mostly constructing the list) and not much better looking, so ; this is a bit ugly. (let fizzbuzz ((pos 1)) (display pos) (display #\space) (display (+ pos 1)) (display " fizz ") (display (+ pos 3)) (display " buzz fizz ") (display (+ pos 6)) (display #\space) (display (+ pos 7)) (display " fizz buzz ") (display (+ pos 10)) (display " fizz ") (display (+ pos 12)) (display #\space) (display (+ pos 13)) (display " fizzbuzz ") (fizzbuzz (+ pos 15)))

Since FizzBuzz is being posted, here's a trigonometric version: a circle with its radius perturbed by a cosine wave with the Fizz pattern, filled with a radial sweep with the color intensity given by a cosine wave with the Buzz pattern. def makePoints (self, cfg): width = cfg ["width" ] height = cfg ["height"] cx, cy = width // 2, height // 2 sin, cos = math.sin, math.cos radians = math.radians ro = round circle = self.makeCircle (cfg) steps, degfun, radfun = circle def p (k): rads = radians (degfun (k)) radi = radfun (k) x = cx + radi * cos (rads) y = cy - radi * sin (rads) return (ro (x), ro (y)) return [ p (k) for k in range (steps) ] def makeCircle (self, cfg): radius = cfg ["radius"] expand = cfg ["expand"] steps = cfg ["steps" ] sides = cfg ["sides" ] sin, cos = math.sin, math.cos radians = math.radians return ( steps, lambda k: 360 * k / steps, lambda k: radius + expand * cos (radians (360 * sides * k / steps)) ) The radial sweep: @ staticmethod def color_radial_slices (ctx, ox, oy, points, closeit, colors): it = iter (points) p = next (it) if closeit: it = itertools.chain (it, [p]) for q, c in zip (it, colors): ctx.move_to (ox, oy) ctx.line_to (* p) ctx.line_to (* q) ctx.close_path () ctx.set_source_rgba (*c) ctx.fill () p = q The color intensity: def inter_colors (colors, n): segments = len (colors) - 1 modf = math.modf for k in range (n): zerotoseg = segments * k / n frac, inte = modf (zerotoseg) inte = int (inte) invfrac = 1 - frac yield tuple (map ( lambda u, v: u * invfrac + v * frac, colors [inte], colors [inte + 1])) class Sources: @ staticmethod def pair_cycle_stops (color1, color2, repeat): ret = [color1, color2] * repeat ret.append (color1) return ret @ staticmethod def inter_buzz (n): colors = Sources.pair_cycle_stops ( (0, 0, 1), (0, 0, 0.1), 3) yield from kit.inter_colors (colors, n)

Monochrome Raster Zigzag Fizzbuzz An image filter that traverses the raster in parallel with the fizzbuzz stream and edits pixels as: fizz -> black, buzz -> white, fizzbuzz -> gray, else -> no-op. To make it slightly less trivial the raster is traversed in zigzag order, starting from the top left and following diagonals of constant sum of coordinates, in alternating direction and in increasing order of sums. The result has 7/15 of pixels encoding the fizzbuzz stream and 8/15 of pixels unmodified from the source, which keeps the original image easily recognizable. class ZigzagFizzbuzz (PrePostEdit): @ staticmethod def lcm2 (a, b): return a * b // fractions.gcd (a, b) @ staticmethod def lcm (ints): return functools.reduce (ZigzagFizzbuzz.lcm2, ints) @ staticmethod def base_loop (data): return tuple (itertools.islice ( map ( lambda tup: functools.reduce (lambda acc, item: acc | (item [1] << item [0]), enumerate (tup), 0), zip (* map (lambda x: itertools.cycle ((1,) + (0,) * (x - 1)), data))), ZigzagFizzbuzz.lcm (data))) def edit_main (self, image): w, h = image.width, image.height setp = image.setpixel data = [3, 5] base = ZigzagFizzbuzz.base_loop (data) fb = itertools.cycle (base) zz = ither.Zigzag.stream (w - 1, h - 1, True) black = ( 0, 0, 0) gray = (128, 128, 128) white = (255, 255, 255) color = (None, black, white, gray) for fbval, pos in zip (fb, zz): if fbval: setp (pos [0], pos [1], color [fbval]) And the zigzag traversal: class Zigzag: @ staticmethod def rangedir (a, b, up): return range (a, b + 1, 1) if up else range (b, a - 1, -1) @ staticmethod def limits4 (maxc, maxr, sum): return ( max (0, sum - maxr), min (maxc, sum), max (0, sum - maxc), min (maxr, sum) ) @ staticmethod def diag (maxc, maxr, sum, up): l4 = Zigzag.limits4 (maxc, maxr, sum) c1, c2, r1, r2 = l4 return ( Zigzag.rangedir (c1, c2, up), Zigzag.rangedir (r1, r2, not up) ) @ staticmethod def stream (maxc, maxr, startdir): dirs = it.cycle ((startdir, not startdir)) sums = range (maxc + maxr + 1) return it.chain.from_iterable (map ( lambda sum, dir: zip (* Zigzag.diag (maxc, maxr, sum, dir)), sums, dirs))

>>2182 Mathlet way in ruby. def print_spiral_like_dumb_coder(rc,cc) m = Array.new(rc) {Array.new(cc,0)} nc = 1 #current north row ec = cc - 1 #current east column sr = rc - 1 #current south row wc = 0 #current west column i = 0 #current row is topmost j = 0 #current column is leftmost d = :e#ast #direction (1..rc*cc).each{|n| #pp "#{n} -> (#{i},#{j})" m[i][j] = n #pp m case(d) when :e#ast j == ec ? (d = :s; i += 1; ec -= 1) : j += 1 when :s#outh i == sr ? (d = :w; j -= 1; sr -= 1) : i += 1 when :w#est j == wc ? (d = :n; i -= 1; wc += 1) : j -= 1 when :n#orth i == nc ? (d = :e; j += 1; nc += 1) : i -= 1 end n += 1 } m.each{|r|r.each{|c| print "%3d " % [c]};puts "\n\n"} nil#hism end print_spiral_like_dumb_coder(15,15)

>>11869 You don't appear to compute the maximum number of digits required by a cell.

>>11871 I left it on purpose since %3d with blank lines creates a more or less square grid for any input up to 37x27=999. If you want you can interpolate in the format string like "#{(rc*cc).digits.count}d " % [c]

Uniform Spiral Card Gradient Fizzbuzz Mapping: fizz -> spade, buzz -> club, fizzbuzz -> diamond, else -> heart. The color interpolation and the fizzbuzz stream are posted above. To make a uniform spiral: def makeCircle (self, cfg): low = cfg ["low" ] high = cfg ["high" ] arc = cfg ["arc" ] turns = cfg ["turns"] pi = math.pi data = [] radius_factor = (high - low) / (turns * 2 * pi) total_angle = 0 while True: radius = low + radius_factor * total_angle if radius > high: # print (len (data), high, radius) break data.append ((total_angle, radius)) angle = arc / radius total_angle += angle return { "mode": "radradlist", "list": data } To consume the spiral points with radial alignment: def draw (self, ctx, width, height): data = self.fosl.radradlist (width, height) count = len (data) fosp = self.fosp fosp.before (ctx) paint = fosp.painter (count) cx = width / 2 cy = height / 2 sin = math.sin cos = math.cos pi = math.pi def p (angle, radius): return ( cx + radius * cos (angle), cy - radius * sin (angle), ) for angle, radius in data: x, y = p (angle, radius) # print (x, y, angle, radius) ctx.save () ctx.translate (x, y) ctx.rotate (pi / 2 - angle) paint (ctx) ctx.restore () fosp.after (ctx) The pointwise card fizzbuzz painter: class CHAR: BLACK_SPADE_SUIT = chr (0x2660) BLACK_CLUB_SUIT = chr (0x2663) BLACK_HEART_SUIT = chr (0x2665) BLACK_DIAMOND_SUIT = chr (0x2666) class FreeOverSpiralPaint: def before (self, ctx): pass def painter (self, count): pass def after (self, ctx): pass @ staticmethod def place (ctx, text, *, color = None, extents = None): if color is not None: ctx.set_source_rgba (* color) if extents is None: extents = ctx.text_extents (text) xb, yb, w, h, xa, ya = extents # print (ext) ctx.move_to (-w / 2 - xb, -h / 2 - yb) ctx.show_text (text) class FOSP_cards (FreeOverSpiralPaint): def before (self, ctx): ctx.set_font_size (24) self.text = [ CHAR.BLACK_HEART_SUIT, CHAR.BLACK_SPADE_SUIT, CHAR.BLACK_CLUB_SUIT, CHAR.BLACK_DIAMOND_SUIT, ] def painter (self, count): colors = Sources.inter_thirds (count) fb = ither.Fizzbuzz.stream_bits ((3, 5)) next (fb) def p (ctx): FreeOverSpiralPaint.place ( ctx, self.text [next (fb)], color = next (colors)) return p

>>2186 You mean you can't write an interpreter without YACC? lmao

>>11886 Different anon, but you don't need a parser generator to parse (YACC). I for one parse bnf and ebnf with parser combinators (the model is a lot easier and it's in language)

>>8799 >>10114 >>11884 I'll assume the person who's been uploading these nice machine-generated images on dithering and whatnot in this thread and others is a single person. I just wanted to let you know that I've been reading your posts (I haven't taken a look at DHALF.TXT yet.) and while I haven't made an effort to understand it much, I like what you're doing.

>>12058 >I'll assume the person who's been uploading these nice machine-generated images on dithering and whatnot in this thread and others is a single person. The first post you linked is by a different Lain. The other two are by this Lain >>9374. >I just wanted to let you know that I've been reading your posts (I haven't taken a look at DHALF.TXT yet.) You can find it in >>8821. >and while I haven't made an effort to understand it much, I like what you're doing. Thanks, Lain.

A cosine loop flower crown. Rather than combining several separate loops, a single cosine wave advances 10/3 phases per turn and loops back on itself after 3 turns. Generating the path: def makeCircle (self, cfg): radius = cfg ["radius" ] expand = cfg ["expand" ] steps = cfg ["steps" ] turns = cfg ["turns" ] turnmul = cfg ["turnmul"] turndiv = cfg ["turndiv"] sin, cos = math.sin, math.cos radians = math.radians return { "mode" : "degradfun", "steps" : steps * turns, "degfun": lambda k: 360 * k / steps, "radfun": lambda k: radius + expand * cos (radians (360 * turnmul * k / (steps * turndiv))) } List of polar coordinates: @ staticmethod def convert_degradfun (degradfun): steps = degradfun ["steps" ] degfun = degradfun ["degfun"] radfun = degradfun ["radfun"] rads = math.radians return [ (rads (degfun (k)), radfun (k)) for k in range (steps) ] The pointwise flower painter uses another cosine wave for the text stamp sizes: class FOSP_flowers2 (FreeOverSpiralPaint): def before (self, ctx): self.text = [ chr (x) for x in Char8.FLOWERS ] def painter (self, count): cos = math.cos radians = math.radians period = 18 def sizes (): for k in range (count): yield 24 + 8 * cos (radians (360 * k / period)) sizeit = sizes () colors = Sources.inter_thirds (count // 3) colors = iter (3 * list (colors)) def p (ctx): FreeOverSpiralPaint.place ( ctx, random.choice (self.text), color = next (colors), size = next (sizeit)) return p And some ancillary bits: class FreeOverSpiralPaint: @ staticmethod def place (ctx, text, *, color = None, size = None, extents = None): if color is not None: ctx.set_source_rgba (* color) if size is not None: ctx.set_font_size (size) if extents is None: extents = ctx.text_extents (text) xb, yb, w, h, xa, ya = extents ctx.move_to (-w / 2 - xb, -h / 2 - yb) ctx.show_text (text) class Char8: FLOWERS = [ 0x273B, 0x273C, 0x273D, 0x273E, 0x273F, 0x2740, 0x2741, 0x2743, ] class Sources: @ staticmethod def inter_thirds (n): colors = [ (1, 0, 0), (1, 1, 0), (0, 1, 0), (0, 1, 1), (0, 0, 1), (1, 0, 1), (1, 0, 0), ] yield from kit.inter_colors (colors, n) The color interpolation is posted in >>11758. The path walk with radial alignment is in >>11884.

While we wait for the answer to >>13856 here is the animated alpha fadeout version of the flower crown to show the triple weave structure. class FOSP_flowers2 (FreeOverSpiralPaint): def before (self, ctx): self.text = [ chr (x) for x in kcha.CHAR8.FLOWERS ] def alphas (self, count): for k in range (count): yield 1 def colors (self, count): col = psou.Source_interthirds.make (count // 3) col = iter (3 * list (col)) return col def sizes (self, count): cos = math.cos radians = math.radians period = 18 for k in range (count): yield 24 + 8 * cos (radians (360 * k / period)) def indices (self, count): nchars = len (self.text) r = random.randrange for k in range (count): yield r (0, nchars) def painter (self, count): colorit = self.colors (count) sizeit = self.sizes (count) indexit = self.indices (count) alphait = self.alphas (count) def p (ctx): FreeOverSpiralPaint.place ( ctx, self.text [next (indexit)], color = next (colorit) + (next (alphait),), size = next (sizeit)) return FreeOverSpiralPaint.adapter_context (p) class FOSP_flowers2b (FOSP_flowers2): def alphas (self, count): time = 0 for k in range (count): yield max (0, (k + time) % 60 - 29) / 30 def indices (self, count): nchars = len (self.text) for k in range (count): yield k % nchars APNG chunk assembly is in >>13347.

This still alive? C?W? write a program which renders a rotating 3D ring in the terminal using ASCII or ANSI art. The program is not allowed to rely on variables external to the inner loop apart from j and tick . uint8_t not_allowed = 0; uint64_t tick = 0; for ( ; ; ) { for (uint16_t j = 0; j<512 ; j++) { for (uint16_t i = 0; i<512 ; i++) { uint8_t allowed = i + j; printf("\x1b[48;2;%d;%d;%dm", allowed, i, j); } puts(""); } tick += 1; } You'll probably need to resize your terminal. C?+1W?: Draw a visual representation of the audio challenge with the same constraints (well, they work the same anyway). Fun task: Implement all of the above with the linux framebuffer device and a higher resolution if you have one. Post a picture along with the code as this will likely not run well on other resolutions. Note, the smaller and simpler the program the better. Using tables is only going to be ugly.

Since this is the challenge thread as well as the de facto fizzbuzz thread, here is a fizzbuzz and a small challenge. A cleartext input is selected, in this case a png image. The standard fizzbuzz [1...100] stream is iterated and for each step a transformation is applied. The input for the first step is the original cleartext, while for every other step the input is the output of the previous step. The transformations are: - Fizz -> every byte is xored with 0xFF. - Buzz -> every byte is decremented by 1 modulo 256, e.g. 0 -> 255, 1 -> 0, 255 -> 254. - FizzBuzz -> every byte is incremented by 1 modulo 256, e.g. 0 -> 1, 254 -> 255, 255 -> 0. - else -> every byte has its low and high nibbles swapped, e.g. 0xC4 -> 0x4C. These transformations are obviously bijections, and their inverses form the same set. The challenge file is the output of step 100. Write a small program to reverse the effect of this chain of transformations and recover the original image. For example, you might apply the inverse transformations in reverse order. Once you have the goal image, try viewing it in a browser as well in addition to a regular image viewer. A backup of the challenge file: http://0x0.st/z6Fv.dat The sha1sum of the challenge file: 9b3a71388720e5fd0cbac43730f453543b70e364 The sha1sum of the goal image: 5703b90dd58a5fe001c69f2896b43e08c8ba22e1

>>2067 This. To memoize you need to do so explicitly. Compare fibonacci :: Int -> Int fibonacci n = fibs!!n where fibs :: [Int] fibs = 1:1:(map (\n -> fibs!!(n-1) + fibs!!(n-2) ) [2..]) Against the naive fibonacci' :: Int -> Int fibonacci' 0 = 1 fibonacci' 1 = 1 fibonacci' n = fibonacci' (n-1) + fibonacci' (n-2)

>>2103 What

>>13899 Not to derail the thread, but why does the png have animations in it? And how did you do that?

>>15984 It's an APNG, as in >>13347. You will find something similar in the challenge >>15979.

Since there seems to be some slight confusion about this >>2103 Lain's method, here are some explanations. Given a set of N points (xi, yi) with no x collisions, we can find an interpolation polynomial P of degree at most N-1 such that P(xi)=yi for all i. For fizzbuzz the x set is {1,2...100} and the y set is at the implementor's discretion. The saner path to take would be to map each input to a 2-bit selector into a 4-element table of printf format strings, but this was too boring for Lain. Instead, each input is mapped to the uint64le reinterpretation of the zero-padded ascii encoding of the final fizzbuzz output string. So Lain's program assumes both ascii and little endian as environment characteristics, thus forfeiting portability. To see the generated values, Lain's polynomial can be extracted: $ grep -e 'polynomial=' 2103.txt | tr -d ' ' | sed -r -e 's/^[^"]*"([^"]+)"[^"]*$/\1/' > 2103-poly.txt parsed: class Parse: @ staticmethod def line (s, *, debug = False): data = [] patt = "([+-]?)([^+-]+)" expos = {} for k, mobj in enumerate (re.finditer (patt, s)): sign, rest = mobj.groups () if not sign: sign = "+" pmult = rest.partition ("*") rest = pmult [0] if pmult [1]: ppower = pmult [2].partition ("^") if ppower [1]: expo = int (ppower [2]) else: expo = 1 else: expo = 0 pdiv = rest.partition ("/") rest = pdiv [0] if pdiv [1]: denom = int (pdiv [2]) else: denom = 1 numer = int (rest) if debug: print (k, sign, expo, numer, denom) if expo in expos: raise ValueError ( "expo clash".format (expo)) expos [expo] = len (data) data.append ((expo, sign, numer, denom)) data.sort (reverse = False) return data @ staticmethod def tofractions (data): ff = fractions.Fraction return [ (expo, ff (numer if sign == "+" else -numer, denom)) for expo, sign, numer, denom in data ] evaluated: class Parse: @ staticmethod def eval (data, x): return sum ( coeff * x ** expo for expo, coeff in data ) and tabulated with a bit of postprocessing: def work_eval (src): with open (src, "r") as f: s = f.readline ().rstrip ("\n") # print (ord (s [-1])) poly = Parse.line (s) poly = Parse.tofractions (poly) for k in range (1, 101): y = Parse.eval (poly, k) print ("{:3d} {} {}".format (k, y, "!" if y.denominator != 1 else struct.pack ("<Q", y.numerator).rstrip (b'\x00'))) yielding: 1 49 b'1' 2 50 b'2' 3 2054842694 b'Fizz' 4 52 b'4' 5 2054845762 b'Buzz' 6 2054842694 b'Fizz' 7 55 b'7' 8 56 b'8' 9 2054842694 b'Fizz' 10 2054845762 b'Buzz' 11 12593 b'11' 12 2054842694 b'Fizz' 13 13105 b'13' 14 13361 b'14' 15 8825495348169042246 b'FizzBuzz' 16 13873 b'16' 17 14129 b'17' 18 2054842694 b'Fizz' 19 14641 b'19' 20 2054845762 b'Buzz' [...] 90 8825495348169042246 b'FizzBuzz' 91 12601 b'91' 92 12857 b'92' 93 2054842694 b'Fizz' 94 13369 b'94' 95 2054845762 b'Buzz' 96 2054842694 b'Fizz' 97 14137 b'97' 98 14393 b'98' 99 2054842694 b'Fizz' 100 2054845762 b'Buzz' The second column is the value of the polynomial when given the first column entry as argument. The third column is the byte view of the uint64le reinterpretation, with trailing null bytes removed. This is how the printf line outputs the fizzbuzz stream.

As you know, the interpolation polynomial is uniquely determined, so we can recover its coefficients from the set of interpolation points. Since for fizzbuzz the number of points and therefore the degree of the polynomial are both relatively small, with a bit of patience we can attempt this in vanilla python. To generate the set of interpolation points from a fizzbuzz stream: class Data: @ staticmethod def fizzbuzz (): data = ["", "Fizz", "Buzz", "FizzBuzz"] return list (itertools.starmap ( lambda k, bits: (k, struct.unpack ("<Q", (data [bits] or str (k)).encode ("ascii").ljust (8, b'\x00')) [0]), itertools.islice ( enumerate (ither.Fizzbuzz.stream_bits ((3, 5))), 1, 101))) Checking against the previous evaluation def work_values (): data = Data.fizzbuzz () print ("\n".join (itertools.starmap ( "{:3d} {}".format, data))) yields 1 49 2 50 3 2054842694 4 52 5 2054845762 6 2054842694 7 55 8 56 9 2054842694 10 2054845762 [...] 95 2054845762 96 2054842694 97 14137 98 14393 99 2054842694 100 2054845762 a match: $ diff <(sed -e "s/ b'.*'\$//" 2103-eval.txt) <(python3 -m flake.poly values) $ Since all methods of finding the coefficients must yield the same result, we are free to choose. For fizzbuzz the set of x values is evenly spaced, so Newton's method readily suggests itself: class Inter: @ staticmethod def polate (data, *, tofractions = False): n = len (data) - 1 ff = fractions.Fraction conv = ff if tofractions else lambda x: x x = [conv (tup [0]) for tup in data] c = [conv (tup [1]) for tup in data] for diff in range (1, n + 1): for index in range (0, n - diff + 1): c [index] = (c [index + 1] - c [index]) / (x [index + diff] - x [index]) # print (c) for diff in range (0, n): for index in range (1, n - diff + 1): c [index] = c [index] - x [index + diff] * c [index - 1] # print (c) return [ (n - k, coeff) for k, coeff in enumerate (c) ] The coefficients can be output in the same format as the debug parse: def work_inter (): data = Data.fizzbuzz () starmap = itertools.starmap poly = Inter.polate (data, tofractions = True) # poly.sort (reverse = False) print ("\n".join ( starmap ("{} {} {} {} {}".format, starmap (lambda k, expo, coeff: (k, "-" if coeff.numerator < 0 else "+", expo, abs (coeff.numerator), coeff.denominator), starmap (lambda k, tup: (k, * tup), enumerate (poly)))))) after which checking for differences 0 - 99 74101463560860539810482394216134472786413399 404009590666424903383979388988167534591844018460526499864038804741201731572423877094984692537474105135297393596654648304117684895744000000000000000000000 1 + 98 1786563401621773217421750502452955853226339781 1943688752347061390850759947022111850270039951356484879070977067483444756705819339975871373032521468004867185688372878439054154137600000000000000000000 2 - 97 27321291157050372775340569532625689973429185264741 12024094960310264981666053243695462339042976739896622019763059664916718201560234437350734896948634081407660523709959770955883479040000000000000000000000 3 + 96 4936870031754926645682423836151042176171669450909 1336493173680525187613977630110369004256312194947800263402124063124652591386915768177479078216982141485276408003996973457735680000000000000000000000 4 - 95 24473118674386691114350902920738421254018653211816783 5509321860394164940053174453010521117545464714729265530246533637991623460050508333264941533538893049900417192994254190586888192000000000000000000000 5 + 94 70955958411301680398144914423031079100489676862292337 16836637264362473501011399007502441445809569586608960873418948299465063205587292020164507233225819167777166294379724181759590400000000000000000000 6 - 93 195447077296046760793596568618244974213518999474673944569 59491157507055852416187694061463793862473023235358410381001603033824121509104033353924740603723449187054576952025621158952960000000000000000000000 7 + 92 3768626467981309940958736490271243883241827877353984433113 1741377967649783015766545405622297430860918840222558009957633541611339279537692732582026457337185788629435048177571945512960000000000000000000000 8 - 91 191509882478118568856359420149388488822510646896295367474063 155750575988552021285951887210938403754020070181396492816086478255921028728837113970069447115872517740769345923956745437184000000000000000000000 9 + 90 3183866784865674882357102610975198211538996328248908077727 5202263802683857887235775650854684650590202417629062186983081540997395662140923677145844788265223211889820833159315456000000000000000000000 10 - 89 122602932492183703177160295045719667527195624761716927797149793 453818694605338057359999671360543134481410460901504932447804423822613720072369213199503051036924585491752173438102405120000000000000000000000 [...] 95 + 4 4353021818054714966071462648760277711153534016644344923943146869615792271408761430298366331628293793025968463884636790147283786917 496033364224659187070987762864106829927030999109098167982161569347924055040000000 96 - 3 17934197811098234184505120115663766567573276672505192992389452402896782102712277769577690491407663682938300855043475172034999 2228626918501922802407296720467981544473059155918725971419226798581248000000 97 + 2 1303108301384808337973620234796701312222070184700563764341592967729999824895692746034116894896125573186790521 247498906860949796715298301379869660653853847832267818240000 98 - 1 151468218170625212862842711493259019240062636816524589253941649490484680567005204702903861 69720375229712477164533808935312303556800 99 + 0 420738786047480639992166722750192850597637924243 1 yields a match: $ diff 2103-parse.txt 2103-inter.txt $ This explains the genesis of Lain's magic constants. The five-star generalship seen in Lain's code is merely a historical artifact of Lain turning up the volume of Pretty-Deep-Parsing to 11.

Here are Lain's >>2103 magic constants >>16020 reobtained with Lagrange interpolation. class LagrangeHelper: @ staticmethod def multiply_xma (poly_n, a): p = poly_n n = len (p) - 1 for k in range (n): p [k] = p [k + 1] - a * p [k] p [n] = - a * p [n] @ staticmethod def add_times (poly_addend, factor, poly_target): for k in range (len (poly_target)): poly_target [k] += factor * poly_addend [k] class Inter: @ staticmethod def split (data, *, tofractions = False): ff = fractions.Fraction conv = ff if tofractions else lambda x: x x = [conv (tup [0]) for tup in data] c = [conv (tup [1]) for tup in data] return (x, c) @ staticmethod def lagrange (data, *, tofractions = False): n = len (data) x, y = Inter.split (data, tofractions = tofractions) c = [1] * n for k in range (n): for other in range (k): c [k] *= x [k] - x [other] for other in range (k + 1, n): c [k] *= x [k] - x [other] c [k] = y [k] / c [k] # print (c) out = [0] * n mulxma = LagrangeHelper.multiply_xma addtim = LagrangeHelper.add_times for k in range (n): term = [0] * n term [-1] = 1 for other in range (k): mulxma (term, x [other]) for other in range (k + 1, n): mulxma (term, x [other]) addtim (term, c [k], out) # print (out) ret = [ (n - 1 - k, coeff) for k, coeff in enumerate (out) ] # return Parse.normalize (ret) return ret The match is just a sanity check: def work_lagrange (): data = Data.fizzbuzz () poly = Inter.lagrange (data, tofractions = True) poly2 = Inter.polate (data, tofractions = True) print (poly == poly2) $ python3 -m flake.poly lagrange True but the fact that this naive implementation is cubic means that it is much slower than the Newton above.

Here is the marginally less naive version of Lagrange >>16071 interpolation. class LagrangeHelper: @ staticmethod def divide_xma (poly_dividend, a, poly_target): p = poly_dividend t = poly_target n = len (p) t [0] = 1 for k in range (1, n): t [k] = p [k - 1] + a * t [k - 1] diff = p [n - 1] + a * t [n - 1] if diff != 0: print ("[divide_xma]", a, diff) class Inter: @ staticmethod def lagrange_basis_coeff (n, x, y): c = [1] * n for k in range (n): for other in range (k): c [k] *= x [k] - x [other] for other in range (k + 1, n): c [k] *= x [k] - x [other] c [k] = y [k] / c [k] # print (c) return c @ staticmethod def lagrange2 (data, *, tofractions = False): n = len (data) x, y = Inter.split (data, tofractions = tofractions) c = Inter.lagrange_basis_coeff (n, x, y) out = [0] * n term = [0] * n product = [0] * n product [-1] = 1 mulxma = LagrangeHelper.multiply_xma addtim = LagrangeHelper.add_times divxma = LagrangeHelper.divide_xma for k in range (n): mulxma (product, x [k]) for k in range (n): divxma (product, x [k], term) addtim (term, c [k], out) # print (out) ret = [ (n - 1 - k, coeff) for k, coeff in enumerate (out) ] # return Parse.normalize (ret) return ret The basis polynomials are computed downward from the common product rather than from the ground up, which brings the order down to quadratic. This means that this version is now competitive with Newton >>16020 and loses on the multiplicative constant and lower terms. def work_lagrange2 (): data = Data.fizzbuzz () names = ["lagrange", "lagrange2", "newton"] n = len (names) t = [None] * (n + 1) poly = [None] * n for k in range (n): t [k] = taim.point () poly [k] = getattr (Inter, names [k]) (data, tofractions = True) t [-1] = taim.point () pairs = ither.pairs pos = max (len (s) for s in names) print ("match:", all (pa == pb for pa, pb in pairs (poly))) print ("\n".join (map ( lambda s, tt: "{:>{pos}}: {}".format ( s, taim.diffstr (* tt), pos = pos), names, pairs (t)))) $ python3 -m flake.poly lagrange2 match: True lagrange: 00:00:03.756 lagrange2: 00:00:00.248 newton: 00:00:00.108

>>16095 >>16071 perfect, can you guys explain to me how Chebyshev Approximation works, in context of GSL library?

>>16098 Function approximation is a quite distant relative of the interpolation posts you link to, so I'm not sure how you arrived at that "perfect". As for "how Chebyshev Approximation works, in context of GSL library", you have here: https://www.gnu.org/software/gsl/doc/html/cheb.html a much better presentation than I could give, including example code.

The scipy version of Lagrange lives in scipy.interpolate at: https://github.com/scipy/scipy/blob/77e5b9b/scipy/interpolate/interpolate.py#L37 def lagrange(x, w): [...] M = len(x) p = poly1d(0.0) for j in xrange(M): pt = poly1d(w[j]) for k in xrange(M): if k == j: continue fac = x[j]-x[k] pt *= poly1d([1.0, -x[k]])/fac p += pt return p It computes the basis polynomials from the ground up. The multiplication with poly1d in the inner loop is linear, so the entire thing is cubic, the equivalent of >>16071. A bit disappointing. Also, while the use of 0.0 and 1.0 makes sense to get the default float dtype, it hampers the use of fractions.Fraction, along with other forced conversions to float throughout numpy.poly1d. >>> import fractions >>> ff = fractions.Fraction >>> ff (1, 3) Fraction(1, 3) >>> ff (1, 3) * 1 Fraction(1, 3) >>> ff (1, 3) * 2 Fraction(2, 3) >>> ff (1, 3) * 1.0 0.3333333333333333 >>> ff (1, 3) + 0.0 0.3333333333333333 >>>

The use of floats in scipy's >>16134 lagrange, which is a problem for fractions.Fraction, is compounded by another issue. When passing in fractions for the interpolation points, the division by 'fac' does not consider 'fac' to be a scalar in the numpy version from the Ubuntu LTS, and attempts to treat it as a polynomial and perform polynomial division with it, which inevitably breaks down: >>> import fractions >>> ff = fractions.Fraction >>> import scipy.interpolate as si >>> si.lagrange ([ff (0), ff (1)], [ff (1), ff (0)]) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/usr/lib/python3/dist-packages/scipy/interpolate/interpolate.py", line 74, in lagrange pt *= poly1d([1.0, -x[k]])/fac File "/usr/lib/python3/dist-packages/numpy/lib/polynomial.py", line 1232, in __div__ return polydiv(self, other) File "/usr/lib/python3/dist-packages/numpy/lib/polynomial.py", line 899, in polydiv q = NX.zeros((max(m - n + 1, 1),), w.dtype) AttributeError: 'float' object has no attribute 'dtype' >>> The decision is made in https://github.com/numpy/numpy/blob/v1.17.0/numpy/lib/polynomial.py @set_module('numpy') class poly1d(object) [...] def __div__(self, other): if isscalar(other): return poly1d(self.coeffs/other) else: other = poly1d(other) return polydiv(self, other) This has been fixed in commit https://github.com/numpy/numpy/commit/cd76e630972c67364dfb12f2fc264becd8e40d12 """PEP 3141 numbers should be considered scalars""" @@ -2162,11 +2163,19 @@ def isscalar(num): [...] - if isinstance(num, generic): - return True - else: - return type(num) in ScalarType + return (isinstance(num, generic) + or type(num) in ScalarType + or isinstance(num, numbers.Number)) but in Ubuntu LTS we have: >>> import fractions >>> ff = fractions.Fraction >>> import numpy >>> numpy.isscalar (ff (1)) False >>> import numbers >>> isinstance (ff (1), numbers.Number) True >>> We have to take another route to get scipy's lagrange to use fractions.

A convenient interpolation point list specification is chosen as: data = "0=1/2 4=3/2" It can be parsed by: class Data: @ staticmethod def pairs (data, convert): return list ( map (lambda p: (convert (p [0]), convert (p [1])), map (lambda s: s.split ("="), data.split ()))) @ staticmethod def split (data): return ( [p [0] for p in data], [p [1] for p in data] ) @ staticmethod def tofloat (s): p = s.split ("/") return ( float (p [0]) if len (p) == 1 else int (p [0]) / int (p [1]) ) By picking either fractions.Fraction or Data.tofloat as the converter, the same specification can be parsed into fractions or floats. Giving the float version to scipy's >>16134 lagrange def work_fake (): ff = fractions.Fraction flo = Data.tofloat data = "0=1/2 4=3/2" xy = Data.pairs (data, flo) x, y = Data.split (xy) print ("data", xy) p = sip.lagrange (x, y) print (p) results in $ python3 -m numsci.interp fake data [(0.0, 0.5), (4.0, 1.5)] 0.25 x + 0.5 which is correct. Even though fractions cause >>16184 to happen, the same result will be possible with fractions after a bit of preparatory work.

One thing to note about scipy's >>16134 lagrange is that it uses the name 'poly1d' brought into its module by an import statement. We can therefore assign something else to 'poly1d' in lagrange's module for the duration of one lagrange run, and the function will use our substitute. This is obviously not threadsafe in principle but for our purposes it will do. To accomplish this: class Patch: @ staticmethod def run (replace, action): old = sip.interpolate.poly1d print ("[patch]", old, "->", replace) sip.interpolate.poly1d = replace try: ret = action () finally: sip.interpolate.poly1d = old return ret The numpy reference page on polynomials states: https://docs.scipy.org/doc/numpy/reference/routines.polynomials.html """Polynomials in NumPy can be created, manipulated, and even fitted using the Using the Convenience Classes of the numpy.polynomial package, introduced in NumPy 1.4. Prior to NumPy 1.4, numpy.poly1d was the class of choice and it is still available in order to maintain backward compatibility. However, the newer Polynomial package is more complete than numpy.poly1d and its convenience classes are better behaved in the numpy environment. Therefore Polynomial is recommended for new coding.""" The bit about >are better behaved in the numpy environment is newspeak for not clobbering your datatype, unlike poly1d. We can therefore try substituting the new Polynomial class https://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.Polynomial.html#numpy.polynomial.polynomial.Polynomial for poly1d as used by lagrange, but since it is not a drop-in replacement we start with a more talkative subclass that shows the arguments it receives. class PolyTalk (numpy.polynomial.polynomial.Polynomial): def __init__ (self, coef, * rest): print ("[talk]", type (coef).__name__, coef, "and", len (rest)) super ().__init__ (coef, * rest) Running this through scipy's >>16134 lagrange with the fractions parse of >>16222 def work_fake (): ff = fractions.Fraction flo = Data.tofloat data = "0=1/2 4=3/2" xy = Data.pairs (data, ff) x, y = Data.split (xy) print ("data", xy) p = Patch.run (PolyTalk, lambda: sip.lagrange (x, y)) print (p) results in $ python3 -m numsci.interp fake data [(Fraction(0, 1), Fraction(1, 2)), (Fraction(4, 1), Fraction(3, 2))] [patch] <class 'numpy.lib.polynomial.poly1d'> -> <class '__main__.PolyTalk'> [talk] float 0.0 and 0 [talk] Fraction 1/2 and 0 [talk] list [1.0, Fraction(-4, 1)] and 0 [talk] ndarray [-0.25 Fraction(1, 1)] and 2 [talk] ndarray [0.0] and 2 [talk] ndarray [-0.125 Fraction(1, 2)] and 2 [talk] ndarray [-0.125 Fraction(1, 2)] and 2 [talk] Fraction 3/2 and 0 [talk] list [1.0, Fraction(0, 1)] and 0 [talk] ndarray [0.25] and 2 [talk] ndarray [0.0] and 2 [talk] ndarray [0.375] and 2 [talk] ndarray [0.25 Fraction(1, 2)] and 2 poly([0.25 Fraction(1, 2)]) This is both floaty and backwards. The former stems from the hardcoded 0.0 and 1.0, while the latter is explained by the same numpy reference page on polynomials: """Transition notice The various routines in the Polynomial package all deal with series whose coefficients go from degree zero upward, which is the reverse order of the Poly1d convention. The easy way to remember this is that indexes correspond to degree, i.e., coef[i] is the coefficient of the term of degree i.""" So at a minimum we have to handle the hardcoded floats and reversed coefficients.

A Polynomial subclass that passes its arguments to a processor, and interprets its decision to proceed with the modified parent call or complain: class Fix: # processor # > coef, rest, what # < handled, decision, coef, rest @ staticmethod def poly (processor, debug): nppp = numpy.polynomial.polynomial.Polynomial class PolyFix (nppp): def __init__ (self, coef, * rest): what = type (coef) if debug: print ("[fix]", what.__name__, coef, "and", len (rest), end = "") ret = processor (coef, rest, what) handled, decision, pcoef, prest = ret if debug: print ("", decision) if handled: nppp.__init__ (self, pcoef, * prest) else: raise NotImplementedError ( "{} {} and {}".format ( what.__name__, coef, len (rest))) return PolyFix A processor that implements the strategy >>16246 of fixing floats and coefficient order: class Fix: @ staticmethod def nofloat (coef, rest, what): if len (rest) > 0: return (True, "leave", coef, rest) if what is list: if (len (coef) == 2) and (type (coef [0]) is float) and (coef [0] == 1.0): return (True, "fix", [coef [1], C.FF_ONE], rest) else: return (False, "what", coef, rest) if what is float: if coef == 0.0: return (True, "fix", [C.FF_ZERO], rest) else: return (False, "what", coef, rest) return (True, "leave", coef, rest) class C: FF_ZERO = fractions.Fraction (0) FF_ONE = fractions.Fraction (1) Running this through scipy's >>16134 lagrange with the fractions parse of >>16222 def work_fake (): ff = fractions.Fraction flo = Data.tofloat data = "0=1/2 4=3/2" xy = Data.pairs (data, ff) x, y = Data.split (xy) print ("data", xy) p = Patch.run ( Fix.poly (Fix.nofloat, True), lambda: sip.lagrange (x, y)) print (p) results in $ python3 -m numsci.interp fake data [(Fraction(0, 1), Fraction(1, 2)), (Fraction(4, 1), Fraction(3, 2))] [patch] <class 'numpy.lib.polynomial.poly1d'> -> <class '__main__.Fix.poly.<locals>.PolyFix'> [fix] float 0.0 and 0 fix [fix] Fraction 1/2 and 0 leave [fix] list [1.0, Fraction(-4, 1)] and 0 fix [fix] ndarray [Fraction(1, 1) Fraction(-1, 4)] and 2 leave [fix] ndarray [Fraction(0, 1)] and 2 leave [fix] ndarray [Fraction(1, 2) Fraction(-1, 8)] and 2 leave [fix] ndarray [Fraction(1, 2) Fraction(-1, 8)] and 2 leave [fix] Fraction 3/2 and 0 leave [fix] list [1.0, Fraction(0, 1)] and 0 fix [fix] ndarray [Fraction(0, 1) Fraction(1, 4)] and 2 leave [fix] ndarray [Fraction(0, 1)] and 2 leave [fix] ndarray [Fraction(0, 1) Fraction(3, 8)] and 2 leave [fix] ndarray [Fraction(1, 2) Fraction(1, 4)] and 2 leave poly([Fraction(1, 2) Fraction(1, 4)]) which is correct and exceedingly fractious. We can in fact make lagrange give us the full precision interpolating polynomial with fractions if we cheat enough, and so are now ready to try for the full >>2103 fizzbuzz polynomial.

The classes from >>16020 can provide inputs usable by scipy's >>16134 lagrange. These are combined with the patched run of >>16279. def work_fb (): data = fpoly.Data.fizzbuzz () t1 = taim.point () lag2 = fpoly.Inter.lagrange2 (data, tofractions = True) t2 = taim.point () x, y = fpoly.Inter.split (data, tofractions = True) t3 = taim.point () scil = Patch.run ( Fix.poly (Fix.nofloat, False), lambda: sip.lagrange (x, y)) t4 = taim.point () scil = list (scil.coef) [::-1] scil = fpoly.Inter.withexpodown (scil) print ("match:", lag2 == scil) print ("lagrange2:", taim.diffstr (t1, t2)) print ("sci.lagra:", taim.diffstr (t3, t4)) The result is a match of the full >>2103 fizzbuzz polynomial. $ python3 -m numsci.interp fb [patch] <class 'numpy.lib.polynomial.poly1d'> -> <class '__main__.Fix.poly.<locals>.PolyFix'> match: True lagrange2: 00:00:00.249 sci.lagra: 00:00:08.037 The reason scipy's cubic is more than twice as slow as the cubic of >>16071 is that scipy uses the generic strategy of dividing early, without any special information about the values. But knowing that our interpolation points are all integral, we can delay the division to keep computing with integers for as long as possible, thus saving on gcd calls. Scipy ends up doing so many extra gcd calls that they outweigh any gain from reducing the numerators. It will be seen in the profiling run that the overhead added by the nofloat processor is negligible.

The cProfile module gained the ability to profile modules from the command line in python 3.7. https://docs.python.org/3/whatsnew/3.7.html#cprofile In previous versions a top-level wrapper script is needed, such as the very basic: import importlib import sys def main (): argv = sys.argv # print (argv) name = argv [1] mod = importlib.import_module (name) del argv [1] mod.main () if __name__ == "__main__": main () A profiling run adds overhead, but still indicates the relative weights that various functions contribute to the overall time. def work_profile (): data = fpoly.Data.fizzbuzz () x, y = fpoly.Inter.split (data, tofractions = True) scil = Patch.run ( Fix.poly (Fix.nofloat, False), lambda: sip.lagrange (x, y)) $ python3 -m cProfile -s time wrapper.py numsci.interp profile | tee numsci/profile.txt | less [patch] <class 'numpy.lib.polynomial.poly1d'> -> <class 'numsci.interp.Fix.poly.<locals>.PolyFix'> 17942930 function calls (17938946 primitive calls) in 11.080 seconds Ordered by: internal time ncalls tottime percall cumtime percall filename:lineno(function) 1534503 3.144 0.000 3.144 0.000 {built-in method math.gcd} 1544602 2.006 0.000 5.335 0.000 fractions.py:84(__new__) 999900 0.987 0.000 4.401 0.000 fractions.py:417(_mul) 495001 0.690 0.000 2.963 0.000 fractions.py:401(_add) 1524601 0.579 0.000 8.230 0.000 fractions.py:374(forward) 9900 0.526 0.000 8.494 0.001 {built-in method numpy.core.multiarray.correlate} 119301 0.404 0.000 1.331 0.000 polyutils.py:124(as_series) 238697 0.273 0.000 0.273 0.000 {built-in method numpy.core.multiarray.array} 3039302 0.227 0.000 0.227 0.000 fractions.py:274(numerator) 3039302 0.217 0.000 0.217 0.000 fractions.py:278(denominator) 1591357 0.190 0.000 0.207 0.000 {built-in method builtins.isinstance} 1544780 0.186 0.000 0.186 0.000 {built-in method __new__ of type object at 0x9d12c0} 119301 0.184 0.000 0.184 0.000 polyutils.py:173(<listcomp>) 9900 0.161 0.000 0.481 0.000 polynomial.py:371(polydiv) 49800 0.110 0.000 0.138 0.000 polyutils.py:94(trimseq) 39801 0.091 0.000 1.055 0.000 _polybase.py:246(__init__) 1 0.082 0.082 10.945 10.945 interpolate.py:43(lagrange) 59600 0.078 0.000 0.140 0.000 type_check.py:571(common_type) 20000 0.058 0.000 0.058 0.000 {method 'reduce' of 'numpy.ufunc' objects} 39801 0.053 0.000 1.161 0.000 interp.py:72(__init__) 119301 0.048 0.000 0.230 0.000 polyutils.py:165(<listcomp>) 19900 0.048 0.000 0.169 0.000 _polybase.py:210(_get_coefficients) 79641 0.041 0.000 0.041 0.000 {method 'copy' of 'numpy.ndarray' objects} 9900 0.033 0.000 1.099 0.000 _polybase.py:364(__divmod__) 39801 0.032 0.000 0.052 0.000 interp.py:91(nofloat) Essentially all of the time is spent on Fraction operations, especially gcd calls. The only user methods that take any time are the PolyFix constructor > 39801 0.053 0.000 1.161 0.000 interp.py:72(__init__) and nofloat > 39801 0.032 0.000 0.052 0.000 interp.py:91(nofloat) The constructor takes 0.053s, with the rest up to 1.161s being the superclass init that would happen anyway, and nofloat. Nofloat itself takes 0.052s with all subcalls. Together the user methods take 0.105s, which out of a total lagrange time of 10.945s is less than one percent. Since the unprofiled scipy >>16305 takes 8.037s, some random Lain's >>16071 cubic takes 3.756s and the quadratic >>16095 takes 0.248s, scipy needs to step up its game with regard to the quality of its lagrange function.

Here is how to produce the >>2103 fizzbuzz polynomial using sympy. https://www.sympy.org/en/index.html Sympy includes both the Lagrange and Newton methods of interpolation. Naturally they will agree on the final result. The Lagrange method is provided by the interpolating_poly function, which lives in the "Polynomials Manipulation Module": https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.specialpolys.interpolating_poly It can produce symbolic Lagrange polynomials: >>> import sympy as s >>> s.init_printing () >>> x = s.symbols ("x") >>> lag = s.polys.specialpolys.interpolating_poly >>> p = lag (4, x) >>> p y₀⋅(x - x₁)⋅(x - x₂)⋅(x - x₃) y₁⋅(x - x₀)⋅(x - x₂)⋅(x - x₃) y₂⋅(x - x₀)⋅(x - x₁)⋅(x - x₃) y₃⋅(x - x₀)⋅(x - x₁)⋅(x - x₂) ───────────────────────────── + ────────────────────────────── + ─────────────────────────────── + ──────────────────────────────── (x₀ - x₁)⋅(x₀ - x₂)⋅(x₀ - x₃) (-x₀ + x₁)⋅(x₁ - x₂)⋅(x₁ - x₃) (-x₀ + x₂)⋅(-x₁ + x₂)⋅(x₂ - x₃) (-x₀ + x₃)⋅(-x₁ + x₃)⋅(-x₂ + x₃) >>> Sympy includes a Poly type that can be queried for its coefficients: https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.Poly https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.Poly.all_coeffs The interpolating_poly function produces a generic expression rather than a Poly: >>> type (p) <class 'sympy.core.add.Add'> >>> There are two conversion routines from generic expressions to Polys, a faster, more specialized one and a slower, more general one: https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.poly https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.poly_from_expr We will use the faster one. Sympy also includes a Rational type compatible with python's built-in Fraction type: https://docs.sympy.org/latest/modules/core.html#sympy.core.numbers.Rational >>> import fractions >>> ff = fractions.Fraction >>> r = s.Rational >>> r (2, 4) 1/2 >>> ff (3, 6) 1/2 >>> r (2, 4) == ff (3, 6) True >>> We can proceed by converting the fizzbuzz data >>16020 to Rationals, generating the Lagrange polynomial, converting it from a generic expression to a Poly and extracting its coefficients. To perform the latter steps: class Inter: @ staticmethod def sympy_lagrange (x, y, *, from_expr = False, domainQ = False): n = len (x) s = sympy z = s.symbols ("x") lag = s.polys.specialpolys.interpolating_poly expr = lag (n, z, x, y) p = Inter.topoly (expr, from_expr = from_expr, domain = 'QQ' if domainQ else None) coef = p.all_coeffs () return coef Poly conversion: class Inter: @ staticmethod def topoly (expr, *, from_expr = False, domain = None): s = sympy if from_expr: p = s.polys.polytools.poly_from_expr (expr) p = p [0] else: fun = s.polys.polytools.poly if domain is not None: p = fun (expr, domain = domain) else: p = fun (expr) return p To coordinate with the cubic Lagrange of >>16071 def work_fb_compare (fpoly_label, fpoly_name, sympy_label, sympy_name): data = fpoly.Data.fizzbuzz () f1 = getattr (fpoly.Inter, fpoly_name) t1 = taim.point () p1 = f1 (data, tofractions = True) t2 = taim.point () x, y = fpoly.Inter.split_convert (data, sympy.Rational) f2 = getattr (Inter, "sympy_" + sympy_name) t3 = taim.point () p2 = f2 (x, y, domainQ = True) t4 = taim.point () p2 = fpoly.Inter.withexpodown (p2) pos = max (map (len, ("match", fpoly_label, sympy_label))) fmt = "{:>{pos}}: {}" print (fmt.format ("match", p1 == p2, pos = pos)) print (fmt.format (fpoly_label, taim.diffstr (t1, t2), pos = pos)) print (fmt.format (sympy_label, taim.diffstr (t3, t4), pos = pos)) def work_fblag (): work_fb_compare ("fplag", "lagrange", "symlag", "lagrange") The match is just a sanity check: $ python3 -m numsci.symterp fblag match: True fplag: 00:00:03.770 symlag: 00:00:02.390 The sympy version has no choice but to be cubic since it has to be able to produce symbolic Lagrange polynomials. The reason it is faster than the cubic Fraction version is that Fractions run in python code while sympy Rationals are backed by the mpq functions from GMP. https://gmplib.org/#FUNCCLASSES The explicit request for the Q domain during Poly conversion is needed because the fizzbuzz polynomial hits a bug in sympy.

Without the explicit request for the Q domain during Poly conversion >>16725 we get: $ python3 -m pdb wrapper.py numsci.symterp fblag > /home/name/src/python/wrapper.py(1)<module>() -> import importlib (Pdb) continue Traceback (most recent call last): File "/usr/lib/python3.6/pdb.py", line 1667, in main pdb._runscript(mainpyfile) File "/usr/lib/python3.6/pdb.py", line 1548, in _runscript self.run(statement) File "/usr/lib/python3.6/bdb.py", line 434, in run exec(cmd, globals, locals) File "<string>", line 1, in <module> File "/home/name/src/python/wrapper.py", line 1, in <module> import importlib File "/home/name/src/python/wrapper.py", line 11, in main mod.main () File "/home/name/src/python/numsci/symterp.py", line 116, in main toyargs.main (globals ()) File "/home/name/src/python/lib/toyargs.py", line 30, in main proc (argv [2:]) File "/home/name/src/python/numsci/symterp.py", line 112, in main_fblag work_fblag () File "/home/name/src/python/numsci/symterp.py", line 109, in work_fblag work_fb_compare ("fplag", "lagrange", "symlag", "lagrange") File "/home/name/src/python/numsci/symterp.py", line 90, in work_fb_compare p2 = f2 (x, y, domainQ = False) File "/home/name/src/python/numsci/symterp.py", line 54, in sympy_lagrange domain = 'QQ' if domainQ else None) File "/home/name/src/python/numsci/symterp.py", line 28, in topoly p = fun (expr) File "/usr/lib/python3/dist-packages/sympy/polys/polytools.py", line 6957, in poly return _poly(expr, opt) File "/usr/lib/python3/dist-packages/sympy/polys/polytools.py", line 6923, in _poly product = product.mul(factor) File "/usr/lib/python3/dist-packages/sympy/polys/polytools.py", line 1389, in mul return f.mul_ground(g) File "/usr/lib/python3/dist-packages/sympy/polys/polytools.py", line 1206, in mul_ground result = f.rep.mul_ground(coeff) File "/usr/lib/python3/dist-packages/sympy/polys/polyclasses.py", line 407, in mul_ground return f.per(dmp_mul_ground(f.rep, f.dom.convert(c), f.lev, f.dom)) File "/usr/lib/python3/dist-packages/sympy/polys/domains/domain.py", line 145, in convert return self.from_sympy(element) File "/usr/lib/python3/dist-packages/sympy/polys/domains/gmpyintegerring.py", line 39, in from_sympy raise CoercionFailed("expected an integer, got %s" % a) sympy.polys.polyerrors.CoercionFailed: expected an integer, got -4952578758792953/18939667865925099543556892782958840305894396259959071067427295492617499220040727470616346835140592271360000000000000000000000 Uncaught exception. Entering post mortem debugging Running 'cont' or 'step' will restart the program > /usr/lib/python3/dist-packages/sympy/polys/domains/gmpyintegerring.py(39)from_sympy() -> raise CoercionFailed("expected an integer, got %s" % a) (Pdb) longlist 32 def from_sympy(self, a): 33 """Convert SymPy's Integer to ``dtype``. """ 34 if a.is_Integer: 35 return GMPYInteger(a.p) 36 elif a.is_Float and int(a) == a: 37 return GMPYInteger(int(a)) 38 else: 39 -> raise CoercionFailed("expected an integer, got %s" % a) (Pdb) whatis a <class 'sympy.core.numbers.Rational'> (Pdb) whatis self <class 'sympy.polys.domains.gmpyintegerring.GMPYIntegerRing'> So the integer domain is asked to convert a Rational to the domain's type, but it can only handle integer types and float types with no fractional part. Walking up the stack will show how we got into this situation. (Pdb) up > /usr/lib/python3/dist-packages/sympy/polys/domains/domain.py(145)convert() -> return self.from_sympy(element) (Pdb) longlist 106 def convert(self, element, base=None): 107 """Convert ``element`` to ``self.dtype``. """ 108 if base is not None: 109 return self.convert_from(element, base) 110 111 if self.of_type(element): 112 return element 113 114 from sympy.polys.domains import PythonIntegerRing, GMPYIntegerRing, GMPYRationalField, RealField, ComplexField 115 116 if isinstance(element, integer_types): 117 return self.convert_from(element, PythonIntegerRing()) 118 119 if HAS_GMPY: 120 integers = GMPYIntegerRing() 121 if isinstance(element, integers.tp): 122 return self.convert_from(element, integers) 123 124 rationals = GMPYRationalField() 125 if isinstance(element, rationals.tp): 126 return self.convert_from(element, rationals) 127 128 if isinstance(element, float): 129 parent = RealField(tol=False) 130 return self.convert_from(parent(element), parent) 131 132 if isinstance(element, complex): 133 parent = ComplexField(tol=False) 134 return self.convert_from(parent(element), parent) 135 136 if isinstance(element, DomainElement): 137 return self.convert_from(element, element.parent()) 138 139 # TODO: implement this in from_ methods 140 if self.is_Numerical and getattr(element, 'is_ground', False): 141 return self.convert(element.LC()) 142 143 if isinstance(element, Basic): 144 try: 145 >> return self.from_sympy(element) 146 except (TypeError, ValueError): 147 -> pass 148 else: # TODO: remove this branch 149 if not is_sequence(element): 150 try: 151 element = sympify(element) 152 153 if isinstance(element, Basic): 154 return self.from_sympy(element) 155 except (TypeError, ValueError): 156 pass 157 158 raise CoercionFailed("can't convert %s of type %s to %s" % (element, type(element), self)) There is a test for rationals in the "if HAS_GMPY:" block but that branch is not taken because a Rational doesn't count as an mpq: (Pdb) p HAS_GMPY 2 (Pdb) p isinstance (element, GMPYRationalField ().tp) False (Pdb) whatis element <class 'sympy.core.numbers.Rational'> (Pdb) p GMPYRationalField ().tp <class 'mpq'>

Further up: (Pdb) up > /usr/lib/python3/dist-packages/sympy/polys/polyclasses.py(407)mul_ground() -> return f.per(dmp_mul_ground(f.rep, f.dom.convert(c), f.lev, f.dom)) (Pdb) longlist 405 def mul_ground(f, c): 406 """Multiply ``f`` by a an element of the ground domain. """ 407 -> return f.per(dmp_mul_ground(f.rep, f.dom.convert(c), f.lev, f.dom)) (Pdb) up > /usr/lib/python3/dist-packages/sympy/polys/polytools.py(1206)mul_ground() -> result = f.rep.mul_ground(coeff) (Pdb) longlist 1191 def mul_ground(f, coeff): 1192 """ 1193 Multiply ``f`` by a an element of the ground domain. [...] 1204 """ 1205 if hasattr(f.rep, 'mul_ground'): 1206 -> result = f.rep.mul_ground(coeff) 1207 else: # pragma: no cover 1208 raise OperationNotSupported(f, 'mul_ground') 1209 1210 return f.per(result) (Pdb) up > /usr/lib/python3/dist-packages/sympy/polys/polytools.py(1389)mul() -> return f.mul_ground(g) (Pdb) longlist 1369 def mul(f, g): 1370 """ 1371 Multiply two polynomials ``f`` and ``g``. [...] 1385 """ 1386 g = sympify(g) 1387 1388 if not g.is_Poly: 1389 -> return f.mul_ground(g) 1390 1391 _, per, F, G = f._unify(g) 1392 1393 if hasattr(f.rep, 'mul'): 1394 result = F.mul(G) 1395 else: # pragma: no cover 1396 raise OperationNotSupported(f, 'mul') 1397 1398 return per(result) So the 'mul' function from Poly considers anything that is not a Poly to be not only a scalar but a scalar from the domain of the Poly's coefficients. If the scalar is not from that domain, multiplication via 'mul' will fail. As such, 'mul' does not function as a general multiplication routine with a scalar, but that is how it is used further up the stack: (Pdb) up > /usr/lib/python3/dist-packages/sympy/polys/polytools.py(6923)_poly() -> product = product.mul(factor) (Pdb) longlist 6896 def _poly(expr, opt): 6897 terms, poly_terms = [], [] 6898 6899 for term in Add.make_args(expr): 6900 factors, poly_factors = [], [] 6901 6902 for factor in Mul.make_args(term): 6903 if factor.is_Add: 6904 poly_factors.append(_poly(factor, opt)) 6905 elif factor.is_Pow and factor.base.is_Add and factor.exp.is_Integer: 6906 poly_factors.append( 6907 _poly(factor.base, opt).pow(factor.exp)) 6908 else: 6909 factors.append(factor) 6910 6911 if not poly_factors: 6912 terms.append(term) 6913 else: 6914 product = poly_factors[0] 6915 6916 for factor in poly_factors[1:]: 6917 product = product.mul(factor) 6918 6919 if factors: 6920 factor = Mul(*factors) 6921 6922 if factor.is_Number: 6923 -> product = product.mul(factor) 6924 else: 6925 product = product.mul(Poly._from_expr(factor, opt)) 6926 6927 poly_terms.append(product) 6928 6929 if not poly_terms: 6930 result = Poly._from_expr(expr, opt) 6931 else: 6932 result = poly_terms[0] 6933 6934 for term in poly_terms[1:]: 6935 result = result.add(term) 6936 6937 if terms: 6938 term = Add(*terms) 6939 6940 if term.is_Number: 6941 result = result.add(term) 6942 else: 6943 result = result.add(Poly._from_expr(term, opt)) 6944 6945 return result.reorder(*opt.get('gens', ()), **args) (Pdb) whatis factor <class 'sympy.core.numbers.Rational'> (Pdb) whatis product <class 'sympy.polys.polytools.Poly'> (Pdb) whatis product.rep <class 'sympy.polys.polyclasses.DMP'> (Pdb) p product.rep.dom ZZ (Pdb) up > /usr/lib/python3/dist-packages/sympy/polys/polytools.py(6957)poly() -> return _poly(expr, opt) (Pdb) longlist 6879 @public 6880 def poly(expr, *gens, **args): 6881 """ 6882 Efficiently transform an expression into a polynomial. [...] 6893 """ 6894 options.allowed_flags(args, []) 6895 6896 def _poly(expr, opt): [...] 6945 return result.reorder(*opt.get('gens', ()), **args) 6946 6947 expr = sympify(expr) 6948 6949 if expr.is_Poly: 6950 return Poly(expr, *gens, **args) 6951 6952 if 'expand' not in args: 6953 args['expand'] = False 6954 6955 opt = options.build_options(gens, args) 6956 6957 -> return _poly(expr, opt) So the '_poly' helper of the faster 'poly' conversion function goes out of its way with the "if factor.is_Number:" block to treat scalars differently, but passes them to a multiplication routine that cannot handle general multiplication with a scalar. This is the source of the bug. The fizzbuzz interpolation >>16020 happens to produce a 'factor' that is a Rational and a 'product' that is a Poly over the integer domain. General multiplication via the * operator works fine. That is backed by the '__mul__' method. The difference between '__mul__' and 'mul' can be demonstrated in a fresh REPL: >>> import sympy as s >>> x = s.symbols ("x") >>> p = s.Poly (x + 1) >>> p Poly(x + 1, x, domain='ZZ') >>> p.mul (s.Rational (1, 2)) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/usr/lib/python3/dist-packages/sympy/polys/polytools.py", line 1389, in mul return f.mul_ground(g) File "/usr/lib/python3/dist-packages/sympy/polys/polytools.py", line 1206, in mul_ground result = f.rep.mul_ground(coeff) File "/usr/lib/python3/dist-packages/sympy/polys/polyclasses.py", line 407, in mul_ground return f.per(dmp_mul_ground(f.rep, f.dom.convert(c), f.lev, f.dom)) File "/usr/lib/python3/dist-packages/sympy/polys/domains/domain.py", line 145, in convert return self.from_sympy(element) File "/usr/lib/python3/dist-packages/sympy/polys/domains/gmpyintegerring.py", line 39, in from_sympy raise CoercionFailed("expected an integer, got %s" % a) sympy.polys.polyerrors.CoercionFailed: expected an integer, got 1/2 >>> p Poly(x + 1, x, domain='ZZ') >>> p * s.Rational (1, 2) Poly(1/2*x + 1/2, x, domain='QQ') >>> The '__mul__' method takes care to try to process non-Polys in a way that includes converting scalars to constant polynomials. This makes general multiplication with scalars work. The price is that '__mul__' and therefore * may produce a general Expr instead of a Poly, whereas 'mul' produces Polys on success. https://github.com/sympy/sympy/blob/7f92c96417d066700d096105282d6fe47ea085ef/sympy/polys/polytools.py#L4038 @_sympifyit('g', NotImplemented) def __mul__(f, g): if not g.is_Poly: try: g = f.__class__(g, *f.gens) except PolynomialError: return f.as_expr()*g return f.mul(g) The sympy people clearly know their math, but they need to step up their programming game.

Sympy's Newton interpolation is provided by the apply_finite_diff function, which lives in the calculus module: https://docs.sympy.org/latest/modules/calculus/index.html#sympy.calculus.finite_diff.apply_finite_diff It can produce symbolic interpolation polynomials, though not quite with the same sign discipline as the Lagrange >>16725 version: >>> import sympy as s >>> s.init_printing () >>> x = s.symbols ("x") >>> newt = s.calculus.finite_diff.apply_finite_diff >>> p = newt (0, s.symbols ("x0:4"), s.symbols ("y0:4"), x) >>> p y₀⋅(-x + x₁)⋅(-x + x₂)⋅(-x + x₃) y₁⋅(-x + x₂)⋅(-x + x₃)⋅(x - x₀) y₂⋅(-x + x₁)⋅(-x + x₃)⋅(x - x₀) y₃⋅(-x + x₁)⋅(-x + x₂)⋅(x - x₀) ──────────────────────────────── + ──────────────────────────────── - ──────────────────────────────── + ──────────────────────────────── (-x₀ + x₁)⋅(-x₀ + x₂)⋅(-x₀ + x₃) (-x₀ + x₁)⋅(-x₁ + x₂)⋅(-x₁ + x₃) (-x₀ + x₂)⋅(-x₁ + x₂)⋅(-x₂ + x₃) (-x₀ + x₃)⋅(-x₁ + x₃)⋅(-x₂ + x₃) >>> It returns a generic expression rather than a Poly, so it needs to be converted as above before it can be queried for its coefficients. class Inter: @ staticmethod def sympy_afd (x, y, *, from_expr = False, domainQ = False): s = sympy z = s.symbols ("x") afd = s.calculus.finite_diff.apply_finite_diff expr = afd (0, x, y, z) p = Inter.topoly (expr, from_expr = from_expr, domain = 'QQ' if domainQ else None) coef = p.all_coeffs () return coef To coordinate with the Newton of >>16020 def work_fbnew (): work_fb_compare ("newton", "newton", "findiff", "afd") The match is just a sanity check: $ python3 -m numsci.symterp fbnew match: True newton: 00:00:00.107 findiff: 00:00:11.716 The sympy version of Newton also has no choice but to be cubic since it has to be able to produce symbolic interpolation polynomials. The way in which it arranges its expression means that it gets a worse constant factor than the Lagrange >>16725 version. Being cubic also means it has no hope of competing with the quadratic >>16020 Newton.

In preparation for parallelization, the lagrange computation of >>16071 can be split into several phases. Some of these can be easily done in parallel while others gather up and combine the partial results of a preceding parallel phase. The multipliers of the basis polynomials can be computed independently from each other, so disjoint chunks of this computation can run at the same time. class LagrangePieces: @ staticmethod def basis_coeff (startstop, nxy): start, stop = startstop n, x, y = nxy todo = stop - start out = [1] * todo for k in range (todo): pos = start + k for other in range (pos): out [k] *= x [pos] - x [other] for other in range (pos + 1, n): out [k] *= x [pos] - x [other] out [k] = y [pos] / out [k] return (start, stop, out) To combine the chunks of the basis multiplier phase: @ staticmethod def basis_coeff_combine (pieces, extra): n = max (p [1] for p in pieces) out = [None] * n for start, stop, coeff in pieces: out [start : stop] = coeff return out The common product for the basis polynomials is kept as a unit, because combining its partial results would be quadratic. @ staticmethod def product (n, x): mulxma = LagrangeHelper.multiply_xma prod = [0] * n prod [-1] = 1 for k in range (n): mulxma (prod, x [k]) return prod Forming the terms of the final sum from the multipliers and the common product can be done in chunks: @ staticmethod def terms (startstop, nxcoefprod): addtim = LagrangeHelper.add_times divxma = LagrangeHelper.divide_xma start, stop = startstop n, x, c, p = nxcoefprod todo = stop - start out = [0] * n term = [0] * n for k in range (todo): pos = start + k divxma (p, x [pos], term) addtim (term, c [pos], out) return out Combining the terms into a plain sum: @ staticmethod def terms_combine (pieces, extra): add = LagrangeHelper.add out = pieces [0] [:] for k in range (1, len (pieces)): add (pieces [k], out) return out These phases will be passed to an Executor. https://docs.python.org/3/library/concurrent.futures.html#executor-objects

An example of a single-phase computation with a ProcessPoolExecutor is given in the documentation: https://docs.python.org/3/library/concurrent.futures.html#processpoolexecutor-example Here is a structure for a multiphase computation. A round function receives the pieces computed by the previous parallel round, along with auxiliary data that does not pass through the Executor and acts as a closure. It returns whether the computation is at its end and if so the result, otherwise the data chunks and processing function for the next parallel round. https://docs.python.org/3/library/concurrent.futures.html#concurrent.futures.Executor.map The Executor and its pool of workers are preserved across rounds because tearing them down and constructing them again would be expensive. class Engine: # roundfun # > count, extra, last # < final, data, itemfun @ staticmethod def rounds (roundfun, extra, initargs, workers): last = initargs count = 0 with cofu.ProcessPoolExecutor (max_workers = workers) as e: while True: final, data, itemfun = roundfun (count, extra, last) if final: ret = data break last = list (e.map (itemfun, * data)) count += 1 return ret For a fixed pipeline computation where the rounds are predetermined, the overall round function can be assembled from individual round functions: class Engine: @ staticmethod def roundfun_list (funlist): limit = len (funlist) - 1 def ret (count, extra, pieces): data, itemfun = funlist [count] (pieces, extra) return (count == limit, data, itemfun) return ret The individual round functions might share a common naming pattern and container: class Engine: @ staticmethod def roundfun_attr (attrof, fmt, count): funlist = [ getattr (attrof, fmt.format (k)) for k in range (count) ] return Engine.roundfun_list (funlist) Simple process count selection: class Workers: @ staticmethod def all (): return os.cpu_count () @ staticmethod def leave1 (): return max (os.cpu_count () - 1, 1) @ staticmethod def half (): return max (os.cpu_count () // 2, 1) @ staticmethod def half_leave1 (): return max (os.cpu_count () // 2 - 1, 1) To partition a range of tasks into disjoint chunks: class Pieces: @ staticmethod def start_stop (start, stop, n): have = stop - start di, mo = divmod (have, n) out = [None] * n pos = start for k in range (n): now = di + (1 if k < mo else 0) out [k] = (pos, pos + now) pos += now return out An example set of individual round functions will be the phases >>17187 of Lagrange interpolation.

The individual round functions >>17263 for Lagrange interpolation will coordinate the flow of arguments between the phases >>17187 of the computation. The first round sets up the basis multiplier phase: class Lagrounds: @ staticmethod def r0 (pieces, extra): n = extra ["n" ] data = extra ["data" ] stst = extra ["chunks"] rep = itertools.repeat x, y = fpoly.Inter.split (data, tofractions = True) extra ["x"] = x extra ["y"] = y return ( [stst, rep ((n, x, y), len (stst))], fpoly.LagrangePieces.basis_coeff ) The next round combines the partial results of the basis multiplier phase, executes the common product phase and sets up the term phase: class Lagrounds: @ staticmethod def r1 (pieces, extra): fplp = fpoly.LagrangePieces c = fplp.basis_coeff_combine (pieces, None) stst = extra ["chunks"] n = extra ["n" ] x = extra ["x" ] p = fplp.product (n, x) rep = itertools.repeat return ( [stst, rep ((n, x, c, p), len (stst))], fplp.terms ) The final round combines the partial results of the term phase and reframes the output into the form produced by the non-parallel >>16095 version, for easy comparison: class Lagrounds: @ staticmethod def r2 (pieces, extra): c = fpoly.LagrangePieces.terms_combine (pieces, None) out = fpoly.Inter.withexpodown (c) return (out, None) These rounds will be consumed by the >>17263 rounds engine.

What on fuarrrking earth is going on in this thread. I don't seem to understand

>>17447 Since this is also the challenge thread, a level 1 challenge for novice adventurers >>15979 was posted before the interpolation interlude.

Comparing a single-threaded run and a parallelized run of Lagrange interpolation: def work_testlag_compare (mononame, parafmt, paracount): data = fpoly.Data.fizzbuzz () t1 = taim.point () co1 = getattr (fpoly.Inter, mononame) (data, tofractions = True) t2 = taim.point () n = len (data) w = Workers.half () stst = Pieces.start_stop (0, n, w) t3 = taim.point () co2 = Engine.rounds ( Engine.roundfun_attr (Lagrounds, parafmt, paracount), {"n": n, "data": data, "chunks": stst}, None, w) t4 = taim.point () print ("match:", co1 == co2) print (" mono:", taim.diffstr (t1, t2)) print (" para:", taim.diffstr (t3, t4)) def work_testlag (): work_testlag_compare ("lagrange2", "r{}", 3) The fizzbuzz interpolation points are generated in >>16020. The quadratic lagrange is in >>16095. The parallel rounds engine as well as the round function construction based on naming pattern and container are in >>17263. Individual round functions >>17358 are based on split >>17187 lagrange phases. Results: $ python3 -m flake.futures testlag match: True mono: 00:00:00.248 para: 00:00:00.112 The match is just a sanity check, whereas the improvement factor will increase for longer runs, such as the slower, cubic lagrange.

>>17491 What should I read up on to even begin to understand this sorta stuff? My maths knowledge pretty much stops at high-school level

>>17494 High school math should be sufficient to be able to follow the method of computation. For this aspect see e.g. the wiki pages: https://en.wikipedia.org/wiki/Polynomial_interpolation https://en.wikipedia.org/wiki/Lagrange_polynomial https://en.wikipedia.org/wiki/Newton_polynomial For the theoretical basis that justifies the method of computation consult a text on numerical analysis. For example, libgen has this copy of Theoretical Numerical Analysis - A Functional Analysis Framework by Kendall Atkinson and Weimin Han, from which section 3.1 and chapter 5 may be of particular interest.

For the cubic >>16071 Lagrange, the common product phase is absent and the term phase looks like: class LagrangePieces: @ staticmethod def cubic_terms (startstop, nxcoef): addtim = LagrangeHelper.add_times mulxma = LagrangeHelper.multiply_xma start, stop = startstop n, x, c = nxcoef todo = stop - start out = [0] * n for k in range (todo): pos = start + k term = [0] * n term [-1] = 1 for other in range (pos): mulxma (term, x [other]) for other in range (pos + 1, n): mulxma (term, x [other]) addtim (term, c [pos], out) return out The second round function changes to: class Lagrounds: @ staticmethod def rc1 (pieces, extra): fplp = fpoly.LagrangePieces c = fplp.basis_coeff_combine (pieces, None) stst = extra ["chunks"] n = extra ["n" ] x = extra ["x" ] rep = itertools.repeat return ( [stst, rep ((n, x, c), len (stst))], fplp.cubic_terms ) rc0 = r0 rc2 = r2 The outer rounds are unchanged >>17358. The comparison is run by def work_testcubic (): work_testlag_compare ("lagrange", "rc{}", 3) which results in $ python3 -m flake.futures testcubic match: True mono: 00:00:03.685 para: 00:00:00.977 As we would expect from our familiarity with orders of growth, for even such a relatively small number of interpolation points the parallelized cubic cannot compete with the single-threaded >>16095 quadratic.

Since this is the challenge thread, Lains might enjoy trying their hand at generating the values of the Hofstadter R sequence: https://en.wikipedia.org/wiki/Hofstadter_sequence#Hofstadter_Figure-Figure_sequences Here is a moderately efficient version that uses stacked self-generation of the S sequence and second-order S-groups to achieve O(log(log(n))) space consumption and O(n^0.25) runtime. def seq2_next (levels, level): if level < len (levels): state = levels [level] start, count, used, S = state if used < count: length = start + (used - 1) state [2] = used + 1 state [3] = S + (length + 1) return (S, length) else: start, count = seq2_next (levels, level + 1) length = start - 1 state [0] = start state [1] = count state [2] = 1 state [3] = S + (length + 1) return (S, length) else: levels.append ([4, 3, 2, 13]) return (8, 4) def seq3_sum1sc (s, c): # start, count return (c - 1 + s + s) * c // 2 def seq3_covered (s, c): # start, count return (c - 3 + s + s) * c // 2 def seq3_inc (s, c, S): # start, count, S # incS = (c - 1 + s + s) * c // 2 # countS = (c - 3 + s + s) * c // 2 cm1 = c - 1 incS = (cm1 + s + s) * c // 2 countS = incS - c overR = (incS - 1) * incS // 2 # backR = (s - 1) * c + s * (c - 1) * c // 2 + (c - 1) * c * (c + 1) // 6 cm1c = cm1 * c backR = (s - 1) * c + s * cm1c // 2 + cm1c * (c + 1) // 6 incR = overR - backR + S * countS return (incR, incS) def work_seq3 (n): seq3covered = seq3_covered seq3inc = seq3_inc seq3next = seq2_next seq3sum1sc = seq3_sum1sc levels = [] R, S = 3, 4 start, count = 4, 3 advance = seq3covered (start, count) covered = 2 + advance while covered <= n: incR, incS = seq3inc (start, count, S) R += incR S += incS start, count = seq3next (levels, 0) advance = seq3covered (start, count) covered += advance print ("levels:", len (levels)) covered -= advance fakelevels = [[start, count, 0, S]] start, count = seq3next (fakelevels, 0) covered += count while covered < n: R += seq3sum1sc (start, count) start, count = seq3next (fakelevels, 0) covered += count R += seq3sum1sc (start, count - (covered - n)) return R The first integer k for which R(10^k) overflows an u64 is k=10. The value of R(10^10) is 50000940106333921296 and it is computed in about 850 microseconds in vanilla python: $ g 10 k 10 index 10000000000 levels: 4 10000000000 -> 50000940106333921296 0x2B5E7061C419DA010 66 $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq3 (10000000000)' 1000 loops, best of 3: 841 usec per loop The first integer k for which R(10^k) overflows an u128 is k=20. The value of R(10^20) is 5000000000942800098290022420982686040347 and it is computed in 300 milliseconds in vanilla python: $ g 20 k 20 index 100000000000000000000 levels: 5 100000000000000000000 -> 5000000000942800098290022420982686040347 0xEB194F8ED94AC565124861C3B7D0C411B 132 $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq3 (100000000000000000000)' 10 loops, best of 3: 300 msec per loop Here is a table showing the exponent k up to 30, the value of R(10^k), the bit count of the R value and the space consumption of the algorithm, to verify other implementations: $ cat levels.txt 1 69 7 0 2 5764 13 1 3 526334 20 2 4 50874761 26 3 5 5028514748 33 3 6 500918449417 39 3 7 50029364025646 46 3 8 5000934571821489 53 4 9 500029664631299282 59 4 10 50000940106333921296 66 4 11 5000029765607020319048 73 4 12 500000941936492050650505 79 4 13 50000029798618763894670256 86 4 14 5000000942529842698007077786 93 4 15 500000029809255627825531266625 99 5 16 50000000942720152238624116795401 106 5 17 5000000029812655507343465281696595 112 5 18 500000000942780823112495107784753816 119 5 19 50000000029813737262126730811322149673 126 5 20 5000000000942800098290022420982686040347 132 5 21 500000000029814080548392288266955229183571 139 5 22 50000000000942806209878293665994446398371544 146 5 23 5000000000029814189323670710814676032031444555 152 5 24 500000000000942808145472258657037814775197247031 159 5 25 50000000000029814223760934912839828327249688721190 166 5 26 5000000000000942808758091081952084125868709915711089 172 5 27 500000000000029814234658067055252766458087914346630413 179 5 28 50000000000000942808951913512756782469229223199209677641 186 5 29 5000000000000029814238105320163714566005337846445873165675 192 5 30 500000000000000942809013222649956573426148951590151576716407 199 6

The OEIS page for the Hofstadter R sequence https://oeis.org/A005228 has C and C++ programs for generating the sequence. Neither of them currently go above an u64, and they run in O(sqrt(n)) and O(n) respectively, so >>17859 is considerably faster than both.

The C++ version from OEIS >>17889 has the same O(log(log(n))) space consumption as >>17859. While the interplay of 'test' and 'generate' confused me for a moment, it becomes much clearer if 'generate' is renamed to 'advance' and 'test' to 'needs_jump'. It would be even better if 'test' did not have two responsibilities. The technique used is almost the same stack of generators as in >>17859, but with a variation. In >>17859 generator k only resorts to generator k+1 when this is needed. But in OEIS/C++ the call to 'test' is always made, and it both returns whether a jump is needed at level k, and advances level k+1 if a check determines this is needed. The generator at level k does not possess enough information to avoid the need for the 'while' check in 'test' on every single invocation, even though this is avoidable. The same extra information would also prevent the need for most 'test' calls. Currently, other than the code golfed increment inside the argument list, the overwhelming majority of 'test' calls are useless no-ops, hence the inline request. In >>17859 this extra information is available so the call pattern is much more disciplined. From a design perspective, the C++ version sends out notifications all the time and the recipient decides whether they were needed, but in my version notifications are only sent out when it is known that they are needed. The OEIS/C++ version for reference: /*From Catalin Francu C++ source code that generates A005228. Runs in (approx.) O( n + sqrt(n) + sqrt(sqrt(n)) + sqrt(sqrt(sqrt(n)))...) to generate the first n elements. Tested up to n = 1,000,000,000. Input is n, output is a(n). Needs quasi-constant memory. %o A005228*/ #include <stdio.h> typedef unsigned long long uint64; typedef struct { uint64 value; uint64 inc; } pair; pair series[100]; void init_all_series() { pair initial_pair = { 7, 5 }; // pregenerate the first three terms for (int i = 0; i < 99; series[i++] = initial_pair); } void generate(int k); inline bool test(int k, unsigned number) { switch (number) { case 1: case 3: case 7: return 1; case 2: case 4: case 5: case 6: return 0; default: while (series[k].value < number) generate(k); return (series[k].value == number); } } inline void generate(int k) { series[k].value += series[k].inc; if (test(k+1, ++series[k].inc)) series[k].inc++; } int main(void) { int x; init_all_series(); scanf("%d", &x); // Enter the index switch (x) { case 1: puts("1"); break; case 2: puts("3"); break; default: for (int k = 4; k <= x; k++) generate(0); printf("%llu", series[0].value); // print the x'th term } return 0; }

Here is the upgraded version of >>17859 with third-order S-groups and a runtime of O(n^(1/8)). def seq4_covered (s, c, S1): # start, count incR, incS = seq3_inc (s, c, S1 - 1) return incR def seq4_inc (s, c, S0, S1, covered): # start, count incS1 = seq3_sum1sc (s, c) countS1 = incS1 - c incS0 = countS1 * S1 + seq3_inc (s, c, 0) [0] countS0 = covered overR = (incS0 - 1) * incS0 // 2 backR = ( countS1 * (countS1 + 1) // 2 * S1 - countS1 + seq4_backR (s, c) ) incR = overR - backR + S0 * countS0 return (incR, incS0, incS1) def seq4_backR (s, c): # start, count return c*(((((5*c + (30*s - 29))*c + ((60*s - 130)*s + 45))*c + (((40*s - 140)*s + 90)*s + 85))*c + ((-60*s + 250)*s - 290))*c + (20*s - 160)*s + 184) // 240 def work_seq4 (n): seq3covered = seq3_covered seq3inc = seq3_inc seq4next = seq2_next seq3sum1sc = seq3_sum1sc seq4covered = seq4_covered seq4inc = seq4_inc levels = [] R, S0, S1 = 3, 4, 4 start, count = 4, 3 advance = seq4covered (start, count, S1) covered = 2 + advance while covered <= n: incR, incS0, incS1 = seq4inc (start, count, S0, S1, advance) R += incR S0 += incS0 S1 += incS1 start, count = seq4next (levels, 0) advance = seq4covered (start, count, S1) covered += advance print ("levels:", len (levels)) covered -= advance fakelevels = [[start, count, 0, S1]] start, count = seq4next (fakelevels, 0) advance = seq3covered (start, count) covered += advance while covered <= n: incR, incS0 = seq3inc (start, count, S0) R += incR S0 += incS0 start, count = seq4next (fakelevels, 0) advance = seq3covered (start, count) covered += advance covered -= advance fakelevels = [[start, count, 0, S0]] start, count = seq4next (fakelevels, 0) covered += count while covered < n: R += seq3sum1sc (start, count) start, count = seq4next (fakelevels, 0) covered += count R += seq3sum1sc (start, count - (covered - n)) return R The timing for R(10^20) is down to about 80 milliseconds in vanilla python: $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq4 (100000000000000000000)' 10 loops, best of 3: 81.2 msec per loop The timing for R(10^30) is about 14 seconds in vanilla python: $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq4 (1000000000000000000000000000000)' 10 loops, best of 3: 14.3 sec per loop The first integer k for which R(10^k) overflows an u256 is k=39. The value of R(10^39) is 500000000000000000029814239694953305682145634502486118475654808671411192582303 $ g4 39 k 39 index 1000000000000000000000000000000000000000 levels: 5 1000000000000000000000000000000000000000 -> 500000000000000000029814239694953305682145634502486118475654808671411192582303 0x4516DF8A16FE63D603196794A678455A328B511811E5818CC6F66DDCD067DD09F 259 Here is a table showing the exponent k up to 40, the value of R(10^k), the bit count of the R value and the space consumption of the algorithm, to verify other implementations: 1 69 7 0 2 5764 13 0 3 526334 20 1 4 50874761 26 2 5 5028514748 33 2 6 500918449417 39 2 7 50029364025646 46 2 8 5000934571821489 53 3 9 500029664631299282 59 3 10 50000940106333921296 66 3 11 5000029765607020319048 73 3 12 500000941936492050650505 79 3 13 50000029798618763894670256 86 3 14 5000000942529842698007077786 93 3 15 500000029809255627825531266625 99 4 16 50000000942720152238624116795401 106 4 17 5000000029812655507343465281696595 112 4 18 500000000942780823112495107784753816 119 4 19 50000000029813737262126730811322149673 126 4 20 5000000000942800098290022420982686040347 132 4 21 500000000029814080548392288266955229183571 139 4 22 50000000000942806209878293665994446398371544 146 4 23 5000000000029814189323670710814676032031444555 152 4 24 500000000000942808145472258657037814775197247031 159 4 25 50000000000029814223760934912839828327249688721190 166 4 26 5000000000000942808758091081952084125868709915711089 172 4 27 500000000000029814234658067055252766458087914346630413 179 4 28 50000000000000942808951913512756782469229223199209677641 186 4 29 5000000000000029814238105320163714566005337846445873165675 192 4 30 500000000000000942809013222649956573426148951590151576716407 199 5 31 50000000000000029814239195666197369782628797813806672178028303 205 5 32 5000000000000000942809032613363227344483520545074111347580106666 212 5 33 500000000000000029814239540504829346639521751526121067855028048469 219 5 34 50000000000000000942809038745792444395013325500575549841178605176222 225 5 35 5000000000000000029814239649559693010895101002840721426848127286038482 232 5 36 500000000000000000942809040685134541874198687080125967207639651015633313 239 5 37 50000000000000000029814239684047173947897406203055196803310606765583116781 245 5 38 5000000000000000000942809041298425788194699108125613517442497904284636150437 252 5 39 500000000000000000029814239694953305682145634502486118475654808671411192582303 259 5 40 50000000000000000000942809041492368615124254362375085779982992292343853681977675 265 5

>>17859 here is my attempt, in rust. Simple naive algorithm. I believe it to take O(log(n)) space and O(nlog(n)) runtime, but I am no good at asymptotic analysis. It is very slow on my laptop, takes about 30 seconds to find R(10^10). Can you explain how yours works? What is "stacked self-generation"?

>>18216 use std::collections::VecDeque; use std::str::FromStr; fn get_RS(iter_to: usize) -> (u128,u128) { if iter_to < 3 { (0,0) } else { //state: //aside for r_prev and s_prev //we only hold pairs of R[k_i], S[k_i] which satisfy the following properties: //k_0 = 3 (i.e. R[k_0] = 7, S[k_0] = 5) //k_i is the smallest number such that R[k_i] > S[k_(i+1)] //k_N (the last element of the state) is the smallest number such that R[k_N] > s_prev //thus the state stays at O(logn) let mut state: VecDeque<(u128,u128)> = VecDeque::with_capacity(iter_to.next_power_of_two().trailing_zeros() as usize + 100); let mut r_prev = 7; let mut s_prev = 5; state.push_back((7,5)); for _ in 3..iter_to { //if it happens that s_prev + 1 == R[k_N], then we have to update R[k_N] if state.back().unwrap().0 == s_prev + 1 { r_prev = r_prev + s_prev; s_prev = s_prev + 2; //N == 0 => R[k_N] = (7,5), and this would cause a segfault. So we use a special //update rule. if state.len() > 1 { let mut got_to = 0; //we proceed by going backwards through the state for j in 0..state.len() - 1 { let this_index = state.len() - (j + 1); let (r_this,s_this) = state[this_index]; let r_this_prev = state[this_index - 1].0; //if it happens that R[k_(i-1)] == S[k_i] + 1, then we must update //R[k_(i-1)] if s_this + 1 == r_this_prev { state[this_index] = (r_this + s_this, s_this + 2); //otherwise, there is no reason to continue updating, and we halt. } else { state[this_index] = (r_this + s_this, s_this + 1); got_to = j; break; } } //if we got all the way to k_0, we increment it and add a (7,5) pair at the //back of the state. if got_to == state.len() - 2 { state[0] = (12,6); state.push_front((7,5)); } } else { state[0] = (12,6); state.push_front((7,5)); } } else { r_prev = r_prev + s_prev; s_prev = s_prev + 1; } } (r_prev,s_prev) } } fn main() { let a = std::env::args().last().unwrap(); let iter_to:usize = usize::from_str(&a).unwrap(); let (R,S) = get_RS(iter_to); println!("{}, {}",R,S); }

>>18217 ah -- I think I understand how this works now! Improved version does 10^10 in a few milliseconds. And it would be simple to go optimize this much more, with a little bit of fancy arithmetic. I might work on it tomorrow. This is a fun challenge! use std::collections::VecDeque; use std::str::FromStr; fn get_RS(iter_to: usize) -> (u128,u128) { if iter_to < 3 { (0,0) } else { //state: //aside for r_prev and s_prev //we only hold pairs of R[k_i], S[k_i] which satisfy the following properties: //k_0 = 3 (i.e. R[k_0] = 7, S[k_0] = 5) //k_i is the smallest number such that R[k_i] > S[k_(i+1)] //k_N (the last element of the state) is the smallest number such that R[k_N] > s_prev //thus the state stays at O(logn) let mut state: VecDeque<(u128,u128)> = VecDeque::with_capacity(iter_to.next_power_of_two().trailing_zeros() as usize + 100); let mut r_prev = 7; let mut s_prev = 5; let mut num_iters:usize = 3; state.push_back((7,5)); while num_iters < iter_to { //get the distance to the next update let next_update = state.back().unwrap().0; let dist_to_update = next_update - s_prev; let dist_to_finish = (iter_to - num_iters) as u128; if dist_to_finish < dist_to_update { r_prev += s_prev*dist_to_finish + ((dist_to_finish - 1)*dist_to_finish)/2; s_prev += dist_to_finish; break; } else { //add to r_prev the result of adding s_prev, s_prev + 1, ... , s_prev + //dist_to_updated - 1 r_prev += s_prev*dist_to_update + ((dist_to_update - 1)*dist_to_update)/2; //add to r_prev dist_update + 1, to increment it one past the update boundary s_prev += dist_to_update + 1; num_iters += dist_to_update as usize; } //now we update the state //N == 0 => R[k_N] = (7,5), and this would cause a segfault. So we use a special //update rule. if state.len() > 1 { let mut got_to = 0; //we proceed by going backwards through the state for j in 0..state.len() - 1 { let this_index = state.len() - (j + 1); let (r_this,s_this) = state[this_index]; let r_this_prev = state[this_index - 1].0; //if it happens that R[k_(i-1)] == S[k_i] + 1, then we must update //R[k_(i-1)] if s_this + 1 == r_this_prev { state[this_index] = (r_this + s_this, s_this + 2); //otherwise, there is no reason to continue updating, and we halt. } else { state[this_index] = (r_this + s_this, s_this + 1); got_to = j; break; } } //if we got all the way to k_0, we increment it and add a (7,5) pair at the //back of the state. if got_to == state.len() - 2 { state[0] = (12,6); state.push_front((7,5)); } } else { state[0] = (12,6); state.push_front((7,5)); } } (r_prev,s_prev) } } fn main() { let a = std::env::args().last().unwrap(); let iter_to:usize = usize::from_str(&a).unwrap(); let (R,S) = get_RS(iter_to); println!("{}, {}",R,S); }

>>18216 >>18217 >>18218 >I believe it to take O(log(n)) space >What is "stacked self-generation"? It appears you might already be using a variation of this technique, although I'm not familiar with rust so I cannot be sure yet. This can easily be resolved if you post a list of final stack lengths for computing R(10^k) for k from 1 to 10. In that case the space consumption would be O(log(log(n))) because the stack length would increase by one roughly when the exponent in the index of R doubles. It would also make more sense to reverse your VecDeque so that new entries come in at the top instead of shifting everything up. The name was chosen based on the development history on SchemeBBS where the space consumption was brought down first to O(n^(1/2)) then to O(n^(1/4)) and finally to O(log(log(n))), by placing a series of generators on top of each other where each generator provides a contiguous run for the previous generator, similar to your reversed state. Self-generation is justified by the observation that R entries are not strictly required for producing the S sequence, similar information can be kept using only the S sequence itself. To see why this is, since you already have (first-order) S-groups in >>18218, write out the S sequence and in parallel the sequence of dist_to_update values, and you will easily see the pattern. >and O(nlog(n)) runtime >And it would be simple to go optimize this much more, with a little bit of fancy arithmetic The stack length itself does not appear in the runtime because of the access pattern of the state array. In >>18217 one check is performed for n items, a second check is performed roughly for n^(1/2) items, a third one roughly for n^(1/4) and so on. The sum can be seen to be in O(n). By a similar argument, the runtime of >>18218 would be O(n^(1/2)), but it's not immediately clear to me that this implementation is correct. You seem to be properly increasing s_prev, but the top of the stack, which you need for next_update, only advances one entry, which would be wrong and result in the wrong final R. Also, if you try your hand at second-order and third-order S-groups and still find the formulas "simple", you must be better at basic algebra than I am. To verify that our indexing of the R sequence lines up, please post the value of R(123456789). Its md5 hash as a decimal string is: 1b9d1bf230f7fa6b4917f61d3c5a2c04

>>18223 (R,S) = 7622071922811773, 123472393 the 'state' vector doesn't keep the top (R,S) value, that is stored in r_prev and s_prev. It just keeps the smallest (R,S) such that R > s_prev

>>18228 >7622071922811773 Thanks, our R indexing lines up. >the 'state' vector doesn't keep the top (R,S) value I realize that. When I referred to top it meant top of the stack, not the running R value that will eventually be returned. Upon further reflection the increments of both s_prev and the top of the stack in >>18218 are correct, and what happens instead is that the previously present check for "state.back().unwrap().0 == s_prev + 1" would now be true on every pass except the last partial one where you break, which is why it was removed. The issue I brought up in >>18223, that of s_prev going too fast for the top of the stack, will manifest for second-order S-groups. You did not post a list of final stack lengths for computing R(10^k) for k from 1 to 10.

>>18229 3, 3, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7 up to 10^15 seems to be the same as yours just shifted by 3 (probably because mine holds some extra information)

(should be four 5's, not five, I mistyped)

>>18230 >>18231 Thanks, your space consumption is O(log(log(n))) as well. The shift comes from the higher S-group order in my version. This means that >>18218 is better than both the C and C++ versions on OEIS >>17889, and equal to their best combination.

Sorry, but could someone please explain what a fuarrrking "s-group" is? I can't find anything online, and either it wasn't explained or I'm too smoothed brained to understand the explanation.

>>18242 it just means skipping over parts of calculating the S sequence

>>18242 S is the established notation for that particular sequence. https://en.wikipedia.org/wiki/Hofstadter_sequence#Hofstadter_Figure-Figure_sequences If you look at the elements of the S sequence https://oeis.org/A030124 you will find that they form contiguous blocks separated by jumps. These blocks are the first-order S-groups, or simply S-groups. The first-order S-groups can themselves be written as a sequence that again forms blocks separated by jumps. These blocks are the second-order S-groups, and their elements are first-order S-groups. Repeating the process yields third-order S-groups, whose elements are second-order S-groups.

[1/2] The stated runtime of >>17859 is correct, but the stated runtime of >>18215 was based on a mistaken assumption about the third loop and is wrong. Here is the corrected version that uses partial S-groups at the order boundaries to eliminate the lower loops and achieve the O(n^(1/8)) runtime. def work_seq3_cheat (start, count, nleft, coverfun): wantlow = coverfun (start, 1) if wantlow > nleft: return (0, 0) low, high = 1, count diff = high - low while diff > 1: now = (low + high) // 2 want = coverfun (start, now) if want <= nleft: low = now wantlow = want else: high = now diff = high - low return (low, wantlow) def work_seq4 (n): seq3covered = seq3_covered seq3inc = seq3_inc seq4next = seq2_next seq3sum1sc = seq3_sum1sc seq4covered = seq4_covered seq4inc = seq4_inc levels = [] R, S0, S1 = 3, 4, 4 start, count = 4, 3 advance = seq4covered (start, count, S1) covered = 2 + advance while covered <= n: incR, incS0, incS1 = seq4inc (start, count, S0, S1, advance) R += incR S0 += incS0 S1 += incS1 start, count = seq4next (levels, 0) advance = seq4covered (start, count, S1) covered += advance print ("levels:", len (levels)) covered -= advance used, advance = work_seq3_cheat ( start, count, n - covered, lambda s, c: seq4covered (s, c, S1)) if used > 0: incR, incS0, incS1 = seq4inc (start, used, S0, S1, advance) R += incR S0 += incS0 S1 += incS1 covered += advance fakelevels = [[start, count, used, S1]] start, count = seq4next (fakelevels, 0) used, advance = work_seq3_cheat ( start, count, n - covered, seq3covered) if used > 0: incR, incS0 = seq3inc (start, used, S0) R += incR S0 += incS0 covered += advance fakelevels = [[start, count, used, S0]] start, count = seq4next (fakelevels, 0) nleft = n - covered if nleft > 0: R += seq3sum1sc (start, nleft) return R Here are the updated timings, all in vanilla python. - The u64 overflow R(10^10) is computed in under 100 microseconds. - The u128 overflow R(10^20) is computed in under 2 milliseconds. - R(10^30) is computed in under 40 milliseconds. - The u256 overflow R(10^39) is computed in just over 500 milliseconds. $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq4 (10000000000)' 10000 loops, best of 3: 96.6 usec per loop $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq4 (100000000000000000000)' 1000 loops, best of 3: 1.81 msec per loop $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq4 (1000000000000000000000000000000)' 10 loops, best of 3: 35.9 msec per loop $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq4 (1000000000000000000000000000000000000000)' 10 loops, best of 3: 527 msec per loop

[2/2] Here is a table showing the exponent k up to 60, the value of R(10^k), the bit count of the R value and the space consumption of the algorithm, to verify other implementations: $ cat levels4.txt 1 69 7 0 2 5764 13 0 3 526334 20 1 4 50874761 26 2 5 5028514748 33 2 6 500918449417 39 2 7 50029364025646 46 2 8 5000934571821489 53 3 9 500029664631299282 59 3 10 50000940106333921296 66 3 11 5000029765607020319048 73 3 12 500000941936492050650505 79 3 13 50000029798618763894670256 86 3 14 5000000942529842698007077786 93 3 15 500000029809255627825531266625 99 4 16 50000000942720152238624116795401 106 4 17 5000000029812655507343465281696595 112 4 18 500000000942780823112495107784753816 119 4 19 50000000029813737262126730811322149673 126 4 20 5000000000942800098290022420982686040347 132 4 21 500000000029814080548392288266955229183571 139 4 22 50000000000942806209878293665994446398371544 146 4 23 5000000000029814189323670710814676032031444555 152 4 24 500000000000942808145472258657037814775197247031 159 4 25 50000000000029814223760934912839828327249688721190 166 4 26 5000000000000942808758091081952084125868709915711089 172 4 27 500000000000029814234658067055252766458087914346630413 179 4 28 50000000000000942808951913512756782469229223199209677641 186 4 29 5000000000000029814238105320163714566005337846445873165675 192 4 30 500000000000000942809013222649956573426148951590151576716407 199 5 31 50000000000000029814239195666197369782628797813806672178028303 205 5 32 5000000000000000942809032613363227344483520545074111347580106666 212 5 33 500000000000000029814239540504829346639521751526121067855028048469 219 5 34 50000000000000000942809038745792444395013325500575549841178605176222 225 5 35 5000000000000000029814239649559693010895101002840721426848127286038482 232 5 36 500000000000000000942809040685134541874198687080125967207639651015633313 239 5 37 50000000000000000029814239684047173947897406203055196803310606765583116781 245 5 38 5000000000000000000942809041298425788194699108125613517442497904284636150437 252 5 39 500000000000000000029814239694953305682145634502486118475654808671411192582303 259 5 40 50000000000000000000942809041492368615124254362375085779982992292343853681977675 265 5 41 5000000000000000000029814239698402168816728284052931912470389462588984951289354475 272 5 42 500000000000000000000942809041553699275352369083555067298879759994519671549809149522 279 5 43 50000000000000000000029814239699492802485126708959655333290571225887358803916579127627 285 5 44 5000000000000000000000942809041573093831553392834637054332803551711977694994007009790372 292 5 45 500000000000000000000029814239699837692448235838620313436065493069882669969298165276420932 298 5 46 50000000000000000000000942809041579226946272526629650233821592037083538450290978659280817924 305 5 47 5000000000000000000000029814239699946756464770327217558477196037097548855222595700238539054865 312 5 48 500000000000000000000000942809041581166410560146744996087637489747926009831507674531860092771246 318 5 49 50000000000000000000000029814239699981245576707388365411525341190531682837962715909950754430162916 325 5 50 5000000000000000000000000942809041581779723574457735429357187680244102657245323510219796435077511331 332 5 51 500000000000000000000000029814239699992151998934816806362938022380717694881944763055115621651771537182 338 5 52 50000000000000000000000000942809041581973670277645147248950274174680459501036705469450283128240202936345 345 5 53 5000000000000000000000000029814239699995600913788820676216348225539327181484799269054153129487836460510789 352 5 54 500000000000000000000000000942809041582035001627896095619939474929936571762744906288012004273501460096218956 358 5 55 50000000000000000000000000029814239699996691556662690868009166942952903517096871912188005052332391102154504541 365 5 56 5000000000000000000000000000942809041582054396306898966297540047317634200313888303569207665197822737010723880453 372 5 57 500000000000000000000000000029814239699997036448263909956827998702688566561312699712649393246034276487069718675623 378 5 58 50000000000000000000000000000942809041582060529443468690126506942010712319519060683290576216509510379646108458690288 385 6 59 5000000000000000000000000000029814239699997145512571896807332586675028913049270348666820028707699723929284676993408671 391 6 60 500000000000000000000000000000942809041582062468911643702033358845623678711109326514636451791089982268036336668374054905 398 6

[1/7] Here is the upgraded version of >>18317 using stacked self-generation of the S sequence and layered, partial, fourth-order S-groups to achieve O(log(log(n))) space consumption and O(n^(1/16)) runtime. class Layer: @ classmethod def count (cls, c, Slist, counts, incs): pass @ classmethod def increment (cls, c, Slist, counts, incs): pass class ChainedLayers: class Layer1 (Layer): @ classmethod def count (cls, c, Slist, counts, incs): counts [0] = c @ classmethod def increment (cls, c, Slist, counts, incs): s = Slist [0] incs [0] = (c - 1 + s + s) * c // 2 def _access (setasup): def fun (cls): cls.up = setasup.increment return cls return fun @ _access (Layer1) class Layer2 (Layer): @ classmethod def count (cls, c, Slist, counts, incs): s = Slist [0] # Slist [0] = s - 1 # cls.up (c, Slist, counts, incs) # Slist [0] = s # counts [1] = incs [0] counts [1] = (c - 3 + s + s) * c // 2 @ classmethod def increment (cls, c, Slist, counts, incs): countS = counts [1] incS = countS + c s = Slist [0] overR = (incS - 1) * incS // 2 cm1c = (c - 1) * c backR = (s - 1) * c + (c + 1 + 3 * s) * cm1c // 6 incR = overR - backR + Slist [1] * countS incs [0] = incS incs [1] = incR @ _access (Layer2) class Layer3 (Layer): @ classmethod def count (cls, c, Slist, counts, incs): S = Slist [1] Slist [1] = S - 1 cls.up (c, Slist, counts, incs) Slist [1] = S counts [2] = incs [1] @ classmethod def increment (cls, c, Slist, counts, incs): # names from seq4_inc countS1 = counts [1] countS0 = counts [2] incS0 = countS0 + countS1 overR = (incS0 - 1) * incS0 // 2 s = Slist [0] s20 = 20 * s s30 = 30 * s s60 = s30 + s30 backR = c*(((((5*c + (s30 - 29))*c + ((s60 - 130)*s + 45))*c + (((s20 + s20 - 140)*s + 90)*s + 85))*c + ((250 - s60)*s - 290))*c + (s20 - 160)*s + 184) // 240 backR = countS1 * (countS1 + 1) // 2 * Slist [1] - countS1 + backR incR = overR - backR + Slist [2] * countS0 incs [1] = incS0 incs [2] = incR @ _access (Layer3) class Layer4 (Layer): @ classmethod def count (cls, c, Slist, counts, incs): S = Slist [2] Slist [2] = S - 1 cls.up (c, Slist, counts, incs) Slist [2] = S counts [3] = incs [2]

[2/7] ... @ classmethod def increment (cls, c, Slist, counts, incs): # s A B C countB = counts [2] countC = counts [3] incC = countC + countB s = Slist [0] A = Slist [1] B = Slist [2] part1 = A*(A*(A*(A*c**3*(c*(c*(c*(3628800*c + 29030400*s - 43545600) + s*(87091200*s - 261273600) + 195955200) + s*(s*(116121600*s - 522547200) + 783820800) - 391910400) + s*(s*(s*(58060800*s - 348364800) + 783820800) - 783820800) + 293932800) + c**2*(c*(c*(c*(c*(c*(c*(3628800*c + 36288000*s - 44755200) + s*(145152000*s - 362880000) + 195955200) + s*(s*(290304000*s - 1103155200) + 1175731200) - 263692800) + s*(s*(s*(290304000*s - 1490227200) + 2264371200) - 667699200) - 572140800) + s*(s*(s*(s*(116121600*s - 754790400) + 1219276800) + 754790400) - 3476390400) + 2340576000) + s*(s*(s*(-348364800*s + 2070835200) - 5080320000) + 5965747200) - 2743372800) + s*(s*(-309657600*s + 1393459200) - 2090188800) + 1045094400)) + c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(1360800*c + 16329600*s - 17236800) + s*(81648000*s - 175996800) + 74239200) + s*(s*(217728000*s - 718502400) + 596332800) - 57093120) + s*(s*(s*(326592000*s - 1466035200) + 1735776000) - 52980480) - 518222880) + s*(s*(s*(s*(261273600*s - 1495065600) + 2036966400) + 1269112320) - 3868421760) + 1589353920) + s*(s*(s*(s*(s*(87091200*s - 609638400) + 462067200) + 3914265600) - 9503222400) + 6943708800) - 754034400) + s*(s*(s*(s*(-503193600*s + 3241728000) - 7576934400) + 6643123200) + 1337817600) - 3952126080) + s*(s*(s*(290304000*s - 2157926400) + 8087385600) - 13536875520) + 7980094080) + s*(s*(s*(-19353600*s + 1645056000) - 7063096320) + 10578677760) - 5430136320) + s*(232243200*s - 696729600) + 503193600) - 38707200*s + 58060800) + c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(226800*c + 3175200*s - 2948400) + s*(19051200*s - 36288000) + 12625200) + s*(s*(63504000*s - 185976000) + 127461600) - 3707760) + s*(s*(s*(127008000*s - 508032000) + 500774400) + 41845440) - 133524720) + s*(s*(s*(s*(152409600*s - 780192000) + 925344000) + 542263680) - 1263376800) + 366997680) + s*(s*(s*(s*(s*(101606400*s - 638668800) + 700358400) + 1967454720) - 4436147520) + 2285216640) + 69970320) + s*(s*(s*(s*(s*(s*(29030400*s - 217728000) - 50803200) + 2999404800) - 6769929600) + 4090443840) + 2021513760) - 1954478160) + s*(s*(s*(s*(s*(-241920000*s + 1678924800) - 3536870400) + 392716800) + 8533123200) - 10621235520) + 3069944640) + s*(s*(s*(s*(464486400*s - 3215923200) + 10204992000) - 15705043200) + 8335474560) + 1164072000) + s*(s*(s*(s*(-19353600*s + 899942400) - 3424619520) + 710277120) + 9379426560) - 8806936320) + s*(s*(s*(58060800*s - 2532096000) + 10668994560) - 16847631360) + 9472619520) + s*(s*(58060800*s - 1403136000) + 3741050880) - 2127928320) + s*(3870720*s + 1033482240) - 1552711680) part2 = B*(A*(A*(A*c**2*(c*(c*(29030400*c + 174182400*s - 261273600) + s*(348364800*s - 1045094400) + 783820800) + s*(s*(232243200*s - 1045094400) + 1567641600) - 783820800) + c*(c*(c*(c*(c*(c*(21772800*c + 174182400*s - 203212800) + s*(522547200*s - 1248307200) + 551577600) + s*(s*(696729600*s - 2554675200) + 2090188800) + 246758400) + s*(s*(s*(348364800*s - 1741824000) + 1393459200) + 2815948800) - 3418329600) + s*(s*(-1161216000*s + 5167411200) - 8360755200) + 4833561600) + s*(-696729600*s + 2090188800) - 1567641600)) + c*(c*(c*(c*(c*(c*(c*(c*(5443200*c + 54432000*s - 52617600) + s*(217728000*s - 435456000) + 131846400) + s*(s*(435456000*s - 1349913600) + 742694400) + 256677120) + s*(s*(s*(435456000*s - 1857945600) + 1117670400) + 2228567040) - 1665740160) + s*(s*(s*(s*(174182400*s - 958003200) - 241920000) + 6178636800) - 7575724800) + 1616630400) + s*(s*(s*(-1122508800*s + 5554483200) - 8157542400) + 1485388800) + 3752179200) + s*(s*(1161216000*s - 5786726400) + 12275020800) - 9252472320) + s*(s*(-38707200*s + 2922393600) - 8387850240) + 5904783360) - 464486400*s + 696729600) + B*(A*(A*c*(c*(58060800*c + 232243200*s - 348364800) + s*(232243200*s - 696729600) + 522547200) + c*(c*(c*(c*(29030400*c + 174182400*s - 183859200) + s*(348364800*s - 774144000) + 145152000) + s*(s*(232243200*s - 812851200) - 58060800) + 1112832000) + s*(-928972800*s + 2748211200) - 2032128000)) + c*(c*(c*(c*(c*(c*(3628800*c + 29030400*s - 24192000) + s*(87091200*s - 154828800) + 4032000) + s*(s*(116121600*s - 329011200) - 67737600) + 290304000) + s*(s*(s*(58060800*s - 232243200) - 377395200) + 1354752000) - 473760000) + s*(s*(-464486400*s + 1606348800) - 774144000) - 846720000) + s*(928972800*s - 2709504000) + 1975680000)) + c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(453600*c + 5443200*s - 4536000) + s*(27216000*s - 47174400) + 10735200) + s*(s*(72576000*s - 195955200) + 81043200) + 36281280) + s*(s*(s*(108864000*s - 406425600) + 197164800) + 371589120) - 204876000) + s*(s*(s*(s*(87091200*s - 420940800) + 91929600) + 1411845120) - 1371565440) + 133156800) + s*(s*(s*(s*(s*(29030400*s - 174182400) - 263692800) + 2364364800) - 2937916800) + 19595520) + 935373600) + s*(s*(s*(s*(-270950400*s + 1475712000) - 1703116800) - 2165184000) + 5123865600) - 2034002880) + s*(s*(s*(677376000*s - 3506227200) + 7190668800) - 5481745920) - 251395200) + s*(s*(s*(-19353600*s + 541900800) - 782853120) - 3935554560) + 5662406400) + s*(s*(77414400*s - 3286886400) + 8872980480) - 5715763200) + 928972800*s - 1354752000) part3 = c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(c*(14175*c + 226800*s - 189000) + s*(1587600*s - 2721600) + 812700) + s*(s*(6350400*s - 16783200) + 9903600) + 52080) + s*(s*(s*(15876000*s - 57456000) + 49215600) + 6347040) - 11050550) + s*(s*(s*(s*(25401600*s - 117936000) + 125193600) + 65677920) - 128034480) + 28664160) + s*(s*(s*(s*(s*(25401600*s - 145152000) + 162388800) + 281756160) - 588183120) + 233335200) + 23121812) + s*(s*(s*(s*(s*(s*(14515200*s - 99187200) + 76204800) + 609799680) - 1329558720) + 628068000) + 363318480) - 231117040) + s*(s*(s*(s*(s*(s*(s*(3628800*s - 29030400) - 37497600) + 662054400) - 1443304800) + 413965440) + 1801674000) - 1744031520) + 303711135) + s*(s*(s*(s*(s*(s*(-38707200*s + 287884800) - 517708800) - 719913600) + 3666835200) - 4386013920) + 1191840000) + 491453160) + s*(s*(s*(s*(s*(122572800*s - 891072000) + 2710512000) - 3678716160) + 179323200) + 3719604000) - 1855381864) + s*(s*(s*(s*(s*(-4838400*s + 75479040) + 11289600) - 2663646720) + 7910091840) - 7185104640) + 1227340800) + s*(s*(s*(s*(25804800*s - 972518400) + 4195215360) - 5737536000) + 766308480) + 2760323440) + s*(s*(s*(-6451200*s + 669312000) - 3513619200) + 7885409280) - 5741832320) + s*(s*(s*(1612800*s - 119347200) + 2311787520) - 5494809600) + 2494340352) + s*(s*(s*(645120*s - 5160960) + 48015360) - 2009825280) + 2832168960 extra = c * (part1 + part2 + part3) // 464486400 incR = extra + Slist [3] * countC incs [2] = incC incs [3] = incR

[3/7] class Group: def counts (self, c, Slist, counts, incs): pass def increments (self, c, Slist, counts, incs): pass class LayeredGroups: class Group1 (Group): def __init__ (self, layer1): self.layer1 = layer1 def counts (self, c, Slist, counts, incs): self.layer1.count (c, Slist, counts, incs) def increments (self, c, Slist, counts, incs): self.layer1.increment (c, Slist, counts, incs) class GroupN (Group): def __init__ (self, order, layers): self.order = order self.layers = layers self.cache_coufun = [l.count for l in layers [:order]] self.cache_incfun = layers [order - 1].increment def counts (self, c, Slist, counts, incs): for fun in self.cache_coufun: fun (c, Slist, counts, incs) def increments (self, c, Slist, counts, incs): self.cache_incfun (c, Slist, counts, incs) def work_layers_partial (c, Slist, counts, incs, nleft, countfun, order): order -= 1 countfun (1, Slist, counts, incs) want = counts [order] if want > nleft: return 0 low, high = 1, c diff = high - low savec = counts [:] savei = incs [:] while diff > 1: now = (low + high) // 2 countfun (now, Slist, counts, incs) want = counts [order] if want <= nleft: low = now savec [:] = counts savei [:] = incs else: high = now diff = high - low counts [:] = savec incs [:] = savei return low def work_layers_loop (n, groups): glen = len (groups) gm1 = glen - 1 R = 3 Slist = [ 4] * glen counts = [None] * glen incs = [None] * glen levels = [] s, c = 4, 3 done = 2 next = seq2_next g = groups [-1] gcounts = g.counts increments = g.increments slistrange = range (gm1) gcounts (c, Slist, counts, incs) advance = counts [gm1] done += advance while done <= n: increments (c, Slist, counts, incs) R += incs [gm1] for k in slistrange: Slist [k + 1] += incs [k] s, c = next (levels, 0) Slist [0] = s gcounts (c, Slist, counts, incs) advance = counts [gm1] done += advance done -= advance print ("levels:", len (levels)) for gorder in range (glen, 1, -1): g = groups [gorder - 1] used = work_layers_partial ( c, Slist, counts, incs, n - done, g.counts, gorder) if used > 0: advance = counts [gorder - 1] g.increments (used, Slist, counts, incs) R += incs [gorder - 1] for k in range (gorder - 1): Slist [k + 1] += incs [k] done += advance levels = [[s, c, used, Slist [1]]] s, c = next (levels, 0) Slist [0] = s del Slist [1] if done < n: g = groups [0] c = n - done g.counts (c, Slist, counts, incs) g.increments (c, Slist, counts, incs) R += incs [0] return R def work_layers (n): layers = [ ChainedLayers.Layer1, ChainedLayers.Layer2, ChainedLayers.Layer3, ChainedLayers.Layer4, ] groups = [ LayeredGroups.Group1 (layers [0]), LayeredGroups.GroupN (2, layers), LayeredGroups.GroupN (3, layers), LayeredGroups.GroupN (4, layers), ] return work_layers_loop (n, groups) Here are the updated timings, all in vanilla python. - The u64 overflow R(10^10) is too small for fourth-order S-groups to pay off and the third-order timing >>18317 of under 100 microseconds remains as reference. - The u128 overflow R(10^20) is computed in under one millisecond. - The u256 overflow R(10^39) is computed in 16 milliseconds. - The first integer k for which R(10^k) overflows an u512 is k=78. The value of R(10^78) is computed in under 5 seconds. $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_layers (100000000000000000000)' 1000 loops, best of 3: 945 usec per loop $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_layers (1000000000000000000000000000000000000000)' 100 loops, best of 3: 15.6 msec per loop $ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_layers (1000000000000000000000000000000000000000000000000000000000000000000000000000000)' 10 loops, best of 3: 4.79 sec per loop Here is a table showing the exponent k up to 120, the value of R(10^k), the bit count of the R value and the space consumption of the algorithm, to verify other implementations. This is most likely the last such table I'll post in full so please excuse the length. Consider it raw research data.

[4/7] $ cat table4.txt 1 69 7 0 2 5764 13 0 3 526334 20 0 4 50874761 26 1 5 5028514748 33 1 6 500918449417 39 1 7 50029364025646 46 1 8 5000934571821489 53 2 9 500029664631299282 59 2 10 50000940106333921296 66 2 11 5000029765607020319048 73 2 12 500000941936492050650505 79 2 13 50000029798618763894670256 86 2 14 5000000942529842698007077786 93 2 15 500000029809255627825531266625 99 3 16 50000000942720152238624116795401 106 3 17 5000000029812655507343465281696595 112 3 18 500000000942780823112495107784753816 119 3 19 50000000029813737262126730811322149673 126 3 20 5000000000942800098290022420982686040347 132 3 21 500000000029814080548392288266955229183571 139 3 22 50000000000942806209878293665994446398371544 146 3 23 5000000000029814189323670710814676032031444555 152 3 24 500000000000942808145472258657037814775197247031 159 3 25 50000000000029814223760934912839828327249688721190 166 3 26 5000000000000942808758091081952084125868709915711089 172 3 27 500000000000029814234658067055252766458087914346630413 179 3 28 50000000000000942808951913512756782469229223199209677641 186 3 29 5000000000000029814238105320163714566005337846445873165675 192 3 30 500000000000000942809013222649956573426148951590151576716407 199 4

[5/7] ... 31 50000000000000029814239195666197369782628797813806672178028303 205 4 32 5000000000000000942809032613363227344483520545074111347580106666 212 4 33 500000000000000029814239540504829346639521751526121067855028048469 219 4 34 50000000000000000942809038745792444395013325500575549841178605176222 225 4 35 5000000000000000029814239649559693010895101002840721426848127286038482 232 4 36 500000000000000000942809040685134541874198687080125967207639651015633313 239 4 37 50000000000000000029814239684047173947897406203055196803310606765583116781 245 4 38 5000000000000000000942809041298425788194699108125613517442497904284636150437 252 4 39 500000000000000000029814239694953305682145634502486118475654808671411192582303 259 4 40 50000000000000000000942809041492368615124254362375085779982992292343853681977675 265 4 41 5000000000000000000029814239698402168816728284052931912470389462588984951289354475 272 4 42 500000000000000000000942809041553699275352369083555067298879759994519671549809149522 279 4 43 50000000000000000000029814239699492802485126708959655333290571225887358803916579127627 285 4 44 5000000000000000000000942809041573093831553392834637054332803551711977694994007009790372 292 4 45 500000000000000000000029814239699837692448235838620313436065493069882669969298165276420932 298 4 46 50000000000000000000000942809041579226946272526629650233821592037083538450290978659280817924 305 4 47 5000000000000000000000029814239699946756464770327217558477196037097548855222595700238539054865 312 4 48 500000000000000000000000942809041581166410560146744996087637489747926009831507674531860092771246 318 4 49 50000000000000000000000029814239699981245576707388365411525341190531682837962715909950754430162916 325 4 50 5000000000000000000000000942809041581779723574457735429357187680244102657245323510219796435077511331 332 4 51 500000000000000000000000029814239699992151998934816806362938022380717694881944763055115621651771537182 338 4 52 50000000000000000000000000942809041581973670277645147248950274174680459501036705469450283128240202936345 345 4 53 5000000000000000000000000029814239699995600913788820676216348225539327181484799269054153129487836460510789 352 4 54 500000000000000000000000000942809041582035001627896095619939474929936571762744906288012004273501460096218956 358 4 55 50000000000000000000000000029814239699996691556662690868009166942952903517096871912188005052332391102154504541 365 4 56 5000000000000000000000000000942809041582054396306898966297540047317634200313888303569207665197822737010723880453 372 4 57 500000000000000000000000000029814239699997036448263909956827998702688566561312699712649393246034276487069718675623 378 4 58 50000000000000000000000000000942809041582060529443468690126506942010712319519060683290576216509510379646108458690288 385 5 59 5000000000000000000000000000029814239699997145512571896807332586675028913049270348666820028707699723929284676993408671 391 5 60 500000000000000000000000000000942809041582062468911643702033358845623678711109326514636451791089982268036336668374054905 398 5

[6/7] ... 61 50000000000000000000000000000029814239699997180001735682481622436253734484378611829704069855253539840822385586766235223775 405 5 62 5000000000000000000000000000000942809041582063082225349534414976370874538365419773853994503610373186092493285764626069077930 411 5 63 500000000000000000000000000000029814239699997190908167132758797213601135850071440613709845156968123373754366585399707913755854 418 5 64 50000000000000000000000000000000942809041582063276172175725114596742736120123079153718814281207807930795334286345090895730008850 425 5 65 5000000000000000000000000000000029814239699997194357083627207181322355525685599913243509599846712885018947248295200321875155512194 431 5 66 500000000000000000000000000000000942809041582063337503547853702069293508176406619791980181270066919789232260704875309077658871313179 438 5 67 50000000000000000000000000000000029814239699997195447726792843063653840657372075719381685614976391722679317994047727999252369922856111 445 5 68 5000000000000000000000000000000000942809041582063356898230747593270774582830958194072544271344029190586861314812326101187327152707466137 451 5 69 500000000000000000000000000000000029814239699997195792618445952782479879199592747532578345617801057490177491249280214020810071936484329248 458 5 70 50000000000000000000000000000000000942809041582063363031368009324911993981653436695081910395214825714129009666207952735693967804641114586360 465 5 71 5000000000000000000000000000000000029814239699997195901682763168111957892127609767027356231264653317643923496918836459245769003479450192260756 471 5 72 500000000000000000000000000000000000942809041582063364970836307405265344147379095849591258863086231065338167605132633027300820966076112761019590 478 5 73 50000000000000000000000000000000000029814239699997195936171928594984731404832410431326232603415497802278529513210446399336722685688471505037581717 485 5 74 5000000000000000000000000000000000000942809041582063365584150035124009253385106350912037609197638702902567591834232001254415139381286688289965434265 491 5 75 500000000000000000000000000000000000029814239699997195947078360337128508382533807826849705743505607000402961948719068672551336010824498525861155901131 498 5 76 50000000000000000000000000000000000000942809041582063365778096865206896625389696631480880925445600815945621310718924565484266894856253334433920035745613 504 5 77 5000000000000000000000000000000000000029814239699997195950527276883481022121537204551816939828037724959179119720005114619789384557978355406667770625673021 511 5 78 500000000000000000000000000000000000000942809041582063365839428238027648024522282176747574937498915190189902751724515553262387695706471431607983773634306949 518 5 79 50000000000000000000000000000000000000029814239699997195951617920058347188612863985188088019716697687439529171713502632907206540254112391538006189341119265333 524 5 80 5000000000000000000000000000000000000000942809041582063365858822921044628544419686861515336911850770576983501847405744828313705187013560455193479696110468330206 531 5 81 500000000000000000000000000000000000000029814239699997195951962811713098347227573413736252403144303445747301758908535925291294897583816428354549868387341178588346 538 5 82 50000000000000000000000000000000000000000942809041582063365864956058328249337495524348465181244040575792683558297760599704116002857011675528582962015410215017827124 544 5 83 5000000000000000000000000000000000000000029814239699997195952071876030605575553722003515432640851210314605052934617877898005396030823635769597183022713568800739692533 551 5 84 500000000000000000000000000000000000000000942809041582063365866895526630222251318501129855514085496196105626102656203814913215260423427828184462929343617690500791646666 558 5 85 50000000000000000000000000000000000000000029814239699997195952106365196084356764765809535763331820943974807290580225482459630513467554281256221042352711338115017463240541 564 5 86 5000000000000000000000000000000000000000000942809041582063365867508840358633208941979850651071615586015009119802497660371744120427648094844027760581745342810140498965106367 571 5 87 500000000000000000000000000000000000000000029814239699997195952117271627835731401282096332807194179800666433767421322710901699216993989552590481278739190818085942770801469939 578 5 88 50000000000000000000000000000000000000000000942809041582063365867702787188839192281512305720005862101628358486348391543656444785856690079066484000745660388932846521717583436672 584 5 89 5000000000000000000000000000000000000000000029814239699997195952120720544383725431643624010804777145056002961829350947459896502695796003933627466058904178569549311922074944075782 591 5 90 500000000000000000000000000000000000000000000942809041582063365867764118561681833720315754323670571454304241020343588011942036131390212670993069425493690411199667029506064382198016 597 5

[7/7] ... 91 50000000000000000000000000000000000000000000029814239699997195952121811187558883507241345558737726545999412237587447898252547907582547510519089410040047788072242400367711147080157727 604 5 92 5000000000000000000000000000000000000000000000942809041582063365867783513244702706919210496847388097907490730561168684789641855837414230081302939593583530252590810192511687081389584649 611 5 93 500000000000000000000000000000000000000000000029814239699997195952122156079213686575662975687566014375880870268129216859868227878333405950277188139310550757934138235810778808652737796639 617 5 94 50000000000000000000000000000000000000000000000942809041582063365867789646381987019941821472362222203714728005221324275936722498476128726389007200479444393049485586814721792217509104079110 624 5 95 5000000000000000000000000000000000000000000000029814239699997195952122265143531203035031439823370678039196191802270607001412643779490569427497770692978007295457906377399559843297890568189866 631 5 96 500000000000000000000000000000000000000000000000942809041582063365867791585850289115953723515409412530453716179723949902714208291933178644568856485101792218509649629380619056180826695356377537 637 5 97 50000000000000000000000000000000000000000000000029814239699997195952122299632696683457783499116765673688757669523265396670671188806266268720804525668981781967144106436746280657940665788439590390 644 5 98 5000000000000000000000000000000000000000000000000942809041582063365867792199164017548801668450709192195528441356626675780278327749842187890684434038130314296479515824504606690865756223764948222340 651 5 99 500000000000000000000000000000000000000000000000029814239699997195952122310539128435124332376888152625270659458987595531648669494610577903474208680016494042339796335956144164046454052238695640354016 657 5 100 50000000000000000000000000000000000000000000000000942809041582063365867792393110847758677724581317563632235727078699914934976135220713319517131197635890992941277256245865422676335917679671598995284082 664 5 101 5000000000000000000000000000000000000000000000000029814239699997195952122313988044983170272979669187564090056204087862225062294241222827751828464949986340519571497298198002598648939074996373483914447397 671 5 102 500000000000000000000000000000000000000000000000000942809041582063365867792454442220602011397936848113874420393537095457965215957158912649545530443213392290485913884416759482555267287764026295136483005354 677 5 103 50000000000000000000000000000000000000000000000000029814239699997195952122315078688158337579677655571138789770963267360003715873732916510918408098298214394658187267587530489642973331025727113006001249379524 684 5 104 5000000000000000000000000000000000000000000000000000942809041582063365867792473836903623007695580001287076860764743889778515149564610140577395467174937323081260061662272697625921374924102431637742255986945708 690 5 105 500000000000000000000000000000000000000000000000000029814239699997195952122315423579813142289648032340230234978179722289107606270016375163200939959606513106143600536988097227759180425522664743968616636021441441 697 5 106 50000000000000000000000000000000000000000000000000000942809041582063365867792479970040907342608601729480939538163751892603941518401843547290072675831349350848171553750716945596883505269063591409273766136095805030 704 5 107 5000000000000000000000000000000000000000000000000000029814239699997195952122315532644130659040929889331216248472550850083269783136807766340829778013041746760196610695120850754652637718380322298674397791085034653961 710 5 108 500000000000000000000000000000000000000000000000000000942809041582063365867792481909509209442513232282866157701883060253714497239934368925370275988749321175503399479440603299642882935414549432956007859008157718549272 717 5 109 50000000000000000000000000000000000000000000000000000029814239699997195952122315567133296139515592327460372291709902029314919549982973759829994860762372140116849311835050087465182281795286741721854664899827352151912566 724 5 110 5000000000000000000000000000000000000000000000000000000942809041582063365867792482522822937876053413359268446786533037870819916002566574758232461266293922758424473117794836218022356203203733610226062472495023950284860065 730 5 111 500000000000000000000000000000000000000000000000000000029814239699997195952122315578039727891191372323848942045187196389956266215002510615796659322779348204022642496544547177964957102389077431946880846604910511822868188331 737 5 112 50000000000000000000000000000000000000000000000000000000942809041582063365867792482716769768086052568450245312452843417669093617146330284626106116978987640409292027276080645020792488674761112283184614545920243865425193780466 744 5 113 5000000000000000000000000000000000000000000000000000000029814239699997195952122315581488644439238954477824572294581975143652415875036756283026191242594450815983121437863375107877624835093736796129957926568426462450366758392478 750 5 114 500000000000000000000000000000000000000000000000000000000942809041582063365867792482778101140929408132244800358157207174012141617105525396954793271772415553258087888277333849857396284671289525172375713654636523221821335286376853 757 5 115 50000000000000000000000000000000000000000000000000000000029814239699997195952122315582579287614406553089530324116033795051581537701215496920388649969502926051838237707358517637569703550538467520187019116517002585257142738052253253 764 5 116 5000000000000000000000000000000000000000000000000000000000942809041582063365867792482797495823950408322620234762526209920619206513480935420807211061249643782036650334666400635855807415064897612541238001972361413653203208834167908773 770 6 117 500000000000000000000000000000000000000000000000000000000029814239699997195952122315582924179269211314970329478642833991142699924678528609685030399323964767533396715839851240530141402298048009514860572457480335347009482712253713858492 777 6 118 50000000000000000000000000000000000000000000000000000000000942809041582063365867792482803628961234743927878607006595156633637303703998869376378658739885940334308940718835321426318790230646440539042552038785423215306810427479795851308403 783 6 119 5000000000000000000000000000000000000000000000000000000000029814239699997195952122315583033243586728075483310895233493405113328031324626690805041965550114607412481285291759737215817417394308157788139557311388604511159666876568051518443324 790 6 120 500000000000000000000000000000000000000000000000000000000000942809041582063365867792482805568429536843955608187840000424522078995935973483241408413592544151956341508175361598028482106450141472450620675293371345306727392622125356063314915792 797 6 I think it might be time for a GMP port.

[1/12] Here is the GMP port of the fourth-order version from >>18647. Care is taken to avoid any allocations in the main loop to get a good constant factor, but the orders of growth are unchanged. #include <stdlib.h> #include <stdio.h> #include <gmp.h> #define ENOUGH_BITS 1000 #define ENOUGH_LEVELS 12 #define ENOUGH_TEMPS 34 struct level_output { mpz_t start; mpz_t count; }; struct level { mpz_t start; mpz_t count; mpz_t used; mpz_t S; }; struct levels { struct level *stack; size_t count; }; struct partial_storage { mpz_t now; mpz_t low; mpz_t high; mpz_t diff; mpz_t *counts; mpz_t *incs; }; struct context { size_t gorder; mpz_t *slist; mpz_t *counts; mpz_t *incs; mpz_t *temps; struct levels levels; struct partial_storage partial; }; typedef void (*layerfp) (mpz_t c, struct context *ctx); struct layer { layerfp count; layerfp increment; }; struct groups { size_t count; struct layer *layers; }; /* zconst */ mpz_t z4386013920; mpz_t z4436147520; mpz_t z4833561600; mpz_t z5080320000; mpz_t z5123865600; mpz_t z5167411200; mpz_t z5430136320; mpz_t z5481745920; mpz_t z5494809600; mpz_t z5554483200; mpz_t z5662406400; mpz_t z5715763200; mpz_t z5737536000; mpz_t z5741832320; mpz_t z5786726400; mpz_t z5904783360; mpz_t z5965747200; mpz_t z6178636800; mpz_t z6643123200; mpz_t z6769929600; mpz_t z6943708800; mpz_t z7063096320; mpz_t z7185104640; mpz_t z7190668800; mpz_t z7575724800; mpz_t z7576934400; mpz_t z7885409280; mpz_t z7910091840; mpz_t z7980094080; mpz_t z8087385600; mpz_t z8157542400; mpz_t z8335474560; mpz_t z8360755200; mpz_t z8387850240; mpz_t z8533123200; mpz_t z8806936320; mpz_t z8872980480; mpz_t z9252472320; mpz_t z9379426560; mpz_t z9472619520; mpz_t z9503222400; mpz_t z10204992000; mpz_t z10578677760; mpz_t z10621235520; mpz_t z10668994560; mpz_t z12275020800; mpz_t z13536875520; mpz_t z15705043200; mpz_t z16847631360; void init_zconst (void) { mpz_init_set_str (z4386013920, "1056d42e0", 16); mpz_init_set_str (z4436147520, "1086a3d40", 16); mpz_init_set_str (z4833561600, "1201a4c00", 16); mpz_init_set_str (z5080320000, "12ecf8800", 16); mpz_init_set_str (z5123865600, "13167fc00", 16); mpz_init_set_str (z5167411200, "134007000", 16); mpz_init_set_str (z5430136320, "143a94e00", 16); mpz_init_set_str (z5481745920, "146bcce00", 16); mpz_init_set_str (z5494809600, "147842400", 16); mpz_init_set_str (z5554483200, "14b12b000", 16); mpz_init_set_str (z5662406400, "151817700", 16); mpz_init_set_str (z5715763200, "154afa000", 16); mpz_init_set_str (z5737536000, "155fbda00", 16); mpz_init_set_str (z5741832320, "1563d6880", 16); mpz_init_set_str (z5786726400, "158ea7000", 16); mpz_init_set_str (z5904783360, "15ff3d800", 16); mpz_init_set_str (z5965747200, "163961400", 16); mpz_init_set_str (z6178636800, "170468400", 16); mpz_init_set_str (z6643123200, "18bf60400", 16); mpz_init_set_str (z6769929600, "19384ed80", 16); mpz_init_set_str (z6943708800, "19de09680", 16); mpz_init_set_str (z7063096320, "1a4fe4c00", 16); mpz_init_set_str (z7185104640, "1ac43ff00", 16); mpz_init_set_str (z7190668800, "1ac98e600", 16); mpz_init_set_str (z7575724800, "1c38c6300", 16); mpz_init_set_str (z7576934400, "1c39ed800", 16); mpz_init_set_str (z7885409280, "1d601cc00", 16); mpz_init_set_str (z7910091840, "1d77a6c40", 16); mpz_init_set_str (z7980094080, "1dba69280", 16); mpz_init_set_str (z8087385600, "1e20bb600", 16); mpz_init_set_str (z8157542400, "1e63a3800", 16); mpz_init_set_str (z8335474560, "1f0d53f80", 16); mpz_init_set_str (z8360755200, "1f2570000", 16); mpz_init_set_str (z8387850240, "1f3f47000", 16); mpz_init_set_str (z8533123200, "1fc9d2080", 16); mpz_init_set_str (z8806936320, "20cef2f00", 16); mpz_init_set_str (z8872980480, "210def000", 16); mpz_init_set_str (z9252472320, "2277d8600", 16); mpz_init_set_str (z9379426560, "22f0eb100", 16); mpz_init_set_str (z9472619520, "2349cb400", 16); mpz_init_set_str (z9503222400, "2366faa80", 16); mpz_init_set_str (z10204992000, "26043d200", 16); mpz_init_set_str (z10578677760, "27689d000", 16); mpz_init_set_str (z10621235520, "279133140", 16); mpz_init_set_str (z10668994560, "27bebf000", 16); mpz_init_set_str (z12275020800, "2dba5f400", 16); mpz_init_set_str (z13536875520, "326dc5400", 16); mpz_init_set_str (z15705043200, "3a817f100", 16); mpz_init_set_str (z16847631360, "3ec327400", 16); } void clear_zconst (void) { mpz_clear (z4386013920); mpz_clear (z4436147520); mpz_clear (z4833561600); mpz_clear (z5080320000); mpz_clear (z5123865600); mpz_clear (z5167411200); mpz_clear (z5430136320); mpz_clear (z5481745920); mpz_clear (z5494809600); mpz_clear (z5554483200); mpz_clear (z5662406400); mpz_clear (z5715763200); mpz_clear (z5737536000); mpz_clear (z5741832320); mpz_clear (z5786726400); mpz_clear (z5904783360); mpz_clear (z5965747200); mpz_clear (z6178636800); mpz_clear (z6643123200); mpz_clear (z6769929600); mpz_clear (z6943708800); mpz_clear (z7063096320); mpz_clear (z7185104640); mpz_clear (z7190668800); mpz_clear (z7575724800); mpz_clear (z7576934400); mpz_clear (z7885409280); mpz_clear (z7910091840); mpz_clear (z7980094080); mpz_clear (z8087385600); mpz_clear (z8157542400); mpz_clear (z8335474560); mpz_clear (z8360755200); mpz_clear (z8387850240); mpz_clear (z8533123200); mpz_clear (z8806936320); mpz_clear (z8872980480); mpz_clear (z9252472320); mpz_clear (z9379426560); mpz_clear (z9472619520); mpz_clear (z9503222400); mpz_clear (z10204992000); mpz_clear (z10578677760); mpz_clear (z10621235520); mpz_clear (z10668994560); mpz_clear (z12275020800); mpz_clear (z13536875520); mpz_clear (z15705043200); mpz_clear (z16847631360); }

[2/12] mpz_t* make_array (size_t count) { size_t k; mpz_t *ret = malloc (count * sizeof (mpz_t)); for (k = 0; k < count; k++) { mpz_init2 (ret [k], ENOUGH_BITS); } return ret; } void init_level (struct level *lev) { mpz_init2 (lev->start, ENOUGH_BITS); mpz_init2 (lev->count, ENOUGH_BITS); mpz_init2 (lev->used, ENOUGH_BITS); mpz_init2 (lev->S, ENOUGH_BITS); } struct level* make_level_stack (size_t count) { size_t k; struct level *ret = malloc (count * sizeof (struct level)); for (k = 0; k < count; k++) { init_level (&ret [k]); } return ret; } void init_partial (struct partial_storage *partial, size_t count) { mpz_init2 (partial->now, ENOUGH_BITS); mpz_init2 (partial->low, ENOUGH_BITS); mpz_init2 (partial->high, ENOUGH_BITS); mpz_init2 (partial->diff, ENOUGH_BITS); partial->counts = make_array (count); partial->incs = make_array (count); } void init_context (struct context *ctx) { size_t n = ctx->gorder; ctx->slist = make_array (n); ctx->counts = make_array (n); ctx->incs = make_array (n); ctx->temps = make_array (ENOUGH_TEMPS); ctx->levels.stack = make_level_stack (ENOUGH_LEVELS); ctx->levels.count = 0; init_partial (&ctx->partial, n); } void init_groups (struct groups *gs) { gs->layers = malloc (gs->count * sizeof (struct layer)); } void clear_array (mpz_t *arr, size_t count) { size_t k; for (k = 0; k < count; k++) { mpz_clear (arr [k]); } free (arr); } void clear_level (struct level *lev) { mpz_clear (lev->start); mpz_clear (lev->count); mpz_clear (lev->used); mpz_clear (lev->S); } void clear_level_stack (struct level *stack, size_t count) { size_t k; for (k = 0; k < count; k++) { clear_level (&stack [k]); } free (stack); } void clear_partial (struct partial_storage *partial, size_t count) { mpz_clear (partial->now); mpz_clear (partial->low); mpz_clear (partial->high); mpz_clear (partial->diff); clear_array (partial->counts, count); clear_array (partial->incs, count); } void clear_context (struct context *ctx) { size_t n = ctx->gorder; clear_array (ctx->slist, n); clear_array (ctx->counts, n); clear_array (ctx->incs, n); clear_array (ctx->temps, ENOUGH_TEMPS); clear_level_stack (ctx->levels.stack, ENOUGH_LEVELS); ctx->levels.count = 0; clear_partial (&ctx->partial, n); } void clear_groups (struct groups *gs) { free (gs->layers); } void next (struct levels *levels, size_t level, struct level_output *output) { struct level *stack = levels->stack; #define SL stack [level] if (level < levels->count) { if (mpz_cmp (SL.used, SL.count) < 0) { mpz_set (output->start, SL.S); mpz_add (output->count, SL.start, SL.used); mpz_add_ui (SL.used, SL.used, 1); mpz_add (SL.S, SL.S, output->count); mpz_sub_ui (output->count, output->count, 1); } else { next (levels, level + 1, output); mpz_set (SL.start, output->start); mpz_set (SL.count, output->count); mpz_set_ui (SL.used, 1); mpz_set (output->start, SL.S); mpz_sub_ui (output->count, SL.start, 1); mpz_add (SL.S, SL.S, SL.start); } } else { levels->count += 1; mpz_set_ui (SL.start, 4); mpz_set_ui (SL.count, 3); mpz_set_ui (SL.used, 2); mpz_set_ui (SL.S, 13); mpz_set_ui (output->start, 8); mpz_set_ui (output->count, 4); } #undef SL } void group_count (size_t order, struct layer *layers, mpz_t c, struct context *ctx) { size_t k; for (k = 0; k < order; k++) { (*layers [k].count) (c, ctx); } } void group_increment (size_t order, struct layer *layers, mpz_t c, struct context *ctx) { (*layers [order - 1].increment) (c, ctx); } void transfer_array (mpz_t *src, mpz_t *dst, size_t count) { size_t k; for (k = 0; k < count; k++) { mpz_set (dst [k], src [k]); } } void partial (mpz_t c, struct context *ctx, mpz_t nleft, struct layer *layers, size_t order, mpz_t output) { size_t ordm1 = order - 1; struct partial_storage *p = &ctx->partial; mpz_set_ui (p->now, 1); group_count (order, layers, p->now, ctx); if (mpz_cmp (ctx->counts [ordm1], nleft) > 0) { mpz_set_ui (output, 0); return; } mpz_set_ui (p->low, 1); mpz_set (p->high, c); mpz_sub (p->diff, p->high, p->low); transfer_array (ctx->counts, p->counts, order); transfer_array (ctx->incs, p->incs, order); while (mpz_cmp_ui (p->diff, 1) > 0) { mpz_add (p->now, p->low, p->high); mpz_tdiv_q_2exp (p->now, p->now, 1); group_count (order, layers, p->now, ctx); if (mpz_cmp (ctx->counts [ordm1], nleft) <= 0) { mpz_set (p->low, p->now); transfer_array (ctx->counts, p->counts, order); transfer_array (ctx->incs, p->incs, order); } else { mpz_set (p->high, p->now); } mpz_sub (p->diff, p->high, p->low); } transfer_array (p->counts, ctx->counts, order); transfer_array (p->incs, ctx->incs, order); mpz_set (output, p->low); }

[3/12] void loop_zg (mpz_t n, struct groups *gs, mpz_t output) { struct context ctx; struct level_output levout; struct layer *layers = gs->layers; size_t k, gorder = gs->count; size_t g, ordm1 = gorder - 1; mpz_t done, nleft, partout; ctx.gorder = gorder; init_context (&ctx); mpz_init2 (levout.start, ENOUGH_BITS); mpz_init2 (levout.count, ENOUGH_BITS); mpz_init2 (done, ENOUGH_BITS); mpz_init2 (nleft, ENOUGH_BITS); mpz_init2 (partout, ENOUGH_BITS); mpz_set_ui (output, 3); mpz_set_ui (levout.count, 3); mpz_set_ui (done, 2); for (k = 0; k < gorder; k++) { mpz_set_ui (ctx.slist [k], 4); } group_count (gorder, layers, levout.count, &ctx); mpz_add (done, done, ctx.counts [ordm1]); while (mpz_cmp (done, n) <= 0) { group_increment (gorder, layers, levout.count, &ctx); mpz_add (output, output, ctx.incs [ordm1]); for (k = 0; k < ordm1; k++) { mpz_add (ctx.slist [k + 1], ctx.slist [k + 1], ctx.incs [k]); } next (&ctx.levels, 0, &levout); mpz_set (ctx.slist [0], levout.start); group_count (gorder, layers, levout.count, &ctx); mpz_add (done, done, ctx.counts [ordm1]); } mpz_sub (done, done, ctx.counts [ordm1]); printf ("levels: %zu\n", ctx.levels.count); for (g = gorder; g >= 2; g--) { mpz_sub (nleft, n, done); partial (levout.count, &ctx, nleft, layers, g, partout); if (mpz_cmp_ui (partout, 0) > 0) { mpz_add (done, done, ctx.counts [g - 1]); group_increment (g, layers, partout, &ctx); mpz_add (output, output, ctx.incs [g - 1]); for (k = 1; k < g; k++) { mpz_add (ctx.slist [k], ctx.slist [k], ctx.incs [k - 1]); } } ctx.levels.count = 1; mpz_set (ctx.levels.stack [0].start, levout.start); mpz_set (ctx.levels.stack [0].count, levout.count); mpz_set (ctx.levels.stack [0].used, partout); mpz_set (ctx.levels.stack [0].S, ctx.slist [1]); next (&ctx.levels, 0, &levout); mpz_set (ctx.slist [0], levout.start); for (k = 2; k < g; k++) { mpz_set (ctx.slist [k - 1], ctx.slist [k]); } } if (mpz_cmp (done, n) < 0) { mpz_sub (nleft, n, done); group_count (1, layers, nleft, &ctx); group_increment (1, layers, nleft, &ctx); mpz_add (output, output, ctx.incs [0]); } clear_context (&ctx); mpz_clear (levout.start); mpz_clear (levout.count); mpz_clear (done); mpz_clear (nleft); mpz_clear (partout); } void l1_count (mpz_t c, struct context *ctx) { mpz_set (ctx->counts [0], c); } void l1_increment (mpz_t c, struct context *ctx) { mpz_t *temps = ctx->temps; mpz_mul_2exp (temps [0], ctx->slist [0], 1); mpz_add (temps [0], temps [0], c); mpz_sub_ui (temps [0], temps [0], 1); mpz_mul (temps [1], temps [0], c); mpz_tdiv_q_2exp (ctx->incs [0], temps [1], 1); } void l2_count (mpz_t c, struct context *ctx) { mpz_t *temps = ctx->temps; mpz_mul_2exp (temps [0], ctx->slist [0], 1); mpz_add (temps [0], temps [0], c); mpz_sub_ui (temps [0], temps [0], 3); mpz_mul (temps [1], temps [0], c); mpz_tdiv_q_2exp (ctx->counts [1], temps [1], 1); } void l2_increment (mpz_t c, struct context *ctx) { mpz_t *temps = ctx->temps; mpz_add (temps [0], ctx->counts [1], c); mpz_sub_ui (temps [1], temps [0], 1); mpz_mul (temps [2], temps [1], temps [0]); mpz_tdiv_q_2exp (temps [1], temps [2], 1); mpz_sub_ui (temps [3], ctx->slist [0], 1); mpz_mul (temps [2], temps [3], c); mpz_sub_ui (temps [4], c, 1); mpz_mul (temps [3], temps [4], c); mpz_mul_ui (temps [5], ctx->slist [0], 3); mpz_add (temps [5], temps [5], c); mpz_add_ui (temps [5], temps [5], 1); mpz_mul (temps [4], temps [3], temps [5]); mpz_divexact_ui (temps [3], temps [4], 6); mpz_add (temps [2], temps [2], temps [3]); mpz_set (ctx->incs [0], temps [0]); mpz_mul (temps [3], ctx->counts [1], ctx->slist [1]); mpz_add (temps [3], temps [3], temps [1]); mpz_sub (ctx->incs [1], temps [3], temps [2]); } void l3_count (mpz_t c, struct context *ctx) { mpz_sub_ui (ctx->slist [1], ctx->slist [1], 1); l2_increment (c, ctx); mpz_add_ui (ctx->slist [1], ctx->slist [1], 1); mpz_set (ctx->counts [2], ctx->incs [1]); } void l3_increment (mpz_t c, struct context *ctx) { mpz_t *temps = ctx->temps; mpz_t *slist = ctx->slist; mpz_add (temps [0], ctx->counts [1], ctx->counts [2]); mpz_sub_ui (temps [1], temps [0], 1); mpz_mul (temps [2], temps [1], temps [0]); mpz_tdiv_q_2exp (temps [1], temps [2], 1); /* backR -> temps [2] */ mpz_mul_ui (temps [11], c, 5); mpz_sub_ui (temps [10], temps [11], 29); mpz_mul_ui (temps [11], slist [0], 30); mpz_add (temps [10], temps [10], temps [11]); mpz_mul (temps [9], c, temps [10]); mpz_add_ui (temps [8], temps [9], 45); mpz_mul_ui (temps [11], slist [0], 60); mpz_sub_ui (temps [10], temps [11], 130); mpz_mul (temps [9], slist [0], temps [10]); mpz_add (temps [8], temps [8], temps [9]); mpz_mul (temps [7], c, temps [8]); mpz_add_ui (temps [6], temps [7], 85); mpz_mul_ui (temps [11], slist [0], 40); mpz_sub_ui (temps [10], temps [11], 140); mpz_mul (temps [9], slist [0], temps [10]); mpz_add_ui (temps [8], temps [9], 90); mpz_mul (temps [7], slist [0], temps [8]); mpz_add (temps [6], temps [6], temps [7]); mpz_mul (temps [5], c, temps [6]); mpz_sub_ui (temps [4], temps [5], 290); mpz_mul_ui (temps [7], slist [0], 60); mpz_sub_ui (temps [6], temps [7], 250); mpz_neg (temps [6], temps [6]); mpz_mul (temps [5], slist [0], temps [6]); mpz_add (temps [4], temps [4], temps [5]); mpz_mul (temps [3], c, temps [4]); mpz_add_ui (temps [2], temps [3], 184); mpz_mul_ui (temps [5], slist [0], 20); mpz_sub_ui (temps [4], temps [5], 160); mpz_mul (temps [3], slist [0], temps [4]); mpz_add (temps [2], temps [2], temps [3]); /* high = 11 */ mpz_mul (temps [3], c, temps [2]); mpz_divexact_ui (temps [2], temps [3], 240); mpz_add_ui (temps [3], ctx->counts [1], 1); mpz_mul (temps [4], temps [3], ctx->counts [1]); mpz_tdiv_q_2exp (temps [3], temps [4], 1); mpz_mul (temps [4], temps [3], slist [1]); mpz_add (temps [2], temps [2], temps [4]); mpz_sub (temps [2], temps [2], ctx->counts [1]); mpz_set (ctx->incs [1], temps [0]); mpz_mul (temps [3], slist [2], ctx->counts [2]); mpz_add (temps [3], temps [3], temps [1]); mpz_sub (ctx->incs [2], temps [3], temps [2]); } void l4_count (mpz_t c, struct context *ctx) { mpz_sub_ui (ctx->slist [2], ctx->slist [2], 1); l3_increment (c, ctx); mpz_add_ui (ctx->slist [2], ctx->slist [2], 1); mpz_set (ctx->counts [3], ctx->incs [2]); }

[4/12] void l4_increment (mpz_t c, struct context *ctx) { mpz_t *temps = ctx->temps; mpz_t *slist = ctx->slist; /* part1 -> temps [2] */ mpz_mul_ui (temps [6], slist [0], 38707200); mpz_sub_ui (temps [5], temps [6], 58060800); mpz_neg (temps [5], temps [5]); mpz_mul (temps [9], c, c); mpz_mul_ui (temps [23], c, 3628800); mpz_sub_ui (temps [22], temps [23], 44755200); mpz_mul_ui (temps [23], slist [0], 36288000); mpz_add (temps [22], temps [22], temps [23]); mpz_mul (temps [21], c, temps [22]); mpz_add_ui (temps [20], temps [21], 195955200); mpz_mul_ui (temps [23], slist [0], 145152000); mpz_sub_ui (temps [22], temps [23], 362880000); mpz_mul (temps [21], slist [0], temps [22]); mpz_add (temps [20], temps [20], temps [21]); mpz_mul (temps [19], c, temps [20]); mpz_sub_ui (temps [18], temps [19], 263692800); mpz_mul_ui (temps [23], slist [0], 290304000); mpz_sub_ui (temps [22], temps [23], 1103155200); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 1175731200); mpz_mul (temps [19], slist [0], temps [20]); mpz_add (temps [18], temps [18], temps [19]); mpz_mul (temps [17], c, temps [18]); mpz_sub_ui (temps [16], temps [17], 572140800); mpz_mul_ui (temps [23], slist [0], 290304000); mpz_sub_ui (temps [22], temps [23], 1490227200); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 2264371200); mpz_mul (temps [19], slist [0], temps [20]); mpz_sub_ui (temps [18], temps [19], 667699200); mpz_mul (temps [17], slist [0], temps [18]); mpz_add (temps [16], temps [16], temps [17]); mpz_mul (temps [15], c, temps [16]); mpz_add_ui (temps [14], temps [15], 2340576000); mpz_mul_ui (temps [23], slist [0], 116121600); mpz_sub_ui (temps [22], temps [23], 754790400); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 1219276800); mpz_mul (temps [19], slist [0], temps [20]); mpz_add_ui (temps [18], temps [19], 754790400); mpz_mul (temps [17], slist [0], temps [18]); mpz_sub_ui (temps [16], temps [17], 3476390400); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], temps [14], temps [15]); mpz_mul (temps [13], c, temps [14]); mpz_sub_ui (temps [12], temps [13], 2743372800); /* zconst 5965747200 163961400 */ /* zconst 5080320000 12ecf8800 */ mpz_mul_ui (temps [19], slist [0], 348364800); mpz_sub_ui (temps [18], temps [19], 2070835200); mpz_neg (temps [18], temps [18]); mpz_mul (temps [17], slist [0], temps [18]); mpz_sub (temps [16], temps [17], z5080320000); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], z5965747200, temps [15]); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], temps [12], temps [13]); mpz_mul (temps [11], c, temps [12]); mpz_add_ui (temps [10], temps [11], 1045094400); mpz_mul_ui (temps [15], slist [0], 309657600); mpz_sub_ui (temps [14], temps [15], 1393459200); mpz_neg (temps [14], temps [14]); mpz_mul (temps [13], slist [0], temps [14]); mpz_sub_ui (temps [12], temps [13], 2090188800); mpz_mul (temps [11], slist [0], temps [12]); mpz_add (temps [10], temps [10], temps [11]); mpz_mul (temps [8], temps [9], temps [10]); mpz_mul (temps [12], c, c); mpz_mul (temps [11], temps [12], c); mpz_mul (temps [10], slist [1], temps [11]); mpz_mul_ui (temps [18], c, 3628800); mpz_sub_ui (temps [17], temps [18], 43545600); mpz_mul_ui (temps [18], slist [0], 29030400); mpz_add (temps [17], temps [17], temps [18]); mpz_mul (temps [16], c, temps [17]); mpz_add_ui (temps [15], temps [16], 195955200); mpz_mul_ui (temps [18], slist [0], 87091200); mpz_sub_ui (temps [17], temps [18], 261273600); mpz_mul (temps [16], slist [0], temps [17]); mpz_add (temps [15], temps [15], temps [16]); mpz_mul (temps [14], c, temps [15]); mpz_sub_ui (temps [13], temps [14], 391910400); mpz_mul_ui (temps [18], slist [0], 116121600); mpz_sub_ui (temps [17], temps [18], 522547200); mpz_mul (temps [16], slist [0], temps [17]); mpz_add_ui (temps [15], temps [16], 783820800); mpz_mul (temps [14], slist [0], temps [15]); mpz_add (temps [13], temps [13], temps [14]); mpz_mul (temps [12], c, temps [13]); mpz_add_ui (temps [11], temps [12], 293932800); mpz_mul_ui (temps [18], slist [0], 58060800); mpz_sub_ui (temps [17], temps [18], 348364800); mpz_mul (temps [16], slist [0], temps [17]); mpz_add_ui (temps [15], temps [16], 783820800); mpz_mul (temps [14], slist [0], temps [15]); mpz_sub_ui (temps [13], temps [14], 783820800); mpz_mul (temps [12], slist [0], temps [13]); mpz_add (temps [11], temps [11], temps [12]); mpz_mul (temps [9], temps [10], temps [11]); mpz_add (temps [7], temps [8], temps [9]); mpz_mul (temps [6], slist [1], temps [7]); mpz_add (temps [5], temps [5], temps [6]);

[5/12] /* zconst 5430136320 143a94e00 */ /* zconst 7980094080 1dba69280 */ mpz_mul_ui (temps [26], c, 1360800); mpz_sub_ui (temps [25], temps [26], 17236800); mpz_mul_ui (temps [26], slist [0], 16329600); mpz_add (temps [25], temps [25], temps [26]); mpz_mul (temps [24], c, temps [25]); mpz_add_ui (temps [23], temps [24], 74239200); mpz_mul_ui (temps [26], slist [0], 81648000); mpz_sub_ui (temps [25], temps [26], 175996800); mpz_mul (temps [24], slist [0], temps [25]); mpz_add (temps [23], temps [23], temps [24]); mpz_mul (temps [22], c, temps [23]); mpz_sub_ui (temps [21], temps [22], 57093120); mpz_mul_ui (temps [26], slist [0], 217728000); mpz_sub_ui (temps [25], temps [26], 718502400); mpz_mul (temps [24], slist [0], temps [25]); mpz_add_ui (temps [23], temps [24], 596332800); mpz_mul (temps [22], slist [0], temps [23]); mpz_add (temps [21], temps [21], temps [22]); mpz_mul (temps [20], c, temps [21]); mpz_sub_ui (temps [19], temps [20], 518222880); mpz_mul_ui (temps [26], slist [0], 326592000); mpz_sub_ui (temps [25], temps [26], 1466035200); mpz_mul (temps [24], slist [0], temps [25]); mpz_add_ui (temps [23], temps [24], 1735776000); mpz_mul (temps [22], slist [0], temps [23]); mpz_sub_ui (temps [21], temps [22], 52980480); mpz_mul (temps [20], slist [0], temps [21]); mpz_add (temps [19], temps [19], temps [20]); mpz_mul (temps [18], c, temps [19]); mpz_add_ui (temps [17], temps [18], 1589353920); mpz_mul_ui (temps [26], slist [0], 261273600); mpz_sub_ui (temps [25], temps [26], 1495065600); mpz_mul (temps [24], slist [0], temps [25]); mpz_add_ui (temps [23], temps [24], 2036966400); mpz_mul (temps [22], slist [0], temps [23]); mpz_add_ui (temps [21], temps [22], 1269112320); mpz_mul (temps [20], slist [0], temps [21]); mpz_sub_ui (temps [19], temps [20], 3868421760); mpz_mul (temps [18], slist [0], temps [19]); mpz_add (temps [17], temps [17], temps [18]); mpz_mul (temps [16], c, temps [17]); mpz_sub_ui (temps [15], temps [16], 754034400); /* zconst 6943708800 19de09680 */ /* zconst 9503222400 2366faa80 */ mpz_mul_ui (temps [26], slist [0], 87091200); mpz_sub_ui (temps [25], temps [26], 609638400); mpz_mul (temps [24], slist [0], temps [25]); mpz_add_ui (temps [23], temps [24], 462067200); mpz_mul (temps [22], slist [0], temps [23]); mpz_add_ui (temps [21], temps [22], 3914265600); mpz_mul (temps [20], slist [0], temps [21]); mpz_sub (temps [19], temps [20], z9503222400); mpz_mul (temps [18], slist [0], temps [19]); mpz_add (temps [17], z6943708800, temps [18]); mpz_mul (temps [16], slist [0], temps [17]); mpz_add (temps [15], temps [15], temps [16]); mpz_mul (temps [14], c, temps [15]); mpz_sub_ui (temps [13], temps [14], 3952126080); /* zconst 6643123200 18bf60400 */ /* zconst 7576934400 1c39ed800 */ mpz_mul_ui (temps [22], slist [0], 503193600); mpz_sub_ui (temps [21], temps [22], 3241728000); mpz_neg (temps [21], temps [21]); mpz_mul (temps [20], slist [0], temps [21]); mpz_sub (temps [19], temps [20], z7576934400); mpz_mul (temps [18], slist [0], temps [19]); mpz_add (temps [17], z6643123200, temps [18]); mpz_mul (temps [16], slist [0], temps [17]); mpz_add_ui (temps [15], temps [16], 1337817600); mpz_mul (temps [14], slist [0], temps [15]); mpz_add (temps [13], temps [13], temps [14]); mpz_mul (temps [12], c, temps [13]); mpz_add (temps [11], z7980094080, temps [12]); /* zconst 13536875520 326dc5400 */ /* zconst 8087385600 1e20bb600 */ mpz_mul_ui (temps [18], slist [0], 290304000); mpz_sub_ui (temps [17], temps [18], 2157926400); mpz_mul (temps [16], slist [0], temps [17]); mpz_add (temps [15], z8087385600, temps [16]); mpz_mul (temps [14], slist [0], temps [15]); mpz_sub (temps [13], temps [14], z13536875520); mpz_mul (temps [12], slist [0], temps [13]); mpz_add (temps [11], temps [11], temps [12]); mpz_mul (temps [10], c, temps [11]); mpz_sub (temps [9], temps [10], z5430136320); /* zconst 10578677760 27689d000 */ /* zconst 7063096320 1a4fe4c00 */ mpz_mul_ui (temps [16], slist [0], 19353600); mpz_sub_ui (temps [15], temps [16], 1645056000); mpz_neg (temps [15], temps [15]); mpz_mul (temps [14], slist [0], temps [15]); mpz_sub (temps [13], temps [14], z7063096320); mpz_mul (temps [12], slist [0], temps [13]); mpz_add (temps [11], z10578677760, temps [12]); mpz_mul (temps [10], slist [0], temps [11]); mpz_add (temps [9], temps [9], temps [10]); mpz_mul (temps [8], c, temps [9]); mpz_add_ui (temps [7], temps [8], 503193600); mpz_mul_ui (temps [10], slist [0], 232243200); mpz_sub_ui (temps [9], temps [10], 696729600); mpz_mul (temps [8], slist [0], temps [9]); mpz_add (temps [7], temps [7], temps [8]); mpz_mul (temps [6], c, temps [7]); mpz_add (temps [5], temps [5], temps [6]); mpz_mul (temps [4], slist [1], temps [5]); mpz_sub_ui (temps [3], temps [4], 1552711680); /* zconst 9472619520 2349cb400 */ /* zconst 8806936320 20cef2f00 */ mpz_mul_ui (temps [28], c, 226800); mpz_sub_ui (temps [27], temps [28], 2948400); mpz_mul_ui (temps [28], slist [0], 3175200); mpz_add (temps [27], temps [27], temps [28]); mpz_mul (temps [26], c, temps [27]); mpz_add_ui (temps [25], temps [26], 12625200); mpz_mul_ui (temps [28], slist [0], 19051200); mpz_sub_ui (temps [27], temps [28], 36288000); mpz_mul (temps [26], slist [0], temps [27]); mpz_add (temps [25], temps [25], temps [26]); mpz_mul (temps [24], c, temps [25]); mpz_sub_ui (temps [23], temps [24], 3707760); mpz_mul_ui (temps [28], slist [0], 63504000); mpz_sub_ui (temps [27], temps [28], 185976000); mpz_mul (temps [26], slist [0], temps [27]); mpz_add_ui (temps [25], temps [26], 127461600); mpz_mul (temps [24], slist [0], temps [25]); mpz_add (temps [23], temps [23], temps [24]); mpz_mul (temps [22], c, temps [23]); mpz_sub_ui (temps [21], temps [22], 133524720); mpz_mul_ui (temps [28], slist [0], 127008000); mpz_sub_ui (temps [27], temps [28], 508032000); mpz_mul (temps [26], slist [0], temps [27]); mpz_add_ui (temps [25], temps [26], 500774400); mpz_mul (temps [24], slist [0], temps [25]); mpz_add_ui (temps [23], temps [24], 41845440); mpz_mul (temps [22], slist [0], temps [23]); mpz_add (temps [21], temps [21], temps [22]); mpz_mul (temps [20], c, temps [21]); mpz_add_ui (temps [19], temps [20], 366997680); mpz_mul_ui (temps [28], slist [0], 152409600); mpz_sub_ui (temps [27], temps [28], 780192000); mpz_mul (temps [26], slist [0], temps [27]); mpz_add_ui (temps [25], temps [26], 925344000); mpz_mul (temps [24], slist [0], temps [25]); mpz_add_ui (temps [23], temps [24], 542263680); mpz_mul (temps [22], slist [0], temps [23]); mpz_sub_ui (temps [21], temps [22], 1263376800); mpz_mul (temps [20], slist [0], temps [21]); mpz_add (temps [19], temps [19], temps [20]); mpz_mul (temps [18], c, temps [19]); mpz_add_ui (temps [17], temps [18], 69970320);

[6/12] /* zconst 4436147520 1086a3d40 */ mpz_mul_ui (temps [28], slist [0], 101606400); mpz_sub_ui (temps [27], temps [28], 638668800); mpz_mul (temps [26], slist [0], temps [27]); mpz_add_ui (temps [25], temps [26], 700358400); mpz_mul (temps [24], slist [0], temps [25]); mpz_add_ui (temps [23], temps [24], 1967454720); mpz_mul (temps [22], slist [0], temps [23]); mpz_sub (temps [21], temps [22], z4436147520); mpz_mul (temps [20], slist [0], temps [21]); mpz_add_ui (temps [19], temps [20], 2285216640); mpz_mul (temps [18], slist [0], temps [19]); mpz_add (temps [17], temps [17], temps [18]); mpz_mul (temps [16], c, temps [17]); mpz_sub_ui (temps [15], temps [16], 1954478160); /* zconst 6769929600 19384ed80 */ mpz_mul_ui (temps [28], slist [0], 29030400); mpz_sub_ui (temps [27], temps [28], 217728000); mpz_mul (temps [26], slist [0], temps [27]); mpz_sub_ui (temps [25], temps [26], 50803200); mpz_mul (temps [24], slist [0], temps [25]); mpz_add_ui (temps [23], temps [24], 2999404800); mpz_mul (temps [22], slist [0], temps [23]); mpz_sub (temps [21], temps [22], z6769929600); mpz_mul (temps [20], slist [0], temps [21]); mpz_add_ui (temps [19], temps [20], 4090443840); mpz_mul (temps [18], slist [0], temps [19]); mpz_add_ui (temps [17], temps [18], 2021513760); mpz_mul (temps [16], slist [0], temps [17]); mpz_add (temps [15], temps [15], temps [16]); mpz_mul (temps [14], c, temps [15]); mpz_add_ui (temps [13], temps [14], 3069944640); /* zconst 10621235520 279133140 */ /* zconst 8533123200 1fc9d2080 */ mpz_mul_ui (temps [24], slist [0], 241920000); mpz_sub_ui (temps [23], temps [24], 1678924800); mpz_neg (temps [23], temps [23]); mpz_mul (temps [22], slist [0], temps [23]); mpz_sub_ui (temps [21], temps [22], 3536870400); mpz_mul (temps [20], slist [0], temps [21]); mpz_add_ui (temps [19], temps [20], 392716800); mpz_mul (temps [18], slist [0], temps [19]); mpz_add (temps [17], z8533123200, temps [18]); mpz_mul (temps [16], slist [0], temps [17]); mpz_sub (temps [15], temps [16], z10621235520); mpz_mul (temps [14], slist [0], temps [15]); mpz_add (temps [13], temps [13], temps [14]); mpz_mul (temps [12], c, temps [13]); mpz_add_ui (temps [11], temps [12], 1164072000); /* zconst 8335474560 1f0d53f80 */ /* zconst 15705043200 3a817f100 */ /* zconst 10204992000 26043d200 */ mpz_mul_ui (temps [20], slist [0], 464486400); mpz_sub_ui (temps [19], temps [20], 3215923200); mpz_mul (temps [18], slist [0], temps [19]); mpz_add (temps [17], z10204992000, temps [18]); mpz_mul (temps [16], slist [0], temps [17]); mpz_sub (temps [15], temps [16], z15705043200); mpz_mul (temps [14], slist [0], temps [15]); mpz_add (temps [13], z8335474560, temps [14]); mpz_mul (temps [12], slist [0], temps [13]); mpz_add (temps [11], temps [11], temps [12]); mpz_mul (temps [10], c, temps [11]); mpz_sub (temps [9], temps [10], z8806936320); /* zconst 9379426560 22f0eb100 */ mpz_mul_ui (temps [18], slist [0], 19353600); mpz_sub_ui (temps [17], temps [18], 899942400); mpz_neg (temps [17], temps [17]); mpz_mul (temps [16], slist [0], temps [17]); mpz_sub_ui (temps [15], temps [16], 3424619520); mpz_mul (temps [14], slist [0], temps [15]); mpz_add_ui (temps [13], temps [14], 710277120); mpz_mul (temps [12], slist [0], temps [13]); mpz_add (temps [11], z9379426560, temps [12]); mpz_mul (temps [10], slist [0], temps [11]); mpz_add (temps [9], temps [9], temps [10]); mpz_mul (temps [8], c, temps [9]); mpz_add (temps [7], z9472619520, temps [8]); /* zconst 16847631360 3ec327400 */ /* zconst 10668994560 27bebf000 */ mpz_mul_ui (temps [14], slist [0], 58060800); mpz_sub_ui (temps [13], temps [14], 2532096000); mpz_mul (temps [12], slist [0], temps [13]); mpz_add (temps [11], z10668994560, temps [12]); mpz_mul (temps [10], slist [0], temps [11]); mpz_sub (temps [9], temps [10], z16847631360); mpz_mul (temps [8], slist [0], temps [9]); mpz_add (temps [7], temps [7], temps [8]); mpz_mul (temps [6], c, temps [7]); mpz_sub_ui (temps [5], temps [6], 2127928320); mpz_mul_ui (temps [10], slist [0], 58060800); mpz_sub_ui (temps [9], temps [10], 1403136000); mpz_mul (temps [8], slist [0], temps [9]); mpz_add_ui (temps [7], temps [8], 3741050880); mpz_mul (temps [6], slist [0], temps [7]); mpz_add (temps [5], temps [5], temps [6]); mpz_mul (temps [4], c, temps [5]); mpz_add (temps [3], temps [3], temps [4]); mpz_mul_ui (temps [6], slist [0], 3870720); mpz_add_ui (temps [5], temps [6], 1033482240); mpz_mul (temps [4], slist [0], temps [5]); mpz_add (temps [3], temps [3], temps [4]); mpz_mul (temps [2], slist [1], temps [3]); /* high = 28 */

[7/12] /* part2 -> temps [3] */ mpz_mul_ui (temps [5], slist [0], 928972800); mpz_sub_ui (temps [4], temps [5], 1354752000); mpz_mul_ui (temps [7], slist [0], 464486400); mpz_sub_ui (temps [6], temps [7], 696729600); mpz_neg (temps [6], temps [6]); /* zconst 4833561600 1201a4c00 */ mpz_mul_ui (temps [21], c, 21772800); mpz_sub_ui (temps [20], temps [21], 203212800); mpz_mul_ui (temps [21], slist [0], 174182400); mpz_add (temps [20], temps [20], temps [21]); mpz_mul (temps [19], c, temps [20]); mpz_add_ui (temps [18], temps [19], 551577600); mpz_mul_ui (temps [21], slist [0], 522547200); mpz_sub_ui (temps [20], temps [21], 1248307200); mpz_mul (temps [19], slist [0], temps [20]); mpz_add (temps [18], temps [18], temps [19]); mpz_mul (temps [17], c, temps [18]); mpz_add_ui (temps [16], temps [17], 246758400); mpz_mul_ui (temps [21], slist [0], 696729600); mpz_sub_ui (temps [20], temps [21], 2554675200); mpz_mul (temps [19], slist [0], temps [20]); mpz_add_ui (temps [18], temps [19], 2090188800); mpz_mul (temps [17], slist [0], temps [18]); mpz_add (temps [16], temps [16], temps [17]); mpz_mul (temps [15], c, temps [16]); mpz_sub_ui (temps [14], temps [15], 3418329600); mpz_mul_ui (temps [21], slist [0], 348364800); mpz_sub_ui (temps [20], temps [21], 1741824000); mpz_mul (temps [19], slist [0], temps [20]); mpz_add_ui (temps [18], temps [19], 1393459200); mpz_mul (temps [17], slist [0], temps [18]); mpz_add_ui (temps [16], temps [17], 2815948800); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], temps [14], temps [15]); mpz_mul (temps [13], c, temps [14]); mpz_add (temps [12], z4833561600, temps [13]); /* zconst 8360755200 1f2570000 */ /* zconst 5167411200 134007000 */ mpz_mul_ui (temps [17], slist [0], 1161216000); mpz_sub (temps [16], z5167411200, temps [17]); mpz_mul (temps [15], slist [0], temps [16]); mpz_sub (temps [14], temps [15], z8360755200); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], temps [12], temps [13]); mpz_mul (temps [11], c, temps [12]); mpz_sub_ui (temps [10], temps [11], 1567641600); mpz_mul_ui (temps [13], slist [0], 696729600); mpz_sub_ui (temps [12], temps [13], 2090188800); mpz_neg (temps [12], temps [12]); mpz_mul (temps [11], slist [0], temps [12]); mpz_add (temps [10], temps [10], temps [11]); mpz_mul (temps [9], c, temps [10]); mpz_mul (temps [12], c, c); mpz_mul (temps [11], slist [1], temps [12]); mpz_mul_ui (temps [17], c, 29030400); mpz_sub_ui (temps [16], temps [17], 261273600); mpz_mul_ui (temps [17], slist [0], 174182400); mpz_add (temps [16], temps [16], temps [17]); mpz_mul (temps [15], c, temps [16]); mpz_add_ui (temps [14], temps [15], 783820800); mpz_mul_ui (temps [17], slist [0], 348364800); mpz_sub_ui (temps [16], temps [17], 1045094400); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], temps [14], temps [15]); mpz_mul (temps [13], c, temps [14]); mpz_sub_ui (temps [12], temps [13], 783820800); mpz_mul_ui (temps [17], slist [0], 232243200); mpz_sub_ui (temps [16], temps [17], 1045094400); mpz_mul (temps [15], slist [0], temps [16]); mpz_add_ui (temps [14], temps [15], 1567641600); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], temps [12], temps [13]); mpz_mul (temps [10], temps [11], temps [12]); mpz_add (temps [8], temps [9], temps [10]); mpz_mul (temps [7], slist [1], temps [8]); mpz_add (temps [6], temps [6], temps [7]); /* zconst 5904783360 15ff3d800 */ /* zconst 9252472320 2277d8600 */ mpz_mul_ui (temps [23], c, 5443200); mpz_sub_ui (temps [22], temps [23], 52617600); mpz_mul_ui (temps [23], slist [0], 54432000); mpz_add (temps [22], temps [22], temps [23]); mpz_mul (temps [21], c, temps [22]); mpz_add_ui (temps [20], temps [21], 131846400); mpz_mul_ui (temps [23], slist [0], 217728000); mpz_sub_ui (temps [22], temps [23], 435456000); mpz_mul (temps [21], slist [0], temps [22]); mpz_add (temps [20], temps [20], temps [21]); mpz_mul (temps [19], c, temps [20]); mpz_add_ui (temps [18], temps [19], 256677120); mpz_mul_ui (temps [23], slist [0], 435456000); mpz_sub_ui (temps [22], temps [23], 1349913600); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 742694400); mpz_mul (temps [19], slist [0], temps [20]); mpz_add (temps [18], temps [18], temps [19]); mpz_mul (temps [17], c, temps [18]); mpz_sub_ui (temps [16], temps [17], 1665740160); mpz_mul_ui (temps [23], slist [0], 435456000); mpz_sub_ui (temps [22], temps [23], 1857945600); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 1117670400); mpz_mul (temps [19], slist [0], temps [20]); mpz_add_ui (temps [18], temps [19], 2228567040); mpz_mul (temps [17], slist [0], temps [18]); mpz_add (temps [16], temps [16], temps [17]); mpz_mul (temps [15], c, temps [16]); mpz_add_ui (temps [14], temps [15], 1616630400); /* zconst 7575724800 1c38c6300 */ /* zconst 6178636800 170468400 */ mpz_mul_ui (temps [23], slist [0], 174182400); mpz_sub_ui (temps [22], temps [23], 958003200); mpz_mul (temps [21], slist [0], temps [22]); mpz_sub_ui (temps [20], temps [21], 241920000); mpz_mul (temps [19], slist [0], temps [20]); mpz_add (temps [18], z6178636800, temps [19]); mpz_mul (temps [17], slist [0], temps [18]); mpz_sub (temps [16], temps [17], z7575724800); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], temps [14], temps [15]); mpz_mul (temps [13], c, temps [14]); mpz_add_ui (temps [12], temps [13], 3752179200);

[8/12] /* zconst 8157542400 1e63a3800 */ /* zconst 5554483200 14b12b000 */ mpz_mul_ui (temps [19], slist [0], 1122508800); mpz_sub (temps [18], z5554483200, temps [19]); mpz_mul (temps [17], slist [0], temps [18]); mpz_sub (temps [16], temps [17], z8157542400); mpz_mul (temps [15], slist [0], temps [16]); mpz_add_ui (temps [14], temps [15], 1485388800); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], temps [12], temps [13]); mpz_mul (temps [11], c, temps [12]); mpz_sub (temps [10], temps [11], z9252472320); /* zconst 12275020800 2dba5f400 */ /* zconst 5786726400 158ea7000 */ mpz_mul_ui (temps [15], slist [0], 1161216000); mpz_sub (temps [14], temps [15], z5786726400); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], z12275020800, temps [13]); mpz_mul (temps [11], slist [0], temps [12]); mpz_add (temps [10], temps [10], temps [11]); mpz_mul (temps [9], c, temps [10]); mpz_add (temps [8], z5904783360, temps [9]); /* zconst 8387850240 1f3f47000 */ mpz_mul_ui (temps [13], slist [0], 38707200); mpz_sub_ui (temps [12], temps [13], 2922393600); mpz_neg (temps [12], temps [12]); mpz_mul (temps [11], slist [0], temps [12]); mpz_sub (temps [10], temps [11], z8387850240); mpz_mul (temps [9], slist [0], temps [10]); mpz_add (temps [8], temps [8], temps [9]); mpz_mul (temps [7], c, temps [8]); mpz_add (temps [6], temps [6], temps [7]); mpz_mul (temps [5], slist [1], temps [6]); mpz_add (temps [4], temps [4], temps [5]); mpz_mul_ui (temps [17], c, 29030400); mpz_sub_ui (temps [16], temps [17], 183859200); mpz_mul_ui (temps [17], slist [0], 174182400); mpz_add (temps [16], temps [16], temps [17]); mpz_mul (temps [15], c, temps [16]); mpz_add_ui (temps [14], temps [15], 145152000); mpz_mul_ui (temps [17], slist [0], 348364800); mpz_sub_ui (temps [16], temps [17], 774144000); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], temps [14], temps [15]); mpz_mul (temps [13], c, temps [14]); mpz_add_ui (temps [12], temps [13], 1112832000); mpz_mul_ui (temps [17], slist [0], 232243200); mpz_sub_ui (temps [16], temps [17], 812851200); mpz_mul (temps [15], slist [0], temps [16]); mpz_sub_ui (temps [14], temps [15], 58060800); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], temps [12], temps [13]); mpz_mul (temps [11], c, temps [12]); mpz_sub_ui (temps [10], temps [11], 2032128000); mpz_mul_ui (temps [13], slist [0], 928972800); mpz_sub_ui (temps [12], temps [13], 2748211200); mpz_neg (temps [12], temps [12]); mpz_mul (temps [11], slist [0], temps [12]); mpz_add (temps [10], temps [10], temps [11]); mpz_mul (temps [9], c, temps [10]); mpz_mul (temps [11], slist [1], c); mpz_mul_ui (temps [15], c, 58060800); mpz_sub_ui (temps [14], temps [15], 348364800); mpz_mul_ui (temps [15], slist [0], 232243200); mpz_add (temps [14], temps [14], temps [15]); mpz_mul (temps [13], c, temps [14]); mpz_add_ui (temps [12], temps [13], 522547200); mpz_mul_ui (temps [15], slist [0], 232243200); mpz_sub_ui (temps [14], temps [15], 696729600); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], temps [12], temps [13]); mpz_mul (temps [10], temps [11], temps [12]); mpz_add (temps [8], temps [9], temps [10]); mpz_mul (temps [7], slist [1], temps [8]); mpz_mul_ui (temps [20], c, 3628800); mpz_sub_ui (temps [19], temps [20], 24192000); mpz_mul_ui (temps [20], slist [0], 29030400); mpz_add (temps [19], temps [19], temps [20]); mpz_mul (temps [18], c, temps [19]); mpz_add_ui (temps [17], temps [18], 4032000); mpz_mul_ui (temps [20], slist [0], 87091200); mpz_sub_ui (temps [19], temps [20], 154828800); mpz_mul (temps [18], slist [0], temps [19]); mpz_add (temps [17], temps [17], temps [18]); mpz_mul (temps [16], c, temps [17]); mpz_add_ui (temps [15], temps [16], 290304000); mpz_mul_ui (temps [20], slist [0], 116121600); mpz_sub_ui (temps [19], temps [20], 329011200); mpz_mul (temps [18], slist [0], temps [19]); mpz_sub_ui (temps [17], temps [18], 67737600); mpz_mul (temps [16], slist [0], temps [17]); mpz_add (temps [15], temps [15], temps [16]); mpz_mul (temps [14], c, temps [15]); mpz_sub_ui (temps [13], temps [14], 473760000); mpz_mul_ui (temps [20], slist [0], 58060800); mpz_sub_ui (temps [19], temps [20], 232243200); mpz_mul (temps [18], slist [0], temps [19]); mpz_sub_ui (temps [17], temps [18], 377395200); mpz_mul (temps [16], slist [0], temps [17]); mpz_add_ui (temps [15], temps [16], 1354752000); mpz_mul (temps [14], slist [0], temps [15]); mpz_add (temps [13], temps [13], temps [14]); mpz_mul (temps [12], c, temps [13]); mpz_sub_ui (temps [11], temps [12], 846720000); mpz_mul_ui (temps [16], slist [0], 464486400); mpz_sub_ui (temps [15], temps [16], 1606348800); mpz_neg (temps [15], temps [15]); mpz_mul (temps [14], slist [0], temps [15]); mpz_sub_ui (temps [13], temps [14], 774144000); mpz_mul (temps [12], slist [0], temps [13]); mpz_add (temps [11], temps [11], temps [12]); mpz_mul (temps [10], c, temps [11]); mpz_add_ui (temps [9], temps [10], 1975680000); mpz_mul_ui (temps [12], slist [0], 928972800); mpz_sub_ui (temps [11], temps [12], 2709504000); mpz_mul (temps [10], slist [0], temps [11]); mpz_add (temps [9], temps [9], temps [10]); mpz_mul (temps [8], c, temps [9]); mpz_add (temps [6], temps [7], temps [8]); mpz_mul (temps [5], slist [2], temps [6]); mpz_add (temps [4], temps [4], temps [5]);

[9/12] /* zconst 5715763200 154afa000 */ /* zconst 5662406400 151817700 */ mpz_mul_ui (temps [25], c, 453600); mpz_sub_ui (temps [24], temps [25], 4536000); mpz_mul_ui (temps [25], slist [0], 5443200); mpz_add (temps [24], temps [24], temps [25]); mpz_mul (temps [23], c, temps [24]); mpz_add_ui (temps [22], temps [23], 10735200); mpz_mul_ui (temps [25], slist [0], 27216000); mpz_sub_ui (temps [24], temps [25], 47174400); mpz_mul (temps [23], slist [0], temps [24]); mpz_add (temps [22], temps [22], temps [23]); mpz_mul (temps [21], c, temps [22]); mpz_add_ui (temps [20], temps [21], 36281280); mpz_mul_ui (temps [25], slist [0], 72576000); mpz_sub_ui (temps [24], temps [25], 195955200); mpz_mul (temps [23], slist [0], temps [24]); mpz_add_ui (temps [22], temps [23], 81043200); mpz_mul (temps [21], slist [0], temps [22]); mpz_add (temps [20], temps [20], temps [21]); mpz_mul (temps [19], c, temps [20]); mpz_sub_ui (temps [18], temps [19], 204876000); mpz_mul_ui (temps [25], slist [0], 108864000); mpz_sub_ui (temps [24], temps [25], 406425600); mpz_mul (temps [23], slist [0], temps [24]); mpz_add_ui (temps [22], temps [23], 197164800); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 371589120); mpz_mul (temps [19], slist [0], temps [20]); mpz_add (temps [18], temps [18], temps [19]); mpz_mul (temps [17], c, temps [18]); mpz_add_ui (temps [16], temps [17], 133156800); mpz_mul_ui (temps [25], slist [0], 87091200); mpz_sub_ui (temps [24], temps [25], 420940800); mpz_mul (temps [23], slist [0], temps [24]); mpz_add_ui (temps [22], temps [23], 91929600); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 1411845120); mpz_mul (temps [19], slist [0], temps [20]); mpz_sub_ui (temps [18], temps [19], 1371565440); mpz_mul (temps [17], slist [0], temps [18]); mpz_add (temps [16], temps [16], temps [17]); mpz_mul (temps [15], c, temps [16]); mpz_add_ui (temps [14], temps [15], 935373600); mpz_mul_ui (temps [25], slist [0], 29030400); mpz_sub_ui (temps [24], temps [25], 174182400); mpz_mul (temps [23], slist [0], temps [24]); mpz_sub_ui (temps [22], temps [23], 263692800); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 2364364800); mpz_mul (temps [19], slist [0], temps [20]); mpz_sub_ui (temps [18], temps [19], 2937916800); mpz_mul (temps [17], slist [0], temps [18]); mpz_add_ui (temps [16], temps [17], 19595520); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], temps [14], temps [15]); mpz_mul (temps [13], c, temps [14]); mpz_sub_ui (temps [12], temps [13], 2034002880); /* zconst 5123865600 13167fc00 */ mpz_mul_ui (temps [21], slist [0], 270950400); mpz_sub_ui (temps [20], temps [21], 1475712000); mpz_neg (temps [20], temps [20]); mpz_mul (temps [19], slist [0], temps [20]); mpz_sub_ui (temps [18], temps [19], 1703116800); mpz_mul (temps [17], slist [0], temps [18]); mpz_sub_ui (temps [16], temps [17], 2165184000); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], z5123865600, temps [15]); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], temps [12], temps [13]); mpz_mul (temps [11], c, temps [12]); mpz_sub_ui (temps [10], temps [11], 251395200); /* zconst 5481745920 146bcce00 */ /* zconst 7190668800 1ac98e600 */ mpz_mul_ui (temps [17], slist [0], 677376000); mpz_sub_ui (temps [16], temps [17], 3506227200); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], z7190668800, temps [15]); mpz_mul (temps [13], slist [0], temps [14]); mpz_sub (temps [12], temps [13], z5481745920); mpz_mul (temps [11], slist [0], temps [12]); mpz_add (temps [10], temps [10], temps [11]); mpz_mul (temps [9], c, temps [10]); mpz_add (temps [8], z5662406400, temps [9]); mpz_mul_ui (temps [15], slist [0], 19353600); mpz_sub_ui (temps [14], temps [15], 541900800); mpz_neg (temps [14], temps [14]); mpz_mul (temps [13], slist [0], temps [14]); mpz_sub_ui (temps [12], temps [13], 782853120); mpz_mul (temps [11], slist [0], temps [12]); mpz_sub_ui (temps [10], temps [11], 3935554560); mpz_mul (temps [9], slist [0], temps [10]); mpz_add (temps [8], temps [8], temps [9]); mpz_mul (temps [7], c, temps [8]); mpz_sub (temps [6], temps [7], z5715763200); /* zconst 8872980480 210def000 */ mpz_mul_ui (temps [11], slist [0], 77414400); mpz_sub_ui (temps [10], temps [11], 3286886400); mpz_mul (temps [9], slist [0], temps [10]); mpz_add (temps [8], z8872980480, temps [9]); mpz_mul (temps [7], slist [0], temps [8]); mpz_add (temps [6], temps [6], temps [7]); mpz_mul (temps [5], c, temps [6]); mpz_add (temps [4], temps [4], temps [5]); mpz_mul (temps [3], slist [2], temps [4]); /* high = 25 */

[10/12] /* part3 -> temps [4] */ /* zconst 5741832320 1563d6880 */ mpz_mul_ui (temps [33], c, 14175); mpz_sub_ui (temps [32], temps [33], 189000); mpz_mul_ui (temps [33], slist [0], 226800); mpz_add (temps [32], temps [32], temps [33]); mpz_mul (temps [31], c, temps [32]); mpz_add_ui (temps [30], temps [31], 812700); mpz_mul_ui (temps [33], slist [0], 1587600); mpz_sub_ui (temps [32], temps [33], 2721600); mpz_mul (temps [31], slist [0], temps [32]); mpz_add (temps [30], temps [30], temps [31]); mpz_mul (temps [29], c, temps [30]); mpz_add_ui (temps [28], temps [29], 52080); mpz_mul_ui (temps [33], slist [0], 6350400); mpz_sub_ui (temps [32], temps [33], 16783200); mpz_mul (temps [31], slist [0], temps [32]); mpz_add_ui (temps [30], temps [31], 9903600); mpz_mul (temps [29], slist [0], temps [30]); mpz_add (temps [28], temps [28], temps [29]); mpz_mul (temps [27], c, temps [28]); mpz_sub_ui (temps [26], temps [27], 11050550); mpz_mul_ui (temps [33], slist [0], 15876000); mpz_sub_ui (temps [32], temps [33], 57456000); mpz_mul (temps [31], slist [0], temps [32]); mpz_add_ui (temps [30], temps [31], 49215600); mpz_mul (temps [29], slist [0], temps [30]); mpz_add_ui (temps [28], temps [29], 6347040); mpz_mul (temps [27], slist [0], temps [28]); mpz_add (temps [26], temps [26], temps [27]); mpz_mul (temps [25], c, temps [26]); mpz_add_ui (temps [24], temps [25], 28664160); mpz_mul_ui (temps [33], slist [0], 25401600); mpz_sub_ui (temps [32], temps [33], 117936000); mpz_mul (temps [31], slist [0], temps [32]); mpz_add_ui (temps [30], temps [31], 125193600); mpz_mul (temps [29], slist [0], temps [30]); mpz_add_ui (temps [28], temps [29], 65677920); mpz_mul (temps [27], slist [0], temps [28]); mpz_sub_ui (temps [26], temps [27], 128034480); mpz_mul (temps [25], slist [0], temps [26]); mpz_add (temps [24], temps [24], temps [25]); mpz_mul (temps [23], c, temps [24]); mpz_add_ui (temps [22], temps [23], 23121812); mpz_mul_ui (temps [33], slist [0], 25401600); mpz_sub_ui (temps [32], temps [33], 145152000); mpz_mul (temps [31], slist [0], temps [32]); mpz_add_ui (temps [30], temps [31], 162388800); mpz_mul (temps [29], slist [0], temps [30]); mpz_add_ui (temps [28], temps [29], 281756160); mpz_mul (temps [27], slist [0], temps [28]); mpz_sub_ui (temps [26], temps [27], 588183120); mpz_mul (temps [25], slist [0], temps [26]); mpz_add_ui (temps [24], temps [25], 233335200); mpz_mul (temps [23], slist [0], temps [24]); mpz_add (temps [22], temps [22], temps [23]); mpz_mul (temps [21], c, temps [22]); mpz_sub_ui (temps [20], temps [21], 231117040); mpz_mul_ui (temps [33], slist [0], 14515200); mpz_sub_ui (temps [32], temps [33], 99187200); mpz_mul (temps [31], slist [0], temps [32]); mpz_add_ui (temps [30], temps [31], 76204800); mpz_mul (temps [29], slist [0], temps [30]); mpz_add_ui (temps [28], temps [29], 609799680); mpz_mul (temps [27], slist [0], temps [28]); mpz_sub_ui (temps [26], temps [27], 1329558720); mpz_mul (temps [25], slist [0], temps [26]); mpz_add_ui (temps [24], temps [25], 628068000); mpz_mul (temps [23], slist [0], temps [24]); mpz_add_ui (temps [22], temps [23], 363318480); mpz_mul (temps [21], slist [0], temps [22]); mpz_add (temps [20], temps [20], temps [21]); mpz_mul (temps [19], c, temps [20]); mpz_add_ui (temps [18], temps [19], 303711135); mpz_mul_ui (temps [33], slist [0], 3628800); mpz_sub_ui (temps [32], temps [33], 29030400); mpz_mul (temps [31], slist [0], temps [32]); mpz_sub_ui (temps [30], temps [31], 37497600); mpz_mul (temps [29], slist [0], temps [30]); mpz_add_ui (temps [28], temps [29], 662054400); mpz_mul (temps [27], slist [0], temps [28]); mpz_sub_ui (temps [26], temps [27], 1443304800); mpz_mul (temps [25], slist [0], temps [26]); mpz_add_ui (temps [24], temps [25], 413965440); mpz_mul (temps [23], slist [0], temps [24]); mpz_add_ui (temps [22], temps [23], 1801674000); mpz_mul (temps [21], slist [0], temps [22]); mpz_sub_ui (temps [20], temps [21], 1744031520); mpz_mul (temps [19], slist [0], temps [20]); mpz_add (temps [18], temps [18], temps [19]); mpz_mul (temps [17], c, temps [18]); mpz_add_ui (temps [16], temps [17], 491453160);

[11/12] /* zconst 4386013920 1056d42e0 */ mpz_mul_ui (temps [29], slist [0], 38707200); mpz_sub_ui (temps [28], temps [29], 287884800); mpz_neg (temps [28], temps [28]); mpz_mul (temps [27], slist [0], temps [28]); mpz_sub_ui (temps [26], temps [27], 517708800); mpz_mul (temps [25], slist [0], temps [26]); mpz_sub_ui (temps [24], temps [25], 719913600); mpz_mul (temps [23], slist [0], temps [24]); mpz_add_ui (temps [22], temps [23], 3666835200); mpz_mul (temps [21], slist [0], temps [22]); mpz_sub (temps [20], temps [21], z4386013920); mpz_mul (temps [19], slist [0], temps [20]); mpz_add_ui (temps [18], temps [19], 1191840000); mpz_mul (temps [17], slist [0], temps [18]); mpz_add (temps [16], temps [16], temps [17]); mpz_mul (temps [15], c, temps [16]); mpz_sub_ui (temps [14], temps [15], 1855381864); mpz_mul_ui (temps [25], slist [0], 122572800); mpz_sub_ui (temps [24], temps [25], 891072000); mpz_mul (temps [23], slist [0], temps [24]); mpz_add_ui (temps [22], temps [23], 2710512000); mpz_mul (temps [21], slist [0], temps [22]); mpz_sub_ui (temps [20], temps [21], 3678716160); mpz_mul (temps [19], slist [0], temps [20]); mpz_add_ui (temps [18], temps [19], 179323200); mpz_mul (temps [17], slist [0], temps [18]); mpz_add_ui (temps [16], temps [17], 3719604000); mpz_mul (temps [15], slist [0], temps [16]); mpz_add (temps [14], temps [14], temps [15]); mpz_mul (temps [13], c, temps [14]); mpz_add_ui (temps [12], temps [13], 1227340800); /* zconst 7185104640 1ac43ff00 */ /* zconst 7910091840 1d77a6c40 */ mpz_mul_ui (temps [23], slist [0], 4838400); mpz_sub_ui (temps [22], temps [23], 75479040); mpz_neg (temps [22], temps [22]); mpz_mul (temps [21], slist [0], temps [22]); mpz_add_ui (temps [20], temps [21], 11289600); mpz_mul (temps [19], slist [0], temps [20]); mpz_sub_ui (temps [18], temps [19], 2663646720); mpz_mul (temps [17], slist [0], temps [18]); mpz_add (temps [16], z7910091840, temps [17]); mpz_mul (temps [15], slist [0], temps [16]); mpz_sub (temps [14], temps [15], z7185104640); mpz_mul (temps [13], slist [0], temps [14]); mpz_add (temps [12], temps [12], temps [13]); mpz_mul (temps [11], c, temps [12]); mpz_add_ui (temps [10], temps [11], 2760323440); /* zconst 5737536000 155fbda00 */ mpz_mul_ui (temps [19], slist [0], 25804800); mpz_sub_ui (temps [18], temps [19], 972518400); mpz_mul (temps [17], slist [0], temps [18]); mpz_add_ui (temps [16], temps [17], 4195215360); mpz_mul (temps [15], slist [0], temps [16]); mpz_sub (temps [14], temps [15], z5737536000); mpz_mul (temps [13], slist [0], temps [14]); mpz_add_ui (temps [12], temps [13], 766308480); mpz_mul (temps [11], slist [0], temps [12]); mpz_add (temps [10], temps [10], temps [11]); mpz_mul (temps [9], c, temps [10]); mpz_sub (temps [8], temps [9], z5741832320); /* zconst 7885409280 1d601cc00 */ mpz_mul_ui (temps [15], slist [0], 6451200); mpz_sub_ui (temps [14], temps [15], 669312000); mpz_neg (temps [14], temps [14]); mpz_mul (temps [13], slist [0], temps [14]); mpz_sub_ui (temps [12], temps [13], 3513619200); mpz_mul (temps [11], slist [0], temps [12]); mpz_add (temps [10], z7885409280, temps [11]); mpz_mul (temps [9], slist [0], temps [10]); mpz_add (temps [8], temps [8], temps [9]); mpz_mul (temps [7], c, temps [8]); mpz_add_ui (temps [6], temps [7], 2494340352); /* zconst 5494809600 147842400 */ mpz_mul_ui (temps [13], slist [0], 1612800); mpz_sub_ui (temps [12], temps [13], 119347200); mpz_mul (temps [11], slist [0], temps [12]); mpz_add_ui (temps [10], temps [11], 2311787520); mpz_mul (temps [9], slist [0], temps [10]); mpz_sub (temps [8], temps [9], z5494809600); mpz_mul (temps [7], slist [0], temps [8]); mpz_add (temps [6], temps [6], temps [7]); mpz_mul (temps [5], c, temps [6]); mpz_add_ui (temps [4], temps [5], 2832168960); mpz_mul_ui (temps [11], slist [0], 645120); mpz_sub_ui (temps [10], temps [11], 5160960); mpz_mul (temps [9], slist [0], temps [10]); mpz_add_ui (temps [8], temps [9], 48015360); mpz_mul (temps [7], slist [0], temps [8]); mpz_sub_ui (temps [6], temps [7], 2009825280); mpz_mul (temps [5], slist [0], temps [6]); mpz_add (temps [4], temps [4], temps [5]); /* high = 33 */ mpz_add (temps [0], temps [2], temps [3]); mpz_add (temps [0], temps [0], temps [4]); mpz_mul (temps [1], temps [0], c); mpz_divexact_ui (temps [0], temps [1], 464486400); mpz_add (ctx->incs [2], ctx->counts [2], ctx->counts [3]); mpz_mul (temps [1], slist [3], ctx->counts [3]); mpz_add (ctx->incs [3], temps [0], temps [1]); }

[12/12] int work (char *s) { mpz_t n, R; int e; struct groups gs; e = mpz_init_set_str (n, s, 10); if (e != 0) { printf ("invalid index: \"%s\"\n", s); mpz_clear (n); return 1; } mpz_init2 (R, ENOUGH_BITS); gs.count = 4; init_groups (&gs); gs.layers [0].count = l1_count; gs.layers [0].increment = l1_increment; gs.layers [1].count = l2_count; gs.layers [1].increment = l2_increment; gs.layers [2].count = l3_count; gs.layers [2].increment = l3_increment; gs.layers [3].count = l4_count; gs.layers [3].increment = l4_increment; loop_zg (n, &gs, R); gmp_printf ("%Zd %zu\n", R, mpz_sizeinbase (R, 2)); mpz_clear (n); mpz_clear (R); clear_groups (&gs); return 0; } int main (int argc, char **argv) { int ret; init_zconst (); ret = work (argc < 2 ? "10000000000" : argv [1]); clear_zconst (); return ret; } The time for the u512 overflow R(10^78) is down from under 5 seconds >>18649 to under one second, yielding a speedup by a factor of 5. $ time bin/gmphofs "1$(printf "%078d" 0)" levels: 5 500000000000000000000000000000000000000942809041582063365839428238027648024522282176747574937498915190189902751724515553262387695706471431607983773634306949 518 real 0m0,921s user 0m0,921s sys 0m0,000s $ python3 [...] >>> 4790 / 921 5.2008686210640604 >>> I think under one second for R(10^78) might be an acceptable reference time, for a while.

>>14060 I believe your pic is the graphics coding equivalent of pic related. It's generated with something like color = sin(x)*cos(y), mapped to an indexed color space. I've never seen it discussed by any academic for example on a pop math channel.