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.