November 2006 - mathfun@mailman.xmission.com

[math-fun] First time harmonic series exceeds 100
by Daniel Asimov 10 Nov '06

10 Nov '06

I've long wondered what is the smallest n for which (*) 1 + 1/2 + 1/3 + . . . + 1/n > 100. Roughly, I thought it was where ln(n) = 100, i.e., n ~ e^100 ~ 2.688 x 10^43. But more exactly, it would be closer to where 100 - ln(n) = gamma (Euler's), i.e. let K = floor( e^(100-gamma) ). Then n should be fairly close to K (which is roughly only 1.509 x 10^43). Mathematica gives K as exactly K = 15092688622113788323693563264538101449859497 and to my surprise it claims that 1 + 1/2 + 1/3 + . . . + 1/K > 100 but 1 + 1/2 + 1/3 + . . . + 1/(K-1) < 100 (!). Can someone confirm that this K is exactly the least n for which (*) holds? Thanks, Dan P.S. Anyone know a useful upper bound U(n) for the error E(n) = |(1+1/2+1/3+...+1/n) - ln(n) - gamma| ? Evidently E(n) gets very small very fast.

7 9

[math-fun] continued subscripts
by Marc LeBrun 10 Nov '06

10 Nov '06

Recent correspondence has inspired me that, analogous to stacking superscripts for iterated exponentiation, we may stack the subscripts that are vulgarly used as "base" notation. For example, we may apprehend, among 239's many sterling properties, the following gratuity: 10 11 12 13 14 15 16 = 239. 17 18 19 20 21 22 23 That is, working up from the bottom: ... ... ... 21 = 21 = 97 and so on = 239. 22 48 23 Question: If stacking superscripts produces a "tower", what are stacked subscripts? A "dungeon"? The sequence following this incrementing pattern of subscripts, namely, 10, 10, 10, 10, ... 11 11 11 12 12 13 has (if I've computed them correctly) the values 10 11 13 16 20 25 31 38 46 55 65 87 135 239 463 943 1967 4143 8751 18479 38959 103471 306223 942127 2932783 9153583 28562479 89028655 277145647 861652015 2675637295 10173443119 41132125231 168836688943 695134284847 2862242032687 11778342480943 48433422101551 199016451894319 817199416306735 3353334654921775 16351027752823855 82924577766468655 423717750455364655 2167310727003204655 11083411175259204655 56654591244123204655 289463880726363204655 1478277273827163204655 7546178967779163204655 38504861079779163204655 227971995605219163204655 1387065053878499163204655 8475364910242019163204655 51807613088766179163204655 ... (which I will submit to the OEIS as soon as I can craft a concise explanation!). Puzzle: the least-significant digits seem to be stabilizing ... why or why not? Of course we can write the numerals in the stack in other bases, such as binary: 10, 10, 10, 10, ... 11 11 11 100 100 101 and these have as their first few values 2 3 5 26 1370 9840770 851566070376026 ... but it grows rapidly--the next two terms have ~144 and ~450 digits respectively (and seem beginning to show a distinct tendency to end with the digits ...26 ?!) Question: What is a good approximate formula for a(n)? And of course other continued subscript patterns besides simply incrementing the terms are possible.

2 1

Re: [math-fun] World's hardest sudoku puzzle?
by Steve Rowley 10 Nov '06

10 Nov '06

