HP-15C Program Collection

1 Prime Number Tester

The primes are the atoms of arithmetic—every integer factors uniquely into primes. Yet they scatter through the number line with maddening irregularity. This program uses trial division, the most ancient primality test: divide by every integer up to √n. If none divide evenly, the number is prime. It's not fast for large numbers, but it's foolproof, and watching the HP-15C methodically test divisors connects you to mathematicians who did this by hand for centuries.

LBL A        ; Entry point - number in X
STO 0        ; Store n in R0
2
STO 1        ; Start divisor at 2
LBL 0        ; Loop start
RCL 1
x²
RCL 0
x≤y?         ; If divisor² > n, it's prime
GTO 1
RCL 0
RCL 1
÷
FRAC
x=0?         ; If evenly divisible, not prime
GTO 2
1
STO+ 1       ; Increment divisor
GTO 0        ; Continue loop
LBL 1        ; Prime exit
1
RTN
LBL 2        ; Composite exit
0
RTN

Usage

Enter number, press GSB A. Returns 1 (prime) or 0 (composite).

Examples: 97 → 1 (prime), 91 → 0 (= 7 × 13)

2 Euclidean GCD Algorithm

Euclid's algorithm for finding the greatest common divisor appears in his Elements around 300 BCE, making it one of the oldest algorithms still in daily use. The insight is beautiful: GCD(a,b) = GCD(b, a mod b). Keep replacing the larger with the remainder until you hit zero. What remains is the GCD. It's fast, elegant, and the foundation for everything from reducing fractions to modern cryptography.

LBL B        ; Y: a, X: b
LBL 0
x=0?         ; If b = 0, done
GTO 1
LastX        ; Recall b
x<>y         ; Stack: b, a
LastX        ; Stack: b, a, b
÷
FRAC
LastX
×            ; a mod b
GTO 0
LBL 1
R↓           ; Result in X
RTN

Usage

Enter a ENTER b, press GSB B. Returns GCD(a,b).

Example: 48 ENTER 18 GSB B → 6

3 Fibonacci Generator

Leonardo of Pisa introduced these numbers to Europe in 1202, describing rabbit population growth. Each number is the sum of the two before it: 1, 1, 2, 3, 5, 8, 13... The ratio of consecutive Fibonacci numbers converges to the golden ratio φ = (1+√5)/2 ≈ 1.618. These numbers appear in pine cones, sunflower seeds, and the branching of trees—nature's favorite sequence.

LBL C        ; Initialize: puts 0, 1 in registers
0
STO 0
1
STO 1
RTN

LBL 1        ; Generate next Fibonacci
RCL 0
RCL 1
STO 0        ; F(n-1) becomes F(n-2)
+
STO 1        ; F(n) = F(n-1) + F(n-2)
RTN

Usage

GSB C to initialize, then GSB 1 repeatedly.

Sequence: 1, 2, 3, 5, 8, 13, 21, 34, 55, 89...

Divide consecutive outputs to watch the ratio converge to φ ≈ 1.6180339887...

4 Pascal's Triangle: Binomial Coefficients

Pascal's triangle contains the binomial coefficients C(n,k)—the number of ways to choose k items from n. Each entry is the sum of the two entries above it. This simple rule generates numbers that appear in probability, combinatorics, and algebra. The triangle also hides Fibonacci numbers, powers of 2, and countless other patterns. This program computes any entry using the formula C(n,k) = n!/(k!(n-k)!), optimized to avoid overflow.

LBL D        ; Y=n, X=k, compute C(n,k)
STO 0        ; Store k
R↓
STO 1        ; Store n
RCL 0
RCL 1
RCL 0
-            ; n - k
x<y?         ; Is k > n-k?
x<>y         ; Use smaller for efficiency
STO 0        ; k = min(k, n-k)
1
STO 2        ; Result = 1
1
STO 3        ; i = 1
LBL 0        ; Main loop
RCL 3
RCL 0
x<y?         ; i > k?
GTO 1        ; Done
RCL 1
RCL 3
-
1
+            ; n - i + 1
STO× 2       ; result *= (n-i+1)
RCL 3
STO÷ 2       ; result /= i
1
STO+ 3       ; i++
GTO 0
LBL 1
RCL 2
RTN

