[Python] sympy.solve never returns for certain coefficients

The form of the expression (corresponding the the variable called 'expression' below) whose critical points I'm trying to solve for is identical in all cases since the only thing that is changing between success and failure is the values of the exponents on each term.

Here is a minimum reproducible example:

from sympy import symbols, Mul, diff, solve
from decimal import Decimal


class Histogram:
def __init__(self, bins, pdf):
self.bins = bins
self.pdf = pdf


histograms = [
#BCHUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.5059288537549406772342308613588102161884307861328125'),
Decimal('0'),
Decimal('0.494071146245059267254617907383362762629985809326171875')]
),
#XRPUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.497461928934010144676136633279384113848209381103515625'),
Decimal('0'),
Decimal('0.50253807106598979981271213546278886497020721435546875')]
),
#EOSUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.494252873563218397823249006250989623367786407470703125'),
Decimal('0'),
Decimal('0.5057471264367816576879022250068373978137969970703125')]
),
#ETHUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.483333333333333337034076748750521801412105560302734375'),
Decimal('0'),
Decimal('0.51666666666666671847707448250730521976947784423828125')]
),
#LINKUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.486166007905138364275643425571615807712078094482421875'),
Decimal('0'),
Decimal('0.51383399209486169123550780568621121346950531005859375')]
),
#XLMUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.507246376811594235078928250004537403583526611328125'),
Decimal('0'),
Decimal('0.492753623188405820432222981253289617598056793212890625')]
),
#TRXUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.489361702127659559113936893481877632439136505126953125'),
Decimal('0'),
Decimal('0.51063829787234038537491187526029534637928009033203125')]
),
#ATOMUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.49829351535836174935667486352031119167804718017578125'),
Decimal('0'),
Decimal('0.50170648464163825064332513647968880832195281982421875')]
),
#CHZUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.49056603773584905869853400872671045362949371337890625'),
Decimal('0'),
Decimal('0.50943396226415094130146599127328954637050628662109375')]
),
#DASHUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.5287356321839080663949062000028789043426513671875'),
Decimal('0'),
Decimal('0.4712643678160919336050937999971210956573486328125')]
),
#ZECUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.52941176470588235947190014485386200249195098876953125'),
Decimal('0'),
Decimal('0.47058823529411764052809985514613799750804901123046875')]
),
#NEOUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.479020979020979009721514785269391722977161407470703125'),
Decimal('0'),
Decimal('0.52097902097902093476733398347278125584125518798828125')]
),
#QTUMUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.5112781954887217761296369644696824252605438232421875'),
Decimal('0'),
Decimal('0.488721804511278168359211804272490553557872772216796875')]
),
#IOSTUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.52280701754385960899895735565223731100559234619140625'),
Decimal('0'),
Decimal('0.477192982456140335489891413089935667812824249267578125')]
),
#COTIUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.5'),
Decimal('0'),
Decimal('0.5')]
),
#IOTAUSDT
Histogram(
bins=[-1, 0, 1],
pdf=[Decimal('0.454861111111111104943205418749130330979824066162109375'),
Decimal('0'),
Decimal('0.545138888888888839545643349993042647838592529296875')]
)
]


def computation(histogram):

A = Decimal('1')
e = Decimal('1')
M_factor = Decimal('1.0828567056280801')
x = symbols('x')
product_terms = []
profit = lambda cp: (e * M_factor ** cp - e) * x
commission = lambda cp: (e * M_factor ** cp + e) * Decimal('0.0005') * x

for candles_profit, freq in zip(histogram.bins, histogram.pdf):
profit_term = profit(candles_profit)
commission_term = commission(candles_profit)
term = ((A + profit_term - commission_term) / A) ** freq
product_terms.append(term)

expression = Mul(*product_terms)
derivative = diff(expression, x)
solutions = solve(derivative, x)
print("Solutions:", solutions)


for histogram in histograms:
computation(histogram)


If you create a python file with the example above and run it, you will see only 15 results in the printed output. The 16th computation for IOTAUSDT never returns.

Each example, top to bottom, turns into one of the expressions below when we inspect the value of the expression variable in the debugger:

(1.0 - 0.0774785191289289x)^0.505928853754941 \* (0.0818152772752661x + 1.0)^0.494071146245059

(1.0 - 0.0774785191289289x)^0.49746192893401 \* (0.0818152772752661x + 1.0)^0.50253807106599

(1.0 - 0.0774785191289289x)^0.494252873563218 \* (0.0818152772752661x + 1.0)^0.505747126436782

(1.0 - 0.0774785191289289x)^0.483333333333333 \* (0.0818152772752661x + 1.0)^0.516666666666667

(1.0 - 0.0774785191289289x)^0.486166007905138 \* (0.0818152772752661x + 1.0)^0.513833992094862

(1.0 - 0.0774785191289289x)^0.507246376811594 \* (0.0818152772752661x + 1.0)^0.492753623188406

(1.0 - 0.0774785191289289x)^0.48936170212766 \* (0.0818152772752661x + 1.0)^0.51063829787234

(1.0 - 0.0774785191289289x)^0.498293515358362 \* (0.0818152772752661x + 1.0)^0.501706484641638

(1.0 - 0.0774785191289289x)^0.490566037735849 \* (0.0818152772752661x + 1.0)^0.509433962264151

(1.0 - 0.0774785191289289x)^0.528735632183908 \* (0.0818152772752661x + 1.0)^0.471264367816092

(1.0 - 0.0774785191289289x)^0.529411764705882 \* (0.0818152772752661x + 1.0)^0.470588235294118

(1.0 - 0.0774785191289289x)^0.479020979020979 \* (0.0818152772752661x + 1.0)^0.520979020979021

(1.0 - 0.0774785191289289x)^0.511278195488722 \* (0.0818152772752661x + 1.0)^0.488721804511278

(1.0 - 0.0774785191289289x)^0.52280701754386 \* (0.0818152772752661x + 1.0)^0.47719298245614

(1.0 - 0.0774785191289289x)^0.5 \* (0.0818152772752661x + 1.0)^0.5


But for some reason, the final expression for IOTAUSDT causes sympy.solve to never return when we attempt to find its critical points:

(1.0 - 0.0774785191289289x)^0.454861111111111 \* (0.0818152772752661x + 1.0)^0.545138888888889


The unique feature I see of the expression that fails is that it has repeating decimals in its exponents, like a rational number like 2/3. To confirm my intuition that this is the source of the error, I revisited the trades that formed the histogram for this symbol: 131 were in the -1 bucket and 157 were in the +1 bucket. 131/288 = 0.454861 with a repeating 1, and 157/288 = 0.545138 with a repeating 8. So it appears that even though the Decimal values keep many places of precision in the histogram we pass in to the computation function, those values are inaccurate since the repetition ceases after a while, but even that doesn't matter because sympy cuts the values off before that point. Since sympy.solve coerces my input python Decimals into sympy Floats, not sympy Rationals, is it reasonable to expect that this is the source of the issue? If yes, how could I modify my code to allow the computation to resolve properly instead of getting stuck forever?

Continue reading...