January 2012 - mathfun@mailman.xmission.com

[math-fun] Diophantine equation a^4+b^4=d*c^2 and the specialness of 17, 68, 82...
by Warren Smith 05 Mar '12

05 Mar '12

Fermat in about 1640 famously proved a^4+-b^4=c^2 has no positive integer solutions. Also a^4-b^4=2*c^2 has no positive integer solutions. These results are shown in T.Nagell: Intro to number theory, Wiley 1951, pp.227-230. However, a^4+b^4=2*c^2 trivially has a ton of solutions a=b, c=a^2. If we exclude the trivial by demanding b>a>0 though, I don't know of any solutions. Let us generalize. Ask about a^4 + b^4 = d*c^2 with a,b,c,d integers and b>a>0. My computer looked for solutions. Search I: a^4+b^4=d*c^2 with 0<a<b<65536 and 0<d<65. Result: a ton of solutions for d=17, but zero solutions for any other d. a=1, b=2, d=17, c=1 a=2, b=4, d=17, c=4 a=2, b=13, d=17, c=41 a=3, b=6, d=17, c=9 a=4, b=8, d=17, c=16 ... Well, that's interesting: 17 is special. Search II: a^4+b^4=d*c^2 with 0<a<b<=1024 and 0<d<=1024. Result: The following d are "special": 17, 68, 82, 97, 113, 153, 193, 257, 272, 274, 328, 337, 388, 425, 433, 452, 514, 577, 593, 612, 626, 641, 673, 706, 738, 772, 833, 873, 881, 914, 1017. [This sequence is not in OEIS.] What is going on here? I do not really know, but an initial speculation is this. When d=17, there is the easy solution a=1, b=2, d=17, c=1. But here is the thing: if there is one solution, then it can usually be used to generate an infinite set of solutions. I got that from Tito Piezas's online book http://sites.google.com/site/tpiezas/018 where he attributes it to A.Desboves (probably about 1880?). So there are "special" d with an infinity of solutions, and there are "usual" d with zero solutions (and conceivably there are a few "crazy" d with a finite nonzero number of solutions because the below Desboves method for generating more solutions happens merely to regenerate old ones, which would be an amazing miracle). How to generate more solutions from one: Piezas says Desboves said (as edited by me): if p^4+q^4 = d*r^2, then x^4+y^4 = d*z^2 where x = p*(4*m^2-3*(m+n)^2), y = q*(4*n^2-3*(m+n)^2), z = r*(4*(m+n)^4-3*(m-n)^4) where m=p^4 and n=q^4. So for example starting from p=1, q=2, d=17, r=1 we get x=863, y=314, z=182209 and sure enough 863^4 + 314^4 = 17*182209^2 is another solution. Questions: Q1. Can you find a parameterized family of d=17 solutions? Q2. Can you understand the set of special d? Q3. Can you prove there are no "crazy" d (or find them all)?

4 3

Re: [math-fun] min(x, y) and sorting without if()
by Phil Carmody 14 Feb '12

14 Feb '12

