New years resolutions

Can't wait.
Try it whilst it's fresh rather than waiting for a new year.

Extensible sieve of Eratosthenes

I was creating a Carmichael triplets generator and needed a function to tell me if a number was prime. It wasn't immediately clear what the largest prime for checking would be so I decided to look into creating an extensible prime generator - If given a number to check for prime-ness that was outside its range, it would attempt to extend its range.

I stopped modifying the class as soon as it was working as the Carmichael generator task then ran in acceptable time.

class Isprime():
    '''
    Extensible sieve of Erastosthenes
    
    >>> isprime.check(11)
    True
    >>> isprime.multiples
    {2, 4, 6, 8, 9, 10}
    >>> isprime.primes
    [2, 3, 5, 7, 11]
    >>> isprime(13)
    True
    >>> isprime.multiples
    {2, 4, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 18, 20, 21, 22}
    >>> isprime.primes
    [2, 3, 5, 7, 11, 13, 17, 19]
    >>> isprime.nmax
    22
    >>> 
    '''
    multiples = {2}
    primes = [2]
    nmax = 2

    def __init__(self, nmax):
        if nmax > self.nmax:
            self.check(nmax)

    def check(self, n):
        if type(n) == float:
            if not n.is_integer(): return False
            n = int(n)
        multiples = self.multiples
        if n <= self.nmax:
            return n not in multiples
        else:
            # Extend the sieve
            primes, nmax = self.primes, self.nmax
            newmax = max(nmax*2, n)
            for p in primes:
                multiples.update(range(p*((nmax + p + 1) // p), newmax+1, p))
            for i in range(nmax+1, newmax+1):
                if i not in multiples:
                    primes.append(i)
                    multiples.update(range(i*2, newmax+1, i))
            self.nmax = newmax
            return n not in multiples

    __call__ = check

Now I wouldn't be surprised if some of the loop indices are not optimum and are maybe recomputing a few values that don't need to be, and this code is only efficient enough for my purposes; (a couple of Rosetta Code entries in other languages state that they are based on the Python but make changes so that is an indication).

That grouping by zip/iter trick explained

In Python 2.7 you can get the following neat chunking of items in a list through application of zip and iter:

Python 2.7.3 (default, Apr 10 2012, 23:24:47) [MSC v.1500 64 bit (AMD64)] on win32
Type "copyright", "credits" or "license()" for more information.
>>> items, chunk = [1,2,3,4,5,6,7,8,9], 3
>>> zip(*[iter(items)]*chunk)
[(1, 2, 3), (4, 5, 6), (7, 8, 9)]
>>> 

It is not as neat in Python 3.X as you have to wrap the expression in list(..) to see the final list.

I thought I had to think hard about what was going on here so here is my annottated Python 3.x session that explains what is going on:

Python 3.3.0 (v3.3.0:bd8afb90ebf2, Sep 29 2012, 10:57:17) [MSC v.1600 64 bit (AMD64)] on win32
Type "copyright", "credits" or "license()" for more information.
>>> from pprint import pprint as pp
>>> items, chunk = [1,2,3,4,5,6,7,8,9], 3
>>> items
[1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> chunk
3
>>> # A quick recap of what an iterator does 
>>> i = iter(items)
>>> i.__next__()
1
>>> i.__next__()
2
>>> i.__next__()
3
>>> i.__next__()
4
>>> i.__next__()
5
>>> i.__next__()
6
>>> i.__next__()
7
>>> i.__next__()
8
>>> i.__next__()
9
>>> i.__next__()
Traceback (most recent call last):
  File "", line 1, in 
    i.__next__()
StopIteration
>>>
>>> # Three copies of the same iterator
>>> threei = [iter(items)] * chunk
>>> pp(threei)
[object at 0x0000000003259D30>,
 object at 0x0000000003259D30>,
 object at 0x0000000003259D30>]
>>> # (notice that the addresses above are the same)
>>>
>>> # Lets break-out the three copies
>>> i, j, k = threei
>>> i
object at 0x0000000003259D30>
>>> j
object at 0x0000000003259D30>
>>> k
object at 0x0000000003259D30>
>>> # i,j and k are the same iterator
>>> i.__next__()
1
>>> j.__next__()
2
>>> k.__next__()
3
>>> # ...
>>>
>>> # Lets try that again with zip:
>>> i, j, k = threei = [iter(items)] * chunk
>>> list(zip(i, j, k))
[(1, 2, 3), (4, 5, 6), (7, 8, 9)]
>>> # zip pulls from i, j, k once to form its first tuple of the list
>>> # then pulls from i, j, k again to form its second tuple of the list
>>> # and so on ...
>>> # Because i,j, and k are names for the same iterator we get the
>>> # given grouping action as output.
>>>
>>> # We could use *threei instead of i,j,k as the arguments to zip:
>>> i, j, k = threei = [iter(items)] * chunk
>>> list(zip(*threei))
[(1, 2, 3), (4, 5, 6), (7, 8, 9)]
>>>
>>> # Or we could use the expression generating threei:
>>> i, j, k = threei = [iter(items)] * chunk
>>> list(zip(*[iter(items)] * chunk))
[(1, 2, 3), (4, 5, 6), (7, 8, 9)]
>>> # But note that we are relying on the unary '*' operator
>>> # binding less strongly than binary '*'
>>>

For me it remains a curiosity as I doubt I will have a use for its functionality. (List-size must be a multiple of chunk-size; its functionality is obscure).

Some identities for Python int.to_bytes and int.from_bytes

I read the documentation for int.to_bytes and int.from_bytes and thought that what was missing was a description of what they did in terms of other Python code, rather like what is done for several of the generators in the docs for the itertools module; (I was completing this RC task at the time).

I came up with the following code snippet that I would have found helpful and that I hope helps you.

>>> n = 2491969579123783355964723219455906992268673266682165637887
>>> length = 25
>>> n2bytesbig    = n.to_bytes(length, 'big')
>>> n2byteslittle = n.to_bytes(length, 'little')
>>> assert n2bytesbig    == bytes( (n >> i*8) & 0xff for i in reversed(range(length)))
>>> assert n2byteslittle == bytes( (n >> i*8) & 0xff for i in range(length))
>>> assert n == sum( n2bytesbig[::-1][i] << i*8 for i in range(length) )
>>> assert n == sum( n2byteslittle[i]    << i*8 for i in range(length) )
>>> assert n == int.from_bytes(n2bytesbig,    byteorder='big')
>>> assert n == int.from_bytes(n2byteslittle, byteorder='little')
>>>

Topswop and a convoluted multiple assignments warning

I was messing around with John Conways Topswop problem:

Given a starting permutation of 1..n items, consider the leftmost value and call it L:

• If L is not 1 then reverse the leftmost L values in the sequence.
• Repeat, counting the number of reversions required until L becomes 1

Repeat the process for all permutations of 1..n items and report the maximum number of reversions required amongst all the permutations.

I initially came up with this doodle:

>>> def f1(p):
    i = 0
    while p[0] != 1:
        i += 1
        print(i,p)
        p[:p[0]] = p[:p[0]][::-1]
    return i

>>> f1([4,2,1,5,3])
1 [4, 2, 1, 5, 3]
2 [5, 1, 2, 4, 3]
3 [3, 4, 2, 1, 5]
4 [2, 4, 3, 1, 5]
5 [4, 2, 3, 1, 5]
5
>>> def f1(p):
    i = 0
    while p[0] != 1:
        i += 1
        #print(i,p)
        p[:p[0]] = p[:p[0]][::-1]
    return i

>>> f1([4,2,1,5,3])
5
>>> from itertools import permutations
>>> def fannkuch(n):
    return max(f1(list(p)) for p in permutations(range(1, n+1)))

>>> for n in range(1, 11): print(n,fannkuch(n))

1 0
2 1
3 2
4 4
5 7
6 10
7 16
8 22
9 30
10 38
>>>

Range

Python usually uses ranges 0..n-1 so I made that modification:

>>> def f1(p):
    i = 0
    while p[0]:
        i += 1
        p[:p[0] + 1] = p[:p[0] + 1][::-1]
    return i

>>> def fannkuch(n):
    return max(f1(list(p)) for p in permutations(range(n)))

Speed

I had to wait for fannkuch(10) so decided to do some quick optimization and created name p0 to mop up some array lookups:

>>> def f1(p):
    i, p0 = 0, p[0]
    while p0:
        i  += 1
        p0 += 1
        p[:p0] = p[:p0][::-1]
        p0  = p[0]
    return i

Those assignments

It was at this stage that I thought "Ooh, I haven't used that new assignment in Python 3 where you can use *name on the left hand side of an assignment", plus, I could add it in to an already complex assignment  just for the hell of it.

I am creating p, so I should be able to assign p[0] to p0 and discard p[1:] for the purposes of the assignment into *_. I tried:

>>> def f1(p):
    i, p0 = 0, p[0]
    while p0:
        i  += 1
        p0 += 1
        p0, *_ = p[:p0] = p[:p0][::-1]
    return i

>>> for n in range(1, 11): print(n,fannkuch(n))

1 0
2 1
3 2
Traceback (most recent call last):
  File "", line 1, in 
    for n in range(1, 11): print(n,fannkuch(n))
  File "", line 2, in fannkuch
    return max(f1(list(p)) for p in permutations(range(n)))
  File "", line 2, in 
    return max(f1(list(p)) for p in permutations(range(n)))
  File "", line 6, in f1
    p0, *_ = p[:p0] = p[:p0][::-1]
KeyboardInterrupt
>>>

I had to interrupt it as it failed to return. Something was wrong!

I thought it might be the order of assignment and tried:

>>> def f1(p):
    i, p0 = 0, p[0]
    while p0:
        i  += 1
        p0 += 1
        p[:p0] = p0, *_ = p[:p0][::-1]
    return i

>>> for n in range(1, 11): print(n,fannkuch(n))

1 0
2 1
3 2
4 4
5 7
6 10
7 16
8 22
9 30
10 38
>>>

Yay Success!!

But I did this so you don't have to :-)

It is a few optimizations too far, and, when you look in the docs they state:

"
WARNING: Although the definition of assignment implies that overlaps between the left-hand side and the right-hand side are ‘safe’ (for example a, b= b, a swaps two variables), overlaps within the collection of assigned-to variables are not safe! For instance, the following program prints [0, 2]:

x = [0, 1]
i = 0
i, x[i] = 1, 2
print(x)

"

RGB spectrum. 3D to 2D with smooth transition?

I was messing about with HTML red, green, blue triplets of 0..255 values each that are used for specifying HTML colours (with a u as I am English). I decided that I wanted to be able to take a spread of rgb values and arrange them in a spectrum.

Easily said, but hard to do. I went down a blind alley of converting to Hue, Saturation, Luminance; but could do nothing with that. I decided to ignore our eyes different sensitivities and any difference in 'purity' of pure red/green/blue from monitors and decided to try and work out mathematically how to create such a spectra.

The values for the three colours are independently variable over their ranges and so I got to thinking of the available colours as a colour cube. Think of a a 3D graph with three perpendicular axes red, green and blue. You can do the physics trick of taking your hand; stick the thumb straight up, the finger directly below the thumb straight out in front and the second finger below the thumb at right angles to them both pointing towards your body. There you have it! three perpendicular axes for red, green and blue, with the origin being towards your palm.  

The available colours fill the cube of integers (0, 0, 0) up to (255, 255, 255). I want to create a linear 'spectrum'.

If you consider the three points R, G, B, = (255, 0, 0), (0, 255, 0), (0, 0, 255) respectively,  they are the points of maximum r,g,b intensity, unsullied by other colours. If you cut the cube to expose the plane linking those three points then I would ideally like a spectra that  traversed R->G->B (->R); but what about the off-plane colour triplets?

If given a list of colour triplets to sort into a spectrum I decided to 

1. Find out which of the three cardinal points/primary colours R, G, or B each triplet was closest-to.
2. Sort all the triplets on their closeness to a primary.
3. Extract three lists of the triplets for each primary colour in such a way that they would be sorted with the primaries near the centre of their list; and those with more of the next primary to the 'right' and to the previous primary to the 'left'.
   For example the Red list would have reddy-greens to the right and reddy-blues to the left; the green list would have greeny-reds to the left and greeny-blues to the right; similarly the blue list would have bluey-greens to the left and bluey-reds to the right.   
4. Join the three lists in the R, G,B order - assuming that B wraps-round to R to form the spectrum.

The Code

 1 
 2 '''
 3 Generate rgb colours in order
 4 '''
 5 
 6 from pprint import pprint as pp
 7 from collections import namedtuple
 8 from itertools import product, groupby
 9 
10 
11 Rgb = namedtuple('Rgb', 'r g b')
12 # red, green, blue, dark, light. 'cardinal' points
13 RGBDL = ( Rgb(255, 0, 0), Rgb(0, 255, 0), Rgb(0, 0, 255), Rgb(0, 0, 0), Rgb(255, 255, 255) )
14 RGB   = ( Rgb(255, 0, 0), Rgb(0, 255, 0), Rgb(0, 0, 255) )
15 
16 def dist(p, q):
17     'dist squared between two points'
18     return sum((pp - qq)**2 for pp,qq in zip(p,q))
19 
20 def cmp(x,y):
21     return 1 if x>y else (0 if x==y else -1)
22 
23 def rgbsort2(rgblist):
24     'reddest to greenest to bluest (... to redest, circularly)'
25     # , , 
26     # cindex = 0 for red; 1 for green; 2 for blue.
27     Closest = namedtuple('Closest', 'distance, cindex, col')
28     # find Closest info for all colours
29     colourdistances = []
30     for col in rgblist:
31         # distance to each cardinal colour in turn
32         rgbdistances = tuple(dist(col, cardinal) for cardinal in RGB)
33         distance = min(rgbdistances)        # closest
34         cindex = rgbdistances.index(distance)
35         colourdistances.append(Closest(distance, cindex, col))
36     # Sort closest to a cardinal point first 
37     colourdistances.sort()
38 
39     # clist accumulates redder, greener, bluer colours separately.
40     # Most intense colour towards the middle of each sublist.
41     clist = [list() for i in range(3)]
42 
43     # Accumulate the sub-spectra for each cardinal colour
44     for _, cindex, col in colourdistances:
45         # add right if closest to next colour cardinal point
46         addright = col[(cindex + 1) % 3] > col[(cindex + 2) % 3]
47         if addright:
48             clist[cindex].append(col)
49         else:
50             clist[cindex].insert(0, col)
51     # Final spectrum
52     return clist[0] + clist[1] + clist[2]

53 
54 def rgb2HTML(colours):
55     f, ht = 4, 50

...

67     return txt
68 
69 
70 if __name__ == '__main__':
71     if 1:
72         # See http://www.lynda.com/resources/webpalette.aspx# for info on
73         # HTML colours without dithering on Windows & Mac.
74         nondithers = (0xFF, 0xCC, 0x99, 0x66, 0x33, 0x00)
75     else:
76         nondithers = tuple(range(0,256,32))
77     rgblist = [Rgb(*rgb) for rgb in product(*(nondithers,)*3)]
78     spectra2 = rgbsort2(rgblist)
79     #assert spectra == spectra2
80     with open('rgbgen.htm', 'w') as f:
81         f.write(rgb2HTML(spectra2))
82 
83 

I have had to miss-out the guts of the rgb2HTML function as Blogger has problems with rendering the HTML tags 

The result

I get some sort of progression, but could do better:

#FFFFFF
(255,255,255)
#FFCCFF
(255,204,255)
#FF99FF
(255,153,255)
#CCCCCC
(204,204,204)
#FFCCCC
(255,204,204)
#FF66FF
(255,102,255)
#FF33FF
(255,51,255)
#CC99CC
(204,153,204)
#FF99CC
(255,153,204)
#FF00FF
(255,0,255)
#000000
(0,0,0)
#999999
(153,153,153)
#CC66CC
(204,102,204)
#FF66CC
(255,102,204)
#CC9999
(204,153,153)
#FF9999
(255,153,153)
#CC33CC
(204,51,204)
#333333
(51,51,51)
#FF33CC
(255,51,204)
#CC00CC
(204,0,204)
#996699
(153,102,153)
#666666
(102,102,102)
#330033
(51,0,51)
#FF00CC
(255,0,204)
#330000
(51,0,0)
#CC6699
(204,102,153)
#993399
(153,51,153)
#663366
(102,51,102)
#FF6699
(255,102,153)
#990099
(153,0,153)
#660066
(102,0,102)
#996666
(153,102,102)
#CC3399
(204,51,153)
#663333
(102,51,51)
#FF3399
(255,51,153)
#CC0099
(204,0,153)
#660033
(102,0,51)
#FF0099
(255,0,153)
#CC6666
(204,102,102)
#993366
(153,51,102)
#660000
(102,0,0)
#FF6666
(255,102,102)
#990066
(153,0,102)
#CC3366
(204,51,102)
#993333
(153,51,51)
#FF3366
(255,51,102)
#CC0066
(204,0,102)
#990033
(153,0,51)
#FF0066
(255,0,102)
#990000
(153,0,0)
#CC3333
(204,51,51)
#FF3333
(255,51,51)
#CC0033
(204,0,51)
#FF0033
(255,0,51)
#CC0000
(204,0,0)
#FF0000
(255,0,0)
#FF3300
(255,51,0)
#CC3300
(204,51,0)
#FF6600
(255,102,0)
#993300
(153,51,0)
#CC6600
(204,102,0)
#FF6633
(255,102,51)
#CC6633
(204,102,51)
#996600
(153,102,0)
#996633
(153,102,51)
#FF9900
(255,153,0)
#663300
(102,51,0)
#CC9900
(204,153,0)
#FF9933
(255,153,51)
#CC9933
(204,153,51)
#666600
(102,102,0)
#999900
(153,153,0)
#FF9966
(255,153,102)
#666633
(102,102,51)
#999933
(153,153,51)
#CC9966
(204,153,102)
#FFCC00
(255,204,0)
#333300
(51,51,0)
#999966
(153,153,102)
#CCCC00
(204,204,0)
#FFCC33
(255,204,51)
#CCCC33
(204,204,51)
#FFCC66
(255,204,102)
#CCCC66
(204,204,102)
#FFCC99
(255,204,153)
#FFFF00
(255,255,0)
#CCCC99
(204,204,153)
#FFFF33
(255,255,51)
#FFFF66
(255,255,102)
#FFFF99
(255,255,153)
#FFFFCC
(255,255,204)
#CCFFCC
(204,255,204)
#CCFF99
(204,255,153)
#CCFF66
(204,255,102)
#99CC99
(153,204,153)
#99FF99
(153,255,153)
#CCFF33
(204,255,51)
#CCFF00
(204,255,0)
#003300
(0,51,0)
#99CC66
(153,204,102)
#99FF66
(153,255,102)
#669966
(102,153,102)
#99CC33
(153,204,51)
#336633
(51,102,51)
#99FF33
(153,255,51)
#99CC00
(153,204,0)
#336600
(51,102,0)
#99FF00
(153,255,0)
#66CC66
(102,204,102)
#669933
(102,153,51)
#006600
(0,102,0)
#66FF66
(102,255,102)
#669900
(102,153,0)
#66CC33
(102,204,51)
#339933
(51,153,51)
#66FF33
(102,255,51)
#66CC00
(102,204,0)
#339900
(51,153,0)
#66FF00
(102,255,0)
#009900
(0,153,0)
#33CC33
(51,204,51)
#33FF33
(51,255,51)
#33CC00
(51,204,0)
#33FF00
(51,255,0)
#00CC00
(0,204,0)
#00FF00
(0,255,0)
#00FF33
(0,255,51)
#00CC33
(0,204,51)
#00FF66
(0,255,102)
#009933
(0,153,51)
#00CC66
(0,204,102)
#33FF66
(51,255,102)
#33CC66
(51,204,102)
#009966
(0,153,102)
#00FF99
(0,255,153)
#339966
(51,153,102)
#006633
(0,102,51)
#00CC99
(0,204,153)
#33FF99
(51,255,153)
#33CC99
(51,204,153)
#006666
(0,102,102)
#009999
(0,153,153)
#66FF99
(102,255,153)
#336666
(51,102,102)
#339999
(51,153,153)
#66CC99
(102,204,153)
#00FFCC
(0,255,204)
#003333
(0,51,51)
#00CCCC
(0,204,204)
#33FFCC
(51,255,204)
#669999
(102,153,153)
#33CCCC
(51,204,204)
#66FFCC
(102,255,204)
#66CCCC
(102,204,204)
#00FFFF
(0,255,255)
#99FFCC
(153,255,204)
#33FFFF
(51,255,255)
#99CCCC
(153,204,204)
#66FFFF
(102,255,255)
#99FFFF
(153,255,255)
#CCFFFF
(204,255,255)
#CCCCFF
(204,204,255)
#99CCFF
(153,204,255)
#66CCFF
(102,204,255)
#9999CC
(153,153,204)
#9999FF
(153,153,255)
#33CCFF
(51,204,255)
#00CCFF
(0,204,255)
#000033
(0,0,51)
#6699CC
(102,153,204)
#6699FF
(102,153,255)
#666699
(102,102,153)
#3399CC
(51,153,204)
#333366
(51,51,102)
#3399FF
(51,153,255)
#0099CC
(0,153,204)
#003366
(0,51,102)
#6666CC
(102,102,204)
#336699
(51,102,153)
#0099FF
(0,153,255)
#000066
(0,0,102)
#6666FF
(102,102,255)
#006699
(0,102,153)
#3366CC
(51,102,204)
#333399
(51,51,153)
#3366FF
(51,102,255)
#0066CC
(0,102,204)
#003399
(0,51,153)
#0066FF
(0,102,255)
#000099
(0,0,153)
#3333CC
(51,51,204)
#3333FF
(51,51,255)
#0033CC
(0,51,204)
#0033FF
(0,51,255)
#0000CC
(0,0,204)
#0000FF
(0,0,255)
#3300FF
(51,0,255)
#3300CC
(51,0,204)
#6600FF
(102,0,255)
#330099
(51,0,153)
#6600CC
(102,0,204)
#6633FF
(102,51,255)
#6633CC
(102,51,204)
#660099
(102,0,153)
#663399
(102,51,153)
#9900FF
(153,0,255)
#330066
(51,0,102)
#9900CC
(153,0,204)
#9933FF
(153,51,255)
#9933CC
(153,51,204)
#9966FF
(153,102,255)
#9966CC
(153,102,204)
#CC00FF
(204,0,255)
#CC33FF
(204,51,255)
#CC66FF
(204,102,255)
#CC99FF
(204,153,255)