Beginners are often startled by (c64) (-1/2)! (d64) sqrt(%pi) but if we are to have (n+1/2)! = (n+1/2)(n-1/2)! for all n, then (c86) (n+1/2)! = h*product(k-1/2,k,1,n+1) n + 1 /===\ 1 | | 1 (d86) (n + -)! = h | | (k - -) 2 | | 2 k = 1 and sqrt(pi) is the only value of h with 1 (n + -)! n! 2 (n + 1)! (d68) -------- < -------- < -------- 1 n! 1 (n - -)! (n + -)! 2 2 as n -> oo. E.g., using the products for (n-1/2)!,...,(n+1)!, (c77) [h*prod(k-1/2,k,1,n),prod(k,k,1,n),h*prod(k-1/2,k,1,n+1),prod(k,k,1,n+1)] n n n + 1 n + 1 /===\ /===\ /===\ /===\ | | 1 | | | | 1 | | (d77) [h | | (k - -), | | k, h | | (k - -), | | k] | | 2 | | | | 2 | | k = 1 k = 1 k = 1 k = 1 the ratios in (d68) become (c78) rest(%)/rest(%,-1) n n + 1 n + 1 /===\ /===\ /===\ | | | | 1 | | | | k h | | (k - -) | | k | | | | 2 | | k = 1 k = 1 k = 1 (d78) [---------------, ---------------, ---------------] n n n + 1 /===\ /===\ /===\ | | 1 | | | | 1 h | | (k - -) | | k h | | (k - -) | | 2 | | | | 2 k = 1 k = 1 k = 1 (c81) sfloat(sqrt(%pi)) (d81) 1.77245 Choosing h too small, n = 22 is large enough to fail: (c79) apply_nouns(subst([n = 22,h = 1.5d0],d78)) (d79) [5.5739412602892d0, 4.03664103177732d0, 5.69780662162896d0] Increasing h not quite enough, (c80) apply_nouns(subst([n = 22,h = 1.75d0],d78)) (d80) [4.77766393739075d0, 4.70941453707354d0, 4.88383424711054d0] With h too big, (c82) apply_nouns(subst([n = 22,h = 1.95d0],d78)) (d82) [4.28764712329939d0, 5.24763334131052d0, 4.38292817048382d0] With h very close, (c83) apply_nouns(subst([n = 22,h = 1.77d0],d78)) (d83) [4.72367903414339d0, 4.76323641749724d0, 4.82864967934658d0] and we'd need to increase n>22 to see which way h is off. This process produces the constant 1.7724..., but doesn't tell us that it's sqrt(pi). But the above arguments lead to the generalization (c4) (product(n*(1+1/n)^z/(n+z),n,1,inf),%% = closedform(%%)) inf 1 z /===\ (- + 1) n | | n (d4) | | ---------- = z! | | z + n n = 1 and for z = -1/2, (c91) subst(-1/2,z,%^2) inf /===\ | | n 2 2 (d91) ( | | -------------------) = (-1/2)! | | 1 1 n = 1 sqrt(- + 1) (n - -) n 2 (after squaring) we can get Wallis's product (c92) factor(intoprod(%)) inf /===\ 3 | | 4 n (d92) | | ------------------ = %pi | | 2 n = 1 (n + 1) (2 n - 1) (c93) factor(prodfudge(%,n/(n-1/2)))/2 inf /===\ 2 | | 4 n %pi (d93) | | ------------------- = --- | | (2 n - 1) (2 n + 1) 2 n = 1 which he presumably got from Euler's infinite product for sin(x)/x. But, as Rich and I mentioned a couple of weeks ago, factorial can be q-extended, inf /===\ i | | 1 - q | | ---------- | | z + i i = 1 1 - q (d97) QF(z, q) = ---------------- z (1 - q) (which becomes (d4) when q -> 1) and when z = -1/2 leads to a pi-like function of q: pi(q) := q^(1/4) QF(-1/2,q^2)^2 . E.g., the q-Wallis product becomes (c1) pi(q)/(1+q) =q^(1/4)*product((1-q^(2*n))/(1-q^(2*n-1)).(1-q^(2*n))/(1-q^(2*n+1)),n,1,inf) inf /===\ 2 n 2 n pi(q) 1/4 | | 1 - q 1 - q (d1) ----- = q | | (------------ . ------------). q + 1 | | 2 n - 1 2 n + 1 n = 1 1 - q 1 - q pi(q) can also be written as a rapidly convergent series: (c3) q^(1/4)*(1-q^2)*sum(q^binomial(n,2),n,1,inf)^2 inf (n - 1) n ==== --------- 1/4 2 \ 2 2 (d3) q (1 - q ) ( > q ) / ==== n = 1 Comparing the Taylor expansions, (c4) taylor(% = (1+q)*rhs(d2),q,0,11) 1/4 5/4 17/4 21/4 25/4 29/4 33/4 (d4)/T/ q + 2 q + q - 2 q + q + 2 q - 3 q 41/4 1/4 5/4 17/4 21/4 25/4 29/4 + 2 q + . . . = q + 2 q + q - 2 q + q + 2 q 33/4 41/4 - 3 q + 2 q + . . . pi(q) can also be written as etas: (c19) pi(q) = (1-q^2)*eta(q^2)^4/eta(q)^2 2 4 2 (1 - q ) eta (q ) (d19) pi(q) = -----------------. 2 eta (q) Recall a week or two back our "least abstruse" eta relation: (c99) 16*eta(q)^8*eta(q^4)^16+eta(q)^16*eta(q^4)^8=eta(q^2)^24 8 16 4 16 8 4 24 2 (d99) 16 eta (q) eta (q ) + eta (q) eta (q ) = eta (q ) This leads to a pi(q) relation: (c5) pi(q) = sqrt(pi(q^2)*(4*pi(q^4)/(q^4+1)+(q^4+1)*pi(q^2)^2/pi(q^4)))/(q^2+1) 4 4 2 2 2 4 pi(q ) (q + 1) pi (q ) sqrt(pi(q ) (-------- + ----------------)) 4 4 q + 1 pi(q ) (d5) pi(q) = ------------------------------------------ 2 q + 1 As Gene Salamin pointed out to me, we can now initialize q to a small number, quickly evaluate the series (d3) for q^2 and q^4, and then use (d5) for q, sqrt q, q^(1/4), ..., 1-epsilon, and pi(1-epsilon) ~~ pi. Let's see how well: To use the recurrence (d5), define a[n] := pi(2^-(2^(4-n))), b[n] := 1+2^-(2^(4-n)): (c23) (bfprecision:25,bfloat(subst(2^-8,q,d16))) (d23) pi(3.90625b-3) = 2.519531250577538086392374b-1 (c24) bfloat(subst(2^-16,q,d16)) (d24) pi(1.52587890625b-5) = 6.250190734863281250338803b-2 (c25) (a[1] : rhs(d23),a[0] : rhs(%), b[0] : 1.0b0+2^-16) (d25) 1.0000152587890625b0 (c26) b[n] := 1+sqrt(b[n-1]-1) (d26) b := 1 + sqrt(b - 1) n n - 1 Then (d5) becomes (c27) a[n] := sqrt(a[n-1]*(b[n-2]*a[n-1]^2/a[n-2]+4*a[n-2]/b[n-2]))/b[n-1] 2 b a 4 a n - 2 n - 1 n - 2 sqrt(a (------------- + --------)) n - 1 a b n - 2 n - 2 (d27) a := --------------------------------------- n b n - 1 and now let's try for pi: (c29) a[2] (d29) 5.625067088994892759372818b-1 (c30) a[3] (d30) 1.06226912564319452424773b0 (c31) a[4] (d31) 1.699635053182282093917339b0 (c32) a[5] (d32) 2.26618007091359690481384b0 . . . (c48) a[21] (d48) 3.141576039984956345244343b0 (c49) a[22] (d49) 3.141584346772731862424423b0 Pretty disappointing! But now look what happens if we divide a[n] by 2^(n-4) (1-2^(5-n)): (c105) makelist(a[n+4]/2^n/(1-2.0b0^-2^(1-n)),n,-4,2) (d105) [1.000030517810962749189647b0, 2.015655756928048047972053b0, 2.258850470247360857097163b0, 2.26617413470548165172849b0, 2.266180070909709458556452b0, 2.26618007091359690481384b0, 2.266180070913596904813842b0] Quadratic convergence! But to what? (c106) %*2*log(2.0b0) (d106) [1.386336667789142000816372b0, 2.794292209788197907696069b0, 3.1314316695169296120584b0, 3.141584424257956628976669b0, 3.141592653584404093636867b0, 3.141592653589793238462641b0, 3.141592653589793238462643b0] I.e., pi/(2 ln 2). --rwg PS, speaking of pi/(2 ln 2) the first identity in http://www.tweedledum.com/rwg/idents.htm is equivalent to %pi limit 1/(n log(---) - log(atan(n)) - log(atan(n + 1))- . . . n -> inf 2 %pi - log(atan(2 n - 1))) = -------- 2 log(2)