The attached version works in the mentioned cases; the bug is fixed
reducing to the case ``A > 0, B > 0``
On Thursday, July 18, 2013 11:28:22 AM UTC+2, Thilina Rathnayake wrote:
>
> I should try making input coefficients to be square free. I think the
> algorithm in general,
> does not work when the coefficients are not square free.
>
>
>
> On Thu, Jul 18, 2013 at 2:53 PM, Thilina Rathnayake
> <thilin...@gmail.com<javascript:>
> > wrote:
>
>> We can add descent(23, 616) to the bug list too. It also return (None,
>> None, None)
>> but (6, 1, 38) is a solution. By finding more and more test cases which
>> fail, I think we can
>> identify a pattern and then the cause for the bug.
>>
>> Regards,
>> Thilina
>>
>
>
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sympy+unsubscr...@googlegroups.com.
To post to this group, send email to sympy@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
For more options, visit https://groups.google.com/groups/opt_out.
from sympy import floor, sqrt, Integer, sign, divisors, S, factorint, sympify
def descent(A, B):
"""
Uses Lagrange's method to finda a solution to $w^2 = Ax^2 + By^2$.
Output a tuple $(x_0, y_0, z_0)$ which is a solution to the above equation.
This solution is used as a base to calculate all the other solutions.
"""
def _coeff(f):
c1 = 1
for b, e in f.items():
if e % 2 == 1:
c1 = c1 * b**e
else:
c1 = c1 * b**(e//2)
return c1
A = sympify(A)
B = sympify(B)
if A > B >= 0:
y, x, w = descent(B, A)
return (x, y, w)
if A > 0 and B < 0:
# x**2 = w**2/A - B*y**2/A
# w = c1*w1; y = c2*y1
# x**2 = c1**2/A*w1**2 - B*c2**2/A*y1**2
f = factorint(A)
c1 = _coeff(f)
b1 = -B/A
f = factorint(b1.q)
c2 = _coeff(f)
w1, y1, x = descent(c1**2/A, -B*c2**2/A)
return (x, c2*y1, c1*w1)
if A < 0 and B > 0:
y, x, w = descent(B, A)
return (x, y, w)
if A == 1:
return (S.One, 0, S.One)
if B == 1:
return (0, S.One, S.One)
start = 0
while 1:
r = quadratic_congruence(A, B, start)
if r is None:
break
start = r + 1
Q = (r**2 - A) // B
if Q == 0:
continue
div = divisors(Q)
B_0 = None
for i in div:
if isinstance(sqrt(abs(Q) // i), Integer):
B_0, d = sign(Q)*i, sqrt(abs(Q) // i)
break
if B_0 != None:
X, Y, W = descent(A, B_0)
if X is None:
break
return ((r*X - W), Y*(B_0*d), (-A*X + r*W))
return None, None, None
def quadratic_congruence(a, m, start):
"""
Solves the quadratic congruence $x^2 \equiv a \ (mod \ m)$. Returns the
first solution in the range $0 .. \lfloor k*m \rfloor$.
Return None if solutions do not exist. Currently uses bruteforce.
Good enough for ``m`` sufficiently small.
TODO: An efficient algorithm should be implemented.
"""
m = abs(m)
for i in range(start, m // 2 + 1 if m%2 == 0 else m // 2 + 2):
if (i**2 - a) % m == 0:
return i
return None
def test1():
u = [(5, 4), (13, 23), (3, -11), (41, -113), (4, -7), (234, -65601),
(-7, 4)]
for a, b in u:
x, y, w = descent(a, b)
assert a*x**2 + b*y**2 == w**2
test1()
