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

***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)
```

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.