From: Joerg Arndt <arndt(a)jjj.de> > * Phil Carmody <thefatphil(a)yahoo.co.uk> [Jan 07. 2012 11:41]: > > From: Joerg Arndt <arndt(a)jjj.de> > > > These (branch-avoiding) techniques are "well known", see (fxtbook) > > > section 1.11 "Avoiding branches", pp.15ff where a warning of the form > > > "Your compiler may be smarter than you thought" is given. > > > > Shouldn't the warning be more Nuddsian? > > I do not parse "Nuddsian" here. Pertaining to Scott Nudds, the troll/loon who managed to force us to turn comp.lang.asm.x86 into a moderated newsgroup back in the mid-90s. His name is certainly not forgotten on comp.lang.c either, as he was a notorious cross-poster (he would deliberately try to get the two camps to fight each other). One of his weapons was the invocation of ... > > In C, such expressions are either unnecessary or UB, pure and > > simple. > > Reading UB as "useless b*ll****": "Undefined behaviour", as per the ISO C specification. In reality, it can equate to the same thing, as a sufficiently smart compiler can do exactly what you *don't* want in cases where UB is encountered. Many of the tricks are far from useless, but sometimes need a little bit of massaging before the C standard guarantees they will work. > There are a couple of tricks that prove useful in practice. > E.g., replacing > if ( (a<b) || (c<d) ) { ... } > by > if ( (a<b) | (c<d) ) { ... } > can be an improvement. Anything's possible. The latter has more re-ordering flexibility in some situations, so might fill bubbles in pipelines better. In reality, the programmer may not give a hoot which of the two comparisons expressions are evaluated first, but the short-circuit operator obliges him to express a preference. > One that I use regularly: > Replace > if ( (a<0) || (a>=n) ) { ... } > by > if ( (unsigned)a>=n ) { ... } Yup, that type of trick is a good one and it's safe. I've even gone as far as subtracting a lower bound before doing that, and then comparing >= range rather than >= max. > Yes, there are assumptions made, but I am not aware of > any _existing_ system where these do no hold. They certainly do exist. I perhaps immerse myself in a wider range of architectures than most, both at work and at home, so I am more aware of them. Shifting sign bits around on some DSPs will cause mayhem. Some Crays have trap representations for certain types which mean you can't hold certain intermediate values that would be useful. Worse, arch's with separate address and integer registers lose some of the efficiency of Knuth's dancing links, for example, in particular if their AGU is slow. > > Your compiler doesn't have to show smart behaviour at all - > > it may kill your cat if it wants to. > > Haven't spotted the part about cat slaughter in the standard yet :-) The problem is that it's not forbidden, so when the C standard says it imposes no restriction on what the implementation may do, that's implicitly giving it permission to kill kittens. The original absurd behaviour example was "nasal demons" (demons may fly out of your nose). That was hackneyed enough to have mutated to having the aforementioned Scott Nudds fly out of your nose. All such outcomes are to be avoided. > And quite a few people work extremely hard to make compilers smart. That's one of the things that makes UB so dangerous. They can be so smart that they realise that certain inputs would lead to UB, and therefore they only have to generate the correct results for the inputs that don't lead to UB, no matter what the consequences are in the UB case. GCC being smart regarding being able to not care about what happens when there's UB has introduced security holes into the linux kernel. It can just throw away a null pointer check, for example: http://lwn.net/Articles/341773/ I seem to have deviated both from mathiness and funness here, I apologise. It's just that I live, eat and breathe strict portable C currently, in particular other people's, and I'm hypersensitive to anything which is UB. Phil

2 1

[math-fun] Warren D Smith video lecture: How I broke AES -- if I did it
by Warren Smith 10 Feb '12

10 Feb '12

hmm. Did not realize I'd been videod+internetted. Hope it came out well. The paper that lecture was based on, can be read here: http://rangevoting.org/WarrenSmithPages/homepage/works.html paper#100 and it described what I called an alleged "nonconstructive break" of the AES cryptosystem. Later an alleged "genuine break" of AES was produced: http://research.microsoft.com/en-us/projects/cryptanalysis/aesbc.pdf but both theirs and mine have no real practical impact. My paper however also cited a "(cheating) break" of AES by Daniel J. Bernstein which I think is very serious practically speaking.

2 1

[math-fun] simple pleasures: subset-lex Gray code
by Joerg Arndt 07 Feb '12

07 Feb '12

