In the attached file there is a new version of ``descent``; I followed your
suggestion about
reducing the input to the square free part; in this way the lowest modular
square root is
always used; there is no need to use others; so the code in Smart seems to
be correct,
provided the input is reduced to the square free part.
It is used the fast ``sqrt_mod`` in https://github.com/sympy/sympy/pull/2307
On Saturday, July 20, 2013 7:40:05 PM UTC+2, Thilina Rathnayake wrote:
>
> Thanks Ondrej for putting the link of the PR. I am sorry I forgot to put
> it.
>
> I am really grateful If other's can play with this and give feedback on
> how this can be improved,
> especially reports any bugs in the programme.
>
> Thanks in advance.
>
>
>
> On Sat, Jul 20, 2013 at 11:00 PM, Ondřej Čertík
> <[email protected]<javascript:>
> > wrote:
>
>> On Sat, Jul 20, 2013 at 11:25 AM, Thilina Rathnayake
>> <[email protected] <javascript:>> wrote:
>> > Hi Ondrej,
>> >
>> > I am pleased to say that I managed to implement the solutions for
>> quadratic
>> > ternary forms.
>> > Only thing we have to verify is that whether they are complete. I found
>> one
>> > incident in the
>> > current implementation where it returned a partial solution. I added it
>> as a
>> > different test. I plan
>> > to implement the algorithm Implemented in MAGMA in the upcoming days,
>> and
>> > also it seems like
>> > MAGMA's algorithm is faster. So we can make it the main algorithm used
>> by
>> > `diop_solve()`.
>> >
>> > There is one major factor affecting the speed of the algorithm, the
>> overhead
>> > of the function
>> > `quadratic_congruence()` which solves a quadratic congruence. This
>> function
>> > was also used in solving
>> > binary quadratic forms also. I found a link to a reference material on
>> the
>> > problem.
>> >
>> > Also the parametric solution returned by the solver would be more
>> elegant if
>> > we could reduce the
>> > basic solution we found using the `descent()` by the Hozler's reduction
>> > algorithm.
>> >
>> > These are the area's I am going to focus in the future.
>>
>> Very cool, I agree. Thanks for the work, you have done excellent
>> progress.
>> I'll try to play with your branch in the next few days.
>>
>> >
>> > (PS - You have already commented in the pull request. Thank you for
>> that.
>> > I'll get to it as soon as
>> > I finish writing this week's blog post)
>>
>> Here is the PR if others would like to comment:
>>
>> https://github.com/sympy/sympy/pull/2303
>>
>> Ondrej
>>
>> --
>> 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] <javascript:>.
>> To post to this group, send email to [email protected] <javascript:>
>> .
>> Visit this group at http://groups.google.com/group/sympy.
>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>>
>>
>
--
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 time import time
from sympy import floor, sqrt, Integer, sign, S, factorint, gcd
def to_squarefree(n):
"""
return ``sf, q``, where ``n = sf * q**2``
"""
f = factorint(n)
q = S(1)
for b, e in f.items():
if e > 1:
q = q * b**(e // 2)
return n // q**2, q
def descent(A, B):
"""
Uses Lagrange's method to find 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)
sf1, q1 = to_squarefree(A)
sf2, q2 = to_squarefree(B)
if q1 != 1 or q2 != 1:
x1, y1, w1 = descent(sf1, sf2)
if x1 is None:
return None, None, None
else:
return (q2*x1, q1*y1, q1*q2*w1)
start = 0
r = quadratic_congruence(A, B)
if r is None:
return None, None, None
start = r + 1
Q = (r**2 - A) // B
B_0, d = to_squarefree(Q)
X, Y, W = descent(A, B_0)
if X is None:
return None, None, None
return ((r*X - W), Y*(B_0*d), (-A*X + r*W))
def quadratic_congruence(a, m):
from sympy.ntheory.residue_ntheory import sqrt_mod
m = abs(m)
rv = sqrt_mod(a, m)
if not rv:
return None
return rv[0]
def quadratic_congruence0(a, m):
m = abs(m)
for i in range(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
def test2():
a = 32101
a = 3210
a = 321010
for b in range(1, 1000):
r = descent(a, b)
if r[0]:
x, y, w = r
assert a*x**2 + b*y**2 == w**2
print 'a=%s b=%s r=%s'%(a, b, r)
t0 = time()
test2()
t1 = time()
print '%.2f' %(t1 - t0)
