Как C вычисляет sin () и другие математические функции?

Question

Я просматривал разборки .NET и исходный код GCC, но, похоже, нигде не могу найти фактическую реализацию sin() и других математических функций ... кажется, что они всегда ссылаются на что-то еще.

Может кто-нибудь помочь мне найти их? Я чувствую, что маловероятно, что ВСЕ оборудование, на котором будет работать C, поддерживает аппаратные функции триггера, поэтому должен быть программный алгоритм где-то , верно?

---

Мне известно о нескольких способах вычисления функций , и я написал свои собственные процедуры для вычисления функций с использованием рядов Тейлора для забавы. Мне любопытно, насколько реальны это производственные языки, так как все мои реализации всегда на несколько порядков медленнее, хотя я думаю, что мои алгоритмы довольно умны (очевидно, это не так).

Answer 1

В GNU libm реализация sin зависит от системы. Поэтому вы можете найти реализацию для каждой платформы где-нибудь в соответствующем подкаталоге sysdeps .

Один каталог включает реализацию на C, предоставленную IBM. С октября 2011 года этот код фактически запускается при вызове sin() в типичной системе Linux x86-64. Это, очевидно, быстрее, чем инструкция по сборке fsin. Исходный код: sysdeps / ieee754 / dbl-64 / s_sin.c , ищите __sin (double x).

Этот код очень сложный. Ни один программный алгоритм не является настолько быстрым, насколько это возможно, а также точным во всем диапазоне значений x , поэтому библиотека реализует множество различных алгоритмов, и ее первой задачей является поиск x и принятие решения. какой алгоритм использовать. В некоторых регионах он использует то, что похоже на знакомую серию Тейлора. Некоторые из алгоритмов сначала вычисляют быстрый результат, затем, если он недостаточно точен, отбрасывают его и прибегают к более медленному алгоритму.

В старых 32-разрядных версиях GCC / glibc использовалась инструкция fsin, что на удивление неточно для некоторых входных данных. Там есть увлекательное сообщение в блоге, иллюстрирующее это всего двумя строками кода .

Реализация sin в чистом C в fdlibm намного проще, чем в glibc, и хорошо прокомментирована. Исходный код: fdlibm / s_sin.c и fdlibm / k_sin.c

Answer 2

ОК, детишки, время для профи .... Это одна из моих самых больших претензий к неопытным разработчикам программного обеспечения. Они приходят к вычислению трансцендентных функций с нуля (используя ряды Тейлора), как будто никто никогда не делал эти вычисления раньше в своей жизни. Не правда. Это хорошо определенная проблема, к которой тысячи раз подходили очень умные разработчики программного и аппаратного обеспечения, и она имеет четко определенное решение. В основном, большинство трансцендентных функций используют полиномы Чебышева для их вычисления. То, какие полиномы используются, зависит от обстоятельств. Во-первых, библией по этому вопросу является книга под названием «Компьютерные приближения» Харта и Чейни. В этой книге вы можете решить, есть ли у вас аппаратный сумматор, множитель, делитель и т. Д., И решить, какие операции выполняются быстрее всего. например Если бы у вас был действительно быстрый делитель, самый быстрый способ вычисления синуса мог бы быть P1 (x) / P2 (x), где P1, P2 - полиномы Чебышева. Без быстрого делителя это может быть просто P (x), где P имеет гораздо больше членов, чем P1 или P2 .... так что это будет медленнее. Итак, первый шаг - определить ваше оборудование и его возможности. Затем вы выбираете подходящую комбинацию полиномов Чебышева (обычно для косинуса, например, cos (ax) = aP (x), опять же, где P полином Чебышева). Затем вы решаете, какую десятичную точность вы хотите. например если вам нужна точность в 7 цифр, вы посмотрите это в соответствующей таблице в упомянутой мной книге, и она даст вам (для точности = 7,33) число N = 4 и полиномиальное число 3502. N - это порядок полинома. (так что это p4.x ^ 4 + p3.x ^ 3 + p2.x ^ 2 + p1.x + p0), потому что N = 4. Затем вы ищите фактическое значение значений p4, p3, p2, p1, p0 в конце книги под 3502 (они будут в плавающей запятой). Затем вы реализуете свой алгоритм в программном обеспечении в виде: (((p4.x + p3) .x + p2) .x + p1) .x + p0 .... и это как вычислить косинус до 7 десятичных знаков на этом оборудовании.

Обратите внимание, что большинство аппаратных реализаций трансцендентных операций в FPU обычно включают в себя некоторый микрокод и подобные операции (зависит от аппаратного обеспечения). Полиномы Чебышева используются для большинства трансцендентных, но не для всех. например Квадратный корень быстрее использовать двойную итерацию метода Ньютона-Рафсона, используя сначала таблицу поиска. Опять же, эта книга «Компьютерные аппроксимации» скажет вам это.

