Длинное выражение приводит к сбою SymPy

Я использую 64-битный Python 3.3.1, pylab и 32 ГБ оперативной памяти. Эта функция:

def sqrt2Expansion(limit):
    x    = Symbol('x')
    term = 1+1/x
    for _ in range(limit):
        term = term.subs({x: (2+1/x)})
    return term.subs({x: 2})

Производит выражения такого вида: 1 + 1/(2 + 1/(2 + 1/(2 + 1/(2 + 1/(2 + 1/(...)))))). При вызове как: sqrt2Expansion(100) возвращает действительный результат, но sqrt2Expansion(200) выдает RuntimeError с большим количеством страниц трассировки и зависает интерпретатор pylab/IPython с большим количеством неиспользуемой системной памяти. Любые идеи, как реализовать это более эффективно? Я хотел бы позвонить sqrt2Expansion(1000) и все равно получить результат.


python-3.x sympy expression
person Paul Jurczak    schedule 06.05.2013    source источник
comment
Замена term = term.subs({x: (2+1/x)}) на term = term.subs({x: (2+1/x)}).factor() должна помочь. Не могли бы вы объяснить цель того, что вы делаете?   -  person Krastanov    schedule 06.05.2013
comment
@Krastanov Спасибо, это действительно работает для sqrt2Expansion(1000). Моя цель состояла в том, чтобы изучить SymPy, решая в этом случае проблему Эйлера 57. Я знаю более элегантные решения, но я хотел посмотреть, сработает ли в этом случае грубая сила. Если вы хотите опубликовать свой комментарий в качестве ответа, я проголосую за него и выберу его в качестве принятого ответа.   -  person Paul Jurczak    schedule 06.05.2013
comment
Я попытался ответить более подробно ниже. Я надеюсь, что вы найдете это полезным.   -  person Krastanov    schedule 06.05.2013

Ответы (2)


arrow_upward
2
arrow_downward

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

Выражения Sympy — это деревья. Каждая операция представляет собой узел, который имеет в качестве ветвей свои операнды. Например, x+y выглядит как Add(x, y), а x*(y+z) как Mul(x, Add(y, z)).

Обычно эти выражения автоматически выравниваются, например, Add(x, Add(y, z)) становится Add(x, y, z), но для более сложных случаев можно получить очень глубокие деревья.

И глубокие деревья могут создавать проблемы, особенно когда либо интерпретатор, либо сама библиотека ограничивают глубину разрешенной рекурсии (в качестве защиты от бесконечной рекурсии и взрывного использования памяти). Скорее всего, это причина вашего RuntimeError: каждый subs делает дерево глубже, и по мере того, как дерево становится глубже, рекурсивный subs должен вызывать себя больше раз, пока не доберется до самого глубокого узла.

Вы можете упростить дерево до вида polynomial/polynomial с постоянной глубиной, используя метод factor. Просто измените term = term.subs({x: (2+1/x)}) на term = term.subs({x: (2+1/x)}).factor().

person Krastanov    schedule 06.05.2013
comment
Я подозревал, что проблема заключается в глубине рекурсии, но было бы неплохо, если бы SymPy выдавал соответствующее сообщение об ошибке, а не просто зависал. В очередной раз благодарим за помощь. - person Paul Jurczak; 07.05.2013
comment
Я думаю (могу ошибаться), что эта ошибка контролируется интерпретатором и что ее вообще сложно отловить в коде (это не SymPy выполняет raise RuntimeError). Таким образом, SymPy не имеет большого контроля над сообщением. С другой стороны, если какой-нибудь добровольец немного рефакторит subs, чтобы не использовать рекурсию, ошибка и ограничение исчезнут. - person Krastanov; 07.05.2013
comment
cancel будет более эффективным, чем фактор. Если вы не беспокоитесь об устранении общих факторов, expr.as_numer_denom является наиболее эффективным, но в некоторых случаях это может привести к очень большому росту выражения. - person asmeurer; 13.05.2013
comment
@asmeurer Я рассчитал время версии с cancel, и она работает в 66% времени, которое занимает версия factor. Спасибо. - person Paul Jurczak; 13.05.2013

arrow_upward
1
arrow_downward

Есть ли причина, по которой вам нужно сделать это символически? Просто начните с 2 и работайте оттуда. Вы никогда не столкнетесь с ошибками рекурсии, потому что дробь будет сглаживаться на каждом этапе.

In [10]: def sqrt2(limit):
   ....:     expr = Rational(1, 2)
   ....:     for i in range(limit):
   ....:         expr = 1/(2 + expr)
   ....:     return 1 + expr
   ....:

In [11]: sqrt2(100)
Out[11]:
552191743651117350907374866615429308899
───────────────────────────────────────
390458526450928779826062879981346977190

In [12]: sqrt2(100).evalf()
Out[12]: 1.41421356237310

In [13]: sqrt(2).evalf()
Out[13]: 1.41421356237310

In [15]: print sqrt2(1000)
173862817361510048113392732063287518809190824104684245763570072944177841306522186007881248757647526155598965224342185265607829530599877063992267115274300302346065892232737657351612082318884085720085755135975481584205200521472790368849847501114423133808690827279667023048950325351004049478273731369644053281603356987998498457434883570613383878936628838144874794543267245536801570068899/122939577152521961584762100253767379068957010866562498780385985503882964809611193975682098617378531179669585936977443997763999765977165585582873799618452910919591841027248057559735534272951945685362378851460989224784933532517808336113862600995844634542449976278852113745996406252046638163909206307472156724372191132577490597501908825040117098606797865940229949194369495682751575387690
person asmeurer    schedule 12.05.2013
comment
Поскольку вы решаете эту задачу проекта Эйлера, вы, вероятно, захотите изменить алгоритм для проверки результата на каждом этапе, а не пересчитывать все это 1000 раз (хотя, честно говоря, даже если вы сделаете это таким образом, это завершится через несколько секунд). - person asmeurer; 13.05.2013
comment
Какой импорт мне нужен для Rational? - person Paul Jurczak; 13.05.2013
comment
Я заменил Rational на Fraction, чтобы он работал в Python 3.3. У вас самое элегантное декларативное решение, которое я когда-либо видел. Я начал в этом направлении, но застрял — я новичок в Python. Я играл с SymPy в основном для будущих приложений — я знаю, что это излишество для этой проблемы. Спасибо за вашу помощь. - person Paul Jurczak; 13.05.2013
comment
Rational исходит от SymPy. Дробь так же хорошо работает и в этой задаче (хотя с ней вы не получите красивую 2D-печать :). - person asmeurer; 13.05.2013
comment
И SymPy отлично работает в python 3 (вам нужно скачать версию python 3). - person asmeurer; 13.05.2013