February 2012 - mathfun@mailman.xmission.com

[math-fun] PsuRandom Generator
by Warren Smith 09 Feb '12

09 Feb '12

> From: Phil Carmody <thefatphil(a)yahoo.co.uk> >> uint64 Rand64(){ >> ? uint i; >> ? uint64 v; >> ? v = statevec[0]; >> ? for(i=NumWordsInState-1; i>0; i--){ >> //can be unrolled and all loop iterations run parallel >> ? ? v ^= statevec[i]; >> ? ? statevec[i] += UpdatorF(statevec[i-1]); >> ? } >> ? statevec[0] += RandStrideConst; >> ? return(v); //returns the XOR of all words in state vector >> } > > I like the idea, but that XOR worries me. You've not demonstrated independence of the bits being XORed together and therefore the XOR may cancel out randomness. > Phil --no, that is ok. The fact that the proposed generator has maximum period 2^B where B is the number of bits in the state, means every bit pattern is equally likely and all state-patterns occur exactly once. If we now XOR various state words together, we therefore get something exactly uniformly distributed. Plus the fact that this has apparently passed the BigCrush tests (it slightly failed one of the 160 tests but that was presumably just because that was supposed to happen if you run that many tests) confirms. Furthermore, ANY combining function FF can be used instead of XOR(a,b) if it has the property that FF(a,b) is exactly uniformly distributed if at least one of its arguments is. For example in C using uints, FF(a,b) = (a|1)*b FF(a,b) = a XOR b FF(a,b) = a-b FF(a,b) = CyclicShift(a,b) all would work. Actually just FF(a,b)=b works too, but that did not seem good enough intuitively; if I just returned the top word in the state I suspect I might fail tests of independence of successive differences, because those are really just the outputs of my same sort of generator with one less stage. The main problem with my RandGen idea is it is not very fast on my machine. To solve that, I came up with the idea of outputting more than one machine word per "count" if the outputs are generated in "sufficiently independent" ways from the state bits. I suggest the following: let the words of the state be (in order from least to most sig) A B C D E... P Q R S ... W X Y Z and output E.S.Z, D.R.Y, C.Q.X, B.P.W where "." is some combining operator, such as XOR, and we make about n outputs for a (3n+1)-word state. That way each output word arises from XORing 3 state words together and has at least as much randomness as a (2n+1)-stage generator of the original type. Anyway, here is the current code which outputs only 1 word per "count" so it is not as fast as it could be. Perhaps you will refine the code further, it clearly can be done. #define uint32 uint32_t #define uint64 uint64_t #define uint128 __uint128_t #define RandStrideConst 14923729446516375051ull //must be odd //The below is written for a 64-bit machine. //The 32-bit version with NumWordsInState=4 passes BigCrush U01 tests. //I do not know whether decreasing to 3 words will still pass BigCrush. //If you have 32-bit machine then text-replace 64-->32 and then 128-->64 throughout //the following: //period will be 2^(64 * NumWordsInState), increase for more randomness but slower runtime: #define NumWordsInState 4 inline uint64 UpdatorF(uint64 x){ uint128 y; y = x; if( x <= (1<<(64-7)) ){ x ^= RandStrideConst; } //if(x<=Z), important that Z be even //the above line was done wrong in my original post so I had failed to get odd-parity F. y *= y; x += y ^ (y>>64); return(x); } uint64 statevec[NumWordsInState]; uint64 Rand64slow(){ uint i; uint64 v; v = statevec[0]; for(i=NumWordsInState-1; i>0; i--){ //can be unrolled and all loop iterations run parallel v ^= statevec[i]; statevec[i] += UpdatorF(statevec[i-1]); } statevec[0] += RandStrideConst; return(v); //currently returns the XOR of all words in state vector }

1 0

Re: [math-fun] Turning discrete games into continuous games ??
by Dan Asimov 09 Feb '12

09 Feb '12

<< Are there any examples of an obviously discrete game -- e.g., tic-tac-toe, checkers, chess, go, etc. -- that can be converted into a continuous game ? I've defined a continuous version of Hex. At any time, each player can either draw a curve at speed S, or else move without drawing at speed 2S. In computer simulation it seems to work well. --Dan ________________________________________________________________________________________ It goes without saying that .

1 0

[math-fun] OK, I programmed my proposed psu-random generator....
by Phil Carmody 09 Feb '12

09 Feb '12

From: Warren Smith > uint64 Rand64(){ > uint i; > uint64 v; > v = statevec[0]; > for(i=NumWordsInState-1; i>0; i--){ > //can be unrolled and all loop iterations run parallel > v ^= statevec[i]; > statevec[i] += UpdatorF(statevec[i-1]); > } > statevec[0] += RandStrideConst; > return(v); //returns the XOR of all words in state vector > } I like the idea, but that XOR worries me. You've not demonstrated independence of the bits being XORed together and therefore the XOR may cancel out randomness. Phil

1 0

[math-fun] Puzzle that Fermat posed to Mersenne
by Dan Asimov 08 Feb '12

08 Feb '12

According to what I've read, Fermat posed this problem to Mersenne in 1943. This puzzle is intended for people who haven't seen the solution already: Find a primitive Pythagorean triple a^2 + b^2 = c^2 such that both a+b and c are each square numbers. (Please e-mail your answers straight to me so they don't keep others from solving this.) --Dan ________________________________________________________________________________________ It goes without saying that .

3 3

[math-fun] Goldbach progress: odd# = sum of five primes
by rcs＠xmission.com 08 Feb '12

08 Feb '12

Date: Tue, 31 Jan 2012 19:37:03 GMT (31kb) Title: Every odd number greater than 1 is the sum of at most five primes Authors: Terence Tao Comments: 44 pages, no figures, submitted, Mathematics of Computation \\ We prove that every odd number $N$ greater than 1 is the sum of at most five primes, improving the result of Ramar\'e that every even natural number is the sum of at most six primes. We follow the circle method of Hardy-Littlewood and Vinogradov, together with Vaughan's identity; our innovations, which may be of interest for other Goldbach-type problems, include the use of smoothed exponential sums and optimisation of the Vaughan identity parameters to save or reduce some logarithmic losses, the use of multiple scales following some ideas of Bourgain, and the use of Montgomery's uncertainty principle and the large sieve to improve the $L^2$ estimates on major arcs. Our argument relies on some previous numerical work, namely the verification of Richstein of the even Goldbach conjecture up to $4 \times 10^{14}$, and the verification of van de Lune and (independently) of Wedeniwski of the Riemann hypothesis up to height $3.29 \times 10^9$. \\ ( http://arxiv.org/abs/1201.6656 , 31kb)

2 1

Re: [math-fun] Goldbach progress: odd# = sum of five primes
by Dan Asimov 08 Feb '12

08 Feb '12

Rich wrote: << . . . Authors: Terence Tao . . . >> Can't say I'm surprised that this listing considers Terence Tao to be plural. --Dan ________________________________________________________________________________________ It goes without saying that .

1 0

[math-fun] Puzzle that Fermat posed to Mersenne
by Dan Asimov 08 Feb '12

08 Feb '12

Make that 1643. ________________________________________________________________________________________ It goes without saying that .

1 0

[math-fun] Fwd: Re: Another puzzle
by rcs＠xmission.com 08 Feb '12

08 Feb '12

This is mainly a correction, pointing out a lucuna. In the problem discussed below, a couple of us wondered (privately) > Why must N be odd? After a bit of back & forth, it appears that N cannot be a multiple of 4, but the N = 4K+2 case is open. Note that even solutions of 2^N = 2 (mod N) exist, but are rare. Pinch's list of base 2 pseudoprimes < 10^13[?] included a small side list of even psps. Rich ------ Quoting Jakob Jonsson <jakobj(a)kth.se>: > On 23 Jan 2012, David Wilson asked: > >> What n satisfy 2^n == 3n+2 (mod n^2 - 2n) > > On 30 Jan 2012, Hans Havermann wrote: > >> At first glance, except for n=5, these appear to be the greater of >> twin primes (A006512). However, on closer inspection, some n (563, >> 645, 647, 1105, 1907, 2467, 2701, 2821, 4371, 4373, 4681, 6601, ...) >> are clearly not. > > (Spoilers ahead) > > The clue is to observe that n^2 - 2n = n(n-2). First note that n is > necessarily odd, which means that n and n-2 are relatively prime. By > the Chinese Remainder Theorem, the equation is equivalent to the > following system of equations: > > 2^n == 2 (mod n) > 2^n == 8 (mod n-2) > > Since n-2 is odd, the second equation is equivalent to > > 2^{n-2} == 2 (mod n-2). > > Thus n is a solution if and only if each of n and n-2 is a prime or > a Sarrus/Poulet pseudoprime. > > Jakob > > _______________________________________________ > math-fun mailing list > math-fun(a)mailman.xmission.com > http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun > ----- End forwarded message -----

1 0

[math-fun] Fwd: Somos++
by Bill Gosper 08 Feb '12

08 Feb '12

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Trilinear Recursion Michael Somos somos(a)cis.csuohio.edu <http://gosper.org/webmail/src/compose.php?send_to=somos%40cis.csuohio.edu> It has been known for some time that the recursion A(n) = (A(n-1) + C1) / A(n-2) (of which C1=1 is the Lyness 5-cycle) has a solution in terms of a general Somos-5 sequence a(n) as follows A(n) = (a(n+2) * a(n-3)) / (a(n) * a(n-1) * C2) 1 + A(n) = (a(n+1) * a(n-2) * C3) / (a(n) * a(n-1) * C2) for some constansts C2, C3. There are similar kinds of examples. One of the simplest is the following. If A(n) = (A(n-1) + A(n-2)) / A(n-3), A(1) = A(2) = A(3) = 1, then A(n) = (a(n) * a(n+7)) / (a(n+3) * a(n+4)) where a(n)=(a(n-1)*a(n-6)*a(n-8) + a(n-2)*a(n-4)*a(n-9)) / (a(n=5)*a(n-10)), a(1) = a(2) = ... = a(10) = 1. which is now sequence A205303 in the OEIS at oeis.org. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ While the idea of an integer trilinear recurrence seems startling, you can make one simply by shifting n in Somos 4, Out[7]= a[-4 + n] a[n] == a[-2 + n]^2 + a[-3 + n] a[-1 + n] In[8]:= % /. n -> n - 1 Out[8]= a[-5 + n] a[-1 + n] == a[-3 + n]^2 + a[-4 + n] a[-2 + n] and eliminating a[n-1] In[9]:= Eliminate[{%, %%}, a[n - 1]] Out[9]= a[-5 + n] a[-4 + n] a[n] == a[-3 + n]^3 + a[-4 + n] a[-3 + n] a[-2 + n] + a[-5 + n] a[-2 + n]^2 In[10]:= Solve[%, a[n]] Out[10]= {{a[n] -> ( a[-3 + n]^3 + a[-4 + n] a[-3 + n] a[-2 + n] + a[-5 + n] a[-2 + n]^2)/(a[-5 + n] a[-4 + n])}} I tend to harp on the mysterious polynomial sequence from initializing Somos 4 to 0,1,1,-1,x,... . (See Somos Sequence Near-Addition Formulas and Modular Theta Functions <http://richard.schroeppel.name:8015/somosadd-intro-070305a.pdf> ~p10.) Besides the mystery of why they're (an elliptic divisibility sequence of) polynomials, we have the following. (Strong) EDS => gcd(p[a](x),p[a*b](x)) = p[a](x). But is p[prime](x) always irreducible? Seems so. The converse, p[n](x) irreducible => n prime is false, despite EDS! Because factors of 2 and 3 only guarantee multiples of p[2](x) = -p[3](x) = 1. But p[n](x) irreducible <=> n prime holds for 9<n<146. And then starts failing like crazy, at almost every n=2*prime and n=3*prime, e.g. 201 and 202. Uh oh, NJAS, what makes you think I said anything about 1,4,6,9,146,158,166,178,183,194,201,202,... ? --rwg

1 3

[math-fun] OK, I programmed my proposed psu-random generator....
by Warren Smith 07 Feb '12

07 Feb '12

========= Summary results of BigCrush ========= Version: TestU01 1.2.3 Generator: warren/32/4 Number of statistics: 160 Total CPU time: 06:19:48.03 The following tests gave p-values outside [0.001, 0.9990]: (eps means a value < 1.0e-300): (eps1 means a value < 1.0e-15): Test p-value ---------------------------------------------- 7 CollisionOver, t = 7 3.1e-4 ---------------------------------------------- All other tests were passed (This is with a revised version of Warren's code -- he found a bug -- and for a 32-bit machine, with 4 words of state.) --does it provide a description of this test? I'm not sure whether to believe this is really a test failure; if BigCrush does 160 tests we'd expect 0.32 tests to have p outside [0.001, 0.999]. E.g. they say in the BigCrush paper http://www.iro.umontreal.ca/~lecuyer/myftp/papers/testu01.pdf that the AES cryptosystem failed I think the same CollisionOver randomness test min did, with even worse p=9.2e-5, but then on re-run with different seed it easily passed it, so the failure was not genuine. They also say a "clear failure" is min(p,1-p)<10^(-10). For example "unix random 32" clearly failed 101 different Crush tests... Among generators listed in the BigCrush paper as passing all tests the fastest included ran2, MRGk5-93, Marsa-LFIB4, SuperDuper64, KISS99, Brent-xor4096s, LFib(2^64 ,55, 24,*), and ISAAC. Of those, the four best were probably Marsa-LFIB4, KISS99, Brent-xor4096s, LFib(2^64 ,55, 24,*). Note that all these 4 best generators were invented by Marsaglia and/or Brent. They are definitely the Kings. MarsaLFIB4 is: x[i] = (x[i-55] + x[i-119] + x[i-179] + x[i-256]) using a 256-element state vector and (mod 2^power) addition (any of the + may be replaced by -). Very fast and simple. I just timed it versus my own rand gen MarsaLFIB4: 2.83 sec WarrenRand: 13.23 sec and I'm 4.7X slower. However... it depends on your computer. I presume they weren't lying in the prizewinning paper when they claimed they had the fastest crushproof generators. Only reason MarsaLFIB4 can be slower than way-more-complicated PRNGs like in http://www.thesalmons.org/john/random123/papers/random123sc11.pdf is because of the memory use. With modern processors it pays to do more computing if you touch less memory. Incidentally if you re-run with only 3 state words (or surely if only 2) my randgen will presumably fail more tests, and it would be interesting to know the descriptions of the failed tests so that one then can attempt to strengthen the randgen.

2 1