TI-84 Plus Fast Fourier Transform

TI-84 Plus Fast Fourier Transform

HAPPY NEW YEAR!

Introduction

The program FFT1 performs the fast Fourier transform of discrete data points named in List 1 (small x, signal at time points) to  List 2 (big X, frequency), using the formula:

X_k = ∑( x_n * e^(-i*2*π*k*m)/n from m = 0 to n - 1)

For the set of n signals.

The program IFFT1, the Inverse Fast Fourier Transform, reverses the process (big X to small x, List 2 to List 1).

x_k = 1/n * ∑( x_n * e^(i*2*π*k*m)/n from m = 0 to n - 1)


TI-84 Plus Program FFT1

"FFT VERSION 1"
"EWS 2018-12-09"
Input "SMALL X LIST:",L₁
"SETUP"
dim(L₁)→N
L₁→L₂
a+bi
Radian
"LOOP"
For(M,0,N-1)
0→T
For(K,0,N-1)
T+L₁(K+1)*e^(-­2*π*i*K*M/N)→T
End
T→L₂(M+1)
End
Pause L₂

TI-84 Plus Program IFFT1

"IFFT VERSION 1"
"EWS 2018-12-09"
Input "BIG X LIST:",L₂
"SETUP"
dim(L₂)→N
L₂→L₁
a+bi
Radian
"LOOP"
For(M,0,N-1)
0→T
For(K,0,N-1)
T+L₂(K+1)*e^(2*π*i*K*M/N)→T
End
T/N→T
T→L₁(M+1)
End
Pause L₁

Example 

FFT Example:
L1 = {2i, 1+i, 3}
Result (Fix 4):   {4.0000+3.0000*i, -1.1340+3.2321*i, -2.8660-0.2321*i}

Inverse FFT Example:
L2 = {0.5, -0.7i, 0.9+0.3i}
Result (Fix 4): {0.4667-0.1333*i, 0.3053-0.1931*i, -0.2720+0.3265*i}

Trigonometric Version

For calculators that do not handle complex numbers, here are the sample code that can handle FFT and IFFT.  The techinque here is to use a pair of lists, one for the real parts, the other for imaginary parts.  Example:  {2, 3+3i, -2i} get split into {2, 3, 0} and {0, 3, -2}.

The following identities and calculation are used:

e^(i*θ) = cos θ + i * sin θ

(a + b*i)*e^(-i*θ)
=  (a + b*i) * (cos θ - i * sin θ)
=  a * cos θ +  i * b * cos θ - i * a * sin θ + b * sin θ
= (a * cos θ + b * sin θ) + i * (b * cos θ - a * sin θ)

(a + b*i)*e^(i*θ)
=  (a + b*i) * (cos θ + i * sin θ)
=  a * cos θ +  i * b * cos θ + i * a * sin θ - b * sin θ
= (a * cos θ - b * sin θ) + i * (b * cos θ + a * sin θ)

where θ = 2*π*m*k/n 

TI-84 Plus Program FFT1

"FFT VERSION 2"

"EWS 2018-12-10"

Disp "SMALL X","L₁ + I*L₂"

Input "REAL: ",L₁

Input "IMAG: ",L₂

"SETUP"

dim(L₁)→N

L₁→L₃

L₂→L₄

Radian

"LOOP"

For(M,0,N-1)

0→S

0→T

For(K,0,N-1)

2*π*M*K/N→θ

S+L₁(K+1)*cos(θ)+L₂(K+1)*sin(θ)→S

T+L₂(K+1)*cos(θ)-L₁(K+1)*sin(θ)→T

End

S→L₃(M+1)

T→L₄(M+1)

End

Disp "REAL = L₃","IMAG = L₄"

TI-84 Plus Program IFFT1

"IFFT VERSION 2"

"EWS 2018-12-10"

Disp "BIG X","L₃ + I*L₄"

Input "REAL: ",L₃

Input "IMAG: ",L₄

"SETUP"

dim(L₃)→N

L₃→L₁

L₄→L₂

Radian

"LOOP"

For(M,0,N-1)

0→S

0→T

For(K,0,N-1)

2*π*M*K/N→θ

S+L₃(K+1)*cos(θ)-L₄(K+1)*sin(θ)→S

T+L₄(K+1)*cos(θ)+L₃(K+1)*sin(θ)→T

End

S/N→S

T/N→T

S→L₁(M+1)

T→L₂(M+1)

End

Disp "REAL = L₁","IMAG = L₂"

Source:  

"An Interactive Guide To The Fourier Transform"  Better Explained.  Retrieved December 9, 2018.  https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

HP 42S/DM 42/Free42: Argument, Real Part, Imaginary Part, and Sign Function

HP 42S/DM 42/Free42: Argument, Real Part, Imaginary Part, and Sign Function

Introduction

The HP 42S* has only one extraction function for complex numbers: ABS (absolute value) by default.  To complete the list, here are some programs for ARG (argument), REAL (real part), and IMAG (imaginary part).  The three functions take a complex number, in either rectangular or polar form. 

