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

248

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

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


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

моток
источник
2
Обратите внимание, что эта реализация зависит. Вы должны указать, какая реализация вас больше всего интересует.
Джейсон
3
Я пометил .NET и C, потому что я смотрел в обоих местах и ​​не мог понять, либо. Хотя, глядя на разборку .NET, она выглядит так, как будто она вызывает неуправляемый C, так что, насколько я знаю, они имеют одинаковую реализацию.
Хэнк

Ответы:

213

В 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

Джейсон Орендорфф
источник
35
Чтобы увидеть, что это действительно код, работающий на x86: скомпилируйте вызывающую программу sin(); типа gdb a.out, тогда break sin, потом run, потом disassemble.
Джейсон Орендорф
5
@ Генри: не делайте ошибку, думая, что это хороший код, хотя. Это действительно ужасно , не учитесь так кодировать!
Томас Бонини
2
@ Андреас Хмм, вы правы, код IBM выглядит ужасно по сравнению с fdlibm. Я отредактировал ответ, добавив ссылки на синус рутины fdlibm.
Джейсон Орендорфф
3
@Henry: __kernel_sinопределен в k_sin.c, и, тем не менее, он чистый C. Нажмите его еще раз - я запустил URL в первый раз.
Джейсон Орендорф
3
Связанный код sysdeps особенно интересен, потому что он правильно округлен. То есть, по-видимому, он дает наилучший ответ для всех входных значений, что стало возможным только в последнее время. В некоторых случаях это может быть медленным, поскольку для обеспечения правильного округления может потребоваться вычисление многих дополнительных цифр. В других случаях это очень быстро - для достаточно маленьких чисел ответ - только угол.
Брюс Доусон
67

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

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

Джон Д. Кук
источник
10
трансцендентные математические функции, такие как синус и косинус, могут быть реализованы в микрокоде или в виде аппаратных инструкций в современных 32-битных процессорах для настольных компьютеров и серверов. Это не всегда имело место до тех пор, пока в i486 (DX) все вычисления с плавающей запятой не были выполнены в программном обеспечении («soft-float») для серии x86 без отдельного сопроцессора. Не все из которых (FPU) включали трансцендентные функции (например, Weitek 3167).
mctylr
1
Можете быть более конкретными? Как «отшлифовать» приближение, используя ряд Тейлора?
Хэнк
4
Что касается «доводки» ответа, предположим, что вы вычисляете как синус, так и косинус. Предположим, вы знаете точное значение обоих в одной точке (например, из CORDIC), но хотите получить значение в соседней точке. Тогда для небольшой разности h можно применить приближения Тейлора f (x + h) = f (x) + h f '(x) или f (x + h) = f (x) + h f' (x) + h ^ 2 f '' (x) / 2.
Джон Д. Кук
6
Чипы x86 / x64 имеют инструкцию по сборке для вычисления синуса (fsin), но эта инструкция иногда довольно неточна и поэтому редко используется больше. Смотрите randomascii.wordpress.com/2014/10/09/… для подробностей. Большинство других процессоров не имеют инструкций для синуса и косинуса, потому что их вычисление в программном обеспечении дает большую гибкость и даже может быть быстрее.
Брюс Доусон
3
Чудесные вещи внутри чипов Intel обычно НЕ используются. Во-первых, точность и разрешение операции чрезвычайно важны для многих приложений. Cordic, как известно, неточно, когда вы добираетесь до 7-й цифры или около того, и непредсказуемо. Во-вторых, я слышал, что в их реализации есть ошибка, которая вызывает еще больше проблем. Я посмотрел на функцию sin для linux gcc, и, конечно же, она использует чебышев. встроенный материал не используется. О, также, сердечный алгоритм в чипе медленнее, чем программное решение.
Дональд Мюррей
63

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

Наслаждаться.

