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

Ответ 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, что на удивление неточно для некоторых входов. Там увлекательное сообщение в блоге, иллюстрирующее это всего двумя строками кода.

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

Ответ 2

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

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

Ответ 3

OK детишки, время для профессионалов.... Это одна из моих самых больших жалоб у неопытных инженеров-программистов. Они приходят в вычисление трансцендентных функций с нуля (с использованием серии Тейлора), как если бы никто никогда не делал эти вычисления раньше в своей жизни. Не правда. Это хорошо определенная проблема, и тысячи раз были подобраны очень умными инженерами по программному обеспечению и аппаратным средствам и имеют четко определенное решение. В основном, большинство трансцендентных функций используют Чебышевские многочлены для их вычисления. Что касается того, какие полиномы используются, зависит от обстоятельств. Во-первых, Библия по этому вопросу - книга под названием "Компьютерные аппроксимации" Харта и Чейни. В этой книге вы можете решить, есть ли у вас аппаратный сумматор, множитель, разделитель и т.д., И решить, какие операции наиболее быстрые. например Если бы у вас был очень быстрый делитель, самым быстрым способом вычисления синуса может быть 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 обычно связаны с некоторыми микрокодами и такими операциями (в зависимости от аппаратного обеспечения). Многочлены Чебышева используются для большинства трансцендентных, но не для всех. например Квадратный корень быстрее использует двойную итерацию метода Ньютона raphson, используя первую таблицу поиска. Опять же, эта книга "Компьютерные аппроксимации" расскажет вам об этом.

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

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

Enjoy.

Ответ 4

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

sin (x): = x - x ^ 3/3! + x ^ 5/5! - x ^ 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;

}

Ответ 5

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

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

Ответ 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)

Ответ 7

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

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

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

Ответ 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 y будет меньше, чем sin 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) с полиномом (b x ^ 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 раз больше последнего бита.

Ответ 9

Что касается тригонометрических функций, таких как 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+ и др.) И предоставляет алгоритм решения, который дает качественные результаты для всех 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()

Ответ 10

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

Я понимаю, что абсолютно никакой помощи.

Ответ 11

Как правило, они реализованы в программном обеспечении и в большинстве случаев не будут использовать соответствующие аппаратные (то есть сборочные) вызовы. Однако, как отметил Джейсон, это зависит от реализации.

Обратите внимание, что эти программные подпрограммы не являются частью исходных текстов компилятора, а скорее находятся в соответствующей библиотеке, такой как clib или glibc для компилятора GNU. См. Http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions.

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

Ответ 12

Я попытаюсь ответить на случай sin() в программе на языке C, скомпилированной с компилятором GCC C на текущем процессоре x86 (скажем, Intel Core 2 Duo).

В языке C в стандартную библиотеку C входят общие математические функции, не включенные в сам язык (например, pow, sin и cos для питания, синуса и косинуса соответственно). Заголовки которых включены в math.h.

Теперь в системе GNU/Linux эти функции библиотеки предоставляются glibc (GNU libc или GNU C Library). Но компилятор GCC хочет, чтобы вы ссылались на математическую библиотеку (libm.so) с использованием флага компилятора -lm, чтобы включить использование этих математических функции. Я не уверен, почему он не является частью стандартной библиотеки C. Это будет программная версия функций с плавающей запятой или "soft-float".

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

Теперь компилятор может оптимизировать стандартную библиотечную функцию C sin() (предоставляемую libm.so), которая будет заменена вызовом собственной инструкции для встроенной функции sin() для CPU/FPU, которая существует как FPU (FSIN для x86/x87) на более новых процессорах, таких как Core 2 series (это довольно точно еще в i486DX). Это будет зависеть от флагов оптимизации, переданных компилятору gcc. Если компилятору было предложено написать код, который будет выполняться на любом i386 или более новом процессоре, он не будет делать такую ​​оптимизацию. Флаг -mcpu=486 должен сообщить компилятору, что безопасно сделать такую ​​оптимизацию.

