1

Доброго дня!

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

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, что вполне реально вывести, вопрос только, каким образом этот алго оптимизировать, чтобы он мог отработать за приемлемое время.

13
  • 1
    если не рекурсия, то итерация. делайте бесконечный цикл, обновляйте переменную, пока она не станет выполнять необходимые условия.
    – etki
    18 мар 2015 в 12:52
  • У вас нет условия конца рекурсии, нужно ограничение на m ввести.
    – Deadkenny
    18 мар 2015 в 13:02
  • задача полная, имеет решение, ограничение не нужно, щапистите для маленьких чисел и убедитесь
    – Isaev
    18 мар 2015 в 13:06
  • 1
    @Isaev, а не могли вы бы озвучить изначальную задачу, из которой получается данная вычислительная схема? Мало приятного разбираться в записях вида expr1 and expr2 or expr3 без скобок.
    – dzhioev
    19 мар 2015 в 1:44
  • 3
    @dzhioev, мне тоже кажется, что вот так сформулировать было бы лучше: def X(m): return 1 if m == 0 else (2 * m - 1) * X(m - 1) def Y(m): return 1 if m == 0 else (4 * m * m - 8 * m + 3) * Y(m - 2) + 7 * X(m - 2) 19 мар 2015 в 3:55

2 ответа 2

6

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

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

4
  • И я писал выше про 20 лет - это были примерные результаты тестов на другом компе. Здесь смотрю - время показывает уже больше 100 лет, так что реально время могло быть и больше 200-300 лет )) А значит исходный код работал бы не одну тысячу лет - быстрее было бы подождать появление кубитов и выполнить код там ) Ну либо провести мозговой штурм и получить результат за пару часов ))
    – BOPOH
    22 мар 2015 в 16:23
  • @Isaev, в той ветке у меня комментарии уже закончились. В ответ я и сам мог преобразовать, но там было кусками, поэтому создал полноценный ответ. Ряд получил когда искал "как найти сумму ряда": нашел статью из ответа, привел к виду an = 1 / (2*n - 3) - 1 / (2*n - 1), а там уже просто подставлял значения n = 2, 4, 6..., в результате и получил этот ряд. Вы мой ответ видно не до конца прочитали, там приведены шаги получения конечного результата ))
    – BOPOH
    23 мар 2015 в 8:35
  • @BOPOH, огромное спасибо, красивое решение! Я бы никогда не догадался наверное, что всё сводится к ряду из формулы Лейбница)
    – Isaev
    23 мар 2015 в 11:53
  • "Если нам необходимо вывести только целую часть исходного выражения, то наш ответ должен будет подойти." Да, судя по начальному коду ответ - целое число python почему-то посчитал выражение ´10**12*(1+7/8*math.pi)´ с огромной погрешностью, но калькулятор windows вполне справился, ответ совпал
    – Isaev
    23 мар 2015 в 13:20
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. Осталось доказать эту зависимость аналитически ;)

9
  • 1
    @Isaev, а вам как срочно ответ получить надо? Я немного оптимизировал ваш код: def Y(m): p = 1 j = 1 o = 1 for i in range(2, m + 1, 2): p *= (4 * (i - 1) * (i - 1) - 1); p += 7 * o; #X(i) o *= (4 * j * j - 1) j += 2 return p Для Y(10**12) первая версия работала бы 33 года, эта - всего 20 лет. Здесь намного быстрее будет "хорошо напрячь математику", без этого вроде никак...
    – BOPOH
    18 мар 2015 в 18:07
  • 1
    Ну да, на сообразительность же задача, а не на 20 лет ожидания) Думаю, вполне логично предположить, что деление одного на второе нужно выполнять сразу во второй функции, чтобы не вычислять оба результата отдельно, тогда изначально избежим очень больших чисел
    – Isaev
    18 мар 2015 в 23:36
  • 2
    @Isaev, выразил x(m) и y(m) через x(m-2) и y(m-2), получил следующее: y(m) / x(m) = y(m-2) / x(m-2) + 7 / ((2m-3) * (2m-1)) Продолжая вычисления получим, что необходимо найти сумму: 1 + 7 / (2 * 2 - 3)(2 * 2 - 1) + 7 / (2 * 4 - 3)(2 * 4 - 1) + ... + 7 / (2 * m - 3)(2 * m - 1) Думаю, с этим работать будет проще )
    – BOPOH
    19 мар 2015 в 17:16
  • 1
    @Isaev, так вам моего комментария хватило? Или вам надо более подробно все расписать? Если грубо, мы получаем ряд (как привести к этому читай здесь): 7/2 * (1/1 - 1/3 + 1/5 - 1/7 + ...) Выражение в скобках равно pi / 4 (док-во). Но поскольку pi мы высчитываем сами до определенного члена, то точность будет ниже. Хотя если print (... * b) / ... означает что мы должны первые 12 цифр вывести, то ответ должен будет подойти.
    – BOPOH
    22 мар 2015 в 9:36
  • 1
    @BOPOH, по вашим выкладкам получается что-то около 2.74, но опытным путем получается ближе к 3.75. И, думаю, стоит преобразовать ваши комментарии в ответ.
    – insolor
    22 мар 2015 в 12:55

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge you have read our privacy policy.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.