[Math-Fun added back, as this was kinda fun. To me, anyway.] >Date: Wed, 08 Nov 2006 17:49:48 -0800 >From: Nick Baxter <nickb(a)baxterweb.com> > >I have two questions for you. What is the "R library"? How long does it >take to solve a "trivial" Sudoku puzzle? There are 3 things to understand here: S, SPlus, and R. * S is a language for statisticians designed in the late 80s. Think Scheme + APL + every statistical library you can think of. (And a few odd things, like the sudoku library contributed by Seeliger, Bengtsson, and Brahm.) * SPlus is a commercial implementation of S. * R is a free, open-source implementation of S. Seems to be more reliable than SPlus: I've had same-day fixes from library authors; you can't buy support like that! http://www.r-project.org/ As for your question re trivial sudokus: if by "trivial", you mean the easy ones published in some newspapers, I have no idea how those are selected. (I'm not really a puzzle guy.) On the other hand, if you mean a "typical" sudoku in the sense that it's randomly selected, then we can MEASURE that! Assuming that "difficulty" translates into computer time, which is of course dubious... I took a sample of 500 random sudokus: the usual 9x9 grids, random solutions, then 58 of 81 cells randomly blanked. (The sudoku library suggests 50/81 blanks, but the "world's most difficult sudoku" has 58, so I went with that. Obviously this parameter could be changed, or even randomly drawn from some distribution.) Solution times were measured in seconds, on an IBM T40 laptop with 2gb of memory running R under Windows2000. Here are the mean, std dev, skewness, kurtosis, range, and quartiles. mean stdDev skewness kurtosis 11.09398 51.72290 18.90855 389.56443 0% 25% 50% 75% 100% 2.67000 5.37750 6.32000 7.21500 1095.12000 Note that the mean is much larger than the median, so there are a lot of outliers. That's confirmed by the hardest one, which took 1095 sec to solve. So our "world's hardest sudoku" is about 5x harder than the median one with similar parameters, but nowhere near the hardest found in this sample! The R script which generated these data was: > library("sudoku") > library("moments") > source("profile.r"); profile("solveSudoku") > times <- unlist(sapply(1:500, function(ignore) { solveSudoku(generateSudoku(Nblank = 58), print.it = FALSE); result <- profileResults(); resetProfiles(); result })["user", ]) > c("mean" = mean(times), "stdDev" = sd(times), "skewness" = skewness(times), "kurtosis" = kurtosis(times), quantile(times, probs = c(0.00, 0.25, 0.50, 0.75, 1.00))) -- Steve Rowley <sgr(a)alum.mit.edu> http://alum.mit.edu/www/sgr/ Skype: sgr000

1 0

RE: [math-fun] First time harmonic series exceeds 100
by Daniel Asimov 10 Nov '06

10 Nov '06

The harmonic number H(n) := 1 + 1/2 + 1/3 + . . . + 1/n has the interesting property that although it -> oo as n does, it's never an integer. So let H(n) = p(n)/q(n) in lowest terms. (I have full confidence that p and q are already in the OEIS.) What is known about the factorizations of the p(n)'s and q(n)'s ? In particular, are their factorizations known or believed to have statistical properties atypical of numbers of the same size (re the number, size, and exponents of their prime factors) ? Do the factors of p(n) (or q(n)) have special number-theoretical properties? (I.e., what is known about the primes -- if any -- that *never* occur as a factor of any p(n) (q(n)) ? --Dan

4 4

[math-fun] Topology puzzle
by Daniel Asimov 09 Nov '06

09 Nov '06

Let S^3 be the unit sphere in R^4 = C^2 given by {(z,w} : |z|^2 +|w|^2 = 1}, and let C and D be the great circles C = {(e^is, 0) : 0 <= s < 2pi}, D = {(0, e^it) : 0 <= t < 2pi}. Let R be a rotation of S^3 that interchanges C and D. PUZZLE: Identify the topology of the manifold X obtained as the quotient space of S^3 by the action of R, i.e. X = S^3 / p ~ R(p). --Dan

1 0

[math-fun] World's hardest sudoku puzzle?
by Christian Boyer 08 Nov '06

08 Nov '06

Agence France-Presse news (Escargot = Snail in French). A FINNISH mathematician claimed he had created the world's hardest sudoku puzzle, a brain-teaser which required three months' work and a billion combinations to produce. "AI Escargot is the most difficult sudoku-puzzle known so far," the puzzle's 37-year-old creator and applied mathematician Arto Inkala said. "I called the puzzle AI Escargot, because it looks like a snail. Solving it is like an intellectual culinary pleasure. AI are my initials," Mr Inkala said. According to a rating published on sudoku website www.sudoku.com, AI Escargot claims the top spot for sudoku's most baffling puzzles. Escargot demands those tackling it to consider eight casual relationships simultaneously, while the most complicated variants attempted by the general public require people to think of only one or two combinations at any one time. AI Escargot has, however, been solved by sudoku experts, but its creator has promised to produce more challenging versions, as with Escargot, with the aid of computers. +-------+-------+-------+ | 1 . . | . . 7 | . 9 . | | . 3 . | . 2 . | . . 8 | | . . 9 | 6 . . | 5 . . | +-------+-------+-------+ | . . 5 | 3 . . | 9 . . | | . 1 . | . 8 . | . . 2 | | 6 . . | . . 4 | . . . | +-------+-------+-------+ | 3 . . | . . . | . 1 . | | . 4 . | . . . | . . 7 | | . . 7 | . . . | 3 . . | +-------+-------+-------+ References: Text above is a copy and paste from: http://www.news.com.au/story/0,23599,20714995-38200,00.html The puzzle was unfortunately not in the news, but I found the grid here: http://www.sudoku.com/forums/viewtopic.php?t=4870 Christian.

