sin 1° на калькуляторе

Важное уточнение — калькулятор обычный, без кнопки sin. Как в бухгалтерии или на рынке.

Под катом три разных варианта решения из разных эпох, от древнего Самарканда до США времён холодной войны.

Простое решение

Первое, что приходит в голову — вот такое заклинание:


http://habrahabr.ru/post/271889/
Переведём эту путаную партитуру для калькулятора на более понятный язык bc. Он часто используется как калькулятор в командной строке UNIX-подобных операционных систем. Увидим примерно следующее:

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

Откуда это взялось

Разлагаем синус в ряд около нуля, берём первые несколько членов этого ряда и подставляем один градус. В данном случае угол маленький, поэтому можно ограничиться многочленом третьей степени:

sin(x) ≅ x — x³/6

Перед подстановкой градус придётся перевести в радианы умножением на π и делением на 180°.

Отдельный приз полагается заметившим странные цифры 355 и 113. Их нашёл в наш китайский товарищ Цзу Чунчжи (祖沖之) ещё во времена династии Ци (479—502). Отношение 355/113 — это единственное приближение числа π рациональной дробью, которое короче десятичного представления аналогичной точности.

Интересное решение

Описанный выше общеизвестный трюк появился только в 1715 году. Тем не менее значения тригонометрических функций были известны намного раньше, и с заметно большей точностью.

Заведующий Самаркандской обсерваторией Гияс-ад-дин Джамшид ибн Масуд аль-Каши (غیاث الدین جمشید کاشانی‎‎) составил таблицы тригонометрических функций с точностью до 16-го знака ещё до 1429 года. В переводе с персидского на bc его заклинание применительно к нашей задаче выглядело примерно так:

  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

Обратите внимание на то, что мы по-прежнему используем только сложение, вычитание, умножение, деление и квадратный корень. При желании все эти операции можно выполнить вообще на бумажке в столбик. Cчитать квадратный корень в столбик раньше даже учили в школе. Это занудно, но не очень сложно.

Что это за шаманство

Разберём магию аль-Каши по шагам.

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

  sin45 = sqrt(2)/2
  cos45 = sin45

Синус и косинус 30° и 45° были известны ещё древним грекам.

  sin75 = sin30*cos45+cos30*sin45

Налицо синус суммы углов 30° и 45°. Ещё до аль-Каши эту формулу вывел другой персидский астроном, Абуль-Вафа Мухаммад ибн Мухаммад ибн Яхья ибн Исмаил ибн Аббас аль-Бузджани.

  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

Считаем синус разности 75° и 72° и получаем синус 3°.

Теперь можно разложить 3° на сумму трёх углов по 1°, но возникает заминка — получаем кубическое уравнение:

sin 3° = 3 x — 4 x³

где x = sin 1°. Решать кубические уравнения аналитически тогда ещё никто не умел.

Мудрый аль-Каши заметил, что можно выразить это уравнение в следующей форме:

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

и потом применить к f(x) метод простой итерации. Напоминаю, что в то время ни Ньютон, ни Рафсон ещё не родились.

  (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

Получаем 16 знаков после пяти итераций.

Как считает сам калькулятор

У пытливого читателя может возникнуть законный вопрос: как же считает значение синуса калькулятор, у которого есть такая кнопка?

Оказывается, что большинство калькуляторов используют совершенно третий способ — «цифра за цифрой», родившийся в недрах военно-промышленного комплекса США во время холодной войны.

Причём тут бомбардировщик Б-58

Придумал этот алгоритм Джек Волдер, который тогда работал в компании Конвэйр над навигационным вычислителем вышеупомянутого бомбардировщика.

Главное преимущество метода «цифра за цифрой» в том, что он использует только операции сложения и деления на два (которое легко реализовать сдвигом вправо).

Кроме того, алгоритм можно заставить работать прямо в двоично-десятичном коде, который используется в большинстве калькуляторов, но в приведённом ниже примере мы в эти дебри не полезем.

Алгоритм итеративный и использует таблицу арктангенсов, по одному на итерацию. Таблицу нужно посчитать заранее:

#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);
}

Заодно считается масштабирующий коэффициент cordic_k.

После этого посчитать пресловутый sin 1° можно так:

#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) 
0x0fadbafc, // 0x40000000 * atan(1/4) 
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));
}

Результат:

CORDIC:   0.01745240
Expected: 0.01745241

Тут 32 итерации, поэтому осталась небольшая ошибка. Калькуляторы обычно используют 40 итераций.