Clean up recursive function

edit

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)