3 2

[math-fun] Question someone posted to sci.math.research
by Daniel Asimov 07 Nov '06

07 Nov '06

Here's a possibly interesting question that someone signing his name as Brendan posted to sci.math.research. I have edited it slightly for clarity. << Consider a set of n points in the plane with diameter = 1, whose average (centroid) is the origin. Choose two of the points whose separation is greatest and replace each of them by their midpoint (so two separated points are replaced by two equal points). The number of points remains at n and their centroid remains at the origin. Again choose two points whose separation is (now) greatest... and so on repeatedly. Assume that for any given such point set, when there are choices to be made, they are made so as to *maximize* the number of steps lrequired to reach diameter = 1/2. Let f(n) be the maximum over all such point sets of this maximum number of steps until the point set reaches diameter = 1/2. What is f(n)? >> --Dan

3 2

[math-fun] Convergence question
by Daniel Asimov 06 Nov '06

06 Nov '06

We know the series for Zeta(s) converges for all s > 1. So, I'm wondering about functions e(n) of n that decrease to 1 slowly enough that A(e) = Sum_{n=1..oo} 1/n^e(n) converges. For example, what about letting e(n) = f_c(n) := 1 + 1/n^c ? Does A(f_c) converge for c close enough to 0, and if so, for which c ? (E.g., I tried using Mathematica to get an idea for c = 1/2, and was unable to tell after 30 minutes whether A(f_(1/2)) converges.) --Dan

1 0

[math-fun] A fast (written in assembly) library for continued fractions
by Thomas Baruchel 03 Nov '06

03 Nov '06

