Так-с, давайте по порядку:
def X(m): return m and (2*m-1)*X(m-1) or 1
def Y(m): return m and (4*m*m-8*m+3)*Y(m-2)+7*X(m-2) or 1
b=10**12
print (Y(b*4)*b)/X(b*4)
Даже если забыть про лимит рекурсии в python`е, то про скорость выполнения забыть не получится: она будет явно больше нескольких сотен (если не тысяч) лет.
Немного оптимизируем наши вычисления.
Для начала заметим, что для вычисления, например, X(m)
требуется знание X(m-1)
, поэтому можно будет ускорить вычисления убрав рекурсию (пример тов. @insolor):
def X1(m):
p = 1
for i in range(1,m+1):
p *= 2*i-1
return p
Далее, поскольку Y(m)
выражается через Y(m - 2)
, то выразив аналогичным образом и X(m)
(т.е. через X(m - 2)
) мы сможем проводить вычисления в одном цикле:
X(m) = (2 * m - 1) * X(m - 1)
X(m - 1) = (2 * (m - 1) - 1) * X(m - 2) = (2 * m - 3) * X(m - 2) =>
=> X(m) = (2 * m - 3) * (2 * m - 1) * X(m - 2)
Перед слиянием функций немного изменим нашу функцию Y(m)
:
Y(m) = (4 * m * m - 8 * m + 3) * Y(m - 2) + 7 * X(m - 2) для m > 0, иначе 1
=> Y(m) = (4 * m * m - 8 * m + 4 - 1)...
=> Y(m) = (4 * (m * m - 2 * m + 1) - 1)...
=> Y(m) = (4 * (m - 1)**2 - 1)...
=> Y(m) = ((2 * (m - 1) - 1) * (2 * (m - 1) + 1))...
=> Y(m) = (2 * m - 3) * (2 * m - 1)...
Таким образом получим нашу ускоренную функцию Y(m)
:
def Y(m):
p = 1
o = 1
for i in range(2, m + 1, 2):
current_mul = (2 * i - 3) * (2 * i - 1)
p *= current_mul
p += 7 * o;
#X(i)
o *= current_mul
return p
Вот только для m = 4 * 10**12
результат будет получен через пару десятков лет, если не позже...
Поэтому надо "хорошо напрячь математику" и оптимизировать не только вычисления, но и формулы.
"Напрягаем математику"
Мы получили следующие выражения:
X(m) = (2 * m - 3) * (2 * m - 1) * X(m - 2)
Y(m) = (2 * m - 3) * (2 * m - 1) * Y(m - 2) + 7 * X(m - 2)
Нам необходимо найти Y(m) / X(m)
, т.е.:
(2*m - 3) * (2*m - 1) * Y(m - 2) + 7 * X(m - 2) Y(m - 2) 7
--------------------------------------------------- = -------- + --------------------
(2*m - 3) * (2*m - 1) * X(m - 2) X(m - 2) (2*m - 3) * (2*m - 1)
Продолжая таким образом вычисления получим, что необходимо найти частичную сумму ряда, n-ый член которого определяется формулой:
an = 7 / ((2*n - 3) * (2*n - 1))
Надо только не забыть, что Y(0)
и X(0)
у нас выбиваются из этого ряда. Таким образом нам надо найти сумму:
Y(m) / X(m) = 1 + (a2 + a4 + a6 + ... + am)
Честно, как работать с рядами уже забыл (давно это было), поэтому воспользовался гуглом и нашел эту статью. Там есть подобный пример (Пример 5).
Сделав аналогичные вычисления для нашего ряда получим, что n-ый член ряда равен 7/2 / (2*n - 3) - 7/2 / (2 * n - 1)
. Поскольку 7/2 общий множитель каждого члена, то его можно вынести за скобку:
Y(m) / X(m) = 1 + 3.5 * (a2 + a4 + a6 + ...)
an = 1 / (2*n - 3) - 1 / (2*n - 1)
таким образом исходная сумма ряда (не всего выражения, а ряда!!!) равна:
1/1 - 1/3 + 1/5 - 1/7 + ... + 1/(2*n - 3) - 1/(2 * n - 1)
Этот ряд показался знакомым, но поскольку ничего уже не помню, то решил опять воспользоваться поиском и нашел, что он равен pi / 4
(доказательство)
Таким образом получаем, что
Y(m) / X(m) = 1 + 3.5 * PI / 4 = 3.748893571891069
Вернемся к задаче
Нам надо найти m * Y(4*m) / X(4*m)
для m = 10**12
.
Т.к. мы нашли точное решение (конечно, без части m * ...
), а в задаче приведено примерное (т.е. решение с определенной точностью), то наше решение может немного не совпадать с ожидаемым. Например, найдите ответ для m = 10
в исходном решении и сравните его с нашим - ответы будут разными.
Поскольку мы знаем как найти n-ый член (an = 7 / ((2*n - 3) * (2*n - 1))
), то мы можем оценить его влияние на конечный результат:
an = 7. / (2 * n - 3) / (2 * n - 1) = 1.093750000000547e-25 для n = 4 * 10**12
Если нам необходимо вывести только целую часть исходного выражения, то наш ответ должен будет подойти. Если же надо будет вывести еще и дробную часть, то у меня пока мыслей нет что и как надо делать, т.к. an гораздо меньше найденной нами точности 10**-15