Функция синуса в Си / Си++


Из-за программного ограничения, я не могу использовать стандартные библиотеки, cmath, алгоритмов, шаблонов, встроенных или увеличить. Я также использую стандарт C (по ISO С99) такие, что массив не зарезервированное ключевое слово, как это в Visual студии.

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

double factorial(double x)
{
    double result = 1;
    for (int i=1;i<=x;i++)
        result *= i;
    return result;
}

double pow(double x,double y)
{
    double result = 1;
    for (int i=0;i<y;i++)
        result *= x;
    return result;
}

double sine_taylor(double x,double n)
{
    double sine = x;
    bool isNeg;
    isNeg = true;
    for (double i=3;i<=n;i+=2)
    {
        if(isNeg)
        {
            sine -= pow(x,i)/factorial(i);
            isNeg = false;
        }
        else
        {
            sine += pow(x,i)/factorial(i);
            isNeg = true;
        }
    }
    return sine;
}

int main()
{
   const double pi = 3.14159265;
   double sineTaylor = sine_taylor(pi/2,7);
}


Комментарии
5 ответов

Обновление: (19-март-2015) этот ответ звучит как ответ эксперта.

Оригинальный ответ:

Есть книги, Полный о том, как эффективно и точно рассчитывать эти функции. Так это слишком много для короткого ответа здесь.

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

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

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

Так четко работа для специалистов. И действительно, решена давным-давно.

Быстрый Google обнаружили, что каждый ученый должен знать о плавающей точкой Arithmic.


Примечание – это приложение представляет собой отредактированное переиздание на бумаге то, что каждый
Ученый должен знать об арифметических операций с плавающей точкой, по
Дэвид Голдберг, опубликованном в марте, выпуск 1991 вычислительной техники
Опросы. Авторское право 1991, Ассоциации вычислительной техники, Инк.
Перепечатано с разрешения.

Я просто погуглил на это звание, пари, что кто-то бы использовал это название в какой-то момент времени. Это объясняет проблемы, которые я упомянул выше и гораздо больше, но я не квалифицирован, чтобы судить о нем. Как это документация SPARC от Sun, я ожидал, что это будет одними из лучших можно найти в интернете, хотя это может быть слишком технической для некоторых.

Это приложение части документации солнца для SPARC. Сканировать весь документ, в ссылках есть этот документ:


Тан, Питер пинг так, некоторые программные реализации функций
Sin и cos, технический отчет АНЛ-90/3, математики и компьютерных
Научный Отдел, Аргоннская Национальная Лаборатория, Аргонн, Иллинойс,
Февраля 1990 года.

Я не могу найти его в интернете с помощью Google, так что удачи охоты, что один!

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

Это может показаться много работы, но когда размер массива фиксируется во время компиляции, можно развернуть цикл (см. Артур Мустафин ответ), и просто использовать 4 местные двухкомнатные (каждая проводит различие между двумя терминами, так что те же 8 условия он рассчитывает).

Основываясь на моих высказываниях, вот моя версия. Встроенный, петля-раскатали, и parallizable. Я обычно не беспокоиться о повторном использовании локальных переменных: компиляторы сейчас достаточно умные, чтобы сделать это для меня. Я также предполагаю, что, например, / (10*11) оптимизирован в * some_constant где some_constant равна 1.0 / 110, предполагая, что процессор быстрее на умножений, чем по отделам.

Но я уверен, что версия, используемая реальные библиотеки-это еще умнее.

double sine_taylor(double x)
{
// useful to pre-calculate
double x2 = x*x;
double x4 = x2*x2;

// Calculate the terms
// As long as abs(x) < sqrt(6), which is 2.45, all terms will be positive.
// Values outside this range should be reduced to [-pi/2, pi/2] anyway for accuracy.
// Some care has to be given to the factorials.
// They can be pre-calculated by the compiler,
// but the value for the higher ones will exceed the storage capacity of int.
// so force the compiler to use unsigned long longs (if available) or doubles.
double t1 = x * (1.0 - x2 / (2*3));
double x5 = x * x4;
double t2 = x5 * (1.0 - x2 / (6*7)) / (1.0* 2*3*4*5);
double x9 = x5 * x4;
double t3 = x9 * (1.0 - x2 / (10*11)) / (1.0* 2*3*4*5*6*7*8*9);
double x13 = x9 * x4;
double t4 = x13 * (1.0 - x2 / (14*15)) / (1.0* 2*3*4*5*6*7*8*9*10*11*12*13);
// add some more if your accuracy requires them.
// But remember that x is smaller than 2, and the factorial grows very fast
// so I doubt that 2^17 / 17! will add anything.
// Even t4 might already be too small to matter when compared with t1.

// Sum backwards
double result = t4;
result += t3;
result += t2;
result += t1;

return result;
}

8
ответ дан 25 октября 2011 в 04:10 Источник Поделиться

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


  • Переводить углы из системы 360 градусов до 256 система степень;

  • Хранить значения синуса для всех 256 возможных степеней в таблице подстановки.

Вы потеряете много точности, но вы получаете много скорости.

