# 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:

1. Clean it up*
2. Use any theory/sage tools I should have started with
3. Derive an equation, upper and lower bounds, etc.
4. 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.

*Examples
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.

Function

# 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[0]:
return 0

p = plist[0]; del plist[0]
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


Test program

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?

Alternate (pseudo-code)

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)

edit retag close merge delete

I've not read your code in detail, but concerning k/p: using ceil versus floor is not an implementation detail, they do not perform the same math operation! Use floor? and ceil? to get the documentation. And if you want the euclidean division of k by p, the right way to get it is p // k.
Yes, that is the point I was trying to make. The function approximates the number of expected composites by dividing the range. But using floor implies the composite is at the end of each section and ceil implies it's at the beginning. Since they are both cumulative I'll probably use round, even though the others seem like a cleaner translation to the underlying density/distribution. Of course floor would give a good lower bound and ceil a good upper - as long as the error doesn't diverge.