I'm looking for advice about how to speed up the attached Sage program. I've been commissioned to write an article for a popular math journal debunking the notion that `seven shuffles suffice'. This article will feature a computation done in Sage of the exact probability of winning `New Age Solitaire' after seven shuffles of a brand new deck. New Age Solitaire is a game between two players named Yin and Yang that would be fair using a perfectly shuffled deck. You can read about it in Grinstead and Snell's `Introduction to Probability', either the FDL version http://math.dartmouth.edu/~prob/prob/prob.pdf or the essentially identical official version http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/book.html With the aid of Sage, we find that the exact probability (numerator, denominator) for Yin to win after 7 shuffles of a brand new deck is (28598403744641703248428714722775364576331627268813061462556779320938157681698209904159169812750526279755985216, 37576681324381331646231689548629392438010920782533117931316655544515344401833735095419183974156299248510959616) The numerical value is 0.761067841456, so just over 3/4. So 7 shuffles are far from sufficient.
But, what if you do a cut after those seven shuffles? New Age Solitaire was designed to be very insensitive to a final cut, and simulation (and actual practice) shows that this final cut doesn't much lower the probability that Yin will win. It has been suggested that I could leave this as a challenge to readers, and that might be the right thing to do. However, I'm impatient to know the answer. The problem is that the only way I know how to compute this is to add a couple more layers of looping, which will amount to running the attached program at least 2500 times - and it already takes hours to run. Specific questions I have are: 1. I would like to speed this computation up by a factor of something like 100. I have a notion that this could be done by recoding in C using the GMP routines. Does this seem plausible? If so, can I hope to achieve something similar within Sage? 2. This program as it stands is pure Python, except for the `time' command. When run from the command line, it took 247 minutes of user time to run under Python, and 323 minutes under Sage. I had thought that pure Python might be slower because unlike Sage it does not use the GMP package. (But, if that's true, why would Python not switch to using GMP??) But maybe these integers aren't long enough for the difference to show up? 3. The program as it stands was tweaked by running in Python and using the cProfile package. What would the natural alternative to this be working within Sage? 4. Should I be thinking about recoding with Cython? If so, what issues do I have to bear in mind? I should say that it is conceivable that this computation could be sped up a hundred-fold by improving the algorithm. But since the program includes little in the way of comments, I don't imagine it will be very clear what that algorithm is. And in any case, I would like to understand what kind of speed-ups are possible within Sage for a computation like this. Please note that the attached program doesn't work when pasted directly into a Sage notebook cell, because the notebook has a problem with things like this: sage: if True: ... print 2+\ ... 2 4 You can find the source (without the output) here: http://www.math.dartmouth.edu/~doyle/docs/newage/newage.sage Cheers, Peter Doyle sage: class Memoize: ... """Memoize(fn) - an instance which acts like fn but memoizes its arguments ... Will only work on functions with non-mutable arguments ... """ ... def __init__(self, fn): ... self.fn = fn ... self.memo = {} ... def __call__(self, *args): ... if args not in self.memo: ... self.memo[args] = self.fn(*args) ... return self.memo[args] ... sage: """ sage: Compute interleavings sage: """ sage: from operator import add, mul sage: def intersperse(a,b): ... if a<0 or b<0: ... return 0 ... elif a==0 or b==0: ... return 1 ... else: ... a,b=max(a,b),min(a,b) ... return reduce(mul,xrange(a+1,a+b+1))//reduce(mul,xrange(1,b +1)) ... sage: intersperse=Memoize(intersperse) sage: def fix_a(a,b): return (a-1,b,a,b) sage: def fix_b(a,b): return (a,b-1,a,b) sage: def fix_top(m,n): return (m,n,m,n) sage: def linear_extensions(constraints): ... constraints.sort() ... a0,b0=0,0 ... prd=1 ... for (a1,b1,a2,b2) in constraints: ... prd*=intersperse(a1-a0,b1-b0) ... a0,b0=a2,b2 ... return prd ... sage: def mixAA((a0,b0),(a1,b1),(m,n)): ... return linear_extensions([fix_a(a0,b0),fix_a(a1,b1),fix_top(m,n)]) ... sage: mixAA=Memoize(mixAA) sage: def mixAAfast((a0,b0),(a1,b1),(m,n)): ... if a0>a1: ... return mixAAfast((a1,b1),(a0,b0),(m,n)) ... else: ... return intersperse(a0-1,b0)*\ ... intersperse(a1-1-a0,b1-b0)*\ ... intersperse(m-a1,n-b1) ... sage: mixAAfast=Memoize(mixAAfast) sage: def meldAA(weight,b0,b1,(m,n)): ... return sum(weight(a0,a1)*mixAA((a0,b0),(a1,b1),(m,n)) ... for a0 in xrange(1,m+1) ... for a1 in xrange(1,m+1) ... if cmp(a0,a1)*cmp(b0,b1)>=0) ... sage: def meldAAfast(weight,b0,b1,(m,n)): ... sm=0 ... if b0<=b1: ... sm+=sum( ... intersperse(a0-1,b0)* ... sum( ... intersperse(a1-1-a0,b1-b0)* ... intersperse(m-a1,n-b1)* ... weight(a0,a1) ... for a1 in xrange(a0+1,m+1)) ... for a0 in xrange(1,m+1)) ... if b1<=b0: ... sm+=sum( ... intersperse(a1-1,b1)* ... sum( ... intersperse(a0-1-a1,b0-b1)* ... intersperse(m-a0,n-b0)* ... weight(a0,a1) ... for a0 in xrange(a1+1,m+1)) ... for a1 in xrange(1,n+1)) ... return sm ... sage: def mixAB((a0,b0),(a1,b1),(m,n)): ... return linear_extensions([fix_a(a0,b0),fix_b(a1,b1),fix_top(m,n)]) ... sage: mixAB=Memoize(mixAB) sage: def mixABB((a0,b0),(a1,b1),(a2,b2),(m,n)): ... return linear_extensions([fix_a(a0,b0),fix_b(a1,b1),fix_b(a2,b2),fix_top(m,n)]) ... sage: def mixAABB((a0,b0),(a1,b1),(a2,b2),(a3,b3),(m,n)): ... return linear_extensions([fix_a(a0,b0),fix_a(a1,b1),fix_b(a2,b2),fix_b(a3,b3),fix_top(m,n)]) ... sage: def nrange(n,Y): return xrange(max(1,Y-n),min(n,Y)+1) sage: """ sage: Arrange cards in a single suit, noting position of first and last cards. sage: """ sage: def chi(test): ... if test: ... return 1 ... else: ... return 0 ... sage: def stack1xy(n,k,x,y): ... if (n==k==x==y==1): ... return 1 ... elif x==y or not(1<=k<=n and 1<=x<=n and 1<=y<=n): ... return 0 ... else: ... return sum(stack1xy(n-1,k-chi(y<=y0),x-chi(y<x),y0) for y0 in xrange(1,n)) ... sage: stack1xy=Memoize(stack1xy) sage: def stack1Xy(n,k,(x0,x1),y): ... return sum(stack1xy(n,k,x,y) for x in xrange(x0,x1+1)) ... sage: stack1Xy=Memoize(stack1Xy) sage: def stack1y(n,k,y): ... return stack1Xy(n,k,(1,n),y) ... sage: """ sage: Combine 2 separate suits sage: """ sage: def stack2xyxy(n,k1,k2,X1,Y1,X2,Y2): ... sm=0 ... for x1 in nrange(n,X1): ... for y1 in nrange(n,Y1): ... for x2 in nrange(n,X2): ... for y2 in nrange(n,Y2): ... sm+=stack1xy(n,k1,x1,y1)*stack1xy(n,k2,x2,y2)*\ ... mixAABB((x1,X1-x1),(y1,Y1-y1),(X2-x2,x2),(Y2-y2,y2), (n,n)) ... return sm ... sage: def stack2yxy(n,k1,k2,Y1,X2,Y2): ... sm=0 ... for y1 in nrange(n,Y1): ... for x2 in nrange(n,X2): ... for y2 in nrange(n,Y2): ... sm+=stack1y(n,k1,y1)*stack1xy(n,k2,x2,y2)*\ ... mixABB((y1,Y1-y1),(X2-x2,x2),(Y2-y2,y2),(n,n)) ... return sm ... sage: def stack2yyd((n1,n2),(k1,k2),Y1,Y2,d12): ... sm=0 ... for y1 in nrange(n1,Y1): ... for y2 in nrange(n2,Y2): ... if d12==0: ... x2range=(Y1-y1+1,n2) ... elif d12==1: ... x2range=(1,Y1-y1) ... sm+=(stack1y(n1,k1,y1)*stack1Xy(n2,k2,x2range,y2)* ... mixAB((y1,Y1-y1),(Y2-y2,y2),(n1,n2))) ... return sm ... sage: stack2yyd=Memoize(stack2yyd) sage: def stack2Yyd((n1,n2),(k1,k2),(Y10,Y11),Y2,d12): ... return sum(stack2yyd((n1,n2),(k1,k2),Y1,Y2,d12) for Y1 in xrange(Y10,Y11+1)) ... sage: """ sage: We analyze newage solitaire using a polynomial-time algorithm. This was sage: a bitch to work out. The tough part was realizing that the suits should sage: be ordered hearts;clubs;spades;diamonds. Well, not the only hard part! sage: It will be a struggle to change this to allow a cut after shuffling. sage: """ sage: def chi(test): ... if test: ... return 1 ... else: ... return 0 ... sage: def yinmix((n1,n2,n3,n4),(k1,k2),y1p,y2p,d12,cmpy1y2): ... return sum(stack2yyd((n1,n2),(k1,k2),y1,y2,d12) ... *mixAA((y1,y1p),(y2,y2p),(n1+n2,n3+n4)) ... for y1 in xrange(1,n1+n2+1) ... for y2 in xrange(1,n1+n2+1) ... if cmp(y1,y2)==cmpy1y2) ... sage: yinmix=Memoize(yinmix) sage: def weakcompat(z1,z3,y1p,y3): ... return not (z3>z1 and y3<=y1p or z3<z1 and y3>y1p) ... sage: def yang((n3,n4),(k3,k4),(z1,z2,z3,z4),y1p,y2p,d34): ... sm=0 ... for y3 in xrange(1,n3+n4+1): ... if weakcompat(z1,z3,y1p,y3) and weakcompat(z2,z3,y2p,y3): ... for y4 in xrange(1,n3+n4+1): ... if cmp(y3,y4)==cmp(z3,z4): ... if weakcompat(z1,z4,y1p,y4) and weakcompat(z2,z4,y2p,y4): ... sm+=stack2yyd((n3,n4),(k3,k4),y3,y4,d34) ... return sm ... sage: yang=Memoize(yang) sage: def stack4((n1,n2,n3,n4),(k1,k2,k3,k4),(z1,z2,z3,z4),d12,d34): ... sm=0 ... for y1p in xrange(0,n3+n4+1): ... for y2p in xrange(0,n3+n4+1): ... yinfactor=yinmix((n1,n2,n3,n4), (k1,k2),y1p,y2p,d12,cmp(z1,z2)) ... yangfactor=yang((n3,n4),(k3,k4),(z1,z2,z3,z4),y1p,y2p,d34) ... summand=yinfactor*yangfactor ... sm+=summand ... return sm ... ... sage: def newagerising(n): ... def beats(i,j): ... return k[i-1]<k[j-1] or k[i-1]==k[j-1] and z[i-1]<z[j-1] ... sm=0 ... yin=[0 for _ in xrange(4*n)] ... yang=[0 for _ in xrange(4*n)] ... perm4= [(0, 1, 2, 3), (1, 0, 2, 3), (1, 2, 0, 3), (1, 2, 3, 0), ... (2, 1, 3, 0), (2, 1, 0, 3), (2, 0, 1, 3), (0, 2, 1, 3), ... (0, 2, 3, 1), (2, 0, 3, 1), (2, 3, 0, 1), (2, 3, 1, 0), ... (3, 2, 1, 0), (3, 2, 0, 1), (3, 0, 2, 1), (0, 3, 2, 1), ... (0, 3, 1, 2), (3, 0, 1, 2), (3, 1, 0, 2), (3, 1, 2, 0), ... (1, 3, 2, 0), (1, 3, 0, 2), (1, 0, 3, 2), (0, 1, 3, 2)] ... for z in perm4: ... for k1 in xrange(1,n+1): ... for k2 in xrange(1,n+1): ... for k3 in xrange(1,n+1): ... for k4 in xrange(1,n+1): ... k=(k1,k2,k3,k4) ... for d12 in xrange(2): ... for d34 in xrange(2): ... num=stack4((n,n,n,n),(k1,k2,k3,k4),z,d12,d34) ... sm+=num ... rising = (k1 + d12-1 + k2 + chi(z[2-1]>z[4-1])-1 + ... (n+1-k4) + (1-d34)-1 + (n+1-k3)) ... if (beats(1,3) and beats(2,3)) or (beats(1,4) and beats(2,4)): ... yin[rising-1]+=num ... else: ... yang[rising-1]+=num ... return yin,yang ... sage: newagerising=Memoize(newagerising) sage: assert newagerising(2)==( sage: [1, 247, 4058, 10852, 4767, 235, 0, 0], sage: [0, 0, 235, 4767, 10852, 4058, 247, 1] sage: ) sage: """ sage: Precomputed value of newagerising(13). sage: """ sage: def newagerisingprecomp(): # We call this when we're in a hurry! ... newagerising.memo[(13,)]=newagerising13 ... sage: newagerising13=( sage: [1, sage: 4503599627370443, sage: 6461081650535893048297331, sage: 20282067166317747370548924397305, sage: 2219371090444690280167825067011163404, sage: 28980470297130316118851707371113927682308, sage: 86585645711842456879259291396042785694734772, sage: 86713283824808603371563209439361605206738793756, sage: 37025109959688438829553523840364680262742546084490, sage: 7911300235037463075597685089436522698036110779652974, sage: 945840628557918451844218451393465611283022070265930318, sage: 68592119011285455655624013113555233530495826611028105002, sage: 3204605114791094679078453140281792404372059654677564605036, sage: 101000927132657645557134474918097219636280102533165819882650, sage: 2226078789301170911354255920866370371383973348465648304170860, sage: 35302220045338220225580622115664740445656875982546609370577770, sage: 412144632644620632385776599476606526541305878691838924369169955, sage: 3608546497413591803007825874543114725410240443825876029297769965, sage: 24054766216581786383422508285608975096737966214618638425793792925, sage: 123592913881250007522349524203788602163050983663352993622075332819, sage: 494322137147841863969033929598152124865942190994877141968645537910, sage: 1550031910420000204327475350201638219761937474122577308611269752748, sage: 3813750647234720923739675547523812833991654884390508732045736588918, sage: 7249685929003799823021465015144309625316439904669010543339100489100, sage: 10157138379957908026364259231467096289826421195548843978279274466084, sage: 9527148524576740523122808540583216291094177476227815070827277915360, sage: 5319367426533290957031609635466949819030410602742172288617328099396, sage: 1708761063962935980537633895317266586170119106407739574178157449160, sage: 321392175387682201555174973511266244770986158887278683415949229420, sage: 33867769676448620279983034127340470376102593103310027644564227378, sage: 1867166053182294763880505284437691405616536599181852310825314844, sage: 47942965943052996618642323398791765740587983618546831080507618, sage: 506932516495199166765198140279691356252737553454998858074700, sage: 1815528418847090136286185809255612475043185532772499623160, sage: 1596943372409915924781915953648558110500694836711171936, sage: 155465066578682557696559182721119863427481794460452, sage: 3935654927798325081534219322126166356196012978, sage: 866959791589056662810630736708606751304, sage: 2404660627670201923936813450, sage: 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], sage: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, sage: 2404660627670201923936813450, sage: 866959791589056662810630736708606751304, sage: 3935654927798325081534219322126166356196012978, sage: 155465066578682557696559182721119863427481794460452, sage: 1596943372409915924781915953648558110500694836711171936, sage: 1815528418847090136286185809255612475043185532772499623160, sage: 506932516495199166765198140279691356252737553454998858074700, sage: 47942965943052996618642323398791765740587983618546831080507618, sage: 1867166053182294763880505284437691405616536599181852310825314844, sage: 33867769676448620279983034127340470376102593103310027644564227378, sage: 321392175387682201555174973511266244770986158887278683415949229420, sage: 1708761063962935980537633895317266586170119106407739574178157449160, sage: 5319367426533290957031609635466949819030410602742172288617328099396, sage: 9527148524576740523122808540583216291094177476227815070827277915360, sage: 10157138379957908026364259231467096289826421195548843978279274466084, sage: 7249685929003799823021465015144309625316439904669010543339100489100, sage: 3813750647234720923739675547523812833991654884390508732045736588918, sage: 1550031910420000204327475350201638219761937474122577308611269752748, sage: 494322137147841863969033929598152124865942190994877141968645537910, sage: 123592913881250007522349524203788602163050983663352993622075332819, sage: 24054766216581786383422508285608975096737966214618638425793792925, sage: 3608546497413591803007825874543114725410240443825876029297769965, sage: 412144632644620632385776599476606526541305878691838924369169955, sage: 35302220045338220225580622115664740445656875982546609370577770, sage: 2226078789301170911354255920866370371383973348465648304170860, sage: 101000927132657645557134474918097219636280102533165819882650, sage: 3204605114791094679078453140281792404372059654677564605036, sage: 68592119011285455655624013113555233530495826611028105002, sage: 945840628557918451844218451393465611283022070265930318, sage: 7911300235037463075597685089436522698036110779652974, sage: 37025109959688438829553523840364680262742546084490, sage: 86713283824808603371563209439361605206738793756, sage: 86585645711842456879259291396042785694734772, sage: 28980470297130316118851707371113927682308, sage: 2219371090444690280167825067011163404, sage: 20282067166317747370548924397305, sage: 6461081650535893048297331, sage: 4503599627370443, sage: 1] sage: ) sage: def yinrising(n): ... return newagerising(n)[0] ... sage: def yangrising(n): ... return newagerising(n)[1] ... sage: """ sage: Now we know how winning configureations for Yin and Yang sage: are distributed according to rising number. sage: That was the hard part. Now we can figure the probability of sage: winning after a designated number of shuffles because sage: the probability of an arrangement depends just on the rising sage: number. sage: """ sage: def dot(alist,blist): ... return sum(a*b for (a,b) in zip(alist,blist)) ... sage: def promote(n,k,a): ... return intersperse(a-k,n) ... sage: def shuffleprob(winsbyrising,hands): ... n=len(winsbyrising) ... num=dot(winsbyrising,[promote(n,k,hands) for k in xrange(1,n +1)]) ... denom=hands**n ... return (num,denom) ... sage: assert shuffleprob(newagerising13[0],2**7)==( sage: 28598403744641703248428714722775364576331627268813061462556779320938157681698209904159169812750526279755985216, sage: 37576681324381331646231689548629392438010920782533117931316655544515344401833735095419183974156299248510959616) sage: def makereal((num,denom)): ... window=10**20 ... return (((window*num)//denom)*1.0)/window ... sage: def test(risinglist,kmax): ... for s in range(kmax+1): ... print s, makereal(shuffleprob(risinglist,2**s)) ... sage: #newagerisingprecomp() sage: ranks=13 sage: time print newagerising(ranks) sage: print sage: test(yinrising(ranks),20) sage: print sage: print shuffleprob(yinrising(ranks),2**7) ([1, 4503599627370443, 6461081650535893048297331, 20282067166317747370548924397305, 2219371090444690280167825067011163404, 28980470297130316118851707371113927682308, 86585645711842456879259291396042785694734772, 86713283824808603371563209439361605206738793756, 37025109959688438829553523840364680262742546084490, 7911300235037463075597685089436522698036110779652974, 945840628557918451844218451393465611283022070265930318, 68592119011285455655624013113555233530495826611028105002, 3204605114791094679078453140281792404372059654677564605036, 101000927132657645557134474918097219636280102533165819882650, 2226078789301170911354255920866370371383973348465648304170860, 35302220045338220225580622115664740445656875982546609370577770, 412144632644620632385776599476606526541305878691838924369169955, 3608546497413591803007825874543114725410240443825876029297769965, 24054766216581786383422508285608975096737966214618638425793792925, 123592913881250007522349524203788602163050983663352993622075332819, 494322137147841863969033929598152124865942190994877141968645537910, 1550031910420000204327475350201638219761937474122577308611269752748, 3813750647234720923739675547523812833991654884390508732045736588918, 7249685929003799823021465015144309625316439904669010543339100489100, 10157138379957908026364259231467096289826421195548843978279274466084, 9527148524576740523122808540583216291094177476227815070827277915360, 5319367426533290957031609635466949819030410602742172288617328099396, 1708761063962935980537633895317266586170119106407739574178157449160, 321392175387682201555174973511266244770986158887278683415949229420, 33867769676448620279983034127340470376102593103310027644564227378, 1867166053182294763880505284437691405616536599181852310825314844, 47942965943052996618642323398791765740587983618546831080507618, 506932516495199166765198140279691356252737553454998858074700, 1815528418847090136286185809255612475043185532772499623160, 1596943372409915924781915953648558110500694836711171936, 155465066578682557696559182721119863427481794460452, 3935654927798325081534219322126166356196012978, 866959791589056662810630736708606751304, 2404660627670201923936813450, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2404660627670201923936813450, 866959791589056662810630736708606751304, 3935654927798325081534219322126166356196012978, 155465066578682557696559182721119863427481794460452, 1596943372409915924781915953648558110500694836711171936, 1815528418847090136286185809255612475043185532772499623160, 506932516495199166765198140279691356252737553454998858074700, 47942965943052996618642323398791765740587983618546831080507618, 1867166053182294763880505284437691405616536599181852310825314844, 33867769676448620279983034127340470376102593103310027644564227378, 321392175387682201555174973511266244770986158887278683415949229420, 1708761063962935980537633895317266586170119106407739574178157449160, 5319367426533290957031609635466949819030410602742172288617328099396, 9527148524576740523122808540583216291094177476227815070827277915360, 10157138379957908026364259231467096289826421195548843978279274466084, 7249685929003799823021465015144309625316439904669010543339100489100, 3813750647234720923739675547523812833991654884390508732045736588918, 1550031910420000204327475350201638219761937474122577308611269752748, 494322137147841863969033929598152124865942190994877141968645537910, 123592913881250007522349524203788602163050983663352993622075332819, 24054766216581786383422508285608975096737966214618638425793792925, 3608546497413591803007825874543114725410240443825876029297769965, 412144632644620632385776599476606526541305878691838924369169955, 35302220045338220225580622115664740445656875982546609370577770, 2226078789301170911354255920866370371383973348465648304170860, 101000927132657645557134474918097219636280102533165819882650, 3204605114791094679078453140281792404372059654677564605036, 68592119011285455655624013113555233530495826611028105002, 945840628557918451844218451393465611283022070265930318, 7911300235037463075597685089436522698036110779652974, 37025109959688438829553523840364680262742546084490, 86713283824808603371563209439361605206738793756, 86585645711842456879259291396042785694734772, 28980470297130316118851707371113927682308, 2219371090444690280167825067011163404, 20282067166317747370548924397305, 6461081650535893048297331, 4503599627370443, 1]) Time: CPU 20900.91 s, Wall: 22396.50 s 0 1.00000000000000 1 1.00000000000000 2 1.00000000000000 3 1.00000000000000 4 1.00000000000000 5 0.998659405061071 6 0.924218572288272 7 0.761067841456394 8 0.638326380209901 9 0.570200634084392 10 0.535232079346741 11 0.517632575837904 12 0.508818357004627 13 0.504409437202608 14 0.502204750940860 15 0.501102379512937 16 0.500551190261784 17 0.500275595194057 18 0.500137797604924 19 0.500068898803449 20 0.500034449401848 (28598403744641703248428714722775364576331627268813061462556779320938157681698209904159169812750526279755985216, 37576681324381331646231689548629392438010920782533117931316655544515344401833735095419183974156299248510959616) --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://sage.math.washington.edu/sage/ and http://sage.scipy.org/sage/ -~----------~----~----~----~------~----~------~--~---
