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

Улучшенная версия кода из ответа 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);
}

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

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

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

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

Для 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;

}

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

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

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

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

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

Я попытаюсь ответить на случай 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. Это будет версия программного обеспечения функций с плавающей запятой или «soft-float».

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

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

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

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

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

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

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

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

Я понимаю, что это совершенно бесполезно.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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 использование расширения Тейлора даст вам:

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) работает из-за приближения 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;

}