Гауссовское интегральное и двойное деление

Для удовольствия я пытался оценить гауссовский интеграл от 0 до 1, используя разложение в ряд. По этой причине, я написал функцию факториала, которая хорошо работает до 20 (я проверил), а затем я написал это:Гауссовское интегральное и двойное деление

int main(){ 
    int n; 
    long double result=0; 
    for(n=0; n<=5; n++){ 
     if(n%2==0){ 
      result+=(((long double) 1/(long double)(factorial(n)*(2*n+1)))); 
     } else { 
      result-=(((long double) 1/(long double)(factorial(n)*(2*n+1)))); 
     } 
    }   
    printf("The Gaussian integral from 0 to 1 is %Lf\n", result); 
}

Это дает мне странное отрицательное число, которое, очевидно, даже не близко. Я подозреваю, что проблема связана с актерским составом, но я не знаю, что это такое. Есть предположения? Это не первое, что я пробовал. Я попытался преобразовать что-либо в выражение и поставить явное действие в начале, но это не сработало.

источник

2014-08-28 Gennaro Marco Devincenzis

Что является результатом 'factorial (0)'? – mch

Работает для меня: http://ideone.com/VLgnoS – interjay

Это 1, я обратился к проблеме, когда писал функцию. –

Вы используете компилятор MinGW (порт gcc для Windows), который имеет проблемы с типом long double. Это связано с конфликтами между реализацией GCC long double и библиотекой Microsoft C. См. Также this question.

Согласно this question, определение __USE_MINGW_ANSI_STDIO может решить эту проблему. Если нет, то вместо этого будет работать double.

источник

2014-08-28 15:04:02 interjay

Впечатляющий подвиг вычета. –

В (long double)(factorial(n)*(2*n+1) умножения являются целыми умножениями, и первый может переполняться, если результат factorial уже близок к пределу используемого целочисленного типа.

Напишите ((long double)(factorial(n))*(2*n+1), чтобы первое умножение было умножением с плавающей запятой.

источник

2014-08-28 14:40:35

'n <= 5' здесь, поэтому он не должен переполняться. – interjay

@interjay Да. Я должен предположить, что код, показанный в вопросе, не совсем тот, который переполняется. Однако, принимая факториальную функцию, которая работает для переданных ей значений и возвращает некоторый целочисленный тип (не заданный в вопросе), целочисленное умножение на '2 * n + 1' является очевидным красным флагом в вопросе и немедленно преобразуется с плавающей точкой - это самое простое средство. –

OP утверждает, что «факториальная функция работает до 20!». Я сомневаюсь, что он переполнит даже 64-битное целое число. Поэтому я согласен с Паскалем - это не SSCCE. Итак, +1 ответ (даже если мой симпатичный). – Bathsheba

Вы почти наверняка переполняете свой целочисленный тип. В C это технически неопределенное поведение.

Для 32-разрядного целого числа без знака, 13! будет переполняться. На 64 бит, 21! будет переполняться.

Ваш алгоритм сохранится немного дольше, если вы используете тип с плавающей точкой или расширение, например __uint128 (дает вам, я думаю, до 34!), Если ваш компилятор его поддерживает.

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

источник

2014-08-28 14:40:43 Bathsheba

Он работает с двойным, хотя и с 6 цифрами точности. Как я могу сделать немного больше? –

предложение об улучшении добавлено в ответ – Bathsheba

Гауссовское интегральное и двойное деление

ответ

Смежные вопросы