Using stack commands and flag checks, I am able to get obtain results without affecting the stack much (although the z and t stacks will have the same value).

Also included is the SIGN (signum) function.

* This applies to the Swiss Micros DM42 (it should work), Free42 from Thomas Oakken, and any other HP 42S emulator apps.  The print out is from Free42 (and I also have a physical HP 42S).

HP 42S Program: ARG

00 { 19-Byte Prgm }
01▸LBL "ARG"
02 COMPLEX
03 X<>Y
04 FC? 73
05 →POL
06 R↑
07 STO ST Y
08 R↓
09 R↓
10 RTN
11 .END.

HP 42S Program:  REAL

00 { 21-Byte Prgm }
01▸LBL "REAL"
02 COMPLEX
03 X<>Y
04 FS? 73
05 →REC
06 X<>Y
07 R↑
08 STO ST Y
09 R↓
10 R↓
11 RTN
12 END

HP 42S Program: IMAG

00 { 20-Byte Prgm }
01▸LBL "IMAG"
02 COMPLEX
03 X<>Y
04 FS? 73
05 →REC
06 R↑
07 STO ST Y
08 R↓
09 R↓
10 RTN
11 END

Example - Rectangular Mode:

7+8i   (7 [ENTER] 8 [ (shift) ] (COMPLEX) )
XEQ ARG:  48.8141 (Degrees)
XEQ REAL: 7.0000
XEQ IMAG: 8.0000

Example - Polar Mode: (Degrees):

4 ∡ 60  (4 [ENTER] 60 [ (shift) ] (COMPLEX) )
XEQ ARG:  60.0000
XEQ REAL:  2.0000
XEQ IMAG:  3.4641

HP 42S Program:  SIGN

00 { 15-Byte Prgm }
01▸LBL "SIGN"
02 X=0?
03 RTN
04 ENTER
05 ABS
06 X<>Y
07 ÷
08 RTN
09 END

Example:

-1.4641  XEQ SIGN returns -1
5.525 XEQ SIGN returns 1
0 XEQ SIGN returns 0

Eddie

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

Casio fx-5800p and HP Prime: American Wire Gauge (AWB)

Casio fx-5800p and HP Prime:  American Wire Gauge (AWB)

Everyone, I wish you all a Happy Holiday season, Merry Christmas and Happy New Year!  I'm thankful for all of you.  Wishing you the best!  - Eddie

Introduction

The following formula calculates the diameter of an American Wire Gauge:

d = 460 / (92^((3 + g)/39))  (in mils)

where g is the wire number (integer).

Define g as:
Wire number:  0000000,  g = -6
Wire number:  000000,  g = -5
Wire number:  00000, g = -4
Wire number:  0000, g = -3
Wire number:  000, g = -2
Wire number:  00, g = -1
Wire number:  0, g = 0
Any other wire number, enter g as the wire number (1, 2, 3, 4, etc.)

Conversion factor:  1000 mils = 1 in

Cross area of the wire:

a = d^2

Casio fx-5800p Program WIRES

Given the wire number, the program WIRES returns the diameter of the wire in inches and the cross sectional area in square inches.

"WIRE NUMBER: "? → G
460 ÷ (92 ^((3 + G) ÷ 39)) → D
D ÷ 1000 → D
D^2 → A
"DIAMETER: (IN)" ⊿
D ⊿
"CROSS AREA:  (IN^2)" ⊿
A 

HP Prime Program Function AWB

Given the wire number, the function AWB returns the diameter of the wire in inches. 

EXPORT AWG(g)
BEGIN
// American Wire Gauge Diameter Function
// 2018-12-23 EWS
// in inches
RETURN (460/(92^((3+g)/39)))/1000;
END;

Examples:

g = -2  (000);  d ≈ 0.4096 in, a ≈ 0.1678 in^2
g = 3; d ≈ 0.2294 in, a ≈ 0.0526 in^2
g = 6; d ≈ 0.1620 in, a ≈ 0.0263 in^2


Sources: 

Ball, John A.  Algorithms For RPN Calculators  John Wiley & Sons: New York  1978.  ISBN 0-471-03070-8

Glover, Thomas J.  Pocket Ref 4th Edition. Sequoia Publishing, Inc.:  Littleton, CO  2012.  ISBN 978-1-885071-62-0

Eddie

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

HP Prime and TI-84 Plus: Fibonacci Triangles (updated)

HP Prime and TI-84 Plus:  Fibonacci Triangles

This program is the request of John Cvetan.  I thank you for your suggestion.

Introduction

The program FIBMAT generates the Fibonacci Triangle in matrix form.  The Fibonacci triangle is a triangle generated where the outer entries of each row contain the Fibonacci sequence.  The Fibonacci sequence is generated by:

f_0 = 1
f_ 1 = 2
f_n = f_n-1 + f_n-2

