Maybeso83's profile - activity

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.

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)

I am working on a counting function involving prime numbers and managed to construct a sage function that is correct:

def sum_X_subs(plist, plen):
    psum = 0
    for k in range(plen-1):
        ptwo = 2^(plen-k) - 2
        psum += sum(mul(psub)*ptwo for psub in Subsets(plist, k))
    return psum
# end sum_X_subs

flist = []
plist = [5]; plen = 1
for p in prime_range(7, 41+1):
    plist.append(p)
    plen += 1
    flist.append(sum_X_subs(plist, plen))
flist
# [2, 52, 1162, 28196, 712250, 20379100, 696733730, 25016026280, 1042543611410, 47439135073960]

Now I need a lower bounds function for this: $$ S = prime~range(5, n)$$ $$f(S) = \sum_{k=0}^{|S|-1} {\left(\left(2^{|S|-k}-2\right)\sum_{T\in \mathcal{P}_{k}(S)} {\prod_{t\in T} {t}}\right)}$$ None of the literature I've found covers anything close to it. (I have a vague memory of an Euler project problem that involved subset products).

Ironically, the LabelledBinaryTree() command only lets you label the root. But you can apply a canonical labeling and then map those labels to the desired ones.

import re
words = re.compile('\w+')

# keep int labels as ints
def safeInt(word):
    try:
        return Integer(word)
    except TypeError:
        return word

# labelled binary tree from string
def lbt_from_str(lbstr):
    blank = words.sub('', lbstr)
    work = LabelledBinaryTree(blank).canonical_labelling()
    wstr = '%s'%work
    canon = [Integer(w) for w in words.findall(wstr)]
    labels = [safeInt(w) for w in words.findall(lbstr)]
    fixit = dict((canon[i], labels[i]) for i in range(len(canon)))
    return work.map_labels(lambda z: fixit[z])
# end lbt_from_str

This doesn't check if the string has a valid structure, and it expects all nodes to have labels.

You could replace ' [' with 'None['

Testing:

t = Permutation([1, 3, 2, 6, 5, 7, 4]).increasing_tree()
print t
show(t)
s = '%s'%t
test = lbt_from_str(s)
print test
show(test)

I think when the png is being constructed it doesn't handle the interaction between the opacity of the implicit plot and the point. Try adjusting your opacity. I also tested this:

distance = 56; x,y,z = var('x y z')
mypic = implicit_plot3d(x^2 + y^2 + z^2 == distance, (x, -10, 10), (y, -10, 10), (z, -10, 10), opacity = 0.5)
for k in range(-9, 10, 2): mypic = mypic + point3d(vector([5, k, 4]), color = 'red', size = 10)
mypic.save('mypic.png', compress = false)
mypic

The inline image shows a line of six red points from up-left to down-right ( \ ), with a center gap where four are completely blocked by the sphere. The saved png shows all ten points in a line from down-left to up-right ( / ), with the center four slightly hazy. So the png is also swapping x with y for the points (the axis labels are correct).

Could someone post a wiki suggestion regarding redirects? Then my question can be de-wiki-fied.

Thank you, I will make a habit of reporting both the old and new links to the site whenever I get redirected.

With the help of a PARI programmer and OEIS regular I cleaned up my code:

wief = DiGraph([prime_range(3600), lambda p, q: power_mod(p, q-1, q^2)==1])
sc = [[0]] + [sorted(c[1:], reverse=1) for c in wief.all_simple_cycles()]
tbl = [sorted(filter(lambda c: len(c)==n, sc))[0] for n in range(1, len(sc[-1]))]
for t in tbl: print 'n=%d:'%(len(t)), ', '.join("%s"%i for i in t)

and posted it to OEIS sequence A271100 Wieferich n-tuples. I used it to increase the number of listed terms from 15 to 300. We'll probably bump that to nearly 10k soon : )
*For the community, I'd like to add mouse-over hints describing what each piece of code does but I'm a noob.

The Sage/Sagemath webmasters need to do the fix described below for all search engines:
Some of the online sage resources that were moved still aren't visible to search engines yet.
(Examples for Google MyWay, links trashed because I'm new)

When I search for documentation: "sage intersection of lists"
Most of the results are similar to:
Sets — Sage Reference Manual v7.2: Basic Structures
http://www.sagemath.org/doc/reference/sage/sets/set.html

When I click the links I get:
Redirecting ...
http://doc.sagemath.org/html/en/reference/structure/sage/sets/set.html

FIX: Google Search Help* says that the sage webmaster needs to:
Create a "Google Webmasters Tools Account" - https://www.google.com/webmasters/tools/
then use the Fetch tool to get the spiders to you site so it can be indexed correctly -
https://support.google.com/webmasters/answer/6066468?hl=en&ref_topic=6066464

*https://productforums.google.com/forum/#!topic/websearch/st2z6P9Q7F8;context-place=topicsearchin/websearch/search returns old site url

I wrote a worksheet that defines a digraph and extracts all of its simple cycles.

pr = prime_range(3600)
wief = DiGraph([pr, lambda p, q: power_mod(p, q - 1, q^2) == 1]) # could take a while!
sc = wief.all_simple_cycles()
sc.append([0, 0])

Then it uses mostly procedural code to:

1. Remove the duplicate endpoints
2. Reverse-sort the cycle
3. Group cycles by length
4. Sort the group of cycles
5. Return the first (lexicographically earliest) cycle of each length

maxn = 0
for c in sc:
    del c[0] # duplicate
    maxn = max(maxn, len(c))
    c.sort(reverse=True) # reverse lexicographic
# group by length
tbl = [[c for c in sc if len(c) == n] for n in range(1,maxn)]
n = 1
for t in tbl: # sort tuples by max element, get smallest tuple
    t.sort()
    print 'n=%d:'%(n),
    # why can't join handle a list of integers?
    print ', '.join(map(lambda m: '%s'%(m), t[0]))
    n = n + 1

The worksheet is adequate - it performs step 5 correctly, but a listing of the code will be posted as the PROG paragraph of an OEIS entry. So I would really like to reduce the numbered steps to as few lines as possible. Are there more graceful ways to tweak the output of one function before passing it on to another?