HP 42S/DM42: Error Function, Bernoulli Numbers, Chebyshev Polynomials

HP 42S/DM42:  Error Function, Bernoulli Numbers, Chebyshev Polynomials

Error Function

The program ERF calculates the error function:

erf(x) = 2 / √π * ∫ e^-(t^2) dt from t = 0 to t = x

The calculation uses the series:

erf(x) = 2 * e^(-(x^2)) / √π * Σ( 2^n * x^(2n+1) / (1 * 3 * ... * (2n+1)) from n=0 to ∞

The error function can be used to calculate the cumulative normal distributive function from 0 to x by:

ncdf(x) = 1/2 + 1/2 * erf(x/√2)

HP 42S/DM42 Program: Error Function:  ERF

00 { 46-Byte Prgm }
01▸LBL "ERF"
02 STO 01
03 ENTER
04 ×
05 2
06 ×
07 STO 02
08 1
09 STO 03
10 RCL 01
11▸LBL 01
12 RCL 02
13 RCL 03
14 2
15 +
16 STO 03
17 ÷
18 RCL 01
19 ×
20 STO 01
21 +
22 X≠Y?
23 GTO 01
24 2
25 ×
26 PI
27 SQRT
28 RCL 02
29 2
30 ÷
31 E↑X
32 ×
33 ÷
34 RTN
35 END

Examples

erf(1.34) ≈ 0.9419
erf(2.43) ≈ 0.9994

Source:

Hewlett Packard  HP-65 Math Pac 2.   1974.

Bernoulli Numbers (Approximation)

The program BERN calculates Bernoulli numbers for positive integers n > 2.   The calculator is set to Radian mode.

The formula used:

B_n = (-cos(n * π/2))^((n+2)/2) * 2 * n! * Σ(2 * j * π)^(-n) from j=1 to ∞

HP 42S/DM42 Program: Bernoulli Numbers:  BERN

00 { 70-Byte Prgm }
01▸LBL "BERN"
02 RAD
03 STO 00
04 CLX
05 STO 02
06 1
07 +
08 STO 01
09▸LBL 00
10 2
11 PI
12 ×
13 RCL× 01
14 RCL 00
15 +/-
16 Y↑X
17 STO+ 02
18 1
19 STO+ 01
20 R↓
21 ABS
22 1ᴇ-10
23 X < Y ? 
24 GTO 00
25 RCL 02
26 RCL 00
27 N!
28 ×
29 2
30 ×
31 RCL 00
32 PI
33 ×
34 2
35 ÷
36 COS
37 +/-
38 2
39 RCL+ 00
40 2
41 ÷
42 Y↑X
43 ×
44 RTN
45 END

Examples

B_2 ≈ 0.1667
B_3 returns 1.1184E-85  (0)
B_4 returns -0.0333
B_6 returns 0.0238

Source:

Keith Oldhma, Jan Myland, and Jerome Spanier An Atlas of Functions 2nd Ed.  Springer: New York.  ISBN: 978-0-387-48806-6

Chebyshev Polynomials

The program CHEBY calculates the value of a first-order Chebyshev polynomial using the recurring formula where:

T_n(x) = 2 * x * T_n(x) - T_(n-1)(x)

where T_0(x) = 1,  T_1(x) = x

HP 42S/DM42 Program:  Chebyshev Polynomials: CHEBY

00 { 51-Byte Prgm }
01▸LBL "CHEBY"
02 STO 02
03 2
04 ×
05 STO 03
06 R↓
07 STO 00
08 X=0?
09 GTO 00
10 1
11 STO 01
12 STO- 00
13 X=Y?
14 GTO 01
15▸LBL 02
16 RCL 03
17 RCL× 02
18 RCL- 01
19 RCL 02
20 STO 01
21 R↓
22 STO 02
23 DSE 00
24 GTO 02
25 RTN
26▸LBL 00
27 1
28 RTN
29▸LBL 01
30 RCL 02
31 RTN
32 .END.

Instructions

Input n,  [ENTER], input x, XEQ CHEBY

Examples
T_3(0.4):  3, ENTER, 0.4, XEQ CHEBY
Result:  -0.9440

T_4(1.5):  4, ENTER, 1.5, XEQ CHEBY
Result:  23.5000

Source

Hewlett Packard  HP-65 Math Pac 2.   1974.



Happy programming,

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Casio fx-5800p Special Functions

Casio fx-5800p Special Functions

Programs

An Alternate Way of Extracting the Fraction and Integer Parts of a Number

The fraction part is stored in F, and the integer part is stored in I.  This algorithm can be used when a calculator or programming language does not have a fractional part or integer part function.

This program assumes the program is in Radians mode.

Casio fx-5800P Program FPIP

Rad

0.5    → F

? → X

cos(πX) = 0 ⇒ Goto 1

Abs(X) → F

tan^-1 (tan (πF)) ÷ π → F

X < 0 ⇒ -F → F

Lbl 1

F /right-triangle ([Shift] [x^2])

X – F → I

Bernoulli Numbers

The program BERNOULLI approximates the Bernoulli number of nth order.

Casio fx-5800P Program BERNOULLI

? → N

If N = 0

Then

1 → B

Goto 1

IfEnd

If Frac(N ÷ 2) ≠ 0

Then

Int( Abs(N – 2) ÷ (N – 2) – 1) ÷ 4 → B

Goto 1

IfEnd

Σ ( (2*J*π)^(-N), J, 1, 250) * (-1)^(N ÷ 2 + 1) * 2 * N! → B

Lbl 1

B

Euler Numbers

The program EULERNUM calculates the Euler number of order n.

Casio fx-5800P Program EULERNUM

? → N

0 → E

Frac(N ÷ 2) ≠ 0 ⇒ Goto 1

(2 ÷ π)^(N + 1) * N! ÷ 5 → E

-1 – 10 * Int(E) → E

Frac(N ÷ 4) ≠ 0 ⇒ Goto 1

4 – E → E

N ≠ 0 ⇒ Goto 1

1   → E

Lbl 1

E

Custom Formulas

How to create custom formulas:  [MODE], 5.  PROG, 1.  NEW, give the name, 3. Formula

To calculate the custom formula, press [CALC], enter the value for T (ignore X because X is a dummy variable).

Sine Integral:   SI:  ∫( sin(X)/X dX,0,T)

Dawson Integral:  DAWSON:  ∫(e^(T^2-X^2) dX,0,T)

Bessel Function, with order N:  BESSEL:  1/π * ∫(cos(T * sin(X)-N*X) dX,0,T)

Error Function:  ERF:  2/√π * ∫(e^(-X^2) dX,0,T)

Fresnel Sine:  FRESSIN:  ∫( sin(π*X^2/2) dX,0,T)

Fresnel Cosine: FRESCOS: ∫( cos(π*X^2/2) dX,0,T)

Source:

Jerome Spainer and Keith B. Oldham  An Atlas of Functions Hemisphere Publication Corporation: Washington  1987  ISBN 0-89116-5738-8

Eddie

This blog is property of Edward Shore, 2017

Bernoulli Numbers and Polynomials

Today is a good day. I went to the library of alma mater, Cal Poly Pomona for library day. Once a month I try to visit a university library and peruse through with Mathematics Section.

I finally got the concept of generating functions. In the many years I studied math, generating functions proved to be an elusive topic for me. Not any more.

What are generating functions?

Generating Functions

Generating functions is a power series. The series is usually does not terminate. The coefficients of the power series can reveal a sequence used in various fields in science.

The Ordinary Generating Function:

G(a_n; x^n) = ∑ a_n * x^n

Sometimes we have an Exponential Generating Function:

E(a_n; x^n) = ∑ a_n * (x^n/n!)

There are other types of generating functions, but I will focus on the basic types described above.

---

n! represents the factorial of n. In this case, n is a positive integer, and:

n! = n * (n-1) * (n-2) * ... * 3 * 2 * 1

By definition, 0! = 1

---

To expand the generating function G(a_n; x^n), calculate a Maclaurin Series of G(x). A Maclaurin Series of a function is a Taylor Series about the point x = 0. Hence the Maclaurin Series:

G(x) = G(0) + G'(0) * x/1! + G''(0) * x^2/2! + G'''(0) * x^3/3! + .... + O(x)

where O(x) is the error term G^(n+1)(z) * x^(n+1)/(n+1)! In practice and calculation, O(x) is sometimes "ignored".

You can find more information about generating functions here.

An Example: A Basic Generating Function

Expand the generating function, let's go five terms:

G(a_n; x^n) = 1/(1 - x)

G(x) = 1/(1 - x)
Then:
G(0) = 1/(1 - 0) = 1

dG/dx = 1/(x-1)^2; dG/dx(0) = 1

d^2G/dx^2 = -2/(x-1)^3; d^2G/dx^2(0) = 2

d^3G/dx^3 = 6/(x-1)^4; d^3G/dx^3(0) = 6

d^4G/dx^4 = -24/(x-1)^5; d^4G/dx^4(0) = 24

To five terms...

1/(1-x) = 1 + 1 * x/1! + 2 * x^2/2! + 6 *x^3/3! + 24 * x^4/4! + O(x)
= 1 + x + x^2 + x^3 + x^4 + O(x)

The sequence of coefficients are: {1, 1, 1, 1, 1}.

The Bernoulli Numbers

The Bernoulli Numbers can be found by using the Exponential Generating Function:

E(a_n; x) = ∑ x/(e^x - 1)

The Bernoulli Numbers, B_n, are the "coefficients" of the expansion of ∑ x/(e^x - 1). Recall in the Exponential Generating Function, the "coefficients" are of the term x^n/n!.

With E(x) = x/(e^x - 1), the first two derivatives are:

dE/dx = - ((x-1)*e^x + 1)/(e^(2x) - 2*e^x + 1)

d^2E/dx^2 =
((x - 2)*e^(2x) + (x + 2)*e^x) / (e^(3x) - 3*e^(2x) + 3*e^(2x) - 1)

Note that calculating E(0) gives us 0/0. Same for dE/dx(0) and d^2E/dx^2(0). However, if I use a calculate with CAS capabilities such as the Hewlett Packard HP 50g, or mathematical software such as MathStudio (an app for smartphones and iPads), I get something like this (first eight terms):

Why is that?

Observe that:

lim x/(e^x - 1) as x → 0 = 0/0

Using the L'Hospital's rule, we can take the derivatives of both numerator and denominator,

lim 1/(e^x) as x → 1 = 1

which implies that:

lim x/(e^x - 1) as x → 0 = 1

You can generate terms by taking the limit as x → 0 for each term.

This how we end up with the series. Now to extract the "coefficients", observe that:

12 = 6 * 2!
(no term contains x^3/3!)
720 = 30 * 4!
(no term contains x^5/5!)
30240 = 42 * 6!
(no term contains x^7/7!)
1,209,600 = 30 * 8!

Our sequence for this generating function (for nine terms) is:

{1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30}

These numbers are the Bernoulli numbers. In fact, definition by generating function is:

x/(e^x - 1) = ∑ B_n * (x^n/n!)

---

Bernoulli Numbers

B_0 = 1
B_1 = -1
B_2 = 1/6
B_3 = 0
B_4 = -1/30
B_6 = 1/42
B_8 = -1/30
B_10 = 5/66
B_12 = -691/2730
B_14 = 7/6

B_n = 0 where n is odd and n > 2

---

Bernoulli Polynomials

Bernoulli Polynomials can be generated by the following formula:

β_n(x) = ∑((B_k * n!/(k!*(n-k)!) * x^k, from k = 0 to n)

where B_k is the kth Bernoulli number.

Source:
Krylov, Vladimir Ivanoch, translated by H. Stroud Approximate Calculations of Integrals McMillian Company: New York, 1962

Until next time, be safe everyone! Eddie



This blog is property of Edward Shore, 2012.