You can quickly calculate the nth Fibonacci number by the formula:

f_n = ( (1 + √5)^n - (1 - √5)^n ) / (2^n * √5)

To generate the Fibonacci triangle,

1.  Let r the row and c be the column with

f_0,0 =  1
f_1,0 = 1
f_1,1 = 1
f_2,1 = 1

2.  Each row will be determined by adding the last two terms going diagonally.  You can use one of two formulas:

f_r,c = f_r-1,c + f_r-2,c

f_r,c = f_r-1,c-1 + f_r-1,c-2

The Program FIBMAT

FIBMAT generates a Fibonacci triangle in matrix form.  It's the result is a triangle that is "tilted".  n will need to be 3 or greater.

HP Prime Program FIBMAT

EXPORT FIBMAT(n)
BEGIN
// Fibonacci "triangle" in
// matrix form
// 2018-12-17 EWS
LOCAL M1,k;
M1:=MAKEMAT(0,n+1,n+1);
M1(1,1):=1;
M1(2,1):=1;
M1(2,2):=1;
FOR k FROM 3 TO n+1 DO
M1(k):=row(M1,k-1)+row(M1,k-2);
M1(k,k):=M1(k-1,k-1)+M1(k-2,k-2);
END;
RETURN M1;
END;

Here's an alternate code for Fibonacci Triangle Matrices.  The row command has been eliminated and I used a second For loop in its place.  (added 12/23/2018)

EXPORT FIBMATALT(n)
BEGIN
// 2018-12-23 EWS
// Fibonacci Matrix Alternate
// matrix form
// This version does not have the
// row function.
LOCAL M1,k,j;
M1:=MAKEMAT(0,n+1,n+1);
M1(1,1):=1;
M1(2,1):=1;
M1(2,2):=1;
FOR k FROM 3 TO n+1 DO
FOR j FROM 1 TO n DO
M1(k,j):=M1(k-2,j)+M1(k-1,j);
END;
M1(k,k):=M1(k-1,k-1)+M1(k-2,k-2);
END;
RETURN M1;
END;


TI-84 Plus Program FIBMAT

"2018-12-18 EWS"
"FIBONACCI MATRIX"
Input "ORDER: ",N
{N+1,N+1}→dim([A])
1→[A](1,1)
1→[A](2,1)
1→[A](2,2)
For(K,3,N+1)
For(J,1,N)
[A](K-2,J)+[A](K-1,J)→[A](K,J)
End
[A](K-1,K-1)+[A](K-2,K-2)→[A](K,K)
End
Pause [A]

The Program FIBTRI

This is a visual program for Fibonacci Triangle.

FIBTRI(n) generates a visual Fibonacci Triangle - although I don't recommend going beyond 12 rows due to the constraints of the screen.  I used the small font for the rows.

HP Prime Program FIBTRI

EXPORT FIBTRI(n)
BEGIN
// Fibonacci triangle
// 2018-12-17 EWS
LOCAL M1,k;
M1:=MAKEMAT(0,n+1,n+1);
M1(1,1):=1;
M1(2,1):=1;
M1(2,2):=1;
FOR k FROM 3 TO n+1 DO
M1(k):=row(M1,k-1)+row(M1,k-2);
M1(k,k):=M1(k-1,k-1)+M1(k-2,k-2);
END;

RECT();
LOCAL s;
FOR k FROM 1 TO n+1 DO
s:=STRING(SUB(row(M1,k),1,k));

IF k≤6 THEN
TEXTOUT_P(s,
140-5.5*(k-1),(k-1)*15,2);
END;

IF k>6 AND k≤11 THEN
TEXTOUT_P(s,
140-8*(k-1),(k-1)*15,2);
END;

IF k>11 THEN
TEXTOUT_P(s,
140-11.5*(k-1),(k-1)*15,2);
END;

END;
WAIT(0);
END;

Source:

Hosoya, Haruo.  "Fibonacci Triangle"  Ochanomizu University, Tokyo, Japan.  1976.  https://www.fq.math.ca/Scanned/14-2/hosoya.pdf

Eddie

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

HP Prime and TI-84 Plus: Sieve of Eratosthenes

HP Prime and TI-84 Plus: Sieve of Eratosthenes

Introduction

The program SIEVEMIN shows a miniature version of the famous Sieve of Eratosthenes.  The Sieve of Eratosthenes is a Greek algorithm that determines prime numbers from 2 to N through eliminating multiplies.

Algorithm:

1.  List all integers from 1 to n.  Of course, you choose what n is. 

2.   1 is not a prime number, hence eliminate 1. 

3.  Start with 2.  Eliminate multiple of 2, but not 2 itself.

4.  Then do step 3, but with multiples of 3, but not 3 itself.

5.  Then do step 5, but with multiples of 5, but not 5 itself. 

6.  Continue on with the next available number until there are no multiples remaining.  At the end, you have all the prime numbers from 1 to n.