Usage

Enter n ENTER k, press GSB D.

Examples: C(10,3) = 120, C(52,5) = 2,598,960 (poker hands!)

5 Digital Root: Casting Out Nines

Add the digits of any number. If the sum has more than one digit, add those digits. Repeat until you get a single digit. This is the digital root, and it equals the number's remainder when divided by 9 (except 9 maps to 9, not 0). Medieval accountants used "casting out nines" to check arithmetic. The technique has deep connections to modular arithmetic and has been used for centuries to catch calculation errors.

LBL E        ; n in X
LBL 0        ; Digit sum loop
STO 0        ; Store current n
0
STO 1        ; Digit sum = 0
LBL 1        ; Sum digits loop
RCL 0
x=0?
GTO 2        ; Done with this round
10
÷
FRAC
10
×            ; Extract last digit
STO+ 1       ; Add to sum
RCL 0
10
÷
INT          ; Remove last digit
STO 0
GTO 1
LBL 2
RCL 1        ; Current digit sum
10
x<y?         ; Sum < 10?
RTN          ; Yes: done
R↓
RCL 1        ; No: continue reducing
GTO 0

Usage

Enter number, press GSB E.

Examples: 12345 → 6, 999 → 9, 19683 → 9

6 The Collatz Conjecture: Mathematics' Simplest Unsolved Problem

Take any positive integer. If it's even, halve it. If it's odd, triple it and add one. Repeat. The Collatz Conjecture claims every starting number eventually reaches 1. It's been verified for numbers up to 2⁶⁸, yet no one can prove it must always be true. Paul Erdős said "Mathematics may not be ready for such problems." Your HP-15C can trace these hailstone sequences—called that because they rise and fall unpredictably like hailstones in a thundercloud before finally dropping to earth.

LBL .0       ; Enter n in X
LBL .1       ; Main loop
x=0?
RTN          ; Safety check
1
x=y?         ; Reached 1?
RTN          ; Done!
R↓
2
÷            ; n/2
FRAC
x=0?         ; Was n even?
GTO .2       ; Yes: go halve it
LastX        ; No: recall n/2
2
×            ; Restore n
3
×
1
+            ; 3n + 1
PSE          ; Pause to display
GTO .1
LBL .2       ; Even case
LastX        ; Recall n/2
2
×            ; n
2
÷            ; n/2 (integer)
PSE          ; Pause to display
GTO .1

Usage

Enter starting number, press GSB .0.

Watch the sequence unfold with pauses. Ends when reaching 1.

Try: 27 → Takes 111 steps, reaches a peak of 9,232!

7 Monte Carlo π Estimator

Imagine throwing darts randomly at a square with an inscribed quarter-circle. The ratio of darts landing inside the circle to total darts approximates π/4. This is Monte Carlo simulation—using randomness to solve deterministic problems. The method was formalized during the Manhattan Project, but the underlying idea is ancient. Each dart throw is an experiment; aggregate enough experiments and order emerges from chaos.

LBL .3       ; n = iterations in X
STO I        ; Save original n
STO 0        ; Loop counter
0
STO 1        ; Hits counter = 0
LBL .4       ; Loop
RAN#         ; Random x ∈ [0,1)
x²
RAN#         ; Random y ∈ [0,1)
x²
+            ; x² + y²
1
x≤y?         ; If x² + y² ≤ 1, inside circle
1
STO+ 1       ; Increment hits (conditional)
1
STO- 0       ; Decrement counter
RCL 0
x≠0?
GTO .4       ; Continue if counter > 0
RCL 1        ; Recall hits
4
×
RCL I        ; Recall original n
÷            ; 4 × hits / n ≈ π
RTN

Usage

Enter iteration count (try 100-500), press GSB .3.

More iterations = better estimate (but slower).

The error decreases as 1/√n, so quadrupling iterations only halves the error. Still, watching π emerge from random noise is magical.

8 Wallis Product: A 17th Century Path to π

In 1655, John Wallis discovered an infinite product that converges to π/2. It's hauntingly simple: (2/1)·(2/3)·(4/3)·(4/5)·(6/5)·(6/7)·... Each fraction pairs an even number with its odd neighbors. Wallis found this while trying to compute the area under a circle—the same problem that would later inspire Newton to develop calculus. The convergence is glacially slow, but there's something meditative about watching your calculator inch toward transcendence.