I wrote for my own purposes the following X86 assembly code, intended to be called from GFORTRAN. Such a thing should only work if you have a similar configuration to mine, but anyway, it may interest you. I wanted a fast library for computing the length in the continued fraction expansion of sqrt(n) for any integer n that can fit in a REAL*8 gfortran type for statistical purposes. I also included a faster function that tells if this length is odd or even. Though I am not an assembly expert, there should be a significant speed improvement if you use this code rather than your own fortran functions, since the whole process is performed on the FPU stack with no load/store useless travels. It may help you, it may not. You know about these kinds of software... Instructions are at the beginning of the file, as comments. Read them carefully. Code is below. Regards. ------------------ FILE quad.s BEGINS BELOW -------------------- /* Compile your piece of code with something like gfortran -fno-underscoring mycode.f90 quad.s and add freely your own flags. Enjoy ! Of course, if you don't want to use the -fno-underscoring option for some reason, you can add an underscore at the end of the global functions below, like that: .globl initfpu_ initfpu_: it should work. Remember this piece of code has the following status: Seems to work for me, maybe it will also seem for you, maybe your computer will explode, maybe it will make the earth explode, maybe God will choose not to have created the world, think about it before using it. (Anyway, if you see some kind of bug, please let me know ;-) Remember that a) this is x86 assembly code; b) it uses the GNU AS (gas) syntax; c) it uses the gfortran calling conventions; (now you understand why you can't expect it will works ;-) Usage. In your program (or subroutine or whatever), declare the function you need like: logical :: oddperiod integer*8 :: period (the two subroutines don't need to be declared). Then use integer*4 :: oldfpustatus call initfpu(oldfpustatus) before any computation; if you don't do it, it won't work. When you don't need any more the library, call restorefpu(oldfpustatus) will restore the control register of the FPU. During all that time, you shouldn't be disturbed too much, for your own mathematical operations, but remember the rounding mode of the FPU has been set to "round down" (rather than the usual "nearest"). How does it work? The whole computation is performed on the stack of the FPU; only integer values are needed for these tasks, but they are stored as REAL*10, 80 bits reals. They are not converted to anything during the whole process, thus it should work for great integers stored in your own REAL*8, 64 bits reals. Of course, it may happen (I think it will never, but didn't check) that an intermediate result can not fit in such a REAL*10 register; in that case, read above what is written about the end of the world... */ /* call initfpu(integer*4) ===================== This subroutine should be called before any other call to the library; its purpose is to set the rounding control in the FPU control register to the "round down" mode. The subroutine puts in the integer variable the initial value of the FPU control register. If you want to use often this subroutine (which you shouldn't: speed being an issue, this subroutine is intended to be used once, then the functions billions of times, not the contrary ;-), maybe removing the "fnclex" can improve the speed, since I don't think it plays a great role here. */ .globl initfpu initfpu: pushl %ebp movl %esp, %ebp movl 8(%ebp), %eax fnstcw (%eax) fwait movw (%eax), %dx andb $0xf3, %dh orb $0x4, %dh movw %dx, 2(%eax) fnclex fldcw 2(%eax) fwait movw $0, 2(%eax) leave ret /* call restorefpu(integer*4) ======================== This subroutine should be used for restoring the initial state of the FPU control register; the integer variable must contain the right value (which should have been previously set by the subroutine initfpu). */ .globl restorefpu restorefpu: pushl %ebp movl %esp, %ebp movl 8(%ebp), %eax fclex fldcw (%eax) leave ret /* sqrtint( real double precision : REAL*8 ) ========================================= This function returns the integer part of the square root of the given real. The return value has itself the real type. */ .global sqrtint sqrtint: pushl %ebp movl %esp, %ebp movl 8(%ebp), %eax fldl (%eax) fsqrt frndint leave ret /* oddperiod( real double precision : REAL*8 ) =========================================== Despite the type of the real argument, the variable has to be some integer value stored as a double precision. The function has the logical type. It returns FALSE if the argument is a perfect square or if the continued fraction expansion of its square root has an even period length; it returns TRUE in the remaining case. oddperiod(dble(25)) --> FALSE oddperiod(dble(13)) --> TRUE oddperiod(dble(12)) --> FALSE If it may help, you can declare the function as an integer*4 type rather than a logical one (output will be 0 or 1). Though it seems to work with integer*8, I am not absolutely sure it *always* will (unless you add the following instruction movl $0, %edx at the beginning of the code (for instance between the current second and third lines of code), but really why would you do it? Maybe something in the calling conventions is such that the register EDX already contains 0, but, again, who cares? If you read the instructions for the period() function, you will learn about some 2^32 limitations, but you don't have to care about that here: the length is not kept in any register during the computation and the code should work as long as the values in the FPU stack are the right one. */ .global oddperiod oddperiod: pushl %ebp movl %esp, %ebp movl 8(%ebp), %eax fldln2 fldl (%eax) fld %st fsqrt frndint fld %st fld %st fmul %st, %st fld %st(3) fsubp %st, %st(1) fcom %st(4) fstsw %ax fwait sahf ja oddperiod0 ffree %st ffree %st(1) ffree %st(2) ffree %st(3) ffree %st(4) oddperiod2: movl $0, %eax leave ret oddperiod0: fld %st(1) fadd %st(3), %st fdiv %st(1), %st frndint fld %st(2) fld %st(1) fmul %st(3), %st fsubp %st, %st(1) fxch %st(3) fsub %st(3), %st fabs fcomp %st(6) fstsw %ax fwait sahf ja oddperiod1 fld %st(3) fadd %st, %st fsubp %st, %st(1) fcomp %st(5) fstsw %ax fwait ffree %st ffree %st(1) ffree %st(2) ffree %st(3) ffree %st(4) sahf ja oddperiod2 movl $1, %eax leave ret oddperiod1: ffree %st fincstp fld %st(1) fmul %st, %st fld %st(4) fsubp %st, %st(1) fdivp %st, %st(1) frndint jmp oddperiod0 /* period( real double precision : REAL*8 ) ======================================== Despite the type of the real argument, the variable has to be some integer value stored as a double precision. The function has integer*8 type, though the output can't be greater than 2^32-1 (this choice has been made for avoiding sign problem when the result is greater (or equal) then 2^31 and smaller then 2^32. The result is the length of period in the continued fraction expansion of the square root of the given integer. If speed is really for you an obsession, remember that it is highly unprobable that you use that library for periods between the two fatidic values. Maybe it will increase a little the speed to use integer*4 which allows to delete twice the line movl $0,%edx (not the first time it appears, but the second and third one). If you try such a thing and see some speed improvement, please let me know. */ .global period period: pushl %ebp movl %esp, %ebp movl 8(%ebp), %eax movl $0, %edx fldln2 fldl (%eax) fld %st fsqrt frndint fld %st fld %st fmul %st, %st fld %st(3) fsubp %st, %st(1) fcom %st(4) fstsw %ax fwait sahf ja period0 ffree %st ffree %st(1) ffree %st(2) ffree %st(3) ffree %st(4) movl %edx, %eax movl $0, %edx leave ret period0: inc %edx fld %st(1) fadd %st(3), %st fdiv %st(1), %st frndint fld %st(2) fld %st(1) fmul %st(3), %st fsubp %st, %st(1) fxch %st(3) fsub %st(3), %st fabs fcomp %st(6) fstsw %ax fwait sahf ja period1 fld %st(3) fadd %st, %st fsubp %st, %st(1) fcomp %st(5) fstsw %ax fwait ffree %st ffree %st(1) ffree %st(2) ffree %st(3) ffree %st(4) sahf jb period2 add %edx, %edx period2: movl %edx, %eax movl $0,%edx leave ret period1: ffree %st fincstp fld %st(1) fmul %st, %st fld %st(4) fsubp %st, %st(1) fdivp %st, %st(1) frndint jmp period0

