(28/48/49/50) Lerch Transcendent and Polylogarithm

[John Keith]

This program computes the Lerch transcendent Phi(z, s, a).The arguments may be real or complex. The program is a straightforward implementation of the function definition and it converges fairly quickly if |z| is not close to 1. The execution time rises and accuracy decreases as |z| approaches 1, so this is not an efficient algorithm for the Zeta function. The number 9999 for the maximum number of iterations may be lowered to improve speed at the cost of accuracy for values of |z| close to 1.

The polylogarithm Li_s(z) is defined as z*Phi(z, s, 1). If the program is named LERCH, with s on level 2 and z on level 1, the polylogarithm can be implemented as DUP ROT 1 LERCH *.

Several other useful identities are listed at the Wikipedia links above, and at the Mathworld pages for the Lerch transcendent and the polylogarithm.

Code:


\<< \-> z s a
  \<< a s ^ INV 1. 9999.
    FOR k DUP k z OVER ^ SWAP a + s ^ / + DUP ROT
      IF SAME
      THEN 9999. 'k' STO
      END
    NEXT
  \>>
\>>

[C.Ret]

Hi,

I confirmed that this code runs on a HP-28S.

Here is a simplification, the long IF SAME THEN ... END sequence is a way to leave the FOR / NEXT loop when convergence is achieved. A more academic way would have been to use a DO /UNTIL or WHILE /REPEAT loop structure. But this also add many instructions for increasing the k-indices.

The smarter way is to kept the FOR/NEXT structure in charge of increasing k indices number but with the flexibility of an DO/UNTIL or WHILE/REPEAT non deterministic loop using a tricky IFTE STEP instance:

Code:

« →  z  s  a
   «  a  s  ^  INV                             // Initialize      L( 0 ) = 1/a^s
     1. 9999. FOR k
          DUP                                  // Leave a copy of L(k-1) in the stack
          z  k  ^  k  a  +  s  ^  /  +         // Compute         L( k )= L(k-1) + k^z/(k+a)^s
          DUP  ROT  SAME  MAXR  1.  IFTE       // Test            L(k-1)==L(k)
     STEP » »                                  // Step k by 1 until  L(k-1) is numerically the SAME as L(k), then jump to infinity

By the way, I also remove OVER or SWAP stack operation in the computation. That why the local variables are for; avoid unnecessary obscure stack manipulations from any formulae.

Num.App.:
From WolframAlpha: \( \varphi (0.6\:,3\:,1)= {\sum_{k=0}^{\infty }\frac{0.6^k}{(k+1)^3}}\simeq 1.0933375227216344720577686546855537153004593302300056675987 \) I found \( 1.0933752271 \) on my HP-28S in 0'03"68.


Please to share with you, I hope this details will help you in your number crunching applications.

Sincerely.
C.Ret


EDIT: Correct a broken english style.

[John Keith]

Thanks for sharing your ideas. I use the stack rather than local variables because stack operations are (usually) faster. Your idea for exiting the FOR loop is interesting, I will check it out.