Clausen's Integral for the HP-41
Overview
-The following program computes f(x) = - §0x Ln ( 2 sin (t/2) ) dt ( § = Integral )
-Integration by parts gives: f(x) = -x Ln ( 2 sin (x/2)
) + 2 §0x/2 ( t / Tan t ) dt
-The series expansion t / Tan t = 1 - x2/3 -
x4/45 - 2 x6/945 - ..... - [ (-1)n-1 22n/(2n)!
] B2n x2n - ..................... is used
where B2n are Bernoulli numbers.
-Bernoulli numbers B2n are calculated through Zeta
function for n > 3
Program Listing
Data Registers: R00-R01 are used by "ZETA"
R02 = x R03 thru R06: temp ( R05 =
§0x/2 ( t / Tan t ) dt )
Flags: /
Subroutine: "ZETA" ( cf "Riemann
Zeta Function for the HP-41" )
-If x is small enough ( for example x = 0.15 ), "ZETA"
is not needed.
| 01 LBL "CLAUS"
02 DEG 03 STO 02 04 2 05 / 06 ENTER^ 07 X^2 08 STO 04 09 LASTX 10 * 11 STO 03 12 9 13 / 14 - 15 RCL 03 16 RCL 04 |
17 *
18 STO 03 19 225 20 / 21 - 22 RCL 03 23 RCL 04 24 * 25 STO 03 26 3307.5 27 / 28 CHS 29 X<>Y 30 + 31 STO 05 32 LASTX |
33 X=Y?
34 GTO 00 35 2 36 PI 37 6 38 STO 06 39 Y^X 40 / 41 CHS 42 ST* 03 43 FIX 9 44 LBL 01 45 RCL 06 46 2 47 + 48 STO 06 |
49 XEQ "ZETA"
50 RCL 03 51 RCL 04 52 * 53 PI 54 X^2 55 / 56 STO 03 57 * 58 RCL 06 59 1 60 + 61 / 62 RCL 05 63 + 64 STO 05 |
65 LASTX
66 X#Y? 67 GTO 01 68 LBL 00 69 ST+ X 70 RCL 02 71 2 72 / 73 R-D 74 SIN 75 ST+ X 76 LN 77 RCL 02 78 * 79 - 80 END |
( 108 bytes / SIZE 007 )
| STACK | INPUTS | OUTPUTS |
| Y | / | g(x) |
| X | x | f(x) |
where f(x) = Clausen's Integral
( 0 < x < 2 PI )
and g(x) = §0x/2
( t / Tan t ) dt
Examples:
PI 3 / XEQ "CLAUS" >>>>
f(PI/3) = 1.014941606
---Execution time = 28s---
X<>Y g(PI/3) = 0.507470803
0.15 R/S >>>>
f(0.15) = 0.434614878
( in 2 seconds )
X<>Y g(0.15) = 0.074953114
6 R/S
>>>> f(6) = -0.640782664
( in 7mn49s )
X<>Y g(6) = -4.115383667
Notes:
-Execution time tends to infinity as x tends to 2 PI
-This routine produces DATA ERROR if x = 0 but since f(0) = 0, simply
add X=0? RTN after line 01
-Clausen's integral is also equal to SUMk=1,2,......
[ Sin (kx) ] / k2 but thousands of terms would be required to
return accurate results if we used this formula.
-However, Clausen's integral Cl(x) is closely related to Lobachevsky
function L(x), we have:
Cl(x) = 2 L(x/2) at least if x >= 0 & x <= 2.PI
-So, we can use the program "LOB" listed in "Lobachevsky Function for
the HP-41"
Reference:
[1] Abramowitz and Stegun - "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4