On Tuesday, July 16, 2013 1:38:02 PM UTC+2, mario wrote:
>
> In the attached file there is my attempt; it seems to work; it seems to me
> that there is a bug
> in the code in the book; in line 4 it should be
> ``Find r in [0,..., ceiling(abs(B)/2] such that r**2 == A mod B``
>
>
>
> On Tuesday, July 16, 2013 2:13:32 AM UTC+2, Thilina Rathnayake wrote:
>>
>> Hi Ondrej,
>>
>> Here is my attempt to debug:
>> http://nbviewer.ipython.org/94ee433a4dca26e3fee1
>>
>> 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 [email protected].
To post to this group, send email to [email protected].
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
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.
"""
if abs(A) > abs(B):
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)]
for a, b in u:
x, y, w = descent(a, b)
assert a*x**2 + b*y**2 == w**2
test1()
