Taxicab count for large values

import ma.primee as PR
import itertools as IT
import operator
import math
import cubes as CC
from time import clock
combinations = IT.combinations
mul = operator.mul

def solveq(a, b, c):
    """ solves ax^2 + bx + c = 0"""
    y = b ** 2  - 4 * a * c
    if y < 0:
        return False
    sqrty = math.sqrt(y)
    return ((-b - sqrty) / (2*a), (-b + sqrty) / (2*a))

def check_sol(sol, p):
    """checks if solution of quadratic equation is integer and in (0,p)"""
    res= []
    if sol:
        a = sol[0]
        inta = int(a)
        if abs(inta - a) > 1e-15:
            return 0
        if 0 < inta < p:
            res.append(inta)
        inta = sol[1]
        if 0 < inta < p:
            res.append(inta)
    return res

def posint(a, x1, x2):
    """check if a is an integer in (x1,x2)"""
    if not x1 < a < x2:
        return False
    if abs(int(a) - a) < 1e-15:
        return int(a)

def combs(n):
    """decompose n in (p,q) with p < q and pq = n
        combs returns all possible values for p
    """
    # search through primefactors
    # make sure that no duplications exist:
    #    go from low to higher factors and keep adding factors until p > sqrt(n)
    sqrtn = math.sqrt(n)
    fac =  [(k, len(list(g))) for k, g in IT.groupby(PR.rho_factors(n))]
    pvalues = []
    Q = [(0, 0, 1)]
    while Q:
        # j: current prime, m: exponent of current prime, p: factor of n
        j, m, p = Q.pop()
        for i, (fi, ni) in enumerate(fac[j:], j):
            mi = m if i == j else 0
            if mi < ni:
                pi = p * fi
                if pi <= sqrtn:
                    pvalues.append(pi)
                    Q.append((i, mi + 1, pi))
    return pvalues

def count_cubes(n):
    """method: https://cs.uwaterloo.ca/journals/JIS/wilson10.html#HW54"""
    sol = set()
    for p in combs(n):
        aa = check_sol(solveq(3, -3*p, p**2 - n // p), p)
        if aa:
            for a in aa:
                b = p - a
                if a > b:
                    a, b = b, a
                sol.add((a, p - a))
    return sol

if __name__ == '__main__':
    t0 = clock()
    print count_cubes(19242016141146624)
    print clock() - t0
    t0 = clock()
    print CC.check_cubes(19242016141146624,4)
    print clock() - t0

    t0 = clock()
    print count_cubes(87539319)
    print clock() - t0
    t0 = clock()
    print CC.check_cubes(87539319 ,3)
    print clock() - t0