def gf_csolve(f, n):
 from sympy.polys.domains import ZZ
 P = factorint(n)
 X = [csolve_prime(f, p, e) for (p, e) in P.items()]
 pools = list(map(tuple, X))
 perms = [[]]
 for pool in pools:
  perms = [(x + [y]) for x in perms for y in pool]
 dist_factors = [pow(p, e) for (p, e) in P.items()]
 return sorted([gf_crt(per, dist_factors, ZZ) for per in perms])

def gf_csolve ( f , n ) : from sympy . polys . domains import ZZ P = factorint ( n ) X = [ csolve_prime ( f , p , e ) for ( p , e ) in P . items ( ) ] pools = list ( map ( tuple , X ) ) perms = [ [ ] ] for pool in pools : perms = [ ( x + [ y ] ) for x in perms for y in pool ] dist_factors = [ pow ( p , e ) for ( p , e ) in P . items ( ) ] return sorted ( [ gf_crt ( per , dist_factors , ZZ ) for per in perms ] )
to solve f(x) congruent 0 mod(n) .