3PCF not working for poles >= 2

[AlanPearl]

Issue

Both SurveyData3PCF and SimulationBox3PCF can be instantiated for example when poles=[0, 1]. However, if I include a pole of 2 or higher in the list of poles, it never instantiates. I've tracked the problem down to the YlmCache class, which seems to enter an infinite loop somewhere in expr.factor() of YlmCache._get_Ylm() whenever m = 2 and above. I'm not super familiar with sympy expressions, so I don't know how to fix this.

Producing this bug

Installation (I'm running Ubuntu 18.04.4 on WSL)

conda create --name nbodykit-env python=3.7 conda activate nbodykit-env conda install -c bccp nbodykit

Anything with m < 2 works fine

python -c "from nbodykit.algorithms.threeptcf import YlmCache; print(YlmCache._get_Ylm(None, 2, 1))"

-0.772548404046379*xpyhat*zhat

Runs forever with m >= 2

python -c "from nbodykit.algorithms.threeptcf import YlmCache; print(YlmCache._get_Ylm(None, 2, 2))" Hit Ctrl-C after 60 seconds

^CTraceback (most recent call last):
File "<string>", line 1, in <module>
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/nbodykit/algorithms/threeptcf.py", line 527, in _get_Ylm
expr = expr.factor().factor(extension=[sp.I]).subs(xhat+sp.I*yhat, xpyhat)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/core/expr.py", line 3587, in factor
return factor(self, *gens, **args)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/polytools.py", line 6361, in factor
return _generic_factor(f, gens, args, method='factor')
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/polytools.py", line 6033, in _generic_factor
return _symbolic_factor(sympify(expr), opt, method)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/polytools.py", line 5973, in _symbolic_factor
coeff, factors = _symbolic_factor_list(together(expr, fraction=opt['fraction']), opt, method)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/polytools.py", line 5938, in _symbolic_factor_list
_coeff, _factors = func()
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/polytools.py", line 3303, in factor_list
coeff, factors = f.rep.factor_list()
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/polyclasses.py", line 796, in factor_list
coeff, factors = dmp_factor_list(f.rep, f.lev, f.dom)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/factortools.py", line 1269, in dmp_factor_list
coeff, factors = dmp_ext_factor(f, u, K0)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/factortools.py", line 1150, in dmp_ext_factor
s, g, r = dmp_sqf_norm(F, u, K)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/sqfreetools.py", line 164, in dmp_sqf_norm
if dmp_sqf_p(r, u, K.dom):
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/sqfreetools.py", line 75, in dmp_sqf_p
return not dmp_degree(dmp_gcd(f, dmp_diff(f, 1, u, K), u, K), u)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/euclidtools.py", line 1626, in dmp_gcd
return dmp_inner_gcd(f, g, u, K)[0]
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/euclidtools.py", line 1585, in dmp_inner_gcd
h, cff, cfg = _dmp_inner_gcd(f, g, u, K)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/euclidtools.py", line 1546, in _dmp_inner_gcd
return dmp_qq_heu_gcd(f, g, u, K)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/euclidtools.py", line 1458, in dmp_qq_heu_gcd
h, cff, cfg = dmp_zz_heu_gcd(f, g, u, K1)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/euclidtools.py", line 1335, in dmp_zz_heu_gcd
h, cff, cfg = dmp_zz_heu_gcd(ff, gg, v, K)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/euclidtools.py", line 1310, in dmp_zz_heu_gcd
return dup_zz_heu_gcd(f, g, K)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/euclidtools.py", line 1247, in dup_zz_heu_gcd
x = 73794*x * K.sqrt(K.sqrt(x)) // 27011
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/domains/pythonintegerring.py", line 85, in sqrt
return python_sqrt(a)
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/sympy/polys/domains/groundtypes.py", line 76, in python_sqrt
return int(mlib.isqrt(n))
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/mpmath/libmp/libintmath.py", line 293, in isqrt_python
return sqrtrem_python(x)[0]
File "/home/alan/local/anaconda3/envs/nbodykit-env/lib/python3.7/site-packages/mpmath/libmp/libintmath.py", line 271, in sqrtrem_python
def sqrtrem_python(x):
KeyboardInterrupt

[AlanPearl]

Update

Everything works now after running pip install sympy==1.5.1 (downgraded from 1.6)! Now I'm not sure if this is a bug with nbodykit or sympy

[rainwoodman]

@eelregit, @nickhand any comments?

[rainwoodman]

The travis builds are stalling in a similar looking location: https://travis-ci.org/github/bccp/nbodykit/jobs/703292960

[eelregit]

Seems like a sympy problem. It'd be nice to submit a minimal example to sympy issues, after figuring out which m it got stuck on, and which factor() of the two on the line.

[rainwoodman]

Yes. I can confirm with the sympy master branch (1.7.0dev0) the three point test case no longer stalls.
I suspect https://github.com/sympy/sympy/pull/19537 is the fix. Not sure about their release cycle -- if they'll cut a release soon, then we probably don't need to do anything on our side.

[AlanPearl]

Sounds good. I looked into it a little bit more and it looks like factoring using the imaginary extension is what causes the issue. I opened an issue here https://github.com/sympy/sympy/issues/19688#issue-650114094

[AlanPearl]

Woops, you're right, this does indeed work on the master branch of sympy!

[rainwoodman]

Thanks for pinning this down to the imaginary extension!