LBL .5       ; n = number of terms in X
STO 0        ; Loop counter
1
STO 1        ; Running product
2
STO 2        ; Current even number
LBL .6       ; Main loop
RCL 2        ; even
ENTER
ENTER
1
-            ; even - 1
÷            ; even/(even-1)
STO× 1       ; Multiply into product
RCL 2
RCL 2
1
+            ; even + 1
÷            ; even/(even+1)
STO× 1       ; Multiply into product
2
STO+ 2       ; Next even number
1
STO- 0       ; Decrement counter
RCL 0
x≠0?
GTO .6
RCL 1
2
×            ; Product × 2 = π approximation
RTN

Usage

Enter number of terms, press GSB .5.

Results: 10 terms → 3.067..., 100 terms → 3.133..., 1000 terms → 3.140...

9 e via Its Limit Definition

Euler's number e = 2.71828... can be defined as the limit of (1 + 1/n)ⁿ as n approaches infinity. This captures the essence of continuous compounding: if you earn 100% interest compounded n times, as n grows the final amount approaches e. This program lets you explore this convergence. Watch how the approximation improves as n increases, approaching but never quite reaching the irrational limit.

LBL .7       ; n in X
STO 0        ; Store n
1
RCL 0
÷            ; 1/n
1
+            ; 1 + 1/n
RCL 0
y^x          ; (1 + 1/n)ⁿ
RTN

Usage

Enter n, press GSB .7.

Convergence: n=10 → 2.5937..., n=100 → 2.7048..., n=10000 → 2.7181...

The series e = 1 + 1/1! + 1/2! + 1/3! + ... converges much faster, but the limit definition reveals e's connection to compound interest.

10 Continued Fraction Expansion

Every real number can be expressed as a continued fraction: a₀ + 1/(a₁ + 1/(a₂ + 1/(a₃ + ...))). For rationals, this terminates; for irrationals, it goes forever. The remarkable thing is that truncating a continued fraction gives you the *best* rational approximation with that denominator size. This is how we know 355/113 is an extraordinary approximation to π—it comes directly from π's continued fraction.

LBL .8       ; Number in X
STO 0        ; Working value
LBL .9       ; Extract next coefficient
RCL 0
INT          ; aₙ = floor(x)
PSE          ; Display coefficient
RCL 0
FRAC         ; Fractional part
x=0?         ; If zero, we're done
RTN
1/x          ; 1 / frac(x)
STO 0        ; This is new x
GTO .9       ; Continue

Usage

Enter number, press GSB .8. Watch coefficients display.

Examples:

π → [3; 7, 15, 1, 292, ...] (that 292 explains why 355/113 is so good)

√2 → [1; 2, 2, 2, 2, ...] (all 2s forever!)

φ → [1; 1, 1, 1, 1, ...] (all 1s—the "most irrational" number)

11 Harmonic Numbers: The Slow March to Infinity

The harmonic series 1 + 1/2 + 1/3 + 1/4 + ... diverges, but with excruciating slowness. To exceed 10, you need 12,367 terms. To exceed 100, you need about 10⁴³ terms—more than atoms in the universe. The partial sums, called harmonic numbers Hₙ, appear everywhere: in probability (the coupon collector problem), in computer science (quicksort analysis), and in number theory.

LBL 2        ; n in X, compute Hₙ
STO 0        ; Store n
0
STO 1        ; Sum = 0
1
STO 2        ; k = 1
LBL 3
RCL 2
1/x          ; 1/k
STO+ 1       ; Add to sum
1
STO+ 2       ; k++
RCL 2
RCL 0
x≤y?         ; k > n?
GTO 4
GTO 3
LBL 4
RCL 1        ; Return Hₙ
RTN

Usage

Enter n, press GSB 2.

Values: H₁₀ = 2.9289..., H₁₀₀ = 5.1873..., H₁₀₀₀ = 7.4854...

For large n, Hₙ ≈ ln(n) + γ, where γ ≈ 0.5772... is the Euler-Mascheroni constant.