Если вы планируете реализовать эти функции, я рекомендую всем, кто получит копию этой книги. Это действительно Библия для таких алгоритмов. Обратите внимание, что есть множество альтернативных средств для вычисления этих значений, таких как кордика и т. Д., Но они, как правило, лучше всего подходят для конкретных алгоритмов, где требуется только низкая точность. Чтобы гарантировать точность каждый раз, полиномы Чебышева - это путь. Как я уже сказал, хорошо определенная проблема. Был решен в течение 50 лет ..... и вот как это делается.

Теперь, как уже говорилось, существуют методы, с помощью которых полиномы Чебышева могут использоваться для получения результата с одинарной точностью и полиномом низкой степени (как в примере с косинусом выше). Затем, есть другие методы для интерполяции между значениями, чтобы увеличить точность без необходимости переходить к многочлену большего размера, например, «метод точных таблиц Гала». Этот последний метод - то, на что ссылается пост, ссылающийся на литературу ACM. Но в конечном итоге полиномы Чебышева - это то, что используется для 90% пути туда.

Наслаждайтесь.

Answer 3

Такие функции, как синус и косинус, реализованы в микрокоде внутри микропроцессора.Например, у чипов Intel есть инструкции по сборке.Компилятор AC сгенерирует код, который вызывает эти инструкции по сборке.(В отличие от этого, компилятор Java не будет. Java оценивает функции триггера в программном, а не аппаратном обеспечении, и поэтому работает намного медленнее.)

Чипы не используют ряды Тейлора для вычисления функций триггераПо крайней мере, не совсем.Прежде всего они используют CORDIC , но они также могут использовать короткий ряд Тейлора для полировки результата CORDIC или для особых случаев, таких как вычисление синуса с высокой относительной точностью для очень малых углов.Дополнительную информацию см. В этом ответе StackOverflow .

.

Answer 4

Для sin, в частности, использование расширения Тейлора даст вам:

грех (х): = х - х ^ 3/3! + х ^ 5/5! - х ^ 7/7! + ... (1)

Вы будете продолжать добавлять термины, пока разница между ними не станет ниже допустимого уровня допуска или не будет достигнута только для конечного количества шагов (быстрее, но менее точно). Примером будет что-то вроде:

float sin(float x)
{
  float res=0, pow=x, fact=1;
  for(int i=0; i<5; ++i)
  {
    res+=pow/fact;
    pow*=-1*x*x;
    fact*=(2*(i+1))*(2*(i+1)+1);
  }

  return res;
}

Примечание: (1) работает из-за апроксимации sin (x) = x для малых углов. Для больших углов вам нужно вычислить все больше и больше терминов, чтобы получить приемлемые результаты. Вы можете использовать аргумент while и продолжить с определенной точностью:

double sin (double x){
    int i = 1;
    double cur = x;
    double acc = 1;
    double fact= 1;
    double pow = x;
    while (fabs(acc) > .00000001 &&   i < 100){
        fact *= ((2*i)*(2*i+1));
        pow *= -1 * x*x; 
        acc =  pow / fact;
        cur += acc;
        i++;
    }
    return cur;

}

Answer 5

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

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

Answer 6

Используйте ряд Тейлора и попытайтесь найти связь между терминами ряда, чтобы не вычислять вещи снова и снова

Вот пример косинуса:

double cosinus(double x, double prec)
{
    double t, s ;
    int p;
    p = 0;
    s = 1.0;
    t = 1.0;
    while(fabs(t/s) > prec)
    {
        p++;
        t = (-t * x * x) / ((2 * p - 1) * (2 * p));
        s += t;
    }
    return s;
}

используя это, мы можем получить новый член суммы, используя уже использованный (мы избегаем факториала и x ^2p )

Answer 7

Это сложный вопрос. Intel-подобные процессоры семейства x86 имеют аппаратную реализацию функции sin(), но они являются частью FPU x87 и больше не используются в 64-битном режиме (где вместо этого используются регистры SSE2). В этом режиме используется программная реализация.

Существует несколько таких реализаций. Один из них находится в fdlibm и используется в Java. Насколько я знаю, реализация glibc содержит части fdlibm и другие части, предоставленные IBM.

Программные реализации трансцендентных функций, такие как sin(), обычно используют приближения полиномами, часто получаемые из ряда Тейлора.

Answer 8

Полиномы Чебышева, как упоминалось в другом ответе, являются полиномами, где наибольшая разница между функцией и полиномом настолько мала, насколько это возможно.Это отличное начало.

В некоторых случаях максимальная ошибка - это не то, что вас интересует, а максимальная относительная ошибка.Например, для функции синуса ошибка около x = 0 должна быть намного меньше, чем для больших значений;Вы хотите небольшую относительную ошибку.Таким образом, вы бы вычислили многочлен Чебышева для sin x / x и умножили этот многочлен на x.

Далее вы должны выяснить, как оценивать полином.Вы хотите оценить его таким образом, чтобы промежуточные значения были небольшими и, следовательно, ошибки округления были небольшими.В противном случае ошибки округления могут стать намного больше, чем ошибки в полиноме.А с такими функциями, как функция синуса, если вы небрежны, то может оказаться, что результат, который вы рассчитываете для sin x, больше, чем результат для sin y, даже если x

