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

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

    Калькулятор Casio

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

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


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

    355 / 113 / 180 = MC M+ * = * MR / 6 +- = + MR =

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

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

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

    sin(x) ≅ x — x3/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 x3

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

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

    f(x) = (sin 3° + 4 x3) / 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 итераций.

    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 44

      +54
      В принципе, в этой постановке задача решается в уме, поскольку синус малого угла равен углу.
      31415/180 = 15707/9.

      То есть в уме надо разделить 157075 на 9, что не вызывает сомнений, равно 174527 и совпадает с Вашим ответом до 7 знака.
        +9
        Верно. Это то же самое разложение в ряд, только ряд короче.
        +8
        Экскурс в античную тригонометрию. Текст до хабраката мне показался слишком скудным для такой интересной статьи.
          0
          Спасибо за обратную связь, добавил параграф до ката. Исходно хотел сделать условие задачи как можно короче, чтобы можно было понять, не вчитываясь.
          0
          За 355/113 спасибо, правда не знаю, что проще запомнить, учитывая то, что 3.14 точно уже у всех отскакивает от зубов.
            +3
            Но вот 3,1415926 отскакивает далеко не у всех.
              +4
              Это я знаю и помню прекрасно, пи многие знаки мне лишни напрасны = 3,14159265358
              Помню класса с пятого
                +12
                How I want a drink, alchoholic of course, after the heavy lectures involving quantuum mechanics = 3,14159265358979
                  +1
                  В копилку (менее точный, зато с рифмой):

                  Кто и шутя и скоро пожелаетъ пи узнать число, ужъ знаетъ
                    0
                    «Три, четырнадцать, пятнадцать, девяносто два и шесть»
                    +3
                    quantum вместо quantuum.
                    +3
                    IMHO запомнить «3, 14, 159, 26, 53, 58» намного проще, чем сначала выучить стишок, а потом ещё и считать количество букв в словах.
                    Для меня, во всяком случае.
                    +2
                    у меня до 6 знака (3.141592) отскакивает. У жены — почти до двадцатого :) Для тренировки она себе пи до какого-то там знака на пароль поставила :)
                      +3
                      Чтобы нам не ошибиться,
                      Надо правильно прочесть:
                      Три, четырнадцать, пятнадцать,
                      Девяносто два и шесть

                      Вполне себе от зубов :)
                      +7
                      Нас вот так учили запоминать эту дробь:
                      image
                        +1
                        Почему дробь «в горку»?
                          +23
                          Потому что 113355. Легче запомнить.
                      +3
                      Есть мнение, что
                      а) иногда синус аппаратно считают через таблицы (аля Брадиса) с поправками до 3-ей степени. Для машинной точности типа double таблицы получаются не такие уж и большие по размеру, зато очень быстро — пара умножений и сложений;
                      б) иногда синус считают аппаратно разложением в ряд, в котором коэффициенты посчитаны заранее. Тоже остаётся только сложение и умножение.
                        +1
                        a) Верно, именно так, например, считается синус в библиотеке C когда нет поддержки плавающей точки. Но на калькуляторах так делают редко. Умножение на калькуляторе — дорогая операция, а в CORDIC используются только сложение и сдвиги.

                        б) Я раньше тоже так думал, но на практике ни разу не встречал. Подозреваю, что это весьма редкое явление.
                        0
                        берём первые члены ряда Фурье для синуса, упрощаем/сокращаем и получаем простую формулу для синуса любого угла
                          +12
                          Вы наверное имели в виду ряд Тейлора. Это как раз и есть первое из описанных решений.

                          Ряд Фурье для синуса будет очень коротким (все коэффициенты, кроме одного, равны нулю).
                            0
                            конечно же Тейлора, не внимателен (ещё их называют рядами Маклорена)
                              +3
                              Ряд Маклорена — это частный случай ряда Тейлора
                                +5
                                А ряд Тейлора частный случай ряда Лорана
                              –4
                              Член какой придётся в пору, чтоб заправить в зад Тейлору?
                              Оба члена хороши: и Лагранжа, и Коши!
                            +1
                            Описанный выше общеизвестный трюк появился только в 1715 году.
                            Ну ряды для конкретных функций были известны и раньше.
                              +2
                              По первой картинке я уж подумал, что вычисление синуса будет с использованием ряда Тейлора с помощью Arduino, который будет сам кнопочки нажимать…
                                +1
                                Я в одной из библиотек видел такой алгоритм — сначала по формулам приведения сводим аргумент к отрезку [0,pi/4], потом применяем к нему синус или косинус (в зависимости от того, какая функция получилась), заданный в виде аппроксимирующего многочлена (8-й или 9-й степени, в зависимости от требуемой точности).
                                  +8
                                  Хабр — торт. Простой математический пример в статье, и множество комментариев с методами, алгоритмами, и т.п.
                                    +1
                                    > Отношение 355/113 — это единственное приближение числа π рациональной дробью, которое короче десятичного представления аналогичной точности.

                                    Ну, не единственное. В десятичной системе счисления ничего уникального нет. Но «единственное такое, что одновременно короткое и с отличной точностью» — это да.
                                      +1
                                      Вот тут товарищ Jon McLoone специально искал и не нашёл ничего лучше: All Rational Approximations of Pi Are Useless.
                                        +1
                                        Если брать относительный выигрыш — то может быть, других и нет. Странно, что в этом тексте нет упоминания о разложении в цепную (непрерывную) дробь. Казалось бы, чем больше член разложения, который мы отбрасываем, тем больше должен быть выигрыш в цифрах. Если 355/113 (отбросили 292) даёт 1 знак выигрыша, то приближение первыми 453293 членами разложения (отброшен член 12996958) должно дать 5, а то и 6 знаков.
                                        Ещё странно сравнение sqrt(2) и «e». У sqrt(2) цепная дробь ограничена, поэтому выигрыша больше, чем в 1 цифру, мы не получим никогда, но у «e» её члены неограниченно растут (хотя и довольно медленно — линейно), и, грубо говоря, удлинение записи в 10 раз всегда будет давать лишнюю цифру выигрыша.
                                          0
                                          В комментариях пишут то же самое, про convergents of continued fractions. Их можно посчитать через WolframAlpha, например, Last[Convergents[Pi, 6]] = 104348/33215. Но выигрыша нет — надо запоминать 11 цифр, чтобы получить π с точностью 10 цифр.

                                          В любом случае много цифр запоминать смысла нет, так что это уже искусство для искусства.
                                        0
                                        Спасибо за интересную статью!
                                          +3
                                          Ждём продолжения — как посчитать это же на камушках и римскими цифрами 8)).
                                            0
                                            Интересно, почему в вычислении синуса «по цифрам» не делалось округления при сдвиге:
                                            вместо
                                            int tx = c - (((s >> k) ^ d) - d);

                                            было бы лучше (при k>0, конечно),
                                            int tx = c - (((((s >> (k - 1)) + 1) >> 1) ^ d) - d);

                                            Да и таблица арктангенсов в конце могла бы быть точнее. Выбросить 2/3 её, к сожалению, нельзя — поправочный коэффициент перестанет быть предсказуемым — но заменить вторую половину на обычное умножение (которое наверняка реализовано очень эффективно) было бы и быстрее, и точнее.
                                            Впрочем, военным виднее. Они высоко летают.
                                              0
                                              почему в вычислении синуса «по цифрам» не делалось округления при сдвиге

                                              Думаю, что для простоты. Попробовал предложенный вариант, получил тот же результат.

                                              заменить вторую половину на обычное умножение

                                              Умножение на калькуляторе — дорогая операция. В частности именно поэтому не используются приближения многочленами или рядами.
                                              +2
                                              На физике нас учили, что
                                              sin a ≈ a, если a -> 0

                                              Пригодилось множество раз.
                                                +1
                                                Ибо, первый замечательный предел
                                                0
                                                А аппроксимация полиномами Чебышева не подойдет?
                                                  0
                                                  Умножение на калькуляторе — дорогая операция, реализуется программно. А сложение и сдвиг — дешёвые, реализованы аппаратно. Поэтому и не используются приближения многочленами или рядами.
                                                    –1
                                                    Умножение реализуется через сложение и сдвиги, так что ваше утверждение не очевидно. Функция умножения в любом калькуляторе уже реализована, скорость вычисления в пределах миллисекунд — удовлетворительно, т.к. операции всё равно вводятся вручную. Не вижу смысла отказываться от рядов, если алгоритм будет проще — аппаратная реализация будет в итоге проще и дешевле.
                                                      +2
                                                      Посмотрите на картинку в начале статьи. Этот калькулятор питается от маленькой солнечной батареи, и даже солнца ему не надо — достаточно настольной лампы. Тактовая частота меряется не в гигагерцах и даже не в мегагерцах, а в килогерцах. Для него одно умножение занимает миллисекунды, а чтобы многочлен посчитать уже секунды уйдут, неприемлемо.

                                                      Разработчики калькуляторов это ещё в семидесятые годы прикинули и вышло, что CORDIC лучше. Разумеется, на более привычном процессоре помощнее всё может оказаться наоборот.
                                                  +3
                                                  В бейсике «Спектрум» и во многих других бейсиках того времени для вычисления синусов использовалось именно разложение в ряд Чебышева. Думаю, дело в том, что в литературе по численным методам 80-х годов (всякие книги типа МакКракен и Дорн и т.д.) разложение в ряд Чебышева было самым продвинутым методом аппроксимации, которые рассматривались. Программисты подобного софта по большей части просто не знали, какие еще бывают методы.

                                                  CORDIC более перспективен для аппаратной реализации, т.к. не требуется умножение. На самом деле умножение — очень затратная операция по времени либо аппаратным ресурсам. Если разрабатывать микросхему, вычисляющую синус — то лучше отказаться от умножителя и реализовать CORDIC, чем использовать умножитель и аппроксимацию, основанную на умножении. То же касается ассемблерного кода для процессоров, не имеющих аппаратного умножителя.

                                                  Если есть умножитель, то несколько лучше, чем разложение в ряд Чебышева, применить аппроксимацию минимакс-полиномом. Коэффициенты такого полинома вычислить сложнее, чем Чебышева, но стоит только освоить алгоритм Ремеза — как трудности исчезнут. Аппроксимация же рядом Тейлора (Маклорена) — это очень неэффективный подход.

                                                  Если имеется много памяти — то лучше всего хранить в ней таблицу значений синуса и применять линейную интерполяцию для промежуточных углов. Линейная интерполяция по сути сводится к одному умножению.

                                                  При вычислении синусов еще важно то, в какой последовательности надо их вычислять. Скажем, если нужно вычислить синусы многих углов, равноотстоящих друг от друга — то есть рекуррентные соотношения. y[i] = -y[i-2] + 2*b*y[i-1], константа b задает шаг. Одно умножение и два сложения — и это при том, что данная формула дает точные ответы. Конечно, при ее использовании будет накапливаться ошибка округления, поэтому представление чисел должно иметь много бит, и самые младшие из них будут постепенно зашумляться. Есть и другие рекуррентные формулы, вычисляющие одновременно синус и косинус многих углов, где влияние ошибок округления уменьшено.

                                                  Only users with full accounts can post comments. Log in, please.