Dear Peter,
I must confess, with apologies, that I am an expert in neither
statistics nor SAGE, and that I have not fully understood your code
(or, for that matter, "New Age Solitaire"). But I suspect there is
an error in your computation showing that "seven shuffles are
insufficient." The claim, as I recall it, is that seven "good but
imperfect" shuffles "adequately" randomize a 52-card deck.
Experimental verification should depend on the definition of "good
but imperfect," presumably based on observations of how people
shuffle. I think that, for each shuffle, you must choose a
probability distribution for the number of cards in the left- and
right-hand packets into which the deck is cut, and then choose a
conditional probability distribution, when there are m cards
remaining in the left-hand packet and n cards remaining in the
right-hand packet, of how many cards the shuffler releases from the
left-hand packet before releasing a card from the right-hand packet
(and vice versa).
It is not clear to me that your program does anything like this.
If, with this model, seven "good but imperfect" shuffles yield no
significant advantage or disadvantage in the game, one must decide
whether "New Age Solitaire" is really an acceptable statistical test
for randomness of the deck.
Concerning your question: I can't believe seven shuffles then a cut
will have any significant randomization advantage over seven
shuffles.
A closing question for the SAGE experts: are there special
considerations for generating a pseudo-random number (or
permutation) in SAGE/Python inside a big loop such as required for
this problem? My understanding is that in C, for example, for
security in cryptographic applications one should use "arc4random"
rather than "rand" or "srand." Are there corresponding commands for
SAGE? I would think that within a big loop one would not want to
use the time to seed a pseudo-random-number generator on each
iteration; on the other hand, might one get better randomness
properties by seeding the generator at the outset and using
something like Mersenne Twister each iteration, rather than using
entropy stored on the computer at each iteration? If so, how does
one reconcile this with the security considerations mentioned
earlier? How about speed considerations?
Sincerely,
Greg
---------------------------------------------------------
| Greg Marks |
| Department of Mathematics and Computer Science |
| St. Louis University |
| St. Louis, MO 63103-2007 |
| U.S.A. |
| |
| Phone: (314)977-7206 |
| Fax: (314)977-1452 |
| Web: http://euler.slu.edu/Dept/Faculty/marks/marks.html |
---------------------------------------------------------
On Dec 22 2007, 6:23 pm, Marshall Hampton <[EMAIL PROTECTED]> wrote:
> 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):
>
> ...
>
> read more ยป
--~--~---------~--~----~------------~-------~--~----~
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/
-~----------~----~----~----~------~----~------~--~---