>From the department of ASCII entertainment: a binary Gray code you have not seen before. Dots for zeros to make the structure more visible. n: Gray code transition 0: [ . . . . . . ] ---- 1: [ 1 . . . . . ] 1..... 2: [ 1 1 . . . . ] .1.... 3: [ 1 1 1 . . . ] ..1... 4: [ 1 1 1 1 . . ] ...1.. 5: [ 1 1 1 1 1 . ] ....1. 6: [ 1 1 1 1 1 1 ] .....1 7: [ 1 1 1 1 . 1 ] ....1. 8: [ 1 1 1 . . 1 ] ...1.. 9: [ 1 1 1 . 1 1 ] ....1. 10: [ 1 1 1 . 1 . ] .....1 11: [ 1 1 . . 1 . ] ..1... 12: [ 1 1 . . 1 1 ] .....1 13: [ 1 1 . . . 1 ] ....1. 14: [ 1 1 . 1 . 1 ] ...1.. 15: [ 1 1 . 1 1 1 ] ....1. 16: [ 1 1 . 1 1 . ] .....1 17: [ 1 1 . 1 . . ] ....1. 18: [ 1 . . 1 . . ] .1.... 19: [ 1 . . 1 1 . ] ....1. 20: [ 1 . . 1 1 1 ] .....1 21: [ 1 . . 1 . 1 ] ....1. 22: [ 1 . . . . 1 ] ...1.. 23: [ 1 . . . 1 1 ] ....1. 24: [ 1 . . . 1 . ] .....1 25: [ 1 . 1 . 1 . ] ..1... 26: [ 1 . 1 . 1 1 ] .....1 27: [ 1 . 1 . . 1 ] ....1. 28: [ 1 . 1 1 . 1 ] ...1.. 29: [ 1 . 1 1 1 1 ] ....1. 30: [ 1 . 1 1 1 . ] .....1 31: [ 1 . 1 1 . . ] ....1. 32: [ 1 . 1 . . . ] ...1.. 33: [ . . 1 . . . ] 1..... 34: [ . . 1 1 . . ] ...1.. 35: [ . . 1 1 1 . ] ....1. 36: [ . . 1 1 1 1 ] .....1 37: [ . . 1 1 . 1 ] ....1. 38: [ . . 1 . . 1 ] ...1.. 39: [ . . 1 . 1 1 ] ....1. 40: [ . . 1 . 1 . ] .....1 41: [ . . . . 1 . ] ..1... 42: [ . . . . 1 1 ] .....1 43: [ . . . . . 1 ] ....1. 44: [ . . . 1 . 1 ] ...1.. 45: [ . . . 1 1 1 ] ....1. 46: [ . . . 1 1 . ] .....1 47: [ . . . 1 . . ] ....1. 48: [ . 1 . 1 . . ] .1.... 49: [ . 1 . 1 1 . ] ....1. 50: [ . 1 . 1 1 1 ] .....1 51: [ . 1 . 1 . 1 ] ....1. 52: [ . 1 . . . 1 ] ...1.. 53: [ . 1 . . 1 1 ] ....1. 54: [ . 1 . . 1 . ] .....1 55: [ . 1 1 . 1 . ] ..1... 56: [ . 1 1 . 1 1 ] .....1 57: [ . 1 1 . . 1 ] ....1. 58: [ . 1 1 1 . 1 ] ...1.. 59: [ . 1 1 1 1 1 ] ....1. 60: [ . 1 1 1 1 . ] .....1 61: [ . 1 1 1 . . ] ....1. 62: [ . 1 1 . . . ] ...1.. 63: [ . 1 . . . . ] ..1... Note successive transitions are either adjacent (one-close) or three-close. Guessing the recursion may be a fun task (there is a twist wrt. to the "well known" recursions). Hint: this Gray code is to (subset-lex)... n: delta-set subset 0: ..... { } 1: 1.... { 0 } 2: 11... { 0, 1 } 3: 111.. { 0, 1, 2 } 4: 1111. { 0, 1, 2, 3 } 5: 11111 { 0, 1, 2, 3, 4 } 6: 111.1 { 0, 1, 2, 4 } 7: 11.1. { 0, 1, 3 } 8: 11.11 { 0, 1, 3, 4 } 9: 11..1 { 0, 1, 4 } 10: 1.1.. { 0, 2 } 11: 1.11. { 0, 2, 3 } 12: 1.111 { 0, 2, 3, 4 } 13: 1.1.1 { 0, 2, 4 } 14: 1..1. { 0, 3 } 15: 1..11 { 0, 3, 4 } 16: 1...1 { 0, 4 } 17: .1... { 1 } 18: .11.. { 1, 2 } 19: .111. { 1, 2, 3 } 20: .1111 { 1, 2, 3, 4 } 21: .11.1 { 1, 2, 4 } 22: .1.1. { 1, 3 } 23: .1.11 { 1, 3, 4 } 24: .1..1 { 1, 4 } 25: ..1.. { 2 } 26: ..11. { 2, 3 } 27: ..111 { 2, 3, 4 } 28: ..1.1 { 2, 4 } 29: ...1. { 3 } 30: ...11 { 3, 4 } 31: ....1 { 4 } .. what the usual (binary reflected) Gray code ist to binary numbers in counting order. Can be generalized to mixed radix numbers (and for even radices >=4 the transitions are either adjacent or two-close) and a loopless algorithm that does _not_ use any sort of focus pointers can be obtained (which is actually very fast). One-close binary Gray codes exist for n<=6 but not for n=7 (and apparently not for any n>=8). Finding a construction for two-close binary Gray code would be _very_ neat. ...and subset-lex gives fast (loopless) methods for generating non-adjacent forms, here is the Fibonacci case: 0: ........ = { } 1: 1....... = { 0 } 2: 1.1..... = { 0, 2 } 3: 1.1.1... = { 0, 2, 4 } 4: 1.1.1.1. = { 0, 2, 4, 6 } 5: 1.1.1..1 = { 0, 2, 4, 7 } 6: 1.1..1.. = { 0, 2, 5 } 7: 1.1..1.1 = { 0, 2, 5, 7 } 8: 1.1...1. = { 0, 2, 6 } 9: 1.1....1 = { 0, 2, 7 } 10: 1..1.... = { 0, 3 } 11: 1..1.1.. = { 0, 3, 5 } 12: 1..1.1.1 = { 0, 3, 5, 7 } 13: 1..1..1. = { 0, 3, 6 } 14: 1..1...1 = { 0, 3, 7 } 15: 1...1... = { 0, 4 } 16: 1...1.1. = { 0, 4, 6 } 17: 1...1..1 = { 0, 4, 7 } 18: 1....1.. = { 0, 5 } 19: 1....1.1 = { 0, 5, 7 } 20: 1.....1. = { 0, 6 } 21: 1......1 = { 0, 7 } 22: .1...... = { 1 } 23: .1.1.... = { 1, 3 } 24: .1.1.1.. = { 1, 3, 5 } 25: .1.1.1.1 = { 1, 3, 5, 7 } 26: .1.1..1. = { 1, 3, 6 } 27: .1.1...1 = { 1, 3, 7 } 28: .1..1... = { 1, 4 } 29: .1..1.1. = { 1, 4, 6 } 30: .1..1..1 = { 1, 4, 7 } 31: .1...1.. = { 1, 5 } 32: .1...1.1 = { 1, 5, 7 } 33: .1....1. = { 1, 6 } 34: .1.....1 = { 1, 7 } 35: ..1..... = { 2 } 36: ..1.1... = { 2, 4 } 37: ..1.1.1. = { 2, 4, 6 } 38: ..1.1..1 = { 2, 4, 7 } 39: ..1..1.. = { 2, 5 } 40: ..1..1.1 = { 2, 5, 7 } 41: ..1...1. = { 2, 6 } 42: ..1....1 = { 2, 7 } 43: ...1.... = { 3 } 44: ...1.1.. = { 3, 5 } 45: ...1.1.1 = { 3, 5, 7 } 46: ...1..1. = { 3, 6 } 47: ...1...1 = { 3, 7 } 48: ....1... = { 4 } 49: ....1.1. = { 4, 6 } 50: ....1..1 = { 4, 7 } 51: .....1.. = { 5 } 52: .....1.1 = { 5, 7 } 53: ......1. = { 6 } 54: .......1 = { 7 } Generating the list in reverse order can be done via... ulong prev_lexrev_fib(ulong x) // Return previous Fibonacci word in subset-lex order. // Start with zero to generate all words of length n. { ulong x0 = x & -x; // lowest bit if ( x & (x0<<2) ) // easy case: next higher bit is set { x ^= x0; // clear lowest bit return x; } else { x += x0; // move lowest bit to the left x |= (x0!=1); // and set rightmost bit unless blocked by next bit return x; } } ...which, ladies and gentlemen, ist werry beautiful! [He bows and leaves the stage. Deafening silence in the audience.]

