Как C вычисляет sin() и другие математические функции?
Я просматривал разборки.NET и исходный код GCC, но, похоже, нигде не нашел реальной реализации sin() и другие математические функции... они всегда, кажется, ссылаются на что-то еще.
Может ли кто-нибудь помочь мне найти их? Я чувствую, что маловероятно, что ВСЕ оборудование, на котором будет работать C, поддерживает аппаратные функции триггера, поэтому где-то должен быть программный алгоритм, верно?
Мне известно о нескольких способах вычисления функций, и я написал свои собственные процедуры для вычисления функций, используя ряды Тейлора для развлечения. Мне любопытно, насколько реальны это производственные языки, так как все мои реализации всегда на несколько порядков медленнее, хотя я думаю, что мои алгоритмы довольно умны (очевидно, что нет).
20 ответов
В 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
Такие функции, как синус и косинус, реализованы в микрокоде внутри микропроцессора. Например, у чипов Intel есть инструкции по их сборке. Компилятор A C сгенерирует код, который вызывает эти инструкции по сборке. (Напротив, компилятор Java не будет. Java оценивает функции триггера в программном, а не аппаратном обеспечении, и поэтому работает намного медленнее.)
Чипы не используют ряды Тейлора для вычисления тригонометрических функций, по крайней мере, не полностью. Прежде всего, они используют CORDIC, но они также могут использовать короткие ряды Тейлора для полировки результата CORDIC или для особых случаев, таких как вычисление синуса с высокой относительной точностью для очень малых углов. Для получения дополнительной информации см. Этот ответ Stackru.
ОК, детки, время для профессионалов... Это одна из моих самых больших претензий к неопытным разработчикам программного обеспечения. Они приходят к вычислению трансцендентных функций с нуля (используя ряды Тейлора), как будто никто никогда не делал эти вычисления раньше в своей жизни. Не правда. Это хорошо определенная проблема, к которой тысячи раз подходили очень умные разработчики программного и аппаратного обеспечения, и она имеет четко определенное решение. В основном, большинство трансцендентных функций используют полиномы Чебышева для их вычисления. То, какие полиномы используются, зависит от обстоятельств. Во-первых, библией по этому вопросу является книга под названием "Компьютерные приближения" Харта и Чейни. В этой книге вы можете решить, есть ли у вас аппаратный сумматор, множитель, делитель и т. Д., И решить, какие операции выполняются быстрее всего. Например, если у вас действительно быстрый делитель, самый быстрый способ вычислить синус может быть 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% пути.
Наслаждаться.
За 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;
}
Относительно тригонометрической функции типа 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 ()
Используйте ряды Тейлора и попытайтесь найти связь между терминами серии, чтобы вы не вычисляли вещи снова и снова
Вот пример косинуса:
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)
Да, есть программные алгоритмы для расчета sin тоже. По сути, вычисление такого рода вещей с помощью цифрового компьютера обычно выполняется с использованием численных методов, таких как аппроксимация ряда Тейлора, представляющего функцию.
Численные методы могут аппроксимировать функции с произвольной степенью точности, и поскольку степень точности, которую вы имеете в плавающем числе, конечна, они вполне подходят для этих задач.
Это сложный вопрос. Intel-подобные процессоры семейства x86 имеют аппаратную реализацию sin() функция, но она является частью x87 FPU и больше не используется в 64-битном режиме (где вместо этого используются регистры SSE2). В этом режиме используется программная реализация.
Существует несколько таких реализаций. Один из них находится в fdlibm и используется в Java. Насколько я знаю, реализация glibc содержит части fdlibm и другие части, предоставленные IBM.
Программные реализации трансцендентных функций, таких как sin() обычно используют приближения полиномами, часто полученными из рядов Тейлора.
Многочлены Чебышева, как уже упоминалось в другом ответе, являются многочленами, где наибольшая разница между функцией и многочленом настолько мала, насколько это возможно. Это отличное начало.
В некоторых случаях максимальная ошибка - это не то, что вас интересует, а максимальная относительная ошибка. Например, для функции синуса ошибка около 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 и 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 по сравнению с последним битом.
Реальная реализация библиотечных функций зависит от конкретного компилятора и / или поставщика библиотеки. Будет ли это сделано в аппаратном или программном обеспечении, будь то расширение Тейлора или нет, и т. Д., Будет отличаться.
Я понимаю, что это абсолютно не поможет.
Нет ничего лучше, чем поразить источник и увидеть, как кто-то на самом деле сделал это в библиотеке общего пользования; давайте посмотрим на одну реализацию библиотеки C в частности. Я выбрал uLibC.
Вот функция греха:
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
Как отмечали многие, это зависит от реализации. Но насколько я понимаю ваш вопрос, вы интересовались реальной программной реализацией математических функций, но просто не смогли найти ее. Если это так, то вот вы здесь:
- Загрузите исходный код glibc с http://ftp.gnu.org/gnu/glibc/
- Посмотри файл
dosincos.cнаходится в распакованной папке glibc root \ sysdeps \ ieee754 \ dbl-64 - Точно так же вы можете найти реализации остальной части математической библиотеки, просто найдите файл с соответствующим именем
Вы также можете посмотреть на файлы с .tbl Расширение, их содержимое - не что иное, как огромные таблицы предварительно вычисленных значений различных функций в двоичном виде. Вот почему реализация такая быстрая: вместо того, чтобы вычислять все коэффициенты какой-либо серии, которую они используют, они просто выполняют быстрый поиск, который намного быстрее. Кстати, они используют ряды Tailor для вычисления синуса и косинуса.
Надеюсь, это поможет.
Я постараюсь ответить по делу sin() в программе на C, скомпилированной с помощью компилятора C GCC на текущем процессоре 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. Это будет программная версия функций с плавающей запятой или "мягкое плавание".
Кроме того: причина, по которой математические функции были отделены друг от друга, является исторической и была просто предназначена для уменьшения размера исполняемых программ в очень старых системах Unix, возможно, до того, как стали доступны разделяемые библиотеки, насколько я знаю.
Теперь компилятор может оптимизировать стандартную функцию библиотеки C sin() (предоставлено libm.so) должен быть заменен вызовом собственной инструкции встроенной функции sin() вашего CPU/FPU, которая существует как инструкция FPU (FSIN для x86/x87) на более новых процессорах, таких как серия Core 2 (это верно почти во времена i486DX). Это будет зависеть от флагов оптимизации, передаваемых компилятору gcc. Если компилятору было приказано написать код, который будет выполняться на любом процессоре i386 или новее, он не будет выполнять такую оптимизацию. -mcpu=486 flag сообщит компилятору, что такая оптимизация безопасна.
Теперь, если бы программа выполняла версию программного обеспечения функции sin(), она делала бы это на основе алгоритма CORDIC (Coordinate Rotation DIgital Computer) или BKM, или, более вероятно, расчета таблицы или ряда степеней, который обычно используется сейчас для вычисления такие трансцендентные функции. [Источник: http://en.wikipedia.org/wiki/Cordic
Любая недавняя (начиная с 2.9x прим.) Версия gcc также предлагает встроенную версию sin, __builtin_sin() что он будет использоваться для замены стандартного вызова версии библиотеки C, в качестве оптимизации.
Я уверен, что это так же ясно, как грязь, но, надеюсь, даст вам больше информации, чем вы ожидали, и много прыжков от точек, чтобы узнать больше самостоятельно.
Если вам нужна реализация в программном обеспечении, а не в аппаратном обеспечении, вам следует найти окончательный ответ на этот вопрос в главе 5 " Числовые рецепты". Моя копия в коробке, поэтому я не могу дать подробности, но короткая версия (если я правильно помню), что вы берете tan(theta/2) как ваша примитивная операция и вычислить другие оттуда. Вычисление выполняется в приближении ряда, но оно сходится гораздо быстрее, чем ряд Тейлора.
Извините, я не могу вспомнить больше, не положив руку на книгу.
Как правило, они реализованы в программном обеспечении и в большинстве случаев не будут использовать соответствующие аппаратные (то есть сборочные) вызовы. Однако, как отметил Джейсон, это зависит от реализации.
Обратите внимание, что эти программные подпрограммы не являются частью исходных текстов компилятора, а скорее находятся в соответствующей библиотеке, такой как clib или g libc для компилятора GNU. См. http://www.gnu.org/software/libc/manual/html_mono/libc.html.
Если вы хотите большего контроля, вы должны тщательно оценить, что именно вам нужно. Некоторые из типичных методов - это интерполяция справочных таблиц, вызов ассемблера (который часто медленный) или другие схемы аппроксимации, такие как Ньютон-Рафсон для квадратных корней.
Не используйте серии Тейлор. Полиномы Чебышева и быстрее, и точнее, как отметили несколько человек выше. Вот реализация (первоначально из ROM ZX Spectrum): https://albertveli.wordpress.com/2015/01/10/zx-sine/
Всякий раз, когда такая функция оценивается, на некотором уровне, скорее всего, либо:
- Таблица значений, которая интерполируется (для быстрых, неточных приложений - например, компьютерной графики)
- Оценка ряда, который сходится к желаемому значению - вероятно, не ряд Тейлора, скорее что-то, основанное на причудливой квадратуре, такой как Кленшоу-Кертис.
Если аппаратная поддержка отсутствует, то компилятор, вероятно, использует последний метод, генерирующий только ассемблерный код (без символов отладки), вместо использования библиотеки ac, что затрудняет отслеживание фактического кода в отладчике.
Если вы хотите взглянуть на фактическую реализацию этих функций в C на GNU, посмотрите последнюю магистраль glibc. Смотрите библиотеку GNU C.
Вычисление синуса / косинуса / тангенса на самом деле очень легко сделать с помощью кода с использованием ряда Тейлора. Написание одного занимает около 5 секунд.
Весь процесс можно суммировать с помощью этого уравнения здесь: http://upload.wikimedia.org/math/5/4/6/546ecab719ce73dfb34a7496c942972b.png
Вот некоторые подпрограммы, которые я написал для 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));
}
Улучшенная версия кода из ответа Блинди
#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);
}
Суть того, как это происходит, заключается в этом отрывке из книги Джеральда Уитли " Прикладной численный анализ ":
Когда ваша программа запрашивает у компьютера значение или
, задавались ли вы вопросом, как он может получить значения, если самые мощные функции, которые он может вычислить, являются полиномами? Он не ищет их в таблицах и не интерполирует! Скорее, компьютер аппроксимирует каждую функцию, кроме полиномов от некоторого полинома, который настроен так, чтобы давать значения очень точно.
Несколько моментов, которые следует упомянуть по поводу вышеизложенного, заключается в том, что некоторые алгоритмы фактически интерполируют из таблицы, хотя и только для первых нескольких итераций. Также обратите внимание на то, как в нем упоминается, что компьютеры используют приближающие многочлены без указания типа приближающего многочлена. Как отмечали другие участники обсуждения, в этом случае полиномы Чебышева более эффективны, чем полиномы Тейлора.
Если вы хотите грешить, то asm volatile("fsin": "= t" (vsin): "0" (xrads)); если вы хотите cos, то asm volatile("fcos": "= t" (vcos): "0" (xrads)); если вы хотите sqrt, то asm volatile("fsqrt": "= t" (vsqrt): "0" (значение)); так зачем использовать неточный код, когда машинные инструкции подойдут.