Automatic differentiation using dual numbers

[Thomas Klemm]

TLDR;
A dual number has the form \(x + {x}' \varepsilon\) where \(x, {x}' \in \mathbb{R}\) and \(\varepsilon \ne 0\), but \(\varepsilon^2 = 0\).
Think of \(\varepsilon\) similar to the imaginary unit \(i\) but with a different rule.
They can be used to calculate derivatives.

Motivation

An older thread about calculators and numerical differentiation, and an even older one about what the 1st CAS pocket calculator was, made me wonder if there was ever a calculator using automatic differentiation using dual numbers.

So I thought I give it a try abusing the complex numbers of the HP-42S.
The advantage is that addition and subtraction as well as multiplication by a constant are already handled.
Complex numbers can also be stored and recalled.

Most of the functions simply unwrap the complex number, evaluate the function \(f\) and its derivative \(f'\), and rewrap it into a complex number:

Code:

01▸LBL "f"          # z
02 COMPLEX          # y x
03 X<>Y             # x y
04 f                # f(x) y
05 X<>Y             # y f(x)
06 LASTX            # x y f(x)
07 f'               # f'(x) y f(x)
08 ×                # f'(x)*y f(x)
09 COMPLEX          # f(z)
10 RTN              #

But then there are more complicated cases like * and / or Y↑X.
In the simplest case, the stack is retained except for the T register.
But in other cases, the stack is lost entirely.
Also LASTX isn't handled.

We end up with a somewhat crippled RPN calculator, but programming a function is still easy.
Instead of saving the entire stack before entering and restoring it on exit, I suggest storing the values in variables only when needed later.

Only Y↑X needs two additional registers. All the other operations use only the stack.
Also remember that this function works only with a constant (i.e. non-complex) exponent.
Otherwise, use LN and E↑X instead.

Only after writing the program for the HP-42S I noticed that Ángel Martin already wrote a module for the HP-41C.
I haven't tried it yet, just read the excellent documentation.
Thus, I strongly recommend to give it a try.

Custom Menu

These are my assignments:

LN    E↑X   SINH  ASINH COSH  ACOSH
*     /     1/X   Y↑X   SQRT  X↑2
SIN   ASIN  COS   ACOS  TAN   ATAN


Program

Code:

00 { 468-Byte Prgm }
01▸LBL "DUAL"       #
02▸LBL "*"          # w z
03 COMPLEX          # v u z
04 X<> ST Z         # z u v
05 COMPLEX          # y x u v
06 RCL× ST Z        # u*y x u v
07 R↓               # x u v u*y
08 STO× ST Z        # x u x*v u*y
09 ×                # x*u x*v u*y
10 X<> ST Z         # u*y x*v x*u
11 +                # x*v+u*y x*u
12 COMPLEX          # z*w
13 RTN              #
14▸LBL "/"          # w z
15 COMPLEX          # v u z
16 X<>Y             # u v z
17 STO÷ ST Z        # u v z/u
18 ÷                # v/u z/u
19 X<>Y             # z/u v/u
20 COMPLEX          # y/u x/u v/u
21 RCL ST Z         # v/u y/u x/u
22 RCL× ST Z        # xv/u^2 y/u x/u
23 -                # y/u-xv/u^2 x/u
24 COMPLEX          # z/w
25 RTN              #
26▸LBL "1/X"        # z
27 COMPLEX          # y x
28 X<>Y             # x y
29 1/X              # 1/x y
30 X<>Y             # y 1/x
31 LASTX            # x y 1/x
32 X↑2              # x^2 y 1/x
33 ÷                # y/x^2 1/x
34 +/-              # -y/x^2 1/x
35 COMPLEX          # 1/z
36 RTN              #
37▸LBL "X↑2"        # z
38 COMPLEX          # y x
39 X<>Y             # x y
40 X↑2              # x^2 y
41 X<>Y             # y x^2
42 RCL× ST L        # x*y x^2
43 STO+ ST X        # 2*x*y x^2
44 COMPLEX          # z^2
45 RTN              #
46▸LBL "SQRT"       # z
47 COMPLEX          # y x
48 X<>Y             # x y
49 SQRT             # sqrt(x) y
50 STO÷ ST Y        # sqrt(x) y/sqrt(x)
51 X<>Y             # y/sqrt(x) sqrt(x)
52 2                # 2 y/sqrt(x) sqrt(x)
53 ÷                # y/(2*sqrt(x)) sqrt(x)
54 COMPLEX          # sqrt(z)
55 RTN              #
56▸LBL "Y↑X"        # n z
57 STO 00           # n -> R 00
58 R↓               # z
59 COMPLEX          # y x
60 X<>Y             # x y
61 STO 01           # x -> R 01
62 RCL 00           # n x y
63 STO× ST Z        # n x n*y
64 1                # 1 n x n*y
65 -                # n-1 x n*y
66 Y↑X              # x^{n-1} n*y
67 STO× ST Y        # x^{n-1} n*x^{n-1}*y
68 RCL× 01          # x^n n*x^{n-1}*y
69 X<>Y             # n*x^{n-1}*y x^n
70 COMPLEX          # z^n
71 RTN              #
72▸LBL "E↑X"        # z
73 COMPLEX          # y x
74 X<>Y             # x y
75 E↑X              # e^x y
76 STO× ST Y        # e^x e^x*y
77 X<>Y             # e^x*y e^x
78 COMPLEX          # e^z
79 RTN              #
80▸LBL "LN"         # z
81 COMPLEX          # y x
82 X<>Y             # x y
83 STO÷ ST Y        # x y/x
84 LN               # ln(x) y/x
85 X<>Y             # y/x ln(x)
86 COMPLEX          # ln(z)
87 RTN              #
88▸LBL "SIN"        # z
89 COMPLEX          # y x
90 X<>Y             # x y
91 SIN              # sin(x) y
92 X<>Y             # y sin(x)
93 LASTX            # x y sin(x)
94 COS              # cos(x) y sin(x)
95 ×                # cos(x)*y sin(x)
96 COMPLEX          # sin(z)
97 RTN              #
98▸LBL "ASIN"       # z
99 COMPLEX          # y x
100 X<>Y            # x y
101 ASIN            # asin(x) y
102 X<>Y            # y asin(x)
103 RCL ST Y        # asin(x) y asin(y)
104 COS             # sqrt(1-x^2) y asin(x)
105 ÷               # y/sqrt(1-x^2) asin(x)
106 COMPLEX         # asin(z)
107 RTN             #
108▸LBL "COS"       # z
109 COMPLEX         # y x
110 X<>Y            # x y
111 COS             # cos(x) y
112 X<>Y            # y cos(x)
113 LASTX           # x y cos(x)
114 SIN             # sin(x) y cos(x)
115 +/-             # -sin(x) y cos(x)
116 ×               # -sin(x)*y cos(x)
117 COMPLEX         # cos(z)
118 RTN             #
119▸LBL "ACOS"      # z
120 COMPLEX         # y x
121 X<>Y            # x y
122 ACOS            # acos(x) y
123 X<>Y            # y acos(x)
124 RCL ST Y        # acos(x) y acos(y)
125 SIN             # sqrt(1-x^2) y acos(x)
126 ÷               # y/sqrt(1-x^2) acos(x)
127 +/-             # -y/sqrt(1-x^2) acos(x)
128 COMPLEX         # acos(z)
129 RTN             #
130▸LBL "TAN"       # z
131 COMPLEX         # y x
132 X<>Y            # x y
133 TAN             # tan(x) y
134 ENTER           # tan(x) tan(x) y
135 X↑2             # tan^2(x) tan(x) y
136 1               # 1 tan^2(x) tan(x) y
137 +               # 1+tan^2(x) tan(x) y y
138 R↑              # y 1+tan^2(x) tan(x) y
139 ×               # (1+tan^2(x))*y tan(x)
140 COMPLEX         # tan(z)
141 RTN             #
142▸LBL "ATAN"      # z
143 COMPLEX         # y x
144 X<>Y            # x y
145 ATAN            # atan(x) y
146 X<>Y            # y atan(x)
147 RCL ST Y        # atan(x) y atan(x)
148 COS             # 1/sqrt(1+x^2) y atan(x)
149 X↑2             # 1/(1+x^2) y atan(x)
150 ×               # y/(1+x^2) atan(x)
151 COMPLEX         # atan(z)
152 RTN             #
153▸LBL "SINH"      # z
154 COMPLEX         # y x
155 X<>Y            # x y
156 SINH            # sinh(x) y
157 X<>Y            # y sinh(x)
158 LASTX           # x y sinh(x)
159 COSH            # cosh(x) y sinh(x)
160 ×               # cosh(x)*y sinh(x)
161 COMPLEX         # sinh(z)
162 RTN             #
163▸LBL "ASINH"     # z
164 COMPLEX         # y x
165 X<>Y            # x y
166 ASINH           # asinh(x) y
167 X<>Y            # y asinh(x)
168 RCL ST Y        # asinh(x) y asinh(y)
169 COSH            # sqrt(x^2-1) y asinh(x)
170 ÷               # y/sqrt(x^2-1) asinh(x)
171 COMPLEX         # asinh(z)
172 RTN             #
173▸LBL "COSH"      # z
174 COMPLEX         # y x
175 X<>Y            # x y
176 COSH            # cosh(x) y
177 X<>Y            # y cosh(x)
178 LASTX           # x y cosh(x)
179 SINH            # sinh(x) y cosh(x)
180 ×               # sinh(x)*y cosh(x)
181 COMPLEX         # cosh(z)
182 RTN             #
183▸LBL "ACOSH"     # z
184 COMPLEX         # y x
185 X<>Y            # x y
186 ACOSH           # acosh(x) y
187 X<>Y            # y acosh(x)
188 RCL ST Y        # acosh(x) y acosh(y)
189 SINH            # sqrt(x^2-1) y acos(x)
190 ÷               # y/sqrt(x^2-1) acos(x)
191 COMPLEX         # acosh(z)
192 RTN             #
193▸LBL "TANH"      # z
194 COMPLEX         # y x
195 X<>Y            # x y
196 TANH            # tanh(x) y
197 ENTER           # tanh(x) tanh(x) y
198 X↑2             # tanh^2(x) tanh(x) y
199 1               # 1 tanh^2(x) tanh(x) y
200 X<>Y            # tanh^2(x) 1 tanh(x) y
201 -               # 1-tan^2(x) tanh(x) y y
202 R↑              # y 1-tan^2(x) tanh(x) y
203 ×               # (1-tan^2(x))*y tanh(x)
204 COMPLEX         # tanh(z)
205 RTN             #
206▸LBL "ATANH"     # z
207 COMPLEX         # y x
208 X<>Y            # x y
209 ATANH           # atanh(x) y
210 X<>Y            # y atanh(x)
211 RCL ST Y        # atanh(x) y atanh(x)
212 COSH            # 1/sqrt(1-x^2) y atanh(x)
213 X↑2             # 1/(1-x^2) y atanh(x)
214 ×               # y/(1-x^2) atanh(x)
215 COMPLEX         # atanh(z)
216 RTN             #
217▸LBL "NROOT"     # n z
218 1/X             # 1/n z
219 STO 00          # 1/n -> R 00
220 R↓              # z
221 COMPLEX         # y x
222 X<>Y            # x y
223 STO 01          # x -> R 01
224 ABS             # |x| y
225 RCL 00          # 1/n |x| y
226 STO× ST Z       # 1/n |x| y/n
227 1               # 1 1/n |x| y/n
228 -               # 1/n-1 |x| y/n
229 Y↑X             # |x|^{1/n-1} y/n
230 STO× ST Y       # |x|^{1/n-1} |x|^{1/n-1}*y/n
231 RCL× 01         # x^{1/n} |x|^{1/n-1}*y/n
232 X<>Y            # |x|^{1/n-1}*y/n x^{1/n}
233 COMPLEX         # z^{1/n}
234 END

Examples

Arithmetic

\(
(3 + 4\varepsilon)(5 + 6\varepsilon)
\)

3 ENTER 4 COMPLEX
5 ENTER 6 COMPLEX
XEQ "*"

15 i38

Evaluation of a function and its derivative

\(
\begin{align}
f(x) = \frac{1}{\sqrt{3 + \sin(x)}}
\end{align}
\)

Evaluate \(f(2)\) and \(f'(2)\).

2 ENTER 1 COMPLEX
XEQ "SIN"
3 +
XEQ "SQRT"
XEQ "1/X"

0.505767 i0.026920

Polynomial

Write a program to evaluate the following polynomial:

\(
f(x) = 2x^3 + 3x^2 + 5x + 7
\)

Code:

00 { 35-Byte Prgm }
01▸LBL "f(x)"
02 RCL "x"
03 2
04 ×
05 3
06 +
07 RCL "x"
08 XEQ "*"
09 5
10 +
11 RCL "x"
12 XEQ "*"
13 7
14 +
15 END

Hint: Please note the use of the ordinary multiplication × with a constant in line 04.

Finding a root with Newton's algorithm

\(
\begin{align}
x \mapsto x - \frac{f(x)}{f'(x)}
\end{align}
\)

The program is straight forward since both the function and its derivative can be calculated with a single call to "f(x)":

Code:

00 { 34-Byte Prgm }
01▸LBL "NEWTON"
02 1
03 COMPLEX
04 STO "x"
05 XEQ "f(x)"
06 COMPLEX
07 ÷
08 RCL "x"
09 COMPLEX
10 X<> ST Z
11 -
12 END

Now we can use it to find the root of the previous polynomial:

-1 XEQ "NEWTON" R/S R/S R/S …

-1.6
-1.459479553903345724907063197026022
-1.445652028784242179107592346048789
-1.445528468451557138940594968587292
-1.445528458679522363979375907703971
-1.445528458679522302862795910270397
-1.445528458679522302862795910270395

-1.44552845867952230286279591027039479019

---

\(e^x - 5 = 0\)

Code:

00 { 19-Byte Prgm }
01▸LBL "f(x)"
02 RCL "x"
03 XEQ "E↑X"
04 5
05 -
06 END

1 XEQ "NEWTON" R/S R/S R/S …

1.839397205857211607977618850807305
1.633963151762016740112346092531373
1.609736212513373931404520946197064
1.609437956921145415593181982892169
1.609437912434101364149332898839731
1.609437912434100374600759333226677
1.609437912434100374600759333226188

1.60943791243410037460075933322618763952

---

\(2x - \sin(x) - 1 = 0\)

Code:

00 { 26-Byte Prgm }
01▸LBL "f(x)"
02 RCL "x"
03 2
04 ×
05 RCL "x"
06 XEQ "SIN"
07 -
08 1
09 -
10 END

1 XEQ "NEWTON" R/S R/S R/S …

0.891395995328754051864551902982202
0.8878657494003519361186921183484038
0.8878622115744122839022415445850894
0.8878622115708660240357050746759603
0.8878622115708660240357015114947114
0.8878622115708660240357015114947116
0.887862211570866024035701511494712

0.88786221157086602403570151149471173497

---

Find a root of : \(x^x = \pi\)

Code:

00 { 28-Byte Prgm }
01▸LBL "f(x)"
02 RCL "x"
03 XEQ "LN"
04 RCL "x"
05 XEQ "*"
06 XEQ "E↑X"
07 PI
08 -
09 END

Hint: As already mentioned we can not use Y↑X here since the exponent is not a constant.

2 XEQ "NEWTON" R/S R/S R/S …

1.873252698249433854127306813867173
1.854459717903199817846714313283606
1.854106089961332624497916868416055
1.854105967921040960682023166413
1.854105967921026432748370718616169
1.854105967921026432748370718410293

1.85410596792102643274837071841029324542

---

Solve Leonardo di Pisa's equation: \(x^3+2x^2+10x−20=0\)

Code:

00 { 34-Byte Prgm }
01▸LBL "f(x)"
02 RCL "x"
03 2
04 +
05 RCL "x"
06 XEQ "*"
07 10
08 +
09 RCL "x"
10 XEQ "*"
11 20
12 -
13 END

1 XEQ "NEWTON" R/S R/S R/S …

1.411764705882352941176470588235294
1.369336470588235294117647058823529
1.368808188617531954512471834907757
1.368808107821374524808003343189187
1.368808107821372635227414330022359
1.368808107821372635227414330021326

1.36880810782137263522741433002132553954


Resources

• Automatic differentiation using dual numbers
• HP-41 DoubleDown Module Released
• What was the 1st CAS pocket calculator?
• Calculators and numerical differentiation
• Calculators and numerical differentiation
• The strange cousin of the complex numbers -- the dual numbers