1 1

[math-fun] Wilson's theorem on primes (no, not that one, the one about decimals)
by Warren Smith 05 Feb '12

05 Feb '12

> From: David Wilson <davidwwilson(a)comcast.net> > Let p >= 10^d, and let 0 <= B < 10^d be any d-digit block of digits. > By the definition of the ceiling function > [1] pB/10^d <= ceil(pB/10^d) < pB/10^d + 1. > Dividing by p gives > [2] B/10^d <= ceil(pB/10^d)/p < B/10^d + 1/p > But p >= 10^d, so 1/p <= 1/10^d, and > [3] B/10^d + 1/p <= B/10^d + 1/10^d = (B+1)/10^d > [2] and [3] together give > [4] B/10^d <= ceil(pB/10^d)/p < (B+1)/10^d > Letting n = ceil(pB/10^d), we have integer n with > [5] B/10^d <= n/p < (B+1)/10^d > [5] says that the decimal fraction n/p starts with block B. > So if p >= 10^d, there is a fraction n/p (0 <= n < p) starting with any > given d-digit block B, namely n/p = ceil(pB/10^d)/p. > > If p > 10^d, then by the previous paragraph, there is a fraction n/p (0 > <= n < p) starting with any given d-digit block B. But because p > 10^d, > since 0/p and 1/p both start with block B = 0, so we can strengthen the > condition on n to 1 <= n < p. In other words: > > If p > 10^d, there is a fraction n/p (1 <= n < p) starting with any > given d-digit block B. > > Now suppose p > 10^d and p has primitive root 10. Because of the > primitive root 10, every n with 1 <= n < p is expressible as n = 10^k > mod p for some 0 <= k < p-1. By our previous paragraph, there is then a > fraction of the form (10^k mod p) / p starting with any desired d-digit > block B. From this we conclude that every d-digit block occurs somewhere > in the decimal expansion of 1/p. --this is a nice theorem, and it might be a new one. It is also of some importance because it proves a randomness property of the Marsaglia type of pseudorandom number generator. To determine if it is new, I suggest looking to see if it is in Hardy & Wright chapter 9 as step one. Frankly though, it looks too easy to have escaped previous notice. > In other words, any prime p > 10^d with primitive root 10 includes every > d-digit block in the decimal expansion of its reciprocal. This is why I > must choose the smallest such prime if I am looking for the first number > with that property. > > I can show that the smallest number n whose decimal reciprocal includes > every d-digit block must satisfy n >= 10^d. --well, that seems trivial since the period<n. > I can also show that the > reciprocal of the smallest prime p > 10^d with primitive root 10 > includes every d-digit block. I can show that if p < 2*10^d, --you do not need this "if" because of the known theorem that there is always a prime between X and 2X if X>1: http://en.wikipedia.org/wiki/Bertrand's_postulate > then there > is no other prime q with 10^d <= q < p whose reciprocal includes every > d-digit block. What I can't show is that there is not some composite q > with 10^d <= q < p whose reciprocal includes every d-digit block, though > I seriously doubt such a q could exist. --if q=a*b*...*c where a,b,...,c are primes with 10 primitive modulo each of them, then every n with 1<=n<q and gcd(n,q)=1 is expressible as n = 10^k mod q for some 0<=k<p-1. However, each n with gcd(n,q)>1 definitely is NOT thus-expressible. This tells me that the decimal expansion of 1/q with q composite, q>10, NEVER includes all k-digit blocks where k=floor(log10(q)), because those blocks are exactly the opening bars of 10^j/q. To prove that, there are two cases: case i: q>10 contains a prime p with p>10. Then the multiples of q/p are unrepresentable and include examples below 10^k case ii: q>10 contains a prime p<10, then the multiples of p do the job. QED.

