I found a system (which varies its length as a function of (p + q-4) / 8) which generates all the numbers N = 8 * G + 3 = p * q so by inserting our number to factor in the underlying system instead of N we will have our p, factor of N.
System example up to x = (p + q-4) / 8 <= 3
var('x y p a A b B c z Z w W N v V')
eq0 = N -27 == 0
eq1 = (2*(3*x+1)-3*(a)+1)/24+3/2*x*(x+1) - A == 0
eq2 = A+b*(b-1)/2 - (2*x*(x+1)) == 0
eq3 = (2*(3*(A)+1)-3*(b)+1)/24+3/2*x*(x+1) - B == 0
eq4 = B+c*(c-1)/2 - (2*x*(x+1)) == 0
eq5 = (2*(3*(B)+1)-3*(c)+1)/24+3/2*x*(x+1) - (2*x*(x+1)) == 0
eq6 = a-(2*x+1) == 0
eq7 = (2*(3*(N-3)/8+1)-3*(z)+1)/24+3/2*x*(x+1) -((N-3)/8+Z) == 0
eq8 = (N-3)/8+z*(z-1)/2 - (2*x*(x+1)) == 0
eq9 = (2*(3*((N-3)/8+Z)+1)-3*(w)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W) == 0
eq10 = (N-3)/8+Z+w*(w-1)/2 - (2*x*(x+1)) == 0
eq11 = (2*(3*((N-3)/8+Z+W)+1)-3*(v)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W+V) == 0
eq12 = (N-3)/8+Z+W+v*(v-1)/2 - (2*x*(x+1)) == 0
eq13 = (N-3)/8+Z+W+V - (2*x*(x+1)) == 0
eq14 = (N-3)/8+y*(y-1)/2 - 2*x*(x+1) == 0
eq15 = 4*x+1-2*(y-1)-p == 0
solutions = solve([eq0,eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12,eq13,eq14,eq15],x,y,p,a,A,b,B,c,z,Z,w,W,N,v,V)
sol = solutions
print(sol)
System example up to x = (p + q-4) / 8 <= 7
var('x y p a A b B c C d z Z w W N v V u U')
eq0 = N-59 == 0
eq1 = (2*(3*x+1)-3*(a)+1)/24+3/2*x*(x+1) - A == 0
eq2 = A+b*(b-1)/2 - (2*x*(x+1)) == 0
eq3 = (2*(3*(A)+1)-3*(b)+1)/24+3/2*x*(x+1) - B == 0
eq4 = B+c*(c-1)/2 - (2*x*(x+1)) == 0
eq5 = (2*(3*(B)+1)-3*(c)+1)/24+3/2*x*(x+1) - C == 0
eq6 = C+d*(d-1)/2 - (2*x*(x+1)) == 0
eq7 = (2*(3*(C)+1)-3*(d)+1)/24+3/2*x*(x+1) - (2*x*(x+1)) == 0
eq8 = a-(2*x+1) == 0
eq9 = (2*(3*(N-3)/8+1)-3*(z)+1)/24+3/2*x*(x+1) -((N-3)/8+Z) == 0
eq10 = (N-3)/8+z*(z-1)/2 - (2*x*(x+1)) == 0
eq11 = (2*(3*((N-3)/8+Z)+1)-3*(w)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W) == 0
eq12 = (N-3)/8+Z+w*(w-1)/2 - (2*x*(x+1)) == 0
eq13 = (2*(3*((N-3)/8+Z+W)+1)-3*(v)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W+V) == 0
eq14 = (N-3)/8+Z+W+v*(v-1)/2 - (2*x*(x+1)) == 0
eq15 = (2*(3*((N-3)/8+Z+W+V)+1)-3*(u)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W+V+U) == 0
eq16 = (N-3)/8+Z+W+V+u*(u-1)/2 - (2*x*(x+1)) == 0
eq17 = (N-3)/8+Z+W+V+U - (2*x*(x+1)) == 0
eq18 = (N-3)/8+y*(y-1)/2 - 2*x*(x+1) == 0
eq19 = 4*x+1-2*(y-1)-p == 0
solutions = solve([eq0,eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12,eq13,eq14,eq15,eq16,eq17,eq18,eq19],x,y,p,a,A,b,B,c,C,d,z,Z,w,W,N,v,V,u,U)
sol = solutions
print(sol)
System example up to x = (p + q-4) / 8 <= 15
var('x y p a A b B c C d D f z Z w W N v V u U t T')
eq0 = N - 123 == 0
eq1 = (2*(3*x+1)-3*(a)+1)/24+3/2*x*(x+1) - A == 0
eq2 = A+b*(b-1)/2 - (2*x*(x+1)) == 0
eq3 = (2*(3*(A)+1)-3*(b)+1)/24+3/2*x*(x+1) - B == 0
eq4 = B+c*(c-1)/2 - (2*x*(x+1)) == 0
eq5 = (2*(3*(B)+1)-3*(c)+1)/24+3/2*x*(x+1) - C == 0
eq6 = C+d*(d-1)/2 - (2*x*(x+1)) == 0
eq20 = (2*(3*(C)+1)-3*(d)+1)/24+3/2*x*(x+1) - D == 0
eq21 = D+f*(f-1)/2 - (2*x*(x+1)) == 0
eq7 = (2*(3*(D)+1)-3*(f)+1)/24+3/2*x*(x+1) - (2*x*(x+1)) == 0
eq8 = a-(2*x+1) == 0
eq9 = (2*(3*(N-3)/8+1)-3*(z)+1)/24+3/2*x*(x+1) -((N-3)/8+Z) == 0
eq10 = (N-3)/8+z*(z-1)/2 - (2*x*(x+1)) == 0
eq11 = (2*(3*((N-3)/8+Z)+1)-3*(w)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W) == 0
eq12 = (N-3)/8+Z+w*(w-1)/2 - (2*x*(x+1)) == 0
eq13 = (2*(3*((N-3)/8+Z+W)+1)-3*(v)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W+V) == 0
eq14 = (N-3)/8+Z+W+v*(v-1)/2 - (2*x*(x+1)) == 0
eq15 = (2*(3*((N-3)/8+Z+W+V)+1)-3*(u)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W+V+U) == 0
eq16 = (N-3)/8+Z+W+V+u*(u-1)/2 - (2*x*(x+1)) == 0
eq22 = (2*(3*((N-3)/8+Z+W+V+U)+1)-3*(t)+1)/24+3/2*x*(x+1) -((N-3)/8+Z+W+V+U+T) == 0
eq23 = (N-3)/8+Z+W+V+U+t*(t-1)/2 - (2*x*(x+1)) == 0
eq17 = (N-3)/8+Z+W+V+U+T - (2*x*(x+1)) == 0
eq18 = (N-3)/8+y*(y-1)/2 - 2*x*(x+1) == 0
eq19 = 4*x+1-2*(y-1)-p == 0
solutions = solve([eq0,eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12,eq13,eq14,eq15,eq16,eq17,eq18,eq19,eq20,eq21,eq22,eq23],x,y,p,a,A,b,B,c,C,d,D,f,z,Z,w,W,N,v,V,u,U,t,T)
sol = solutions
print(sol)
What mathematical procedure does sagemath use to solve these types of systems?
Is there a better process?