Most sieves will demonstrate the algorithm from 1 to 100.  However, to fit the calculator screens, SIEVEMIN finds all the prime numbers from 1 to 49. 

HP Prime Program: SIEVEMIN 

sub1(); // subroutine

EXPORT SIEVEMIN()
BEGIN
// 2018-12-18 EWS
// Mini Sieve
HFormat:=0; // standard
PRINT();
PRINT("Mini Sieve");
WAIT(1);
LOCAL M0,T,R,C,K;
T:=1;
M0:=MAKEMAT(0,7,7);

FOR R FROM 1 TO 7 DO
FOR C FROM 1 TO 7 DO
M0(R,C):=T;
T:=1+T;
END;
END;

sub1(M0);

M0(1,1):=0;
FOR K FROM 2 TO 7 DO
FOR R FROM 1 TO 7 DO
FOR C FROM 1 TO 7 DO
IF FP(M0(R,C)/K)==0 AND 
M0(R,C)>K THEN
M0(R,C):=0;
sub1(M0);
END;
END;
END;
END;

// end of program

END;

sub1(M1)
BEGIN
PRINT();
LOCAL I;
FOR I FROM 1 TO 7 DO
PRINT(row(M1,I));
END;
WAIT(0.5);
END;

TI-84 Plus CE Program:  SIEVEMIN

"2018-12-11 EWS"
Disp "MINI SIEVE OF","ERATOSHENES","TI-84 PLUS CE"
Float
Wait 1
{7,7}→dim([J])
1→T
For(R,1,7)
For(C,1,7)
T→[J](R,C)
1+T→T
End
End
ClrHome
Disp [J]
Wait 1

0→[J](1,1)
For(K,2,7)
For(R,1,7)
For(C,1,7)
If fPart([J](R,C)/K)=0 and [J](R,C)>K
Then
0→[J](R,C)
ClrHome
Disp [J]
Wait .5
End
End
End
End

Eddie

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

Retro Review: Texas Instruments Explorer Plus

Retro Review:  Texas Instruments Explorer Plus

General Information

Company:  Texas Instruments
Type:  Scientific - Algebraic
Display:  10 digits or 5 digits with 2 digit exponents
Power:  Solar
Memory:  1
Years:  1995 - 1997 (two editions).  I have a the 2nd edition.
Original Cost: $24.95
Documentation:  Manual
Datamath page:  http://www.datamath.org/Sci/Modern/TI-32EXPLUS_2.htm

Features
* Trigonometry, Logarithms
* Degree and Radian mode
* 1 Variable Statistics
* Backspace key
* Fractions with simplification
* 2 Stored Operations (Macros)
* Random Numbers
* Factorial, Combinations (nCr), Permutations (nPr)

Fractions

The Explorer Plus dedicates six keys to fractions. 

[Unit]:  Enter the unit portion of a mixed fraction.  A small "u" separates the whole unit and its fraction.

[ / ]:  Use this to separate a fraction's numerator and denominator.  This is not the same as the division key [ ÷ ].

[ F <> D]:  Change a result from fraction to decimal and vice versa.  Changing to a fraction is does not always change the decimal to the fraction

[ Simp ]:  Use this key to simplify fractions.  I'll talk about this in a bit.

[ A b/c ]:  Change an improper fraction to a mixed fraction.

[ x<> y ]:  Exchange entries.  This works for operations outside the fraction commands. 


Changing to a mixed fraction with the [ A b/c ] key

Example:  Change the improper fraction 74/11 to a mixed fraction.

Keystrokes:  74 [ / ] 11 [ A b/c ]

Result:  6 u 8/11

Simplifying a Fraction

If a result has the display "N/D → n/d" on the left side, this means the fraction can be simplified.  Fractions can be fully simplified, or you can provide a factor. 

Example:  Fully simplify 88/55.

Keystrokes: 88 [ / ]  55 [ = ] [Simp] [ = ]

Result: 8 / 5

Press [ x<>y ] to get the factor used to reduce the fraction, in this case, it would be 11.

Example:  Reduce the fraction 284/404 by a factor of 2.

Keystrokes:  284 [ / ] 404 [ = ] [Simp] 2 [ = ]

Result:  142 / 202

π and Fractions

In Radians Mode, you can work with fractions involving π.  The display will show "Pi". 

Example:  π/4 + 3π/6

Keystrokes:  [ 2nd ] [ π ] (DR>) until RAD shows in the display

[ π ] [ / ] 4        Display:  Pi/4
[ + ] 3 [ π ] [ / ] 6      Display:  3Pi/6
[ = ]     Display:  9Pi/12

Result: 9π/12

For the decimal answer, press [ F<>D ] to get 2.35619449.

Other Keys of Interest

Integer Division

The integer division [ INT÷ ] will display the quotient and remainder of the division problems of two integers.  The displays the quotient by Q and remainder by R.  If you use the result in future calculations, only the integer portion is carried over. 