10 13

[math-fun] Another puzzle
by David Wilson 02 Feb '12

02 Feb '12

What n satisfy 2^n == 3n+2 (mod n^2 - 2n)

3 2

Re: [math-fun] math-fun Digest, Vol 107, Issue 42
by Warren Smith 31 Jan '12

31 Jan '12

> From: Andy Latto <andy.latto(a)pobox.com> >>> I can also show that the >>> reciprocal of the smallest prime p > 10^d with primitive root 10 >>> includes every d-digit block. I can show that if p < 2*10^d, >> >> --you do not need this "if" because of the known >> theorem that there is always a prime between X and 2X if X>1: >> http://en.wikipedia.org/wiki/Bertrand's_postulate > > Bertrand's postulate tells us there is a prime between 10^d and > 2*10^d. But what I think David Wilson needs as an assumption is the > stronger statement that there is a prime p between 10^d and 2*10^d > such that 10 is a primitive root mod p. > Andy --Ah. Andy is right. And that is devastating because (pathetically!) it is not even known that 10 is a prim-root modulo more than a finite set of primes: http://en.wikipedia.org/wiki/Artin's_conjecture_on_primitive_roots So while I have no doubt whatever that the "10ified Bertrand conjecture" is true, a proof seems out of reach.

1 0

[math-fun] outlandish 2F1[4/5]
by Bill Gosper 31 Jan '12