11
ответ дан 6 октября 2011 в 07:10 Источник Поделиться

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

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

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

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

Попробуйте это такой:

static int factorial(int n)
{
int i;
int result = 1.0;
for (i = 2; i <= n; i++)
result *= i;
return result;
}

static double pow_int(double x, int y)
{
int i;
double result = 1.0;
for (i = 0; i < y; i++)
result *= x;
return result;
}

double sine_taylor(double x, int n)
{
int i;
double result = x;
for (i = 3; i <= n; i += 4)
result += pow_int(x, i) / factorial(i);
for (i = 5; i <= n; i += 4)
result -= pow_int(x, i) / factorial(i);
return result;
}

6
ответ дан 6 октября 2011 в 08:10 Источник Поделиться

Это правильный ответ. Предоставляя число Пи с точностью 1е-8, вы даете мне возможность оптимизировать функцию в тексте:

inline double _sin8(double radians)
{
double x = radians;
double x2 = x*x;
double f = 1;
double s = 0;
int i = 1;

s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;
s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;
s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;
s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;

return s;
}

Я надеюсь, вы сделали домашнее задание так, как я сделал в моем прошлом:

#include <math.h>

double _sin8(double radians);

inline double _sin(double radians, double epsilon, int &count )
{
double x = radians;
double x2 = x*x;
double f = 1;
double s = 0;
int i = 1;

while (x/f >= epsilon)
{
s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;
count++;
}

return s;
}

int main()
{
double e_x = 0.00000001; // 1e-8, because pi differs from real pi at 8 symbol after decimal point

double e_pi = 3.14159265;
double m_pi = M_PI;

int n = 10;

double e_dx = e_pi/(2.0*n);
double m_dx = m_pi/(2.0*n);

for (double e_i=e_pi/2.0, m_i=m_pi; n>=0; e_i-=e_dx, m_i-=m_dx, n--)
{
int c = 0;
double e_f = _sin(e_i, e_x, c);
double m_f = sin(m_i);
printf("e_pi=%.16e e_f=%.16e c2=%d\n", e_pi, e_f, c);
printf("m_pi=%.16e m_f=%.16e\n", m_pi, m_f);
}
// by the experiment, for the given epsion, we can easily take double _sin8(double) function, because:
// e(f(x)) >= e(x), for any given f(x) = a0*x^0 + a1*x^1 + ... that's called "accumuilation of errors"
// in my opinion, you'd better go to scool library, i do not think it is good for you to
// ask community to do your homework. as you see, i got bachelor degree in CS years ago,
// so we know something about it, and you'd better to discover this knowledge yourself
}

Некоторые mathiematical константы, чтобы помнить:

/* Definitions of useful mathematical constants
* M_E - e
* M_LOG2E - log2(e)
* M_LOG10E - log10(e)
* M_LN2 - ln(2)
* M_LN10 - ln(10)
* M_PI - pi
* M_PI_2 - pi/2
* M_PI_4 - pi/4
* M_1_PI - 1/pi
* M_2_PI - 2/pi
* M_2_SQRTPI - 2/sqrt(pi)
* M_SQRT2 - sqrt(2)
* M_SQRT1_2 - 1/sqrt(2)
*/

#define M_E 2.71828182845904523536
#define M_LOG2E 1.44269504088896340736
#define M_LOG10E 0.434294481903251827651
#define M_LN2 0.693147180559945309417
#define M_LN10 2.30258509299404568402
#define M_PI 3.14159265358979323846
#define M_PI_2 1.57079632679489661923
#define M_PI_4 0.785398163397448309616
#define M_1_PI 0.318309886183790671538
#define M_2_PI 0.636619772367581343076
#define M_2_SQRTPI 1.12837916709551257390
#define M_SQRT2 1.41421356237309504880
#define M_SQRT1_2 0.707106781186547524401

4
ответ дан 7 октября 2011 в 07:10 Источник Поделиться

факторный умоляет, чтобы быть массивом (1,2,6,24,120...)
тяп можно оптимизировать немного... 3^101= 3^64*3^32*3^4*3^1= 3^1*3^32*...*
с помощью быстрого возведения в степень.

Проблема в том, что у вас нет инт как мощность, но по-прежнему можно на "Инт специализация" таким образом. О грехе можно вручную развернуть на петли и убедиться, что у вас есть какой-то параллельной обработки (не несколько потоков, но, что максимальное количество инструкция уровень параллельности можно). Это очень спекулятивный, но допустим, что компилятор может быть достаточно умен, чтобы оптимизировать

 x=(a+b)+(c  +d)+(e+f);`   

таким образом, что это быстрее, чем

x=a+b+c+d+e+f;  

так что вы хотите, что явно для Тейлор сумму.

Но я предполагаю, что компилятор будет раскатать и сделать автоматическое преобразование. Также это может быть так, что с плавающей точкой умножение выполняется быстрее, чем деление, так что вы, возможно, захотите сделать что-то как 1/1,1/2,1/6... блок для расчета Тейлор.

Если вы хотите узнать больше об оптимизации:
http://www.agner.org/optimize/

0
ответ дан 9 ноября 2011 в 05:11 Источник Поделиться