12 Newton-Raphson Square Root

Newton's method finds roots by iterating: xₙ₊₁ = xₙ - f(xₙ)/f'(xₙ). For square roots, f(x) = x² - n, giving xₙ₊₁ = ½(xₙ + n/xₙ). The convergence is quadratic—each iteration roughly doubles the correct digits. This is the algorithm inside every CPU's floating-point unit, and it runs beautifully on the HP-15C.

LBL 5        ; n in X, find √n
STO 0        ; Store n
2
÷
STO 1        ; Initial guess = n/2
LBL 6        ; Iteration loop
RCL 0
RCL 1
÷            ; n/x
RCL 1
+            ; x + n/x
2
÷            ; New estimate
STO 2        ; Store temporarily
RCL 1
-
ABS
EEX
8
CHS          ; 1e-8 tolerance
x<y?         ; Converged?
GTO 7
RCL 2
STO 1        ; Update x
GTO 6
LBL 7
RCL 2
RTN

Usage

Enter number, press GSB 5. Converges rapidly.

Example: 2 GSB 5 → 1.414213562

13 Babylonian Square Root: Ancient Wisdom

Over 3,500 years ago, Babylonian mathematicians discovered an algorithm for square roots that we now recognize as Newton's method applied to f(x) = x² - n. Start with a guess, average it with n/guess, repeat. The convergence is quadratic—each iteration roughly doubles the correct digits. This algorithm was carved into clay tablets in cuneiform, and it runs beautifully on silicon.

LBL 8        ; n in X, find √n
STO 0        ; Store n
2
÷
STO 1        ; Initial guess = n/2
LBL 9        ; Iteration loop
RCL 0
RCL 1
÷            ; n/guess
RCL 1
+            ; guess + n/guess
2
÷            ; (guess + n/guess)/2
STO 2        ; New guess
RCL 1
-
ABS
1
EEX
9
CHS          ; 10⁻⁹ tolerance
x>y?         ; Converged?
GTO .0
RCL 2
STO 1        ; Update guess
GTO 9
LBL .0
RCL 2        ; Return √n
RTN

Usage

Enter number, press GSB 8.

Example: 2 GSB 8 → 1.414213562 (√2)

14 Numerical Derivative: Calculus Without Limits

Every calculus student learns that the derivative is the limit of a difference quotient. But what if you can't find the limit symbolically? The HP-15C's elegant solution uses the central difference formula: sample the function at x+h and x-h, then divide by 2h. With a small enough step size, you get remarkable accuracy. This transforms the calculator into a numerical differentiating machine.

; First, define your function in LBL .1
; Example: f(x) = x³ - 2x + 1
LBL .1
ENTER
ENTER
×
×            ; x³
LastX
2
×
-            ; x³ - 2x
1
+            ; x³ - 2x + 1
RTN

; Derivative routine
LBL .2       ; x in X register
STO 0        ; Store x
1
EEX
5
CHS          ; h = 0.00001
STO 1
RCL 0
RCL 1
+            ; x + h
GSB .1       ; f(x + h)
RCL 0
RCL 1
-            ; x - h
GSB .1       ; f(x - h)
-            ; f(x+h) - f(x-h)
RCL 1
2
×            ; 2h
÷            ; [f(x+h) - f(x-h)] / 2h
RTN

Usage

Modify LBL .1 with your function. Enter x value, press GSB .2.

Example: For f(x) = x³ - 2x + 1, derivative at x=2:

2 GSB .2 → 10 (exact answer: 3(2)² - 2 = 10 ✓)

15 Quadratic Solver with Complex Roots

The HP-15C was one of the first handheld calculators with true complex number support. This quadratic solver exploits that capability: when the discriminant is negative, it doesn't give up—it returns the complex roots. The formula is the same quadratic formula from high school, but now it works for every equation, real or complex.

LBL .3       ; Enter: Z=a, Y=b, X=c
STO 0        ; c in R0
R↓
STO 1        ; b in R1
R↓
STO 2        ; a in R2
RCL 1
x²           ; b²
RCL 2
RCL 0
×
4
×            ; 4ac
-            ; discriminant = b² - 4ac
SF 8         ; Set complex mode
SQRT         ; √(discriminant) - complex if negative!
STO 3        ; Store √discriminant
RCL 1
CHS          ; -b
RCL 3
+            ; -b + √d
RCL 2
2
×
÷            ; First root: (-b + √d)/(2a)
RCL 1
CHS
RCL 3
-            ; -b - √d
RCL 2
2
×
÷            ; Second root in X, first in Y
RTN

