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

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

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

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

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

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

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

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

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

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

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

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

вот пример для cosinus:

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)

объяснение http://img514.imageshack.us/img514/1959/82098830.jpg

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

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

Далее вы должны выяснить, как оценить полинома. Вы хотите оценить его таким образом, чтобы промежуточные значения были малы, и поэтому ошибки округления малы. В противном случае ошибки округления могут стать намного больше, чем ошибки в многочлене. И с такими функциями, как функция синуса, если вы небрежны, то возможно, что результат, который вы вычисляете для греха 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, х^5, х^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) является функция, П (Х) является полиномом, и Е (Х) - ошибка округления. Вы не хотите оптимизировать | f (x) - P (x) |, вы хотите оптимизировать | f (x) - P (x) +/- E (x) |. Вы получите немного другой полином, который пытается сохранить полиномиальные ошибки там, где ошибка округления велика, и немного ослабляет полиномиальные ошибки там, где ошибка округления мала.

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

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

sin (x): = x - x^3/3! + x^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) работает из-за aproximation 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, PC, HP, 30+ other) и предоставляет алгоритм решения, который дает качественные результаты для все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 ()

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

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

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

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

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

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

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

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

Я надеюсь, что это помогает.

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

в языке C стандартная библиотека C включает в себя общие математические функции, не включенные в сам язык (например,pow,sin и cos для мощности, синуса и Косинуса соответственно). Заголовки которых включены в математика.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.9 X прибл.) версия gcc также предлагает встроенную версию sin,__builtin_sin() что он будет использоваться для замены стандартного вызова версии библиотеки C в качестве оптимизации.

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

нет ничего лучше, чем попасть в источник и увидеть, как кто-то действительно сделал это в библиотеке общего пользования; давайте рассмотрим одну реализацию библиотеки 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

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

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

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

Если вы хотите посмотреть на фактическую реализацию GNU этих функций в C, проверьте последнюю магистраль 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));  
}

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

Если ты хочешь греха, то asm volatile("fsin": "=t"(vsin) : "0"(xrads)); если вы хотите, потому что asm volatile("fcos": "=t"(VCO) : "0"(xrads)); если вы хотите корень после asm volatile("fsqrt": "=t "(vsqrt) : "0"(значение)); так зачем использовать неточный код, когда машина инструкции будут делать.