Например, sin x = x - x ^ 3/6 + x ^ 5/120 - x ^ 7/5040 ... Если вычислить наивно, sin x = x * (1 - x ^ 2 /6 + x ^ 4/120 - x ^ 6/5040 ...), то эта функция в скобках уменьшается, и будет , что если y будет следующим большим числом после x, то иногда греху будет меньше греха х.Вместо этого вычислите sin x = x - x ^ 3 * (1/6 - x ^ 2/120 + x ^ 4/5040 ...), где это не может произойти.

При вычислении полиномов Чебышева, например, обычно необходимо округлять коэффициенты для удвоения точности.Но хотя полином Чебышева является оптимальным, полином Чебышева с коэффициентами, округленными до двойной точности, не является оптимальным полиномом с коэффициентами двойной точности!

Например, для sin (x), где вам нужны коэффициенты для x, x ^ 3, x ^ 5, x ^ 7 и т. Д., Вы делаете следующее: Рассчитываете наилучшее приближение sin x с помощью полинома (ax + bx ^ 3 + cx ^ 5 + dx ^ 7) с точностью выше двойной, затем округляем до двойной точности, получая A. Разница между a и A будет довольно большой.Теперь вычислите наилучшее приближение (sin x - Ax) с полиномом (bx ^ 3 + cx ^ 5 + dx ^ 7).Вы получаете разные коэффициенты, потому что они адаптируются к разности между a и A. Округлите b до двойной точности B. Затем аппроксимируйте (sin x - Ax - Bx ^ 3) полиномом cx ^ 5 + dx ^ 7 и так далее.Вы получите многочлен, который почти так же хорош, как и исходный полином Чебышева, но намного лучше, чем Чебышев, округленный до двойной точности.

Далее следует учесть ошибки округления при выборе полинома.Вы нашли полином с минимальной ошибкой в ​​полиноме, игнорирующем ошибку округления, но вы хотите оптимизировать полином плюс ошибка округления.Получив полином Чебышева, вы можете вычислить оценки для ошибки округления.Скажем, f (x) - ваша функция, P (x) - полином, а E (x) - ошибка округления.Вы не хотите оптимизировать |f (x) - P (x) |, вы хотите оптимизировать |f (x) - P (x) +/- E (x) |.Вы получите немного другой полином, который попытается сохранить полиномиальные ошибки там, где ошибка округления велика, и немного ослабит полиномиальные ошибки там, где ошибка округления мала.

Все это позволит вам легко округлять ошибки не более чем в 0,55 раз больше последнего бита, где +, -, *, / имеют ошибки округления не более чем в 0,50 раз больше последнего разряда.

Answer 9

Реальная реализация библиотечных функций зависит от конкретного компилятора и / или поставщика библиотеки. Будет ли это сделано в аппаратном или программном обеспечении, будь то расширение Тейлора или нет, и т. Д., Будет отличаться.

Я понимаю, что это абсолютно не поможет.

Answer 10

Относительно тригонометрических функций, таких как sin(), cos(), tan(), после 5 лет не упоминалось о важном аспекте высококачественных функций триггера: Уменьшение диапазона .

Первым шагом в любой из этих функций является уменьшение угла в радианах до диапазона 2 * π-интервала. Но π иррационально, поэтому простые сокращения, такие как x = remainder(x, 2*M_PI), вносят ошибку, поскольку M_PI, или машина pi, является приближением π. Итак, как сделать x = remainder(x, 2*π)?

Ранние библиотеки использовали расширенную точность или специально разработанное программирование для получения качественных результатов, но все еще в ограниченном диапазоне double. Когда запрашивалось большое значение, например sin(pow(2,30)), результаты были бессмысленными или 0.0 и, возможно, с флагом ошибки , установленным на что-то вроде TLOSS полная потеря точности или PLOSS частичная потеря точности .

Хорошее сокращение диапазона больших значений до интервала, подобного -π до π, является сложной задачей, которая конкурирует с вызовами базовой функции триггера, такой как sin(), сама по себе.

Хороший отчет - Сокращение аргументов для огромных аргументов: от хорошего до последнего бита (1992). Он хорошо освещает проблему: обсуждает необходимость и состояние дел на различных платформах (SPARC, ПК, HP, 30+ другие) и предоставляет алгоритм решения, который дает качественные результаты для all double с -DBL_MAX до DBL_MAX.

---

Если исходные аргументы выражены в градусах, но могут иметь большое значение, сначала используйте fmod() для повышения точности. Хороший fmod() введет без ошибок и, таким образом, обеспечит превосходное уменьшение диапазона.

// sin(degrees2radians(x))
sin(degrees2radians(fmod(x, 360.0))); // -360.0 <= fmod(x,360) <= +360.0

Различные триггеры и remquo() предлагают еще больше улучшений. Образец: sind ()