Usage

Enter a ENTER b ENTER c, press GSB .3.

Press Re<>Im to see imaginary parts.

Example: 1 ENTER 2 ENTER 5 GSB .3

Roots of x² + 2x + 5 = 0 are -1 ± 2i

16 Matrix Determinant: The Power of the HP-15C

The HP-15C was revolutionary for including matrix operations in a handheld calculator. With 64 registers configurable as matrices, it could handle serious linear algebra. The determinant tells you if a system of equations has a unique solution (non-zero) or is degenerate (zero). It's the first question to ask about any linear system.

; First, configure and enter your matrix
; Example: 2×2 matrix

f DIM        ; Enter matrix dimension mode
2
ENTER
2
f MATRIX 1  ; Dimension Matrix A as 2×2

f USER       ; Enter User mode (optional)
f MATRIX 1  ; Select Matrix A for input

; Enter elements row by row:
1 STO A     ; a₁₁
2 STO A     ; a₁₂
3 STO A     ; a₂₁
4 STO A     ; a₂₂

; Now compute determinant:
RCL MATRIX A
f MATRIX 7  ; Compute determinant
; Result: -2 (for this example)

Usage

This is more of a workflow than a program—the HP-15C's matrix functions are built in.

Key matrix commands:

f MATRIX 7 → Determinant, f MATRIX 9 → Transpose

The HP-15C can invert matrices up to 8×8 directly. For a 3×3 system of equations, enter the coefficient matrix, invert it, and multiply by the constants vector.

17 The Logistic Map: Chaos in Your Pocket

In 1976, biologist Robert May showed that a simple population model could produce chaos. The equation xₙ₊₁ = r·xₙ·(1-xₙ) looks innocent enough—multiply, subtract, scale. But for certain values of r, the system never settles down. It bounces forever, never repeating, yet never escaping a bounded region. The HP-15C can explore this mathematical wilderness one iteration at a time, letting you watch order dissolve into beautiful chaos.

LBL .4       ; Initialize: X=x₀, Y=r
STO 0        ; Store x in R0
R↓
STO 1        ; Store r in R1
RCL 0        ; Return initial x
RTN

LBL .5       ; Iterate: compute next x
RCL 0        ; x
1
RCL 0
-            ; (1 - x)
×            ; x(1-x)
RCL 1
×            ; r·x·(1-x)
STO 0        ; Store new x
RTN

Usage

Enter r ENTER x₀ then GSB .4 to initialize.

Press GSB .5 repeatedly to iterate.

Try it: r=3.2, x₀=0.5 → Period-2 oscillation

Then try: r=3.57 → Edge of chaos, r=4.0 → Full chaos

For r between 3.57 and 4, the orbit becomes chaotic—tiny changes in x₀ lead to completely different trajectories. This is the famous "sensitive dependence on initial conditions."

18 Compound Interest with Regular Deposits

Einstein allegedly called compound interest the eighth wonder of the world. Whether or not he said it, the math is undeniably powerful. This program answers the practical question: if I start with some principal, add a fixed amount regularly, and earn compound interest, what will I have after n periods? It's the equation underlying every retirement calculator.

LBL .6       ; R0=principal, R1=rate, R2=deposit, R3=periods
; Setup: Store values first, then GSB .6
RCL 3        ; n periods
STO 4        ; Loop counter
RCL 0        ; Starting principal
LBL .7       ; Compound loop
RCL 1
1
+            ; 1 + rate
×            ; balance × (1 + rate)
RCL 2
+            ; + deposit
1
STO- 4
RCL 4
x≠0?
GTO .7
RTN          ; Final balance in X

Usage

R0: Principal R1: Rate (decimal) R2: Deposit/period R3: # Periods

Example: $10,000 initial, 0.5% monthly rate, $500/month, 30 years:

