Clean up recursive function
I've written a recursive function to approximate the count of k < K with 6*k +-1 both composite. It's accurate enough that now I want to:
- Clean it up*
- Use any theory/sage tools I should have started with
- Derive an equation, upper and lower bounds, etc.
- Write a paper on it
I confess that I'm stumped by the recursion to equation/theory part. All my attempts at combinatoric, iterative, etc. failed. Can you help me clean it up and put it in it's best form? I may also have to do a print trace to help me see the math. I'm worried a code to math theory question won't fit anywhere.
The line ki = ceil(k/p) If I use:
- ki = k/p Treats k's as rationals not integers, range divisions always 'overlap'.
- ki = ceil(k/p) k's are integers, but still divisions overlap too often.
- ki = floor(k/p) Never overlap, but 'gaps' mean some k not counted.
This is a coding detail that affects the underlying and any derived math.
# twin composite counting function, recursive def twinc(k, flag, plist): # no more primes to try OR range < prime if len(plist) == 0: return 0 if k < plist: return 0 p = plist; del plist ki = ceil(k/p) # half of twin is composite # ki = count of (k + pq): p | other half # try other primes where p misses if flag > 0: return ki + twinc(k - ki, flag, plist) # previous primes all missed, divide range into: # p | -1 or +1, p missed return twinc(ki, 1, copy(plist))*2 + twinc(k - ki*2, 0, plist) # end twinc
pmax = 201 plist = prime_range(5, pmax) kmax = 5*max(plist) + 1 tlist = [(k, twinc(k, 0, copy(plist))) for k in range(5, kmax, 5)] P = list_plot(tlist, plotjoined=True) # direct plot fails, use stepping? P.show()
Update: Originally the flag was two bits for 0..3. While I was writing pseudo-code I realized that a boolean parameter can be eliminated with a pair of recursive functions. This looks closer to something that can be turned into an equation, yes?
def twinCa(k, [p: plist]): ki = round(k/p) return k < p? 0: twinCa(k - ki*2, plist) + twinCb(ki, plist)*2 def twinCb(k, [p: plist]): ki = round(k/p) return k < p? 0: ki + twinCb(k - ki, plist)