Example:  Use the integer division command on 128 ÷ 35

Keystrokes:  128 [ INT÷ ] 35 [ = ]

Result:  Q:  3,  R:  23


Random Numbers

The RAND command on the Explorer Plus has two functions:

RAND by itself:  a three digit number between 0.000 and 0.999.  Press [ 2nd ] [ = ] (RAND)

Integer Random Number:  You can get a random integer from 1 to n (up to 100) by the this keystroke sequence:  1 [ 2nd ] [ = ] (RAND) n [ = ].

Example:  Generate a random integer from 1 to 50.

Keystrokes:  1 [ 2nd ] [ = ]  (RAND) 50 [ = ]

Results may vary.


 Stored Operations

You can store simple operations in one of two OP registers, OP1 and OP2.  The operations that can be included are +, -, *, ÷, y^x, INT÷, y^(1/x), nCr, nPr, and random numbers.

Example:  Store x + 6 in OP1. 

[CE/C] (or suitable number) [ + ] 6 [ OP1 ]

11 + 6, 15 + 6, 28 + 6

11 [ OP1 ]  Display:  1        17
15 [ OP1 ]  Display:  1       21
28  [ OP1 ] Display:  1      34

You can repeat operations. 

Clear stored operations by either pressing [ON/AC], which will clear everything, or pressing [2nd] (

Verdict

I really like the the Explorer Plus and I think it would have a made a great base for the current TI-30X generation.  Texas Instruments could have used the TI-40 College (France), which added the rectangular/polar and degrees/degrees-minute-seconds conversion, along with the DIV command (don't know what that does). 

This is a solid calculator for a basic scientific calculator.  I like the emphasis on fractions and it contains the basic scientific calculator functions. 

Eddie

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

DM 41L and HP 41C: Generating a Polynomial Given Its Roots

DM 41L and HP 41C:  Generating a Polynomial Given Its Roots

Introduction

Generate the coefficients of a polynomial (up to the order 4) with the roots a_0, a_1, a_2, and a_3.   The resulting polynomial is:

p(x) = (x - a_0) * (x - a_1) * (x - a_2) * (x - a_3) * (x - a_4)

p(x) = r_4 * x^4 + r_5 * x^3 + r_6 * x^2 + r_7 * x + r_8

The default is a polynomial where the lead coefficient is positive.   If you want a polynomial where the lead coefficient is negative, multiply every coefficient by -1.

Instructions

Store the four roots in registers R00, R01, R02, and R03 respectively.  Run POLY4.  Coefficients are shown briefly as they are calculated.  They are can be recalled by the registers in decreasing order of x:  R04, R05, R06, R07, and R08.

DM 41L and HP 41C Program: POLY4

01  LBL^T POLY4
02 1
03 STO 04
04 PSE 
05 RCL 00
06 CHS
07 RCL 01
08 -
09 RCL 02
10 -
11 RCL 03
12 -
13 STO 05
14 PSE
15 RCL 01
16 RCL 02
17 +
18 RCL 03
19 +
20 RCL 00
21 * 
22 RCL 02
23 RCL 03
24 +
25 RCL 01
26 *
27 +
28 RCL 02
29 RCL 03
30 *
31 +
32 STO 06
33 PSE
34 RCL 01
35 RCL 02
36 *
37 RCL 01
38 RCL 03
39 *
40 +
41 RCL 02
42 RCL 03
43 *
44 +
45 RCL 00
46 *
47 CHS
48 RCL 01
49 RCL 02
50 *
51 RCL 03
52 *
53 -
54 STO 07
55 PSE
56 RCL 00
57 RCL 01
58 *
59 RCL 02
60 *
61 RCL 03
62 *
63 STO 08
64 RTN

Example

Roots x = -3, x = 3, x= 4, and x= 6

Coefficients: 
R04 = 1
R05 = -10
R06 = 15
R07 = 90
R08 = -216

Polynomial:  p(x) = x^4 - 10 * x^3 + 15 * x^2 + 90 * x - 216

Eddie

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

DM 41L and HP 41C: Schur-Cohn Algorithm

DM 41L and HP 41C:  Schur-Cohn Algorithm

Introduction

The Schur-Cohn Algorithm tests whether the roots of a polynomial p(x) lies with in the unit circle.  That is for the polynomial p(x):

p(x) = a_0 * x^n + a_1 * x^(n - 1) + a_2 * x^(n - 2) + ... + a_n

For all the roots of p(x), r_0, r_1, r_2, ... , r_n ,   |r_k| < 1.    The test covers both real and complex roots.  Keep in mind that this test doesn't tell us what the roots are, only a clue to whether the roots lie in the unit circle or not.

The Schur-Cohn Algorithm returns the test answers b_0, b_1, etc.  If |b_k| < 1 holds true for each b_k for k = 0 to n, then we can conclude that |r_k| <  1.

The following program is adopted from Peter Henrici's book Computational Analysis With the HP-25 Pocket Calculator [Henrici, 111] for the HP 41C and the Swiss Micros DM 41L.  The program deals with polynomials up to the fourth order, but can be expanded on an RPN calculator with more program space and memory registers.

Instructions:

For the fourth order polynomial:

p(x) = a_0 * x^4 + a_1 * x^3 + a_2 * x^2 + a_1 * x + a_4

To the coefficients in the following registers:

R00 = a_0
R01 = a_1
R02 = a_2
R03 = a_3
R04 = a_4

The test results will pop momentarily during execution.  The program is done when after four test values, 4 is in the display.  The last test value is stored in R05.

DM 41L and HP 41C Program:  SCHUR (77 bytes, 11 registers of memory) 

01  LBL^T SCHUR
02  CLX
03  STO 07
04  LBL 03
05  STO 06
06  RCL 04
07  RCL 00 
08  /
09  STO 05    
10 PSE   // show test values
11 RCL 04
12  *
13 ST- 00
14 RCL 01
15 RCL 02
16 RCL 03
17 LBL 15
18 RCL 06
19 X=0?
20 GTO 24
21 RDN   // R↓
22 RCL 04
23 1
24 ST- 06
25 RDN
26 GTO 15
27 LBL 24
28 RDN
29 RCL 05
30 *
31 ST- 01
32 RDN
33 RCL 05
34 *
35 ST- 02
36 RDN
37 RCL 05
38  *
39 ST- 03
40 RCL 03
41 STO 04
42 RCL 02
43 STO 03
44 RCL 01
45 STO 02
46 RCL 00
47 STO 01
48 1
49 ST+ 07
50 4
51 RCL 07
52 X
53 GTO 03

Example

p(x) = 20 * x^4 - 16 * x^3 - 2 * x^2 + 2.08 * x + 0.21

Store the following values:
R0 = 20
R1 = -16
R2 = -2
R3 = 2.08
R4 = 0.21

Test Values:
0.0105
0.1124
-0.0090
-0.8074

For the record the roots of p(x) are 0.5, -0.1, 0.7, and -0.3.

This blog entry is not for use for commercial purposes.

Source:

Henrici, Peter.  Computational Analysis With the HP-25 Calculator  A Wiley-Interscience Publication. John Wiley & Sons: New York 1977 .  ISBN 0-471-02938-6

Eddie

TI-60: Triangle Numbers

TI-60:  Triangle Numbers

Introduction

The following program calculates T_n, the nth triangle number.  Triangle numbers are determined by the sum:

T_n = ∑k, from k = 0 to n

The TI-60 does not have either comparison commands or a goto command.  Luckily for the TI-60 (not the case for the TI-55 III), the RST sends the program pointer back to the first step (Step 00) and continues to execute.

We simulate a loop by counting down, when the counter reaches -1, the square root function will cause an error, forcing execution to stop.  You'll see this trick towards the end of the program.

Instructions:

1.  In RUN mode (outside of LRN), store the following values:

R0 = 0 
R1 = n

2.  Execute the program by pressing [RST], [R/S].  The program is done when you see "Error".

3.  Clear the error by pressing [CE/C].

4.  Recall R0.  This is your triangle number.  (R1 will have -1).

TI-60 Program: Triangle Numbers

PC00  OP71:  RCL 
PC01  OP01:  1
PC02  OP61:  STO 
PC03  OP85:  +
PC04  OP00:  0     // sum R1 in R0
PC05  OP75:  -
PC06  OP01:  1
PC07  OP95:  =
PC08  OP61:  STO
PC09  OP01:  1    // store R1 -1 in R1
PC10  OP86:  √  // square root of R1, if less than 0, program stops<0 ending="" error="" font="" loop="" occurs="" the="">
PC11  OP22:  RST  // reset back to step 00

Examples:

T_5 = 15  (0 [STO] 0,  5 [ST0] 1, [RST] [R/S], [CE/C]  [RCL] 0)

T_35 = 630  (in about 45 seconds)

Eddie

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

HP Prime and TI-84 Plus CE: Double Integral by Monte Carlo Method

HP Prime and TI-84 Plus CE:  Double Integral by Monte Carlo Method

Introduction

There are several ways to calculate integrals.  The most direct way is to find the anti derivative and evaluate it.  Another way is to use an approximation technique.   One technique to find a double integral is to evaluate random points on a given area (circle, rectangle, etc) and use the average in its calculations.

With the Monte Carlo method, we are at mercy of the random numbers selected.  Hence to improve our chances of getting a more accurate result, we have to use a large amount of random points for the calculation.  The flip side is that the Monte Carlo method, similar to the Romberg sum and Simpsons Method, finding the antiderivatve is not required. 

About MONTEINT

The program MONTEINT calculates the double integral

∫ ∫ f(x,y) dx dy

over a circular area with center (X, Y) and radius R.   Please be sure the function is continuous and exists over the entire circular region.   In order to help with the calculation time is that the program runs the calculations several times and averages the results.  You enter the number of batches (B, number of times the calculation is made) and the number of sample points of each batch (N). 

A first for me:  this is the first diagram I drew and posted using Microsoft Paint 3D! 

HP Prime Monte Carlo Double Integral Program: MONTEINT

EXPORT MONTEINT(f,C,D,R,B,N)
BEGIN
// 2018-12-11 EWS
// 'Z(X,Y)', center X,
// center Y, radius,
// number of batches,
// batch size

LOCAL V,S,J,K,I;

V:=0;
FOR K FROM 1 TO B DO
S:=0;
J:=0;
FOR I FROM 1 TO N DO
X:=RANDOM(C-R,C+R);
Y:=RANDOM(D-R,D+R);
IF (X-C)^2+(Y-D)^2≤R^2 THEN

// for EVAL() to work,
// X and Y must be global
S:=S+EVAL(f);
J:=J+1;
END;
END;
V:=V+(π*R^2/J)*S;
END;
V:=V/B;
RETURN V;
END;

Note:  Enter f(X,Y) in single quotes. 

TI 84 Plus Monte Carlo Double Integral Program: MONTEINT

Program MONTEINT

"2018-12-08 EWS"
Input "Z(X,Y) AS A STRING:",Y₀
FnOff 0
Disp "DISK"
Input "CENTER X: ",C
Input "CENTER Y: ",D
Input "RADIUS: ",R
"SETUP VOLUME"
0→V
Input "BATCHES: ",B
Input "SAMPLE SIZE: ",N
"LOOPS"
For(K,1,B)
"SUM, TRUE COUNT"
0→S
0→J
For(I,1,N)
2*R*rand+(C-R)→X
2*R*rand+(D-R)→Y
If (X-C)²+(Y-D)²≤R²
Then
S+Y₀→S
J+1→J
End
End
V+(π*R²/J)*S→V
Disp "BATCH",K,"/",B
End
V/B→V
Disp "APPROX DOUBLE","INTEGRAL:",V

Note:  I put code that shows progress since the TI-84 Plus CE has a slower processor than the HP Prime.

Examples

Full decimal answers will vary. 

Example 1:

f(X,Y) = √(LN(X + Y)) over a circle with center (5,5) of radius 1.   20 batches of 50 points each.

The value is approximately 4.76.

Example 2 - A Cautionary Tale:

f(X,Y) = X^2 + X*Y + Y^3 over a circle with center (0,0) of radius 2.  20 batches of 50 points each.

To illustrate how the Monte Carlo method is at the mercy of random points, the results vary between 11.3 to 13.1. 

I am going to increase the number of batches to 30 with 100 points each.  Thankfully with the HP Prime's processing speed, this won't increase the calculation by a large amount.  The TI-84 Plus CE will require a longer wait. 

After five trials, I get values between 12.07 and 13.15.

Increasing the number of points per batch to 200, I get results of 11.3 to 12.6.  The function f(X,Y) = X^2 + X*Y + Y^3  may not be a good function for the Monte Carlo method.

CAUTION:  The Monte Carlo is a method to calculate double integrals but it's only for approximate values only.  Several evaluations may be necessary to get a satisfactory approximation. 

Eddie

Source:

Ward Cheney and David Kincaid.  Numerical Mathematics and Computing Sixth Edition Thomson Brooks/Cole: Belmont, CA 2008  ISBN-13 978-0-495-11475-8

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

Fun with the Infinite Series 1 + x + x^2 + x^3 + x^4 + x^5 + ...

Fun with the Infinite Series 1 + x + x^2 + x^3 + x^4 + x^5 + ...


The Series and Its Derivatives

Let F be the infinite series:

F = 1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + ... = ∑ x^k from k = 0 to ∞.

Working with derivatives:

----

First Derivative of F:  (F' = dF/dx)

F' = 1 + 2*x + 3*x^2 + 4*x^3 + 5*x^4 + 6*x^5 + 7*x^6 + ... 
= ∑ (k+1)*x^k from k = 0 to ∞

----

Second Derivative of F:  (F'' = d^2F/dx^2)

F'' = 2 + 6*x + 12*x^2 + 20*x^3 + 30*x^4 + 42*x^5 + 56*x^6 + ....

Factor out a 2:
= 2 * (1 + 3*x + 6*x^2 + 10*x^3 + 15*x^4 + 21*x^5 + 28*x^6 + .... )

Note the sequence 1, 3, 6, 10, 15, 21, 28...  These are triangle numbers, denoted as T_n.  

T_1 = 1
T_2 = 1 + 2 = 3
T_3 = 1 + 2 + 3 = 6
T_4 = 1 + 2 + 3 + 4 = 10 
and so on.

Using summation notation,  T_n = ∑ k from k = 1 to n

Going back to the series:

F'' = 2 + 6*x + 12*x^2 + 20*x^3 + 30*x^4 + 42*x^5 + 56*x^6 + ....
= 2 * (1 + 3*x + 6*x^2 + 10*x^3 + 15*x^4 + 21*x^5 + 28*x^6 + .... )
= 2 * (∑ x^k * T_k+1 from k = 0 to ∞)

In nested summation notation:

= 2 * (∑ x^k * (∑ m from m = 1 to k+1) from k = 0 to ∞)


Addition with F, F', and F"


F + F' = 2 + 3*x + 4*x^2 + 5*x^3 + 6*x^4 + 7*x^5 + 8*x^6 + ...
= ∑ ((k + 2) * x^k from k = 0 to ∞)

----

F + F' + F'' = 4 + 9*x + 16*x^2 + 25*x^3 + 36*x^4 + 49*x^5 + 64*x^6 + ...
= ∑ ((k + 2)^2 * x^k from k = 0 to ∞)

----

F' + F" = 3 + 8*x + 15*x^2 + 24*x^3 + 35*x^4 + 48*x^5 + 63*x^6 + ...

Note the sequence 3, 8, 15, 24, 35, 48, 63... where
3 = 4 - 1 = 2^2 -1
8 = 9 - 1 = 3^2 - 1
15 = 16 - 1 = 4^2 - 1
24 = 25 - 1 = 5^2 - 1
35 = 36 - 1 = 6^2 - 1
48 = 49 - 1 = 7^2 - 1
63 = 64 - 1 = 8^2 - 1
and so on...

This can be summarized as ∑( (k + 2)^2 - 1 from k = 0 to ∞)

Hence:
F' + F" = 3 + 8*x + 15*x^2 + 24*x^3 + 35*x^4 + 48*x^5 + 63*x^6 + ...
= ∑  ((k + 2)^2 - 1) * x^k from k = 0 to ∞)

Multiplying F and F' by x and x^2 

F = 1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + ... = ∑ x^k from k = 0 to ∞.
x * F = x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + ... = ∑ x^(k+1) from k = 0 to ∞.
x^2 * F =  x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + ... = ∑ x^(k+2) from k = 0 to ∞.

-----

F + x * F = 1 + 2*x + 2*x^2 + 2*x^3 + 2*x^4 + 2*x^5 + 2*x^6 + ...
= 2 - 1 + 2*x + 2*x^2 + 2*x^3 + 2*x^4 + 2*x^5 + 2*x^6 + ...
= 2 + 2*x + 2*x^2 + 2*x^3 + 2*x^4 + 2*x^5 + 2*x^6 + ... - 1
= 2 * F - 1

This is one way to dervie the formula for the Infinite Geometric Series (for |x| < 1), to solve for F:

F + x * F  = 2 * F - 1
F + x * F - 2 * F = -1
F * (1 + x - 2) = -1
F * (x - 1) = -1
F = -1 / (x - 1)
F = 1 / (1 - x)   (keep this mind, this is true only when |x| < 1)

----

F - x * F = (1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + ... ) - (x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + ... )
= F * (1 - x)

For |x| < 1,

F * (1 - x) = 1 / (1 - x) * (1 - x) = 1

In general:
F - x * F = (1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + ... ) - (x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + ... )
= 1 + (x - x) + (x^2 - x^2) + (x^3 - x^3) + (x^4 - x^4) + (x^5 - x^5) + (x^6 - x^6) + ...
= 1

F - x * F = 1

----

F' = 1 + 2*x + 3*x^2 + 4*x^3 + 5*x^4 + 6*x^5 + 7*x^6 + ... 
= ∑ ((k + 1) * x^k from k = 0 to ∞)

x * F' = x + 2*x^2 + 3*x^3 + 4*x^4 + 5*x^5 + 6*x^6 + 7*x^7 + ... 
= ∑ ((k + 1) * x^(k + 1) from k = 0 to ∞)


x^2 * F' = x^2 + 2*x^3+ 3*x^4 + 4*x^5 + 5*x^6 + 6*x^7 + 7*x^8 + ... 
= ∑ ((k + 1) * x^(k + 1) from k = 0 to ∞)

----

F' + x * F' = 1 + 3*x + 5*x^2 + 7*x^3 + 9*x^4  + 11*x^5 + 13*x^6 + ...

The sequence of 1, 3, 5, 7, 9, 11, 13, ... is the sequence of odd numbers which can be summarized as:

∑(2 * k + 1 from k = 0 to ∞)

Then:

F' + x * F'  =  F' * (1 + x) = ∑( (2*k + 1) * x^k from k = 0 to ∞)

----

F' + x * F' + x^2 * F' = 1 + 3*x + 6*x^2 + 9*x^3 + 12*x^4 + 15*x^5 + 18*x^6 + ...
= 1 + ( 3*x + 6*x^2 + 9*x^3 + 12*x^4 + 15*x^5 + 18*x^6 + ... )
= 1 + ∑(3 * k * x^k from k = 0 to ∞)

Eddie

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.