Дональд Мюррей
источник
6
Я не мог согласиться с первыми несколькими предложениями. Также стоит напомнить, что вычисление специальных функций с гарантированной точностью является сложной задачей . Умные люди, о которых вы говорите, проводят большую часть своей жизни, занимаясь этим. Кроме того, в более техническом примечании полиномы min-max являются искомыми граалями, а полиномы Чебышева являются для них более простыми прокси.
Александр С.
161
-1 за непрофессиональный и бессвязный (и слегка грубый) тон, а также за то, что фактическое не лишнее содержание этого ответа, лишенное бессвязных и снисходительных слов, в основном сводится к: «Они часто используют полиномы Чебышева; см. Эту книгу для более подробной информации, это действительно хорошо! Что, вы знаете, вполне может быть абсолютно правильным, но на самом деле это не тот самостоятельный ответ, который мы хотим здесь, на SO. Сжатая таким образом, она сделала бы достойный комментарий по этому вопросу.
Илмари Каронен
2
Еще в первые годы разработки игр это обычно делалось с помощью справочных таблиц, критически важных для скорости). Обычно мы не использовали стандартные функции lib для этих целей.
Topspin
4
Я часто использую таблицы поиска во встроенных системах и в битах (вместо радианов), но это для специализированного приложения (такого как ваши игры). Я думаю, что парню интересно, как компилятор c вычисляет грех для чисел с плавающей запятой ....
Дональд Мюррей
1
Ах, 50 лет назад. Я начал играть с такими на Burroughs B220 с серией McLaren. Позже аппаратное обеспечение CDC, а затем Motorola 68000. Arcsin был грязным - я выбрал частное двух полиномов и разработал код, чтобы найти оптимальные коэффициенты.
Рик Джеймс
15

В 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;

}
Blindy
источник
1
Если вы немного подкорректируете коэффициенты (и жестко закодируете их в полином), вы можете остановить примерно на 2 итерации раньше.
Рик Джеймс
14

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

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

Мехрдад Афшари
источник
12
Реальная реализация, вероятно, не будет использовать серию Тейлора, поскольку есть более эффективные способы. Вам нужно только правильно аппроксимировать в области [0 ... pi / 2], и есть функции, которые обеспечат хорошее приближение более эффективно, чем ряд Тейлора.
Дэвид Торнли
2
@ Давид: я согласен. Я был достаточно осторожен, чтобы упомянуть слово «как» в моем ответе. Но разложение Тейлора простое, чтобы объяснить идею методов, аппроксимирующих функции. Тем не менее, я видел программные реализации (не уверен, что они были оптимизированы), которые использовали серии Тейлор.
Мехрдад Афшари
1
На самом деле, полиномиальные аппроксимации являются одним из наиболее эффективных способов вычисления тригонометрических функций.
Джереми Сальвен
13

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

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

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 )

объяснение

Ханнун Ясир
источник
2
Знаете ли вы, что вы можете использовать API Google Chart для создания таких формул, используя TeX? code.google.com/apis/chart/docs/gallery/formulas.html
Габ Ройер,
11

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

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

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

Томас Порнин
источник
3
Регистры SSE2 не используются для вычисления sin () ни в режиме x86, ни в режиме x64, и, конечно, sin рассчитывается аппаратно независимо от режима. Эй, это 2010 год, в котором мы живем :)
Игорь Корхов
7
@Igor: это зависит от того, на какую математическую библиотеку вы смотрите. Оказывается, что наиболее оптимизированные математические библиотеки в x86 используют программные реализации SSE sinи cosработают быстрее, чем аппаратные инструкции на FPU. Проще, наивные библиотеки , как правило, используют fsinи fcosинструкцию.
Стивен Кэнон
@Stephen Canon: Эти быстрые библиотеки имеют 80-битную точность, как регистры FPU? У меня очень подлое подозрение, что они предпочитают скорость, а не точность, что, конечно, разумно во многих сценариях, например в играх. И я считаю, что вычисление синуса с 32-битной точностью с использованием SSE и предварительно вычисленных промежуточных таблиц может быть быстрее, чем с использованием FSINс полной точностью. Буду очень признателен, если вы скажете названия этих быстрых библиотек, интересно посмотреть.
Игорь Корхов
@Igor: на x86 в 64-битном режиме, по крайней мере во всех Unix-подобных системах, о которых я знаю, точность ограничена 64 битами, а не 79 битами x87 FPU. Программная реализация sin()оказывается в два раза быстрее, чем fsinвычисляется (именно потому, что это делается с меньшей точностью). Обратите внимание, что x87, как известно, имеет немного меньшую фактическую точность, чем заявленные 79 бит.
Томас Порнин
1
Действительно, как 32-битные, так и 64-битные реализации sin () в библиотеках времени выполнения msvc не используют инструкцию FSIN . Фактически, они дают разные результаты, например, грех (0.70444454416678126). Это приведет к 0,64761068800896837 (правильно с допуском 0,5 * (eps / 2)) в 32-битной программе и к 0,64761068800896848 (неверно) в 64-битной.
e.tadeu
9

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

