Как уйти от рекурсии?

Рейтинг: 2Ответов: 2Опубликовано: 18.03.2015

Доброго дня!

Есть программка:

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)

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

1: 3
2: 7
3: 11
4: 14
5: 18
...
10: 37
...
100: 374
...
200: 749

Т.е., грубо говоря, ответ около b*3,745, что вполне реально вывести, вопрос только, каким образом этот алго оптимизировать, чтобы он мог отработать за приемлемое время.

Ответы

▲ 7Принят

Так-с, давайте по порядку:

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

▲ 3

Первая функция разворачивается так:

def X1(m):
    p = 1
    for i in range(1,m+1):
        p *= 2*i-1
    return p

Вторая функция разворачивается похожим образом.

UPD. Со второй функцией все оказалось не так уж просто.

Неоптимизированный вариант:

def Y1(m):
    p = 1

    if m>0 and m%2==0:
        k=2
    else:
        k=1

    for i in range(k, m+1, 2):
        p = (4*i*i-8*i+3)*p+7*X1(i-2)

    return p

Оптимизированный вариант:

def Y2(m):
    p = 1

    if m>0 and m%2==0:
        k=2
    else:
        k=1

    x = 1
    for i in range(k, m+1, 2):
        p = (4*i*i-8*i+3)*p+7*x

        if i-1>0:
            x *= 2*(i-1)-1

        x *= 2*i-1
    return p

Замеры производительности:

In [8]:

%timeit Y(800)

1 loops, best of 3: 221 ms per loop

In [9]:

%timeit Y1(800)

10 loops, best of 3: 141 ms per loop

In [10]:

%timeit Y2(800)

100 loops, best of 3: 3.28 ms per loop

UPD2. Если считать, что m всегда четное, то можно избавиться от лишних условий:

def Y3(m):
    assert(m%2==0)
    p = 1
    x = 1
    for i in range(2, m+1, 2):
        p = (4*i*i-8*i+3)*p+7*x
        x *= (2*(i-1)-1)*(2*i-1)
    return p

Производительность выросла незначительно:

In [31]:

%timeit Y2(10000)

1 loops, best of 3: 517 ms per loop

In [30]:

%timeit Y3(10000)

1 loops, best of 3: 486 ms per loop

Теперь вернемся к выражению (Y(b*4)*b)/X(b*4) и предположим, что при больших b выражение приближается к K*b. Это действительно похоже на правду:

for b in range(100, 5000, 100):
    print(b, (Y3(b*4)*b)/X1(b*4) / b)

100 3.746706075309011
200 3.747799822318314
300 3.7481644053509937
400 3.7483466969444748
500 3.748456071918413
600 3.7485289885735598
700 3.7485810719010337
800 3.748620134397745
900 3.7486505163402017
1000 3.748674821894487
...
4700 3.7488470293379104
4800 3.748847998974433
4900 3.7488489290339553

b = 12500
print(b, (Y3(b*4)*b)/X1(b*4) / b)

12500 3.748876071891071

Коэффициент приближается к чему-то вроде 3.74888. Осталось доказать эту зависимость аналитически ;)