Bürgi's Kunstweg to Calculate Sines
05-07-2022, 11:31 AM
Post: #1
 Thomas Klemm Senior Member Posts: 2,141 Joined: Dec 2013
Bürgi's Kunstweg to Calculate Sines
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.
05-08-2022, 01:19 PM
Post: #2
 Thomas Klemm Senior Member Posts: 2,141 Joined: Dec 2013
RE: Bürgi's Kunstweg to Calculate Sines
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)$$

\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$$

\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
05-09-2022, 01:48 AM
Post: #3
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
RE: Bürgi's Kunstweg to Calculate Sines
(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
05-09-2022, 04:55 PM
Post: #4
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
RE: Bürgi's Kunstweg to Calculate Sines
(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.
 « Next Oldest | Next Newest »