### Developers Club geek daily blog

2 years, 11 months ago

Important refining — the calculator normal, without the sin button. As in accounts department or in the market.

Under a cat three different candidate solutions from different eras, from ancient Samarkand to the USA of times of cold war.

### Simple solution

The first that comes to mind — here such spell:

Let's translate this confused score for the calculator into more clear bc language. It is often used as the calculator in the command line of UNIX-like operating systems. Let's see approximately following:

``````  scale = 7
x = 355/113/180
x-x^3/6
.0174524
``````

From where it undertook
We decompose a sine in a row about zero, we take the first several members of this row and we substitute one degree. In this case a corner small therefore it is possible to be limited to a polynom of the third degree:

sin(x) ≅ x — x3/6

Before substitution the degree should be transferred in radians multiplication on `π` and division on 180 °.

The separate prize is necessary noticed strange digits 355 and 113. Found them in our Chinese companion Tszu Chunchzhi (祖沖之) at the time of Qi's dynasty (479 — 502). The relation 355/113 is the only approach of number `π` in rational fraction which is shorter than decimal representation of similar accuracy.

### Interesting solution

The well-known trick described above appeared only in 1715. Nevertheless values of trigonometric functions were known much earlier, and with considerably bigger accuracy.

The manager of the Samarkand observatory of Giyas-ad-din Jamsheed ibn Masoud of al Kashi (غیاث   ) made tables of trigonometric functions to within the 16th sign till 1429. In translation from Persian on bc its spell in relation to our task looked approximately so:

``````  scale = 16

sin30 = .5
cos30 = sqrt(3)/2

sin45 = sqrt(2)/2
cos45 = sin45

sin75 = sin30*cos45+cos30*sin45
cos75 = sqrt(1-sin75^2)

cos36 = (1+sqrt(5))/4
sin36 = sqrt(1-cos36^2)

sin72 = 2*sin36*cos36
cos72 = sqrt(1-sin72^2)

(sin3 = sin75*cos72-cos75*sin72)
.0523359562429430

(x = sin3/3)
.0174453187476476
(x = (sin3+4*x^3)/3)
.0174523978055315
(x = (sin3+4*x^3)/3)
.0174524064267667
(x = (sin3+4*x^3)/3)
.0174524064372703
(x = (sin3+4*x^3)/3)
.0174524064372831
(x = (sin3+4*x^3)/3)
.0174524064372831
``````

Pay attention that we still use only addition, subtraction, multiplication, division and the square root. At desire all these operations can be executed in general on a piece of paper in a column. Learned to consider the square root in a column even at school earlier. It is boring, but it is not really difficult.

What shamanism is
Let's sort magic of al Kashi on steps.

``````  sin30 = .5
cos30 = sqrt(3)/2

sin45 = sqrt(2)/2
cos45 = sin45
``````

The sine and cosine 30 ° and 45 ° were known still to ancient Greeks.

``````  sin75 = sin30*cos45+cos30*sin45
``````

The sine of the sum of corners 30 ° and 45 ° is available. To al Kashi this formula was displayed by other Persian astronomer, Abul-Vafa Mahomed ibn Mahomed ibn Yaha ibn Ismail ibn Abbas al-Buzdzhani.

``````  cos75 = sqrt(1-sin75^2)
``````

Pythagorean trousers are extensively equal.

``````  cos36 = (1+sqrt(5))/4
sin36 = sqrt(1-cos36^2)
``````

It from the regular pentagon known still to ancient Greeks.

``````  sin72 = 2*sin36*cos36
cos72 = sqrt(1-sin72^2)
``````

Again sine of the sum and Pythagorean theorem.

``````  (sin3 = sin75*cos72-cos75*sin72)
.0523359562429430
``````

We consider a sine of a difference 75 ° and 72 ° and we receive a sine 3 °.

Now it is possible to spread out 3 ° to the sum of three corners of 1 °, but there is a hitch — receive cubic equation:

sin 3 ° = 3 x — 4 x3

where x = sin 1 °. Nobody was able to solve cubic equations analytically then still.

Wise al Kashi noticed that it is possible to express this equation in the following form:

f(x) = (sin 3 ° + 4 x3) / 3

and then to apply to f (x) a method of simple iteration. I remind that at that time neither Newton, nor Rafson were born yet.

``````  (x = sin3/3)
``````

First approximation.

``````  .0174453187476476
(x = (sin3+4*x^3)/3)
.0174523978055315
(x = (sin3+4*x^3)/3)
.0174524064267667
(x = (sin3+4*x^3)/3)
.0174524064372703
(x = (sin3+4*x^3)/3)
.0174524064372831
(x = (sin3+4*x^3)/3)
.0174524064372831
``````

We receive 16 signs after five iterations.

### As considers the calculator

