Решение менее типичной задачи

[ahiin]

Условие задачи здесь.

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

Так же не могу не отметить упорство son_0f_morning, который подошел к задаче с другой стороны и, таки, прорвался в итоге к решению.

Под катом само решение (теория) и несколько вариантов реализации.


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


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


Зная, что сумма по 1/k2 даёт π2/6 и вынося x за знак суммы, получим:



эта сумма сходится существенно лучше, поскольку "хвост" пропорционален 1/k2, т.е возьмём в 10 раз больше слагаемых - в 100 раз выше точность.

"Прикинуть" точность можно по x = 1, поскольку нетрудно показать, что результатом должна стать единица:


При суммировании сократятся все слагаемые, кроме самого первого (=1) и самого последнего, стремящегося к нулю. "

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

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

Подобный фокус: выделить некую "главную часть", посчитать аналитически, а оклонения от нее уже считать численно - исключительно полезен и встречается повсеместно, особенно в приложениях. Например, таким способом прекрасно берутся численно некоторые интегралы и т.п.

Очевидно, что с рядом
можно проделать тот же фокус (и так много раз).

Пока мне не надоело, я успел получить следующее выражение:


Здесь
- дзета-функция Римана, которая определяется следующим образом:


Нужные нам значения хорошо известны, есть в том числе в Википедии. Лично я воспользовался Вольфрамом, как заведомо более надежным источником.

Моя версия решения под спойлером.
[Код на C++]

#include
#include
#include

constexpr double zeta2 = 1.644934066848226436472415166646025189218949901206798437735;
constexpr double zeta3 = 1.202056903159594285399738161511449990764986292340498881792;
constexpr double zeta4 = 1.082323233711138191516003696541167902774750951918726907682;
constexpr double zeta5 = 1.036927755143369926331365486457034168057080919501912811974;
constexpr double zeta6 = 1.017343061984449139714517929790920527901817490032853561842;
constexpr double zeta7 = 1.008349277381922826839797549849796759599863560565238706417;

constexpr double epsilon = 1e-14;
constexpr int half_range = 100;

double Sum(double x)
{
double sum1 = 0.0;
double sum2 = 0.0;

for (int i = half_range; i > 0; --i)
{
double index = static_cast(i);
double inv1 = index*index*index*index*index*index*index*(index + x);
sum1 += 1.0 / inv1;

index += half_range;
double inv2 = index*index*index*index*index*index*index*(index + x);
sum2 += 1.0 / inv2;
}

return (sum2 < epsilon)*(sum1 + sum2); //отладочный контроль
}

double f(double x)
{
return zeta2 - x*(zeta3 - x*(zeta4 - x*(zeta5 - x*(zeta6 - x*(zeta7 - x*Sum(x))))));
}

int main()
{
clock_t start = clock();

double fn[10000]; //дабы компилятор не выкинул "лишние" проходы цикла
int cntr = 0;
for (int k = 0; k < 1000; ++k) //чтобы хоть как-то замерить время работы, вычисления делаеются 1000 раз
for (double x = 0.1; x < 1.05; x += 0.1)
{
fn[cntr++] = f(x);
}

std::cout << "Run time: " << clock() - start << " ms" << "\n\n";

std::cout << std::setprecision(17) << zeta2 << '\n';
for (cntr = 0; cntr < 10; ++cntr)
{
std::cout << std::setprecision(17) << fn[cntr] << '\n';
}
}



На моей машине (Intel® Core™ i7-4720HQ) время счета составляет 8 микросекунд. До лимита в 10 секунд достаточно далеко, как видим.

Можно пойти и другим путем. son_0f_morning предложил несколько вариантов, о которых он сам расскажет при желании, здесь рассмотрим его решение, которое я зачел первым.

Оно базируется на следующей идее: при достаточно большом k, члены рядов
и
становятся практически неразличимы.

Таким образом, если взять n достаточно большим, получим:


Идея вполне годная, авторская реализация под спойлером.
[Код на C++]

#include
#include
#include

double calc_summ(double x, double delta)
{
int num = 2 * (int)(1 / sqrt(delta)) + 1;
double eiler_real = M_PI * M_PI / 6;

double sum = 0;
double eiler_sum = 0;

for (int i = num; i >= 1; 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;
}

return sum + (eiler_real - eiler_sum);
}

double delta = 0.00000000000001;
int main()
{
double x;

int i;
for (i = 0; i < 11; i++) {
double x = (double)(i) / 10;
double sum = calc_summ(x, delta);
printf("%.14lf(x = %f)\n", sum, x);
}
}



Для языкового разнообразия, привожу так же вариант уважаемого alexis_m.
[Код на Python]

from __future__ import print_function

import math
import timeit

# https://ahiin.livejournal.com/271654.html

N = 1000000

truth = [
1.6449340668482264,
1.5346072449045607,
1.4408788415467229,
1.360082586782444,
1.2895778007910417,
1.2274112777602189,
1.1721051961250153,
1.1225193425357527,
1.0777588727442429,
1.0371109178506583,
0.99999999999999989
]

def f(x):
a = math.pi ** 2 / 6.0
b = 0
for k in xrange(N, 0, -1):
b += 1.0 / (k + x) / k / k

return a - x * b

def calc(i):
p = truth[i]
x = i / 10.0
q = f(x)
print('{:.1f} {:.16f} {:g}'.format(x, q, abs(q - p)))

def main():
t1 = timeit.default_timer()
print(' x ', '{:18}'.format('f(x)'), 'err')
for n in xrange(11):
calc(n)
t2 = timeit.default_timer()

print(t2 - t1, 'sec')

if __name__ == '__main__':
main()