в точках x=0,0.1, ... , 0.9,1. Абсолютная погрешность вычисленных сумм должна быть
Полное время счета не должно превышать 10 секунд.
Комментарии в этот раз не скрываю, обсуждение открытое. Соответственно, могут быть спойлеры. Vorsicht!
PS. Не мучайте Гугл, суть задачи суметь ряд просуммировать, а не числа найти:
где
UPD. Уменьшил погрешность, дабы "лобовое" решение манило меньше.
Решение с пояснениями.
const = pi*pi/6;
for(x = 0; x < 1.1; x+=0.1)
{
sum1 = 0;
sum2 = 0;
k = 1;
do
{
sum1 += 1/(k*(k+x));
sum2 += 1/(k*k);
k++;
} while(const-sum2 > 0.0000000001);
printf(" %f ",sum1);
}
Сейчас влом проверять, но в большинстве случаев этого достаточно
Так что если цель именно "ряд просуммировать", то задача не выглядит разрешимой, если не использовать необычные компьютеры (ну или хотя бы видеокарты): многопоточность мне только эти первые сто миллионов и позволит сложить в отведённое время (причём для более простой задачи с заведомо неучтёнными накопленными ошибками).
Использование опыта человечества в вычислении специально придуманной и изученной "дигаммы-функции и постоянной Эйлера-Маскерони" выглядит более разумным, а также экономия на совместном вычислении значений в разных точках, но это же получается уже не просто "суммированием ряда".
Чего я не понял про ограничения?
Использование дигамма-функции будет по факту означать, что вы переложили суммирование ряда на автора математической библиотеки, это неинтересно.
Задача состоит в том, чтобы поиграться с рядом и привести его к виду, который будет допускать суммирование.
Стартовые дефиниции
Погрешностюь IEEE 754 (чисел двойной точности) пренебрежём.
1. Первое что приходит в голову -- разложить бесконечный ряд в конечный с остаточным членом и суммировать, пока оценка ост. члена не будет меньше погрешности. Но я уже забыл всё это (по условиям гугл не мучаем) поэтому зафиксируем, что решение есть и пойдём рассуждать.
2. Второе, что приходит в голову -- заключить ряд в две функции (оценка сверху-снизу) и сводить обе, пока точность нас не удовлетворит. Это по сути тот же остаточный член, только с переизобретением велосипеда.
Второй частью и займёмся, когда появятся идеи.
ПС
Спасибо за разминку мозгов.
Сравнение с остатком ряда эйлера.
Обозначения:
остаток нашего ряда через G(k=k0, x),
Остаток ряда эйлера через E(k==k1).
1.
Тогда:
E(k) < G(k, x) < E(k+x) <= E(k + [x+1])
2.
Решение
просуммируем k0 первых членов ряда, пока
E(k0+2) не станет меньше погрешности => G(k+1, x) тоже меньше погрешности.
3. Улучшение оценка остаточного члена как среднего от E(k) и E(k+[x+1])
Мы можем считать остаточные члены эйлеровского ряда:
E(k) и E(k + [x+1])
и считать что наш остаток G(k, x) -- лежит между ними. Как-то пропорционально тому, к какому краю он ближе.
Я бы предложил считать пропорционально квадрату, но могу ошибаться.
4. Код
double calc_summ(double x, double delta)
{
double eiler_real = pi * pi /6;
int extra_iterations = trunc(x + 1);
double eiler_find = 0;
double summ_find = 0;
for (i = 1; i <= extra_iterations; i++)
eiler += 1/(i * i);
while((eiler - eiler_find) > delta) {
summ_find += 1 / (i * (i + x));
eiler_find += 1 / ((i + extra_iterations) * (i + extra_iterations))
}
return summ_find;
}
Re: Стартовые дефиниции
А вот игнорировать погрешность как раз не надо, в этом суть задачи.
В заданной форме этот ряд тупо не суммируется с требуемой точностью, вообще (даже на 80-разрядном флоате). Надо извернуться.
Я же советовал уже попробовать реально посуммировать и посмотреть, что получается.
Re: Стартовые дефиниции
Ну есть решение в лоб -- с увеличивающими точность (квадрат точности ЕМНИП) преобразования при использовании двух float / double вместо одного.
Но речь явно не про них, а про "подумать".
Сейчас посчитаю "в лоб" и пойду думать хД
Re: Стартовые дефиниции
Даже если предположить, что нет погрешности совсем, лобовой ряд сходится так медленно, что ни о каких 10 секундах речи не идет.
Re: Стартовые дефиниции
Рассмотрим ряд эйлера и ряд-эйлера+1.
1/(n^2) < f(n) < 1/((n+1)^2)
Что мы можем сказать?
1. сумма нашей ф-ии лежит между суммами рядов.
2. для n_i -- i-й итерации:
остаточный_ряд_эйлера - остаточный_ряд_эйлера_1 = 1/(n+1)^2.
3. Значит, если:
-а) посчитать 1/sqrt(delta) членов последовательности
-б) взять за ост. член -- ост. член эйлера -- всё сойдётся.
x = 0.100000; sum = 1.534607
x = 0.200000; sum = 1.440879
x = 0.300000; sum = 1.360083
x = 0.400000; sum = 1.289578
x = 0.500000; sum = 1.227411
x = 0.600000; sum = 1.172105
x = 0.700000; sum = 1.122519
x = 0.800000; sum = 1.077759
x = 0.900000; sum = 1.037111
x = 1.000000; sum = 1.000000
Re: Стартовые дефиниции
Теперь по поводу double:
1. Достаточно просто считать с конца, чтобы не было ошибки, т.е.
sum = 1/(n*(n+x)) + 1/((n-1)*(n-1+x) + ...
2. Но по хорошему это надо проверить.
Проверяется это так: вместо суммирования по одному члену ряда будем суммировать группу членов ряда, а в конце -- суммировать между собой частичные суммы.
По-хорошему надо контролировать чтобы при "суммировании частичных сумм" не вылезла ошибка.
Но по-плохому если: суммирование-с-конца-в-лоб == суммирование-частичных-сумм -- Значит всё сошлось.
Код проверки (через частичные суммы). Грязноват, но всё-же:
#define ACCUR 100
double calc_summ(double x, double delta, double *eiler_ret)
{
int num = (int)(1/sqrt(delta)) + 1;
double eiler_real = M_PI * M_PI /6;
double sum = 0;
double eiler_sum = 0;
int i;
int n = 0;
for (i = 1; i <= num; i++) {
double i_e = i;
double curr = 1/(i_e * (i_e + x));
sum += curr;
double eiler = 1 / (i_e * i_e);
eiler_sum += eiler;
if (curr * ACCUR < sum) {
C[n] = sum;
E[n] = eiler_sum;
sum = 0;
eiler_sum = 0;
n++;
}
}
int j;
for (j = n-1; j >= 0; j--) {
sum += C[j];
eiler_sum += E[j];
}
return sum + (eiler_real - eiler_sum);
}
Edited at 2018-03-27 05:16 pm (UTC)
Re: Стартовые дефиниции
Ну тут решение понятно.
Запускаем рекурсию, в которой храним текущую сумму.
Как только очередной член меньше суммы во сколько-то раз -- перезапускаем рекурсивную функцию.
ПС
Решение заведомо некрасивое -- и отношение точностей наколенное
Первые два слагаемых приведём к общему знаменателю, получая:
Зная, что сумма по 1/k2 даёт π2/6 и вынося x за знак суммы, получим:
эта сумма сходится существенно лучше, поскольку "хвост" пропорционален 1/k2, т.е возьмём в 10 раз больше слагаемых - в 100 раз выше точность.
"Прикинуть" точность можно по x = 1, поскольку нетрудно показать, что результатом должна стать единица:
При суммировании сократятся все слагаемые, кроме самого первого (=1) и самого последнего, стремящегося к нулю.
Я взял 6 млн слагаемых и должен был получить точность на несколько порядков выше, чем требовалось. Использовал 80-битные числа с плавающей точкой, суммировал ряд "справа налево", от самых маленьких к самым большим, вроде так учили :)
Такая простыня получается:
x=0; f(x)=1,64493406684823
x=0,1; f(x)=1,53460724490456
x=0,2; f(x)=1,44087884154673
x=0,3; f(x)=1,36008258678245
x=0,4; f(x)=1,28957780079105
x=0,5; f(x)=1,22741127776023
x=0,6; f(x)=1,17210519612502
x=0,7; f(x)=1,12251934253576
x=0,8; f(x)=1,07775887274425
x=0,9; f(x)=1,03711091785067
x=1; f(x)=1,00000000000001
посчитано за 8954 мс
Использовался древний компьютер с Intel Atom D510, задействовано только одно ядро.
PS. Попробовал и слева направо посчитать, разница в 14-м знаке после запятой, причем "справа налево" явно точнее. f(1), посчитанный слева направо вышел 1,00000000000004 - в 4 раза больше ошибка :)
Edited at 2018-03-27 05:35 pm (UTC)