?

Log in

No account? Create an account

Популярно о науке

Previous Entry Share Next Entry
Менее типичная задача
ahiin
На языке программирования общего назначения (С/С++/Python и т.п.), напишите программу, которая численно просуммирует ряд:



в точках x=0,0.1, ... , 0.9,1.  Абсолютная погрешность вычисленных сумм должна быть
Полное время счета не должно превышать 10 секунд.

Комментарии в этот раз не скрываю, обсуждение открытое. Соответственно, могут быть спойлеры. Vorsicht!


PS. Не мучайте Гугл, суть задачи суметь ряд просуммировать, а не числа найти:



где — дигамма-функция, а — постоянная Эйлера—Маскерони.

UPD. Уменьшил погрешность, дабы "лобовое" решение манило меньше.

Решение с пояснениями.

  • 1
Топорное решение, но. Считать две суммы: для 1/k(k+x) и для 1/k*k. Хвост первой суммы - искомая погрешность, хвост второй суммы больше хвоста первой суммы и через него можно оценивать погрешность. В итоге как-то так:

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);
}

Советую все же поэкспериментировать с реальным суммированием, хехе. Хотя бы и 1/k^2.

Ну, прежде всего, суммировать надо задом наперед (сначала меньшие члены), т.е. определить k максимальное, от которого спускаться к единице.

Сейчас влом проверять, но в большинстве случаев этого достаточно


У меня самым тупым способом приближённая сумма первых ста миллионов членов 1/k^2 занимает 30 секунд (Питон). Это не учитывает накопления погрешностей и того, что сумма "хвоста" довольно велика (порядка 1/n, если хвост начинается с n-ного члена, так что считать надо далеко за миллиардный член).

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

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

Чего я не понял про ограничения?

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

Задача состоит в том, чтобы поиграться с рядом и привести его к виду, который будет допускать суммирование.

Стартовые дефиниции

Без математики (в смысле аналитической математики) задачку всё равно не решить. Нас интересует лишь СТЕПЕНЬ, в которой мы привлечём аналитику.

Погрешностюь 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: Стартовые дефиниции

update:
Ну тут решение понятно.
Запускаем рекурсию, в которой храним текущую сумму.

Как только очередной член меньше суммы во сколько-то раз -- перезапускаем рекурсивную функцию.

ПС
Решение заведомо некрасивое -- и отношение точностей наколенное

Внутри суммы прибавим и отнимем 1/k2:



Первые два слагаемых приведём к общему знаменателю, получая:



Зная, что сумма по 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)

Прикрою пока) Все норм, да.

  • 1