You can rewrite some of that in Cython and it will probably speed up by a factor of 10 to 200 (hard to say in advance). Someone more qualified to help with Cython will hopefully write in next - btw for other readers, I think there should be some more advanced/varied examples in the manual for cython.
Cheers, Marshall Hampton On Dec 22, 12:52 pm, pgdoyle <[EMAIL PROTECTED]> wrote: > 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 > versionhttp://math.dartmouth.edu/~prob/prob/prob.pdf > or the essentially identical official > versionhttp://www.dartmouth.edu/~chance/teaching_aids/books_articles/probabi... > 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/ -~----------~----~----~----~------~----~------~--~---