Теперь, если программа выполнила версию программного обеспечения функции sin(), она будет делать это на основе CORDIC (компьютер с коррекцией по кругообороту) или алгоритм BKM, или, более вероятно, вычисление таблицы или мощности, которое обычно используется для расчета таких трансцендентных функций. [Src: http://en.wikipedia.org/wiki/Cordic#Application]

В любой недавней (начиная с версии 2.9x) версии gcc также имеется встроенная версия sin, __builtin_sin(), которая будет использоваться для замены стандартного вызова в библиотеке C-библиотеки в качестве оптимизации.

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

Ответ 13

Если вам нужна реализация в программном обеспечении, а не в аппаратном обеспечении, вам нужно найти окончательный ответ на этот вопрос в главе 5 " Числовые рецепты". Моя копия находится в коробке, поэтому я не могу дать подробности, но короткая версия (если я правильно помню это) состоит в том, что вы берете tan(theta/2) качестве своей примитивной операции и вычисляете другие оттуда. Вычисление выполняется в приближении ряда, но оно сходится гораздо быстрее, чем ряд Тейлора.

Извините, я не могу вспомнить больше, не положив руку на книгу.

Ответ 14

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

• Загрузите исходный код glibc с http://ftp.gnu.org/gnu/glibc/
• Посмотрите на файл dosincos.c находящийся в распакованной папке glibc root\sysdeps\ieee754\dbl-64
• Точно так же вы можете найти реализации остальной части математической библиотеки, просто найдите файл с соответствующим именем

Вы также можете взглянуть на файлы с расширением .tbl, их содержимое - не что иное, как огромные таблицы предварительно вычисленных значений различных функций в двоичном виде. Вот почему реализация такая быстрая: вместо того, чтобы вычислять все коэффициенты какой-либо серии, которую они используют, они просто выполняют быстрый поиск, который намного быстрее. Кстати, они используют ряды Tailor для вычисления синуса и косинуса.

Надеюсь, это поможет.

Ответ 15

Нет ничего похожего на то, чтобы поразить источник и увидеть, как кто-то действительно сделал это в библиотеке общего пользования; рассмотрим, в частности, одну реализацию библиотеки C. Я выбрал uLibC.

Здесь функция sin:

http://git.uclibc.org/uClibc/tree/libm/s_sin.c

который выглядит так, как будто он обрабатывает несколько особых случаев, а затем выполняет некоторую редукцию аргумента, чтобы отобразить вход в диапазон [-pi/4, pi/4] (разделение аргумента на две части, большую часть и хвост) перед вызовом

http://git.uclibc.org/uClibc/tree/libm/k_sin.c

который затем работает на этих двух частях. Если хвост отсутствует, то приближенный ответ генерируется с использованием полинома степени 13. Если есть хвост, вы получаете небольшое корректирующее дополнение, основанное на принципе, что sin(x+y) = sin(x) + sin'(x')y

Ответ 16

Всякий раз, когда такая функция оценивается, то на каком-то уровне есть, скорее всего, либо:

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

Если нет поддержки аппаратного обеспечения, тогда компилятор, вероятно, использует последний метод, испускающий только код ассемблера (без символов отладки), вместо использования библиотеки ac, что затрудняет отслеживание фактического кода в вашем отладчик.

Ответ 17

Если вы хотите посмотреть фактическую реализацию GNU этих функций на C, посмотрите последнюю строку glibc. См. Библиотека GNU C.

Ответ 18

Не используйте серии Тейлора. Многочлены Чебышева являются более быстрыми и точными, о чем указывает несколько человек выше. Вот реализация (первоначально из ROM ZX Spectrum): https://albertveli.wordpress.com/2015/01/10/zx-sine/

Ответ 19

Вычисление синуса/косинуса/тангенса на самом деле очень легко сделать с помощью кода с использованием ряда Тейлора. Написание одного занимает около 5 секунд.

Весь процесс можно суммировать с помощью этого уравнения здесь:

Вот некоторые подпрограммы, которые я написал для C:

double _pow(double a, double b) {
    double c = 1;
    for (int i=0; i<b; i++)
        c *= a;
    return c;
}

double _fact(double x) {
    double ret = 1;
    for (int i=1; i<=x; i++) 
        ret *= i;
    return ret;
}

double _sin(double x) {
    double y = x;
    double s = -1;
    for (int i=3; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _cos(double x) {
    double y = 1;
    double s = -1;
    for (int i=2; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _tan(double x) {
     return (_sin(x)/_cos(x));  
}

Ответ 20

Улучшенная версия кода от ответа Blindy

#define EPSILON .0000000000001
// this is smallest effective threshold, at least on my OS (WSL ubuntu 18)
// possibly because factorial part turns 0 at some point
// and it happens faster then series element turns 0;
// validation was made against sin() from <math.h>
double ft_sin(double x)
{
    int k = 2;
    double r = x;
    double acc = 1;
    double den = 1;
    double num = x;

//  precision drops rapidly when x is not close to 0
//  so move x to 0 as close as possible
    while (x > PI)
        x -= PI;
    while (x < -PI)
        x += PI;
    if (x > PI / 2)
        return (ft_sin(PI - x));
    if (x < -PI / 2)
        return (ft_sin(-PI - x));
//  not using fabs for performance reasons
    while (acc > EPSILON || acc < -EPSILON)
    {
        num *= -x * x;
        den *= k * (k + 1);
        acc = num / den;
        r += acc;
        k += 2;
    }
    return (r);
}

Ответ 21

если хочешь sin то

 __asm__ __volatile__("fsin" : "=t"(vsin) : "0"(xrads));

если вы хотите, cos тогда

 __asm__ __volatile__("fcos" : "=t"(vcos) : "0"(xrads));

если хочешь sqrt то

 __asm__ __volatile__("fsqrt" : "=t"(vsqrt) : "0"(value));

так зачем использовать неточный код, когда машинные инструкции подойдут?