10000 STO 0, .005 STO 1, 500 STO 2, 360 STO 3, GSB .6

Result: ≈ $502,257

19 Day of Week: Zeller's Congruence

Christian Zeller devised this formula in 1887 to compute the day of week for any date. It's a beautiful piece of mathematical engineering, encoding calendar quirks (leap years, the irregular months) into a single modular arithmetic expression. Enter any date and know instantly whether it was a Monday or Friday.

LBL .8       ; Z=year, Y=month, X=day
STO 0        ; day
R↓
STO 1        ; month
R↓
STO 2        ; year
RCL 1
2
x≤y?
GTO .9
RCL 1
12
+
STO 1        ; Adjust month (Jan=13, Feb=14)
1
STO- 2       ; Adjust year
LBL .9
RCL 0        ; d
RCL 1
1
+
13
×
5
÷
INT
+            ; + floor(13(m+1)/5)
RCL 2
+            ; + K
RCL 2
4
÷
INT
+            ; + floor(K/4)
RCL 2
100
÷
INT
STO 3        ; J = century
4
÷
INT
+            ; + floor(J/4)
RCL 3
2
×
-            ; - 2J
7
÷
FRAC
7
×
INT          ; mod 7
RTN

Usage

Enter year ENTER month ENTER day, press GSB .8.

Result: 0=Sat, 1=Sun, 2=Mon, 3=Tue, 4=Wed, 5=Thu, 6=Fri

Example: 2024 ENTER 7 ENTER 4 → July 4, 2024 was a Thursday (5)

20 Gamma Function via Integration

The gamma function extends factorial to non-integers: Γ(n+1) = n! for positive integers, but it's defined for all complex numbers (except non-positive integers). The integral definition is Γ(n+1) = ∫₀^∞ tⁿe⁻ᵗ dt. The HP-15C's built-in numerical integration can evaluate this directly, giving you factorials of fractions.

LBL A        ; n in X (setup for gamma)
STO 4        ; Store n
RTN          ; Return to setup

LBL B        ; Integrand: t^n × e^(-t)
; t is in X from integrate
RCL 4        ; n
y^x          ; t^n
LastX
CHS
e^x          ; e^(-t)
×
RTN

; To use:
; 1. Enter n, GSB A (stores n)
; 2. Enter: 0 ENTER 20
; 3. f ∫ B (integrate label B)
; Result: Γ(n+1) = n!

Usage

Enter n, GSB A (stores n)

Enter 0 ENTER 20 (integration limits)

Press f ∫ B (integrate label B)

Example: n=0.5 gives Γ(1.5) = 0.5! ≈ 0.886 = √π/2

The famous result Γ(1/2) = √π means that (-1/2)! = √π. The gamma function connects factorials to the circle in unexpected ways.

Table of Contents

Part I: Number Theory & Fundamentals

Part II: Approximating Constants

Part III: Numerical Methods

Part IV: Dynamics & Applications

PART I: NUMBER THEORY & FUNDAMENTALS

1 Prime Number Tester

Usage

2 Euclidean GCD Algorithm

Usage

3 Fibonacci Generator

Usage

4 Pascal's Triangle: Binomial Coefficients

Usage

5 Digital Root: Casting Out Nines

Usage

6 The Collatz Conjecture: Mathematics' Simplest Unsolved Problem

Usage

PART II: APPROXIMATING CONSTANTS

7 Monte Carlo π Estimator

Usage

8 Wallis Product: A 17th Century Path to π

Usage

9 e via Its Limit Definition

Usage

10 Continued Fraction Expansion

Usage

11 Harmonic Numbers: The Slow March to Infinity

Usage

PART III: NUMERICAL METHODS

12 Newton-Raphson Square Root

Usage

13 Babylonian Square Root: Ancient Wisdom

Usage

14 Numerical Derivative: Calculus Without Limits

Usage

15 Quadratic Solver with Complex Roots

Usage

16 Matrix Determinant: The Power of the HP-15C

Usage

PART IV: DYNAMICS & APPLICATIONS

17 The Logistic Map: Chaos in Your Pocket

Usage

18 Compound Interest with Regular Deposits

Usage

19 Day of Week: Zeller's Congruence

Usage

20 Gamma Function via Integration

Usage