31 Jan '12

Hypergeometric2F1[1/8, 5/8, 3/4, 4/5] == (Pi*((Sqrt[1 + Sqrt[5]]*(10*(1 + Sqrt[5])*Pi + Sqrt[5]*Log[5] + 5*Log[9 - 4*Sqrt[5]]))/ (Gamma[11/20]*Gamma[19/20]) - (Sqrt[5 - Sqrt[5]]*(2*(-5 + 3*Sqrt[5])*Pi + Sqrt[5]*Log[5] + 5*Log[9 - 4*Sqrt[5]]))/ (Gamma[13/20]*Gamma[17/20])))/(20*Gamma[1/4]^2) Mma 8.04 knows squat about 2F1[4/5]. Macsyma knows most, but divides by 0 here. But it could take the requisite limit, which Mma couldn't. But Macsyma couldn't simplify it, which Mma could. (You can never have too many CASes.) DLMF 15.4iii just quotes A&S and still doesn't even know the 2F1[1/5]s. (And rejected my services, back in the day.) --rwg

1 0

[math-fun] decimal expansions of P/Q in radix R and their uses
by Warren Smith 30 Jan '12

30 Jan '12

Hardy & Wright "Intro to the theory of numbers" chapter 9 is on this topic. Assume gcd(P,Q)=1. thm 135: P/Q has terminating or eventually-periodic expansion. If gcd(Q,R)=1 then periodic, otherwise terminating. The period length is the order of R mod Q. So, the conjecture that 1/Q where Q is "the first prime greater than 10^n" will have period exceeding 10^n is undoubtably false eventually, in fact 10 should presumably be a square mod Q eventually about half the time so it should be false at least half the time. You can check that P/Q with Q prime will have maximum period in radix R by computing the order of R mod Q and verifying it is Q-1, which you need to do by computing R^((Q-1)/d) mod Q for each divisor d of Q-1, and verifying none of these are 1. The idea that we can identify numbers like 1.004009016025036... by use of "generating functions" is silly in comparison with the much better and simpler idea that you can find the optimum rational approximations to ANY decimal X by use of (Euclid's) continued fraction expansion of X: 1.004009016025036049064081100 = [1, 249, 2, 3, 1, 1, 14, 7, 4, 2, 4, 7, 7, 1, large] then 1 + 1/(249+1/(2+1/(3+1/(1+1/(1+1/(14+1/(7+1/(4+1/(2+1/(4+1/(7+1/(7+1)))))))))))); = 1001000000 / 997002999 = 1.004009016025036049064081100121144169196225256289324361... Marsaglia suggested an actual use for all this, which is that the radix R expansion of some rational P/Q with huge period could provide a good and efficient pseudorandom number generator. It is possible by clever choice of P,Q,R to obtain huge periods (e.g. of order 2^500) but with no explicit use of multiprecision arithmetic. These go under the names "Marsaglia/Zaman add-with-carry, subtract-with-borrow, and multiply-with-carry" generators and they are implemented with extremely simple recurrences that effectively implement division of a huge multiprecision number by Q but because of the ultrasparse binary structure of the particular Q he chooses (and in some versions one can generate the quotient "backwards" too) this is easy. As the radix R, Marsaglia uses either 2^wordsize or 2^wordsize-1. As Q, Marsaglia uses a prime of the form Q=a*R^k-1 chosen so that R is a primitive root mod Q or anyhow having as large multiplicative order as possible. As P, pretty much any gunk filling an array of k machine words, will do. Actually, now that I am looking at this, http://en.wikipedia.org/wiki/Multiply-with-carry it occurs to me there ought to be a more general idea: Instead let Q be a prime of form a*R^k-b for integers a,b that fit well within one machine word, such that R has multiplicative order as large as possible mod Q. I think the recurrence-system instead will then be x[n] = (a * x[n-k] + c[n-1]) mod R c[n] = floor( (a*x[n-k] + b*c[n-1]) / R ) for n=k,k+1,k+2,... When implementing, x of course does not need to be an infinite-length array, but only with k entries since you can "wrap" it mod k, and c does not need to be an array at all, only a single machine word. As the math-funners have pointed out, though, certain Q lead to rather ridiculously obviously-nonrandom sequences (expansions of 1/Q mod R), albeit that is disguised to a considerable degree by the fact we are computing P/Q not 1/Q and P is made of truly-random "seed" bits. Any such generator is going to fail a "continued fraction test" of randomness almost instantly. This objection is defeated by: I also once invented a way to do the usual linear congruential generator x <-- (A*x+B) mod M to get huge periods without explicit use of multiprecision arithmetic, but it is considerably more complicated than Marsaglia's idea, although also quite efficient and probably "more random." My idea was to carefully choose M to be a prime of the form 2^a+-2^b+-1 where a and b are multiples of the computer wordsize; A is chosen to be a "sparse" multiprecision (a words long) number which happens to be a primitive root mod M. The sparsity structure of A is carefully selected to try to get both computational efficiency and fast graph-theoretical "expansion."

4 3

[math-fun] a fraction for .000100040009001600250036... (which?)
by Robert Munafo 30 Jan '12

30 Jan '12

The question (from Pacher Christoph): What rational fraction is equal to .000100040009001600250036...? can be answered easily using (ordinary) generating functions*. First look up the sequence 1,4,9,16,25,36 on OEIS: http://oeis.org/A000290 In the "FORMULA" section we find the (odinary) generating function: G.f.: x * (1 + x) / (1 - x)^3. Since we want each term to be 4 digits of the decimal expansion, let x equal 0.0001. x = 0.0001 x*(1+x) / (1-x)^3 = 0.00010001 / 0.999700029999 Then turn it into a reduced fraction of integers: 0.00010001 / 0.999700029999 = 100010000/999700029999 = (2^4*5^4*73*137)/(3^6*11^3*101^3) The prime factorization of 100010000 and 999700029999 shows that the fraction is in reduced form. Verifying the result, to 300 digits using the Unx tool bc: scale=300 100010000/999700029999 .0001000400090016002500360049006400810100012101440169019602250256028\ 90324036104000441048405290576062506760729078408410900096110241089115\ 61225129613691444152116001681176418491936202521162209230424012500260\ 12704280929163025313632493364348136003721384439694096422543564489462\ 44761490050415184532954765625 We get a two-digit version by letting x=0.01, which gives 10100/970299=0.0104091625364964... what I could not find there but it also quite well-known: > > The decimal expansions of the multiples of 1/7 (and of inverses of > so-called full reptend primes) are cyclic shifted versions of each > other. There exists a simple card trick (which I found in a book of Martin > Gardner) based on this fact. > > Christoph > Integers whose reciprocal has N-1 digits are "the primes with primitive root 10", http://oeis.org/A001913 The card trick is here: http://www.ams.org/samplings/feature-column/fcarc-mulcahy2 (found with Google "142857 card trick") - Robert * Finding the ordinary generating function for an integer sequence is the inverse of the process of finding the a Taylor/Maclaurin series expansion for a function. It is a bit surprising that this is how you would answer this type of decimal fraction question, but there it is. -- Robert Munafo -- mrob.com Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 - mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com

2 2