Как C вычисляет sin () и другие математические функции? - PullRequest
226 голосов
/ 18 февраля 2010

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

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


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

Ответы [ 20 ]

5 голосов
/ 18 февраля 2010

Я попытаюсь ответить для случая 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) для замены вызовом встроенной инструкции встроенной функции sin () вашего CPU / FPU, которая существует как Инструкция FPU (FSIN для x86 / x87) на более новых процессорах, таких как серия Core 2 (это верно в значительной степени еще в i486DX). Это будет зависеть от флагов оптимизации, передаваемых компилятору gcc. Если компилятору было приказано написать код, который будет выполняться на любом процессоре i386 или новее, он не будет выполнять такую ​​оптимизацию. Флаг -mcpu=486 сообщит компилятору, что такая оптимизация безопасна.

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

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

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

5 голосов
/ 18 июля 2015

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

5 голосов
/ 18 февраля 2010

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

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

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

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

5 голосов
/ 18 февраля 2010

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

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

5 голосов
/ 18 февраля 2010

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

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

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

4 голосов
/ 18 февраля 2010

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

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

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

3 голосов
/ 18 февраля 2010

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

2 голосов
/ 18 июля 2015

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

1 голос
/ 03 апреля 2014

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

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

sin and cost expansion

Вот некоторые подпрограммы, которые я написал для 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));  
}
0 голосов
/ 20 марта 2015

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

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

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

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

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

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

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

...