[math-fun] atan(tan(x)/2)
For real x, tan(x) 3 x x 1 x 1 x atan(------) = atan(tan (-)) - pi floor(-- + -) + pi floor(---- + -) + -, 2 2 pi 2 2 pi 2 2 implying that 3 2 atan(x) 1 x atan(x ) = pi floor(--------- + -) + atan(------) - atan(x), pi 2 2 1 - x and 2 x atan(x) sqrt(x + 1) - 1 3 atan(-) = ------- + atan((----------------) ) 2 2 x I was fiddling with GFee's Fourier query. A while ago I mentioned gene_salamin's observation that you could get rational function nth harmonics to shrink like n^-k, for any positive integer k, by *composing* the approximated function with one that k-smoothly (actually (k-2)-smoothly) stops and resumes at each discontinuity. The catch was large numerical coefficients delaying the onset of the faster convergence with increasing k. E.g., suppose you draw a square with epicyclic rotors. If you jump from corner to corner, you get k=1 convergence (and Gibbs ringing) and need thousands of rotors (harmonics) for screen resolution. If instead you sweep the square at constant speed, you get k=2 and need only dozens. If you come to a smooth stop at each corner, you get k=3, etc., needing fewer and fewer rotors, leading to the absurd limit that eventually you would need only one or zero! But what actually happens when you get too greedy (k too big) is the postponement of the amplitude dropoff as a function of n. So each resolution target (error bound) has its optimal k, which presumably varies with the function being approximated. But GFee didn't mention a speed-smoother option, so we're left with matching k derivatives at 0 and -1 == 1, which, I think, exhibits the same behavior. So I was looking at patch functions like 2/pi*atan(tan(t*pi/2)^(2*k+1)) for sneaking from 0 to 1 with 2k vanishing derivatives at the endpoints. Or even tanh(tan(t*pi/2)), with *all* derivatives vanishing. I don't really have time to investigate this, but the temptation may be overwhelming. Oh, hell, I gave in and tried f(t) := t 2 2 %e ((%e - 1) tanh(cot(%pi t)) + %e + 1) ------------------------------------------ 2 on [-1,0], which matches all derivatives at both ends, and got g(t) = 1.97473445780695d0 - 0.8981270448759d0 cos(pi t) - 0.41753308527721d0 sin(pi t) - 0.1247188401558d0 cos(2 pi t) + 0.23434279646825d0 sin(2 pi t) + 0.00330568635158d0 cos(3 pi t) + 0.03901681282034d0 sin(3 pi t) + 0.00107335808438d0 cos(4 pi t) + 0.04714152469382d0 sin(4 pi t) + 0.02622378383461d0 cos(5 pi t) + 0.00949607729559d0 sin(5 pi t) + 0.00589694508514d0 cos(6 pi t) + 7.34937084126405d-4 sin(6 pi t) + 0.00925330770208d0 cos(7 pi t) + 7.6772973803263d-6 sin(7 pi t) + 0.00221574684403d0 cos(8 pi t) - 0.00414748855076d0 sin(8 pi t) + 0.00175756423966d0 cos(9 pi t) - 0.00109985016356d0 sin(9 pi t) + 4.19777476400997d-4 cos(10 pi t) - 0.00251568862979d0 sin(10 pi t) - 2.97177678883817d-4 cos(11 pi t) - 6.82066527883218d-4 sin(11 pi t) which is rather disappointing: g(0) = 1.0017, g(1) = 2.7175 . Possibly the numerical integration is poor, but Maple 5.4 gave 13 place agreement. Strangely, I couldn't coax more than six places out of Mma 5.2. Maybe an atan(tan^k) hack would be better. --rwg Input forms: atan(tan(x)/2) = atan(tan(x/2)^3)-%pi*floor(x/%pi+1/2)+%pi*floor(x/(2*%pi)+1/2)+x/2 atan(x^3) = %pi*floor(2*atan(x)/%pi+1/2)+atan(x/(1-x^2))-atan(x) atan(x/2) = atan(x)/2+atan(((sqrt(x^2+1)-1)/x)^3) PS: Mathematica knows I from (-I: In[1]:= N[DedekindEta[I]] Out[1]= 0.7682254223260567 - 1.7568472086707517*^-19*I In[2]:= N[DedekindEta[-I]] Out[2]= DedekindEta[0. - 1.*I] (* Undefined. Diverges. :*) Clear thinking is much too hard. Ask any reporter. --------------------------------- Get easy, one-click access to your favorites. Make Yahoo! your homepage.
rwg>Maybe an atan(tan^k) hack would be better. Weird(?): using 2 3 %pi t 2 (%e - 1) atan(tan (-----)) t 2 %e (1 - -----------------------------) %pi instead of t 2 2 %e ((%e - 1) tanh(cot(%pi t)) + %e + 1) ------------------------------------------ 2 (matching only three derivatives) on [-1,0] gives a whole 'nother digit: f(t)=1.17652167574325d-4*sin(11*%pi*t)+6.08113365897434d-4*cos(11*%pi*t) +0.00100036999827d0*sin(10*%pi*t)-1.50249524904942d-4*cos(10*%pi*t)-3.81283020669115d-4 *sin(9*%pi*t)-0.00172079920634d0*cos(9*%pi*t)-0.00369552568607d0*sin(8*%pi*t) +9.47550861840517d-4*cos(8*%pi*t)+0.00173215817055d0*sin(7*%pi*t)+0.00784143653234d0 *cos(7*%pi*t)+0.0148972048649d0*sin(6*%pi*t)-0.00302598613176d0*cos(6*%pi*t) -0.00704134867511d0*sin(5*%pi*t)-0.02895637122193d0*cos(5*%pi*t)-0.0655859716264d0 *sin(4*%pi*t)+0.01897353641646d0*cos(4*%pi*t)+0.04756611866171d0*sin(3*%pi*t) +0.1616146418401d0*cos(3*%pi*t)+0.39334833081999d0*sin(2*%pi*t)-0.11547127986474d0 *cos(2*%pi*t)-0.39081697791208d0*sin(%pi*t)-0.99853697078804d0*cos(%pi*t) +1.95774972920477d0 [f(0.0d0) = 0.99987335148369d0, f(1.0d0) = 2.71817325043964d0] . For screen resolution, you could discard the five "d-4" terms. This seems in line with 11^-4 = 1/N^(derivs+1), but tan^5 is worse, and tan^7 worse yet. Presumably, with large enough N, tan^5 will win, and then tan^7, and eventually tanh(cot). ? --rwg --------------------------------- Be a better friend, newshound, and know-it-all with Yahoo! Mobile. Try it now.
rwg>[Four+ digit accuracy] seems in line with 11^-4 = 1/N^(derivs+1), but tan^5 is worse, and tan^7 worse yet. Presumably, with large enough N, tan^5 will win, ... Seems to be about N=50(!), where accuracy is about 10^-9. This means tan^3 is almost always what we want, especially since atan(tan(%pi*t/2)^3) =atan(tan(%pi*t)/2)+%pi*floor(t+1/2)-%pi*t/2-%pi*floor(t/2+1/2) 3 pi t tan(pi t) 1 pi t atan(tan (----)) = atan(---------) + pi floor(t + -) - ---- 2 2 2 2 t 1 - pi floor(- + -). 2 2 But we can do much better. If we only demand finitely many matching derivatives, we don't need fancy invocations of atan that only integrate numerically. For f(t) on [0,1] we can match six(!) derivatives with g(t) on [-1,0] := f(t) (f(t + 2) - f(t)) (3 cos(5pi t) - 25 cos(3pi t) + 150 cos(pi t) - 128) - ----------------------------------------------------------------------. 256 Then for f = exp, inf ==== t 6 \ k e = 225 sinh(1)pi ( > ((- 1) + e) (k pi / ==== k = 1 6 6 6 2 2 4 2 4 4 4 2 (k pi - 225pi + 7 k pi (37pi + 50pi + 5) - 777pi - 7 k pi (5pi + 3) 2 2 2 4 4 4 2 2 2 - 175pi - 7) sin(k pi t) - (7 k pi (k pi + 111pi - 5 k pi (5pi + 1) 2 2 2 2 + 50pi + 3) - (pi + 1) (9pi + 1) (25pi + 1)) cos(k pi t)) 2 2 2 2 2 2 2 2 2 2 /((k pi + 1) (k pi - 10 k pi + 25pi + 1) (k pi - 6 k pi + 9pi + 1) 2 2 2 2 2 2 2 2 (k pi - 2 k pi +pi + 1) (k pi + 2 k pi +pi + 1) 2 2 2 2 2 2 2 2 (k pi + 6 k pi + 9pi + 1) (k pi + 10 k pi + 25pi + 1)) e + 1 + ------------------------------------) 2 2 2 2 (pi + 1) (9pi + 1) (25pi + 1) Approximating 11 sin and 12 cos terms, (%e)^t ~ 9.89950946692436d-6*sin(11*pi*t)-2.24141986582817d-6*cos(11*pi*t) +4.46390665096906d-5*sin(10*pi*t)-1.14371353488594d-5*cos(10*pi*t)+4.74200179416212d-5 *sin(9*pi*t)-1.40787399330621d-5*cos(9*pi*t)+2.69674929205567d-4*sin(8*pi*t) -9.63830363878749d-5*cos(8*pi*t)+3.9887744385691d-4*sin(7*pi*t)-1.84767667695833d-4 *cos(7*pi*t)+0.00379228989443d0*sin(6*pi*t)-0.00282457989254d0*cos(6*pi*t) -0.00956883953678d0*sin(5*pi*t)-0.02111980808042d0*cos(5*pi*t)-0.07025349499746d0 *sin(4*pi*t)+0.02402343916233d0*cos(4*pi*t)+0.05218583563153d0*sin(3*pi*t) +0.18293572772707d0*cos(3*pi*t)+0.42360609603345d0*sin(2*pi*t)-0.11580061884067d0 *cos(2*pi*t)-0.38558037682327d0*sin(pi*t)-1.02075497561669d0*cos(pi*t) +1.95385367517543d0. Testing: makelist(''%,t,dfloat([0,log(3/2),log(2),1])) [1.0d0 = 1.00000395163528d0, 1.5d0 = 1.49999731997367d0, 2.0d0 = 1.99999386979907d0, 2.71828182845905d0 = 2.71828423923035d0] or ~4ppm. For screen resolution, you could probably neglect all ten "d-" terms, leaving 13. Then, it might even win to use only a 3rd harmonic "slide function" instead of the 5th we used here. Input form: %E^T = 225*%PI^6*SINH(1)*('SUM(((-1)^K+%E)*(%PI*K*(%PI^6*K^6-7*(5*%PI^2+3)*%PI^4*K^4+7*%PI^2*(37*%PI^4+50*%PI^2+5)*K^2-225*%PI^6-777*%PI^4-175*%PI^2-7)*SIN(%PI*K*T)-(7*%PI^2*K^2*(%PI^4*K^4-5*%PI^2*(5*%PI^2+1)*K^2+111*%PI^4+50*%PI^2+3)-(%PI^2+1)*(9*%PI^2+1)*(25*%PI^2+1))*COS(%PI*K*T))/((%PI^2*K^2+1)*(%PI^2*K^2-10*%PI^2*K+25*%PI^2+1)*(%PI^2*K^2-6*%PI^2*K+9*%PI^2+1)*(%PI^2*K^2-2*%PI^2*K+%PI^2+1)*(%PI^2*K^2+2*%PI^2*K+%PI^2+1)*(%PI^2*K^2+6*%PI^2*K+9*%PI^2+1)*(%PI^2*K^2+10*%PI^2*K+25*%PI^2+1)),K,1,INF)+(%E+1)/(2*(%PI^2+1)*(9*%PI^2+1)*(25*%PI^2+1))) --rwg PARASELENIC <-> LINEAR SPACE . --------------------------------- Looking for last minute shopping deals? Find them fast with Yahoo! Search.
Nine years ago, in the era of the reboot-only-for-major-quakes Lisp Machine environment (note Macsyma command line #), I sent Date: Wed, 14 Oct 1998 11:52-0700 From: Bill Gosper <rwg@[192.156.175.1]> Reply-To: rwg@NEWTON.Macsyma.COM Subject: D^k trigs To: math-fun@optima.CS.Arizona.EDU Message-ID: <19981014185209.7.RWG@[192.233.166.105]> d^k trig(x)/dx^k is just a shift by k quarter-periods for trig in {sin,cos,sinh,cosh}. For k>0, the other eight cases come out as explicit sums of <= k terms via the two cases (c3777) diff(tan(x),x,k) = (-2*sum((-1)^j*eulerian(k,j)*sin((k-2*j-1)*x-%pi*k/2), j,0,floor((k-1)/2)) -nummod(k,2)*eulerian(k,(k-1)/2))/cos(x)^(k+1); [...] Far more tasteful: 'DIFF(COT(X),X,K)=-('SUM(EULERIAN(K,J)*COS((K-2*J-1)*X),J,0,INF))/(-SIN(X))^(K+1) inf ==== \ > eulerian(k, j) cos((k - 2 j - 1) x) / k ==== d j = 0 --- (cot(x)) = - ----------------------------------------- k k + 1 dx (- sin(x)) where the eulerian terminates the series at j=k. Test: makelist(trigsimp(trigexpand(?meqhk(apply_nouns(''(subst(k,inf,%)))))),k,0,4) [0, 0, 0, 0, 0] (Now including k=0. ?meqhk subtracts halves of equations.) For eulerian(x,integer k), Macsyma currently uses the finite series k ==== \ j x eulerian(x, k) = > (- 1) (k - j + 1) binomial(x + 1, j), / ==== j = 0 E.g., eulerian(x,4) = x x x x x (x + 1) 3 (x - 1) x (x + 1) 2 5 - (x + 1) 4 + ------------ - -------------------- 2 6 (x - 2) (x - 1) x (x + 1) + ------------------------- 24 Note that term j=k+1 vanishes and for fractional x, subsequent terms are complex. If we apply the derivative formula to d^s (sum(1/(n+x)) = cot(mumble))/dx^s, for fractional s, but use the finite, real definition of Eulerian, we get complex = real. But if we take real parts, we get (within small factors) things like inf 5 ==== 3 (8 - sqrt(2)) zeta(-) \ 3 j 2 (d122) > eulerian(-, j) (- 1) = ----------------------- / 2 2 ==== 4 %pi j = 0 (c123) dfloat(apply_nouns(subst(69,inf,%))) (d123) 0.67136038306057d0 = 0.67136038786498d0 inf 1/3 1 4 ==== (2 2 - 1) (-)! zeta(-) \ 1 j 3 3 (d113) > eulerian(-, j) (- 1) = ------------------------- / 3 4/3 ==== %pi j = 0 (c114) dfloat(apply_nouns(subst(69,inf,%))) (d114) 1.06216017259015d0 = 1.06215789555839d0 This my first recollection of an algebraic series for (1/3)!. I have yet to reconcile the imaginary parts from the nonterminating complex eulerians, but something seriously interesting is going on. --rwg --------------------------------- Looking for last minute shopping deals? Find them fast with Yahoo! Search.
Bill Gosper <rwmgosper@yahoo.com> wrote: This is my first recollection of an algebraic series for (1/3)!. inf inf j ==== ==== ==== \ \ \ s n
eulerian(s, j) = > > (- n + j + 1) (- 1) binomial(s + 1, n) = s! / / / ==== ==== ==== j = 0 j = 0 n = 0
SUM(EULERIAN(S,J),J,0,INF) = S! EULERIAN(S,J) := SUM((-N+J+1)^S*(-1)^N*BINOMIAL(S+1,N),N,0,J) This is just the interpolation of n! = nth Eulerian row sum, so is unlikely to be new. It's also the b=0 and b=1 limits of the empirical (and unpolished) result, for 0<b<1, inf ==== \
eulerian(s - 1, j) cos(%pi b (s - 2 j - 2)) / ==== j = 0
sin(%pi b) s = (----------) Gamma(s) (hurwitz_zeta(s, 1 - b) cos(%pi s) + hurwitz_zeta(s, b)) %pi SUM(EULERIAN(S-1,J)*COS(%PI*B*(S-2*J-2)),J,0,INF) = (SIN(%PI*B)/%PI)^S*GAMMA(S)*(HURWITZ_ZETA(S,1-B)*COS(%PI*S)+HURWITZ_ZETA(S,B)) This lets us multisect the Eulerian row sums, e.g., n floor(-) 2 n n + 1 ==== B (1) 2 (2 - 1) \ n! n
eulerian(n, 2 j + 1) = -- - ---------------------,
/ 2 n + 1 ==== j = 0 n floor(-) 2 n n + 1 ==== B (1) 2 (2 - 1) \ n! n
eulerian(n, 2 j) = -- + ---------------------,
/ 2 n + 1 ==== j = 0 whose difference is the alternating row sum: n n n + 1 ==== 2 B (1) 2 (2 - 1) \ j n
(- 1) eulerian(n, j) = ----------------------- / n + 1 ==== j = 0
(More evidence that Bern should have been Bernpoly(1).) These multisection formulae should be of the form integer = not obviously integer. I suspect there's an analogous StirlingS2 Gamma(s) series based on the Adamchik formula per Oleg. (Extending S2 to noninteger first arg, e.g., x x x 4 - 4 3 + 6 2 + Stirling_s2(0, x) - 4 Stirling_s2(x, 4) = ---------------------------------------- 24 .) Using the infinite versions of the Eulerian and StirlingS2 series, both Gamma(noninteger) formulas ought to have interesting imaginary parts. --rwg --------------------------------- Be a better friend, newshound, and know-it-all with Yahoo! Mobile. Try it now.
Bill Gosper <rwmgosper@yahoo.com> wrote: >rwg>This my first recollection of an algebraic series for (1/3)!. Foo, the nth difference of a monic nth degree polynomial anywhere is n!, so just forget n is an integer and run the sum to oo. For oo=9, (c540) ('sum((-1)^k*(x-k)^n*binomial(n,k),k,0,9),%% = apply_nouns(%%)) 9 ==== \ k n (d540) > (- 1) binomial(n, k) (x - k) = / ==== k = 0 n n n n (n - 1) n (x - 2) (n - 2) (n - 1) n (x - 3) x - n (x - 1) + ------------------ - -------------------------- 2 6 n (n - 3) (n - 2) (n - 1) n (x - 4) + ---------------------------------- 24 n (n - 4) (n - 3) (n - 2) (n - 1) n (x - 5) - ------------------------------------------ 120 n (n - 5) (n - 4) (n - 3) (n - 2) (n - 1) n (x - 6) + -------------------------------------------------- 720 n (n - 6) (n - 5) (n - 4) (n - 3) (n - 2) (n - 1) n (x - 7) - ---------------------------------------------------------- 5040 n (n - 7) (n - 6) (n - 5) (n - 4) (n - 3) (n - 2) (n - 1) n (x - 8) + ------------------------------------------------------------------ 40320 n (n - 8) (n - 7) (n - 6) (n - 5) (n - 4) (n - 3) (n - 2) (n - 1) n (x - 9) - -------------------------------------------------------------------------- 362880 (c541) block([ratfac : true],makelist(ratsimp(''(rhs(%))),n,0,10)) (d541) [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 10 9 8 7 6 5 4 - x + 100 x - 4500 x + 120000 x - 2100000 x + 25200000 x - 210000000 x 3 2 + 1200000000 x - 4500000000 x + 10000000000 x - 9996371200] (Why not to rely *too* heavily on EIS.) (c543) subst([x = 9,n = 9/2],n! = rhs(d540)) 945 sqrt(%pi) 151263 sqrt(7) 196875 sqrt(5) (d543) ------------- = -------------- - 8505 sqrt(6) + -------------- 32 8 128 1701 sqrt(3) 4718601 sqrt(2) 640843731 - ------------ - --------------- + --------- 1024 128 32768 (c544) dfloat(%) (d544) 52.3427777845535d0 = 52.3427798776902d0 Even easier: Just construct a telescoping series from a product limit for x!: inf x + 1 x x - 1 x x ==== (----- + n) - (----- + n) (- + 1) \ 2 2 n x + 1 x (d488) > ----------------------------------- + (-----) / binomial(x + n, n) 2 ==== n = 1 (c489) apply_nouns(subst([x = 1/2.0d0,inf = 333],d488))^2*4 (d489) 3.14159177222301d0 Finally, a really simple formula for eulerian(n,j): | n + 1 n | (d547) eulerian(n, j) = delta (max (k - j, 0))| k | |k = n | n + 1 p | = delta (max (n - j, 0))|, n | |p = n where delta^j := the jth backward difference operator. (The 'at notation just is to avoid differencing the n in the exponent. There's a subtle difference between the evaluation semantics of the two forms which necessitates quotes in one or 'eval in the other, depending on your implementation of 'delta.) With nth difference vs n+1, you get the eulerian running row sum. --rwg --------------------------------- Looking for last minute shopping deals? Find them fast with Yahoo! Search.
participants (1)
-
Bill Gosper