1 0

[math-fun] 2-adic pi: a null result
by Schroeppel, Richard 01 Nov '06

01 Nov '06

I gave a talk last week on the topic of Numbers, and included a slide about P-adic numbers. I mentioned in passing that a P-adic field includes all rationals, some algebraics, and some of the usual transcendental functions. In particular, those defined by convergent power series, when they converge P-adically. For example, e^4 = 1 + 4 + 16/2! + 64/3! + ... converges 2-adically. Some of the usual identities for exp, log, sin, atan, etc. will hold, because the identities are valid as formal power series. Thus e^log(1+x) = 1+x is 2-adically valid when x is a multiple of 4, since the power series for log(1+x) converges (to a value that's a multiple of 4), and the power series for e^(the log power series) also converges (since the exponent is a multiple of 4). Note that 4/3 qualifies as a multiple of 4 here. I later began wondering about defining 2-adic pi. Ordinary pi can be computed using atan(1/2) + atan(1/3) = pi/4, and I wondered about extending this to P-adics. For my test case, I needed an arctan formula for pi that converges 2-adically. The atan series is atan x = x - x^3/3 + x^5/5 - x^7/7 + ... so any rational x with an even numerator (like 2, or -2/3) will converge 2-adically. A bit of experimenting turned up the complex number product (1+2i)(1+8i)(3+2i) = -65, which gives the arctan identity atan 2 + atan 8 + atan(2/3) = pi. This is valid over the reals; each angle is in the first quadrant. The power series for atan won't converge (on the reals) for 2 and 8, but all three series converge 2-adically. The arithmetic looks like this: atan 2, 2-adic 2 = 1o -8/3 = ... 1010 1010 1ooo 32/5 = ... 1001 101o oooo -128/7 = ... 0100 1ooo oooo 512/9 = ... 001o oooo oooo -2048/11 = ... 1ooo oooo oooo ------------------ ... 0010 1100 1010 atan 8, 2-adic 8 = 1ooo -512/3 = ... 101o oooo oooo ------------------ ... 1010 0000 1000 atan(2/3), 2-adic 2/3 = ... 0101 0101 011o -8/81 = ... 1010 0111 1ooo 32/1215 = ... 0111 111o oooo -128/15309 = ... 1101 1ooo oooo 512/177147 = ... 011o oooo oooo -2048/1948617 = ... 1ooo oooo oooo ------------------ ... 0011 0010 1110 [Lower case o is the same as 0; I've used o for low order 0s for readability. To confirm a fraction, just multiply the binary by the denominator, and verify the numerator: For 2/3, ... 0101 0101 011o * 11 = ... 0000 0000 001o = 2.] However, when I totted up the sum of the three arctangents to compute 2-adic pi, I got an unexpected result: atan 2 = ... 0010 1100 1010 atan 8 = ... 1010 0000 1000 atan 2/3 = ... 0011 0010 1110 --------- ------------------ 2-adic pi = ... 0000 0000 0000 So the conjectured value for 2-adic pi is 0. I thought I might have some kind of wraparound problem, but I can convert my arctangent identity to one that's convergent for both reals and 2-adics, and is definitely valid on the reals, although it's a mediocre pi formula. Use the substitutions atan 2 = atan 2/3 + atan 4/7 atan 8 = atan 2 + atan 6/17 to get pi = 3 atan 2/3 + 2 atan 4/7 + atan 6/17. This converges in both reals and 2-adics, and all the angles are in the first quadrant and less than 45 degrees, so it's a plain vanilla pi formula. But the 2-adic version will also give 0 for pi, since the atan substitutions are 2-adically valid. (All the numerators of the atans are even, so the power series all converge, as does the power series proof of the arctan addition formula atan x + atan y = atan (x+y)/(1-xy). ) This hardly counts as enough numerical work to justify any broad conjecture about P-adic pi, but the result is puzzling. Rich

2 1