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

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

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

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

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

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

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

• Когда х очень очень близко к 0, sin(x) == xэто правильный ответ.

• Чуть дальше sin(x)используется знакомая серия Тейлора. Тем не менее, это точно только около 0, так что ...

• Когда угол больше, чем приблизительно 7 °, используется другой алгоритм, вычисляющий приближения ряда Тейлора для sin (x) и cos (x), а затем используя значения из предварительно вычисленной таблицы для уточнения аппроксимации.

• Когда | х | > 2, ни один из вышеперечисленных алгоритмов не будет работать, поэтому код начинается с вычисления некоторого значения, близкого к 0, которое может быть передано sinили cosвместо него.

• Есть еще одна ветвь, чтобы иметь дело с х, являющимся NaN или бесконечностью.

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

double t = (x * hpinv + toint);
double xn = t - toint;

используются (иногда) при уменьшении x до значения, близкого к 0, которое отличается от x кратным π / 2, в частности xn× π / 2. То, как это делается без деления или ветвления, довольно умно. Но нет комментариев вообще!

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

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

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

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

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

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

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

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

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 )

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

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

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

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

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

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

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

Что касается тригонометрических функций , как 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 ()

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

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

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

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

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

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

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

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

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

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

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

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

• Загрузите исходный код 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 хочет, чтобы вы связались с библиотекой math ( libm.so), используя -lmфлаг компилятора, чтобы разрешить использование этих математических функций.Я не уверен, почему это не является частью стандартной библиотеки C. Это будет программная версия функций с плавающей запятой или «мягкое плавание».

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

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

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

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

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

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

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

Вычисление синуса / косинуса / тангенса на самом деле очень легко сделать с помощью кода с использованием ряда Тейлора. Написание одного занимает около 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));  
}

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

#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);
}

Суть того, как это происходит, заключается в этом отрывке из прикладного численного анализа. Джеральда Уитли:

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

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

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

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

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

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

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

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

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