В некоторых случаях максимальная ошибка - это не то, что вас интересует, а максимальная относительная ошибка. Например, для функции синуса ошибка около 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 по сравнению с последним битом.

gnasher729
источник
1
Это хорошее объяснение того , как один может вычислить грех (х) эффективно, но это на самом деле не кажется , ответить на вопрос OP, который является конкретно о том , как общие библиотеки C / Составителей этого расчета.
Илмари Каронен
Полиномы Чебышева минимизируют максимальное абсолютное значение за интервал, но они не сводят к минимуму наибольшую разницу между целевой функцией и полиномом. Минимаксные полиномы делают это.
Эрик Постпишил
9

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

Chux - Восстановить Монику
источник
6

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

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

Джон Боде
источник
5

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

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

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

mnemosyn
источник
5

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

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

Норман Рэмси
источник
5

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

Moschops
источник
4

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

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

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

Джеймс
источник
4

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

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

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

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

Игорь Корхов
источник
4

Я попытаюсь ответить для случая 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, в качестве оптимизации.

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

mctylr
источник
3

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

Крис Тонкинсон
источник
3

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

Альберт Вели
источник
2
Похоже, это не отвечает на вопрос, который задают. OP спрашивает , как Триг функция будет вычислена общими компиляторы C / библиотеками (и я уверен , что ZX Spectrum не отвечает), а не как они должны быть рассчитаны. Это могло быть полезным комментарием к некоторым из более ранних ответов, все же.
Илмари Каронен
1
Ах, ты прав. Это должен был быть комментарий, а не ответ. Я давно не пользовался SO и забыл, как работает система. В любом случае, я думаю, что реализация Spectrum актуальна, потому что у нее был очень медленный процессор, и скорость была существенной. Тогда лучший алгоритм, безусловно, все еще довольно хорош, поэтому для библиотек C было бы неплохо реализовать тригонометрические функции с использованием полиномов Чебышева.
Альберт Вели
2

Вычисление синуса / косинуса / тангенса на самом деле очень легко сделать с помощью кода с использованием ряда Тейлора. Написание одного занимает около 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));  
}
user1432532
источник
4
Это довольно плохая реализация, так как она не использует то, что последовательные члены ряда синуса и косинуса имеют очень простые отношения. Это означает, что можно уменьшить количество умножений и делений с O (n ^ 2) здесь до O (n). Дальнейшее сокращение достигается путем деления пополам и возведения в квадрат, как, например, это делается в математической библиотеке bc (POSIX multiprecision calculator).
Лутц Леманн
2
Это также, кажется, не отвечает на вопрос как задано; OP спрашивает, как функции trig рассчитываются общими компиляторами / библиотеками C, а не для пользовательских переопределений.
Илмари Каронен
2
Я думаю, что это хороший ответ, поскольку он отвечает духу вопроса, который (и я могу только догадываться, конечно) любопытство по поводу функции «черного ящика», такой как sin (). Это единственный ответ, который дает здесь возможность быстро понять, что происходит, замаскировав его за несколько секунд, а не читая оптимизированный исходный код на Си.
Майк М
на самом деле библиотеки используют гораздо более оптимизированную версию, понимая, что, получив термин, вы можете получить следующий, умножив некоторые значения. Смотрите пример в ответе Блинди . Вы снова и снова вычисляете мощность и факториалы, что намного медленнее
phuclv
0

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

#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);
}
Rugnar
источник
0

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

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

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

Дин П
источник
-1

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

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

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

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

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

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

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

user80998
источник