Bürgi's Kunstweg to Calculate Sines

[Thomas Klemm]

Reference
From: Jost Bürgi’s Method for Calculating Sines

Abstract

Quote:From various sources we know that the Swiss instrument maker and mathematician Jost Bürgi (1552-1632) found a new way of calculating any sine value.
Many mathematicians and historians of mathematics have tried to reconstruct his so-called “Kunstweg” (“skillful / artful method”), but they did not succeed.
Now a manuscript by Bürgi himself has been found which enables us to understand his procedure.
The main purpose of this article is to explain Bürgi’s method. It is totally different from the conventional way to calculate sine values which was used until the 17th century.

Division of a right angle into nine parts

Program to initialize the registers:

Code:

00 { 39-Byte Prgm }
01▸LBL "INIT"
02 CLRG
03 2
04 STO 01
05 4
06 STO 02
07 6
08 STO 03
09 7
10 STO 04
11 8
12 STO 05
13 9
14 STO 06
15 10
16 STO 07
17 11
18 STO 08
19 12
20 STO 09
21 END

Quote:The other nine cells are filled with an arbitrary series of natural numbers, for which Bürgi selects 2, 4, 6, 7, 8, 9, 10, 11, 12.
Bürgi does not mention that the numbers in column 1 are in principle his first approximations for sin 0° to sin 90° in relation to the value 12 for sin 90°.
Nor does he mention that the calculation will be shorter if the numbers correspond approximately to the ratio of the sine values of 10° to 90°.

Code:

00 { 42-Byte Prgm }
01▸LBL "STEP"
02 8
03 STO 00
04 2
05 STO÷ 09
06 RCL 09
07▸LBL 00
08 STO+ IND 00
09 RCL IND 00
10 DSE 00
11 GTO 00
12 2.009
13 STO 00
14 RCL 01
15▸LBL 01
16 STO+ IND 00
17 RCL IND 00
18 ISG 00
19 GTO 01
20 END

After the initialisation we perform 4 steps to get the values of column 5:

XEQ INIT
XEQ STEP
R/S
R/S
R/S

RCL 04
8273441

RCL 09
12871192

÷
0.64278747

Compare this to the correct value:

40
SIN
0.64278761


Division of a right angle into 90 parts

It's easy to adjust the program to do the same for \(n = 90\) to calculate the sines of all degrees from sin 1° to sin 90°.
To initialize the registers I use an approximation of the sines similar to before: \(\text{round}(12 \sin(x))\):

Code:

00 { 34-Byte Prgm }
01▸LBL "INIT-90"
02 CLRG
03 90
04 STO 00
05▸LBL 00
06 RCL 00
07 SIN
08 12
09 ×
10 .5
11 +
12 IP
13 STO IND 00
14 DSE 00
15 GTO 00
16 END

Code:

00 { 46-Byte Prgm }
01▸LBL "STEP-90"
02 89
03 STO 00
04 2
05 STO÷ 90
06 RCL 90
07▸LBL 00
08 STO+ IND 00
09 RCL IND 00
10 DSE 00
11 GTO 00
12 2.09
13 STO 00
14 RCL 01
15▸LBL 01
16 STO+ IND 00
17 RCL IND 00
18 ISG 00
19 GTO 01
20 END

SIZE 91

XEQ INIT-90
XEQ STEP-90
R/S
R/S

RCL 47
310796277956

RCL 90
424964111250

÷
0.73134712

Compare this to the correct value:

47
SIN
0.73135370

Here we performed only 3 steps instead of 4.
That's why the accuracy is not as good.
But we can also notice that the numbers get bigger than before.

[Thomas Klemm]

Double Difference

We can use the addition formualas to simplify the double difference:

\(
\sin(\alpha \pm \beta) = \sin(\alpha) \cos(\beta) \pm \cos(\alpha) \sin(\beta)
\)

This leads to:

\(
\begin{align}
\Delta^2 \sin(\alpha)
&:= \left[\sin(\alpha + \beta) - \sin(\alpha)\right] - \left[\sin(\alpha) - \sin(\alpha - \beta)\right] \\
&= \sin(\alpha + \beta) - 2 \sin(\alpha) + \sin(\alpha - \beta) \\
&= \sin(\alpha) \cos(\beta) + \cos(\alpha) \sin(\beta) - 2 \sin(\alpha) + \sin(\alpha) \cos(\beta) - \cos(\alpha) \sin(\beta) \\
&= 2 \sin(\alpha) \left(\cos(\beta) - 1\right) \\
\end{align}
\)

We notice that it is proportional to \( \sin(\alpha) \) independent of \( \alpha \).

In the limit of \( \beta \to 0 \) we get:
\(
\lim_{\beta \to 0} \frac{\cos(\beta) - 1}{\beta^2} = - \frac{1}{2}
\)

This leads to the 2nd derivative of the sine function:
\(
\sin(\alpha){''} = \lim_{\beta \to 0} \frac{\Delta^2 \sin(\alpha)}{\beta^2} = - \sin(\alpha)
\)

Conclusion

If we apply the double difference on a sequence \( \{a_j\} \) where \( j \in \{1, \cdots, n\} \) and the values are sines of angles in arithmetic progression, we get a sequence that is proportional to the original.

The proportional factor is:

\(
2 \left(\cos(\beta) - 1\right)
\)

Here \( \beta \) is the difference between consecutive angles.

However we have to consider two cases at the boundary.

Lower bound: j = 1

In this case we have:

\(
\alpha = \beta
\)

This leads to:

\(
\begin{align}
\sin(\alpha + \beta) - 2 \sin(\alpha) + \sin(\alpha - \beta)
&= \sin(2 \alpha) - 2 \sin(\alpha) + \sin(0) \\
&= 2 \sin(\alpha) \cos(\alpha) - 2 \sin(\alpha) \\
&= 2 \sin(\alpha) \left(\cos(\alpha) - 1 \right) \\
&= 2 \sin(\alpha) \left(\cos(\beta) - 1 \right) \\
\end{align}
\)

We end up with the same result as before.

Upper bound: j = n

In this case we have:

\(
\alpha = 90^\circ
\)

Here we calculate only the single difference:

\(
\sin(\alpha) - \sin(\alpha - \beta) = 1 - \cos(\beta)
\)

We notice that apart from the factor \( -2 \) we get the same result since \( \sin(\alpha) = 1 \).

Matrix Notation

We can therefore describe the double difference operation with the matrix \( \Delta \):

\(
\Delta = \left[\begin{matrix}
2 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
-1 & 2 & -1 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & -1 & 2 & -1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & -1 & 2 & -1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & -1 & 2 & -1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & -1 & 2 & -1 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & -1 & 2 & -1 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & -1 & 2 & -1\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & -2 & 2
\end{matrix}\right]
\)

Bürgi seems to have noticed that the result was less precise than before.
So he reversed the process.

\(
\Sigma = \Delta^{-1} = \left[\begin{matrix}
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & \frac{1}{2}\\
1 & 2 & 2 & 2 & 2 & 2 & 2 & 2 & 1\\
1 & 2 & 3 & 3 & 3 & 3 & 3 & 3 & \frac{3}{2}\\
1 & 2 & 3 & 4 & 4 & 4 & 4 & 4 & 2\\
1 & 2 & 3 & 4 & 5 & 5 & 5 & 5 & \frac{5}{2}\\
1 & 2 & 3 & 4 & 5 & 6 & 6 & 6 & 3\\
1 & 2 & 3 & 4 & 5 & 6 & 7 & 7 & \frac{7}{2}\\
1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 4\\
1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & \frac{9}{2}
\end{matrix}\right]
\)

Calculating the Sines

We can simply copy and paste the following matrices into Free42:

1 1 1 1 1 1 1 1 0.5
1 2 2 2 2 2 2 2 1
1 2 3 3 3 3 3 3 1.5
1 2 3 4 4 4 4 4 2
1 2 3 4 5 5 5 5 2.5
1 2 3 4 5 6 6 6 3
1 2 3 4 5 6 7 7 3.5
1 2 3 4 5 6 7 8 4
1 2 3 4 5 6 7 8 4.5

ENTER
ENTER
ENTER

2
4
6
7
8
9
10
11
12

×

63
124
181
232
276
312
339
356
362

×

2064
4065
5942
7638
9102
10290
11166
11703
11884

×

67912
133760
195543
251384
299587
338688
367499
385144
391086

×

2235060
4402208
6435596
8273441
9859902
11146776
12094962
12675649
12871192


Note: What appeared to be simple turned out to be a problem.
The separator in the matrix has to be a tabulator.
But I couldn't figure out how to do that when you copy it from this page.
It is always replaced by a blank.
Thus my recommendation for now is to copy the matrix into a text-editor, replace the blanks by tabs and copy the result into Free42.

References

• The mathematics behind Jost Bürgi’s method for calculating sine tables
• Jost Bürgi’s skillful computation of sines
  

[Albert Chan]

(05-08-2022 01:19 PM)Thomas Klemm Wrote:  We can therefore describe the double difference operation with the matrix \( \Delta \):

\(
\Delta = \left[\begin{matrix}
2 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
-1 & 2 & -1 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & -1 & 2 & -1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & -1 & 2 & -1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & -1 & 2 & -1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & -1 & 2 & -1 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & -1 & 2 & -1 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & -1 & 2 & -1\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & -2 & 2
\end{matrix}\right]
\)

Asymmetry of matrix, [2,-1] top, [-2,2] bottom , we should match edge cases coefficients.

Lower bound, α = β, sin(α-β) = sin(0°) = 0

Δ² = sin(α+β) - 2 sin(α) + sin(α-β) = -dot([sin(α), sin(α+β)], [2, -1])

Upper bound, α = 90°, sin(α+β) = sin(α-β) = cos(β)

Δ² = sin(α+β) - 2 sin(α) + sin(α-β) = -dot([sin(α-β), sin(α)], [-2,2])

Quote:Note: What appeared to be simple turned out to be a problem.
The separator in the matrix has to be a tabulator.
But I couldn't figure out how to do that when you copy it from this page.
It is always replaced by a blank.
Thus my recommendation for now is to copy the matrix into a text-editor, replace the blanks by tabs and copy the result into Free42.

Another way is to click QUOTE, copy raw matrix data (with tabs), then paste to Free42

[Albert Chan]

(05-09-2022 01:48 AM)Albert Chan Wrote:  Asymmetry of matrix, [2,-1] top, [-2,2] bottom , we should match edge cases coefficients.

Another way is to split central difference to 2 steps, forward diff, then backward diff.
Below, 1 ≡ 10°, (0-1) ≡ sin(0°) - sin(10°) = -sin(10°) ≡ -1

Code:

0     F        B F
1   (0-1)   (0-1)-(1-2)
2   (1-2)   (1-2)-(2-3)
3   (2-3)   (2-3)-(3-4)
4   (3-4)   (3-4)-(4-5)
5   (4-5)   (4-5)-(5-6)
6   (5-6)   (5-6)-(6-7)
7   (6-7)   (6-7)-(7-8)
8   (7-8)   (7-8)-(8-9)
9   (8-9)   (10-9)-(9-8) = (8-9)*2

>>> from numpy import *
>>> F = [[-1,0,0,0,0,0,0,0,0],
...            [1,-1,0,0,0,0,0,0,0],
...            [0,1,-1,0,0,0,0,0,0],
...            [0,0,1,-1,0,0,0,0,0],
...            [0,0,0,1,-1,0,0,0,0],
...            [0,0,0,0,1,-1,0,0,0],
...            [0,0,0,0,0,1,-1,0,0],
...            [0,0,0,0,0,0,1,-1,0],
...            [0,0,0,0,0,0,0,1,-1]]
>>>
>>> B = [[1,-1,0,0,0,0,0,0,0],
...            [0,1,-1,0,0,0,0,0,0],
...            [0,0,1,-1,0,0,0,0,0],
...            [0,0,0,1,-1,0,0,0,0],
...            [0,0,0,0,1,-1,0,0,0],
...            [0,0,0,0,0,1,-1,0,0],
...            [0,0,0,0,0,0,1,-1,0],
...            [0,0,0,0,0,0,0,1,-1],
...            [0,0,0,0,0,0,0,0, 2]]
>>>
>>> M = -matrix(B) * matrix(F)
>>> print M
[[ 2 -1 0 0 0 0 0 0 0]
 [-1 2 -1 0 0 0 0 0 0]
 [ 0 -1 2 -1 0 0 0 0 0]
 [ 0 0 -1 2 -1 0 0 0 0]
 [ 0 0 0 -1 2 -1 0 0 0]
 [ 0 0 0 0 -1 2 -1 0 0]
 [ 0 0 0 0 0 -1 2 -1 0]
 [ 0 0 0 0 0 0 -1 2 -1]
 [ 0 0 0 0 0 0 0 -2 2]]

If we apply this M operator, it generated massive cancellation errors.
Applying its inverse turn all entries non-negative, without cancellation errors.

>>> print linalg.inv(M)
[[ 1. 1. 1. 1. 1. 1. 1. 1. 0.5]
 [ 1. 2. 2. 2. 2. 2. 2. 2. 1. ]
 [ 1. 2. 3. 3. 3. 3. 3. 3. 1.5]
 [ 1. 2. 3. 4. 4. 4. 4. 4. 2. ]
 [ 1. 2. 3. 4. 5. 5. 5. 5. 2.5]
 [ 1. 2. 3. 4. 5. 6. 6. 6. 3. ]
 [ 1. 2. 3. 4. 5. 6. 7. 7. 3.5]
 [ 1. 2. 3. 4. 5. 6. 7. 8. 4. ]
 [ 1. 2. 3. 4. 5. 6. 7. 8. 4.5]]

We can also think of inv(M) as "going back in time", turning divergence into convergence.