The inquisitive reader can have a legal question: how the calculator which has such button considers value of a sine?

It turns out that the majority of calculators use absolutely the third method — "digit behind digit", been born in a subsoil of military industrial complex of the USA during cold war.

And here B-58 bomber
Jack Volder who then worked in the company Konveyr on the navigation calculator of an above-mentioned bomber thought up this algorithm.

The main benefit of the "digit behind digit" method is that it uses only addition operations and divisions into two (which is easy for implementing the right shift).

Besides, algorithm it is possible to force to work directly in a binary-coded decimal code which is used in the majority of calculators, but in the example given below we will not get into this jungle.

The iterative algorithm also uses the table of arctangents, on one on iteration. The table needs to be counted in advance:

``````#include <stdio.h>
#include <math.h>

int main(int argc, char **argv)
{
int bits = 32;
int cordic_one = 1 << (bits - 2);
printf("// Число с фиксированной точкой, соответствующее единице с плавающей точкой\n");
printf("static const int cordic_one = 0x%08x;\n", cordic_one);
printf("static const int cordic_table[] = {\n");
double k = 1;
for (int i = 0; i < bits; i++) {
printf("0x%08x, // 0x%08x * atan(1/%.0f) \n", (int)(atan(pow(2, -i)) * cordic_one), cordic_one, pow(2, i));
k /= sqrt(1 + pow(2, -2 * i));
}
printf("};\n");
printf("static const int cordic_k = 0x%08x; // %.16f * 0x%08x\n", (int)(k * cordic_one), k, cordic_one);
}
``````

At the same time the scaling coefficient is considered `cordic_k`.

After that it is possible to count notorious sin 1 ° so:

``````#include <stdio.h>
#include <math.h>

// Число с фиксированной точкой, соответствующее единице с плавающей точкой
static const int cordic_one = 0x40000000;
static const int cordic_table[] = {
0x3243f6a8, // 0x40000000 * atan(1/1)
0x1dac6705, // 0x40000000 * atan(1/2)
0x07f56ea6, // 0x40000000 * atan(1/8)
0x03feab76, // 0x40000000 * atan(1/16)
0x01ffd55b, // 0x40000000 * atan(1/32)
0x00fffaaa, // 0x40000000 * atan(1/64)
0x007fff55, // 0x40000000 * atan(1/128)
0x003fffea, // 0x40000000 * atan(1/256)
0x001ffffd, // 0x40000000 * atan(1/512)
0x000fffff, // 0x40000000 * atan(1/1024)
0x0007ffff, // 0x40000000 * atan(1/2048)
0x0003ffff, // 0x40000000 * atan(1/4096)
0x0001ffff, // 0x40000000 * atan(1/8192)
0x0000ffff, // 0x40000000 * atan(1/16384)
0x00007fff, // 0x40000000 * atan(1/32768)
0x00003fff, // 0x40000000 * atan(1/65536)
0x00001fff, // 0x40000000 * atan(1/131072)
0x00000fff, // 0x40000000 * atan(1/262144)
0x000007ff, // 0x40000000 * atan(1/524288)
0x000003ff, // 0x40000000 * atan(1/1048576)
0x000001ff, // 0x40000000 * atan(1/2097152)
0x000000ff, // 0x40000000 * atan(1/4194304)
0x0000007f, // 0x40000000 * atan(1/8388608)
0x0000003f, // 0x40000000 * atan(1/16777216)
0x0000001f, // 0x40000000 * atan(1/33554432)
0x0000000f, // 0x40000000 * atan(1/67108864)
0x00000008, // 0x40000000 * atan(1/134217728)
0x00000004, // 0x40000000 * atan(1/268435456)
0x00000002, // 0x40000000 * atan(1/536870912)
0x00000001, // 0x40000000 * atan(1/1073741824)
0x00000000, // 0x40000000 * atan(1/2147483648)
};
static const int cordic_k = 0x26dd3b6a; // 0.6072529350088813 * 0x40000000

void cordic(int theta, int&s, int&c)
{
c = cordic_k;
s = 0;
for (int k = 0; k < 32; ++k) {
int d = (theta >= 0) ? 0 : -1;
int tx = c - (((s >> k) ^ d) - d);
int ty = s + (((c >> k) ^ d) - d);
c = tx; s = ty;
theta -= ((cordic_table[k] ^ d) - d);
}
}

int main(void)
{
double alpha = M_PI / 180;
int sine, cosine;
cordic(alpha * cordic_one, sine, cosine);
printf("CORDIC:   %.8f\nExpected: %.8f\n", (double)sine / cordic_one, sin(alpha));
}
``````

Result:
``````CORDIC:   0.01745240
Expected: 0.01745241
``````

There are 32 iterations therefore there was a small error. Calculators usually use 40 iterations.