[math-fun] an unreasonably handy tool
[In my integrate atan^2 msg, I should have made it clearer that Mma was refusing on license-enforcement grounds. Ed Pegg suggested integrals.wolfram.com, which gives a somewhat messy, complex indefinite integral, and apparently doesn't offer the definite option. Wouter, the upper limit wants to be Tan[n*Pi/d], not Pi*n/d, thanks.] Back in the 70's, Rich Schroeppel informally described his integer lattice algorithm to Robert Maas, who called it LINREL for LINear RELation finder. Almost trivial by comparison is my TLINREL, which seeks linear relations among Taylor expansions using polynomials, instead of integers, as coefficients. Suppose you need the formula for inf ==== 2 \ n > q , / ==== n = 1 but all you remember is that the bilateral series had something to do with the infinite products inf inf /===\ /===\ | | n | | n | | (1 - q ), | | (q + 1), and maybe | | | | n = 1 n = 1 inf /===\ | | 2 n | | (q + 1). | | n = 1 Equate the logs of all four quantities with expansions to a dozen or so terms: (c197) (log(1+2*sum(q^n^2,n,1,inf)),%% = taylor(%%,q,0,9)) Time= 10 msec. inf ==== 2 3 5 \ n 2 8 q 4 12 q (d197)/T/ log(2 > q + 1) = 2 q - 2 q + ---- - 2 q + ----- / 3 5 ==== n = 1 6 7 9 8 q 16 q 8 26 q - ---- + ----- - 2 q + ----- + . . . 3 7 9 (c198) (log(product(1-q^n,n,1,inf)),%% = taylor(%%,q,0,9)) Time= 10 msec. inf /===\ 2 3 4 5 | | n 3 q 4 q 7 q 6 q 6 (d198)/T/ log( | | (1 - q )) = - q - ---- - ---- - ---- - ---- - 2 q | | 2 3 4 5 n = 1 7 8 9 8 q 15 q 13 q - ---- - ----- - ----- + . . . 7 8 9 (c199) (log(product(1+q^n,n,1,inf)),%% = taylor(%%,q,0,9)) Time= 10 msec. inf /===\ 2 3 4 5 6 | | n q 4 q q 6 q 2 q (d199)/T/ log( | | (q + 1)) = q + -- + ---- + -- + ---- + ---- | | 2 3 4 5 3 n = 1 7 8 9 8 q q 13 q + ---- + -- + ----- + . . . 7 8 9 (c200) (log(product(1+q^(2*n),n,1,inf)),%% = taylor(%%,q,0,9)) Time= 10 msec. inf /===\ 4 6 8 | | 2 n 2 q 4 q q (d200)/T/ log( | | (q + 1)) = q + -- + ---- + -- + . . . | | 2 3 4 n = 1 Then simply (c201) tlinrel(d197,d198,d199,d200) Time= 602 msec. inf ==== 2 inf \ n /===\ log(2 > q + 1) | | 2 n / log( | | (q + 1)) ==== | | n = 1 n = 1 (d201)/T/ - -------------------- - --------------------- 4 4 2 q q inf inf /===\ /===\ | | n | | n 3 log( | | (q + 1)) log( | | (1 - q )) | | | | n = 1 n = 1 + --------------------- + ------------------- = 0 + . . . 4 4 2 q 2 q Solve for the desired sum: (c202) solve(lhs(%),sum(q^n^2,n,1,inf)) Time= 50 msec. inf ==== 2 \ n (d202) [ > q = / ==== n = 1 inf inf inf /===\ /===\ /===\ | | 2 n 2 | | n | | n 3 ( | | (q + 1)) - ( | | (1 - q )) ( | | (q + 1)) | | | | | | n = 1 n = 1 n = 1 - --------------------------------------------------------] inf /===\ | | 2 n 2 2 ( | | (q + 1)) | | n = 1 For brevity, combine the products: (c203) prodcontract(intoprod(expand(%[1]))) Time= 30 msec. inf /===\ n n 3 | | (1 - q ) (q + 1) inf | | ------------------ ==== 2 | | 2 n 2 \ n n = 1 (q + 1) 1 (d203) > q = ------------------------ - - / 2 2 ==== n = 1 Test: (c204) taylor(rhs(%),q,0,100) Time= 2283 msec. 4 9 16 25 36 49 64 81 100 (d204)/T/ q + q + q + q + q + q + q + q + q + q + . . . One more example. The series for log(1-x)^2 is "triangular hypergeometric", i.e., it requires a double sum, though only a single matrix product: [ 2 x ] inf [ x --- 0 ] /===\ [ n ] [ 2 ] | | [ ] [ 0 0 log (1 - x) ] | | [ x ] = [ ] . | | [ 0 x - ] [ 0 0 - log(1 - x) ] n = 1 [ n ] [ ] [ ] [ 0 0 1 ] [ 0 0 1 ] But I wanted [ n x 1 ] inf [ ----- - 0 ] /===\ [ n + 2 2 ] [ 0 0 z(x) ] | | [ ] [ ] | | [ n x ] = [ 0 0 y(x) ], | | [ 0 ----- 1 ] [ ] n = 1 [ n + 2 ] [ 0 0 1 ] [ ] [ 0 0 1 ] which is sort of "contiguous" to the previous product, so maybe z(x) and y(x) likewise to come out in terms of log(1-x) and log(1-x)^2. Taking the right column of the above product to 15 terms: (c209) taylor(col(prud(%,n,1,15),3),x,0,14) Time= 481 msec. 2 3 4 5 6 1 5 x x 49 x 41 x 109 x 853 x (d209)/T/ matrix([- + --- + -- + ----- + ----- + ------ + ------ 2 12 3 180 180 560 5040 7 8 9 10 11 12 209 x 1679 x 19981 x 18167 x 252341 x 3488333 x + ------ + ------- + -------- + --------- + ---------- + ----------- 1400 12600 166320 166320 2522520 37837800 13 2 3 4 5 6 7 8 3694253 x x x x x x x x x + ----------- + . . .], [1 + - + -- + -- + -- + -- + -- + -- + -- 43243200 3 6 10 15 21 28 36 45 9 10 11 12 13 14 x x x x x x + -- + --- + --- + --- + --- + --- + . . .], [1 + . . .]) 55 66 78 91 105 120 Expanding the log (over x, a non-essential micro-optimization): (c210) (log(1-x)/x,%% = taylor(%%,x,0,14)) Time= 0 msec. 2 3 4 5 6 7 8 9 log(1 - x) x x x x x x x x x (d210)/T/ ---------- = - 1 - - - -- - -- - -- - -- - -- - -- - -- - -- x 2 3 4 5 6 7 8 9 10 10 11 12 13 14 x x x x x - --- - --- - --- - --- - --- + . . . 11 12 13 14 15 Finally, as a stunt, solve for both z znd y simultaneously: (c211) tlinrel(%,%^2,a*(z = d209[1,1])+b*(y = d209[2,1])) Time= 27139 msec. 3 2 2 3 2520 (128 b - 216 a b + 3 a b + 4 a ) (a z + b y) (d211)/T/ ---------------------------------------------------- 2 4 b (2 b + a) (16 b - 27 a) x 3 2 2 3 + 2520 (128 b - 216 a b + 3 a b + 4 a ) log(1 - x) 2 6 (2 b x - 2 b + a)/(b (2 b + a) (16 b - 27 a) x ) 3 2 2 3 5040 (128 b - 216 a b + 3 a b + 4 a ) - ---------------------------------------- 2 5 (2 b + a) (16 b - 27 a) x 3 2 2 3 2 2520 a (128 b - 216 a b + 3 a b + 4 a ) log (1 - x) (x - 1) - -------------------------------------------------------------- = 2 7 b (2 b + a) (16 b - 27 a) x 0 + . . . (c212) solve(lhs(%),a*z+b*y)[1] Time= 130 msec. 2 (d212) a z + b y = - ((2 b log(1 - x) - 2 b) x 2 2 3 + ((a - 2 b) log(1 - x) - a log (1 - x)) x + a log (1 - x))/x The function is strikingly short: /* -*- Mode: Macsyma -*- */ tlinrel([eqns]):=catch(block([len:length(eqns),mat,rot], push(1=1,eqns), mat:addcol(transpose(map('rhs,eqns)), apply('taylor,cons(ident(1+len),first(taylorinfo(rhs(eqns[2])))))), rot:addrow(ematrix(1,len+1,1,1,len+1),addcol(ident(len),ematrix(len,1,0,1,1)-1)), do (for k:2 thru len+1 do (mat[k], if %%[1]=0 then throw(eqns.sqfr(rest(%%))) else mat[k]:%%/first(%%[1])), mat:rot.mat)))$ --rwg Any insufficiently advanced technology is indistinguishable from Microsoft.
I said,
But I wanted [ n x 1 ] inf [ ----- - 0 ] /===\ [ n + 2 2 ] [ 0 0 z(x) ] | | [ ] [ ] | | [ n x ] = [ 0 0 y(x) ], | | [ 0 ----- 1 ] [ ] n = 1 [ n + 2 ] [ 0 0 1 ] [ ] [ 0 0 1 ]
and we used the TLINREL function to determine log(1 - x) (log(1 - x) (x - 1) - x) z(x) = -----------------------------------, 3 x 2 (log(1 - x) (x - 1) - x) y(x) = - --------------------------. 2 x I thought you'd never ask! This product is interesting because it is prod(N(2,n),n,1,inf), where [ n x 1 ] [ ----- - 0 ] [ n + k k ] [ ] N(k, n) := [ n x ], [ 0 ----- 1 ] [ n + k ] [ ] [ 0 0 1 ] which with K(k, n) := [ k x (n + 1) x + k (x - 1) - n 1 ] [ --------------- - ------------------------- ------------------ ] [ (n + k) (x - 1) 2 2 ] [ (k + 1) (n + k) (x - 1) k (k + 1) (x - 1) ] [ ] [ k x 1 ] [ 0 --------------- - ----- ], [ (n + k) (x - 1) x - 1 ] [ ] [ 0 0 1 ] forms a path-invariant pair: N(k,n) K(k,n+1) = K(k,n) N(k+1,n). Equating alternate paths from k=a, n=b around two sides of the infinite rectangle with opposite corner approaching k=oo, n=oo gives an identity relating "triangular" sums, one in powers of x, and the other in powers of x/(x-1). This corresponds to the "linear" transformation of ordinary hypergeometrics, A&S 15.3.4 More interestingly, we can take a diagonal path from (a,b), N(a,b) K(a,b+1) N(a+1,b+1) K(a+1,b+2) = J(0) J(1) J(2) ..., where J(j) = N(j+a,j+b) K(j+a,j+b+1), which with a=1, b=1 is [ 2 ] [ j x x (2 x - 1) j (2 x - 2) - 1 ] [ ----------------- - -------------------------- - -------------------- ] [ 1 1 2 1 2 ] [ 4 (j + -) (x - 1) 4 (j + -) (j + 1) (x - 1) 2 j (j + -) (x - 1) ] [ 2 2 2 ] [ ] [ 2 ] [ j x j (x - 2) - 1 ] [ 0 ----------------- ----------------- ] [ 1 1 ] [ 4 (j + -) (x - 1) 2 (j + -) (x - 1) ] [ 2 2 ] [ ] [ 0 0 1 ] This is analogous the the "quadratic transformation" (A&S 15.3.15), and computes [ 1 log(1 - x) (log(1 - x) + 2) 1 ] [ 0 0 - + --------------------------- - ----- ] [ x 2 x - 1 ] [ 2 x ] [ ] [ log(1 - x) ] [ 0 0 - ---------- ] . [ x ] [ ] [ 0 0 0 ] The symbolic expansion at x=0 takes half as many terms as prod(N(1,n),n,1,inf) or prod(K(k,1),k,1,inf). Numerically it wins big at x=-1 and breaks even at x=4/5, and extends the lower limit of convergence from -1 to -2-2 sqrt(2). My previous 3by3 systems came from grafting together contiguous 2by2s, e.g. the "3F2 Rosetta Stone". What's new here is targeting functions that require 3x3 to begin with. Presumably, th1s technique extends to functions more interesting than log(1-x)^2. --rwg
participants (1)
-
R. William Gosper