Hi Ondrej, Here is my attempt to debug: http://nbviewer.ipython.org/94ee433a4dca26e3fee1
Regards, Thilina. On Tue, Jul 16, 2013 at 3:01 AM, Thilina Rathnayake <[email protected]>wrote: > Hi Ondrej, > > Here is the notebook I created: > http://nbviewer.ipython.org/3dfe97025f06876820b4 > > Please let me know if anything needs to be changed. > > Regards, > Thilina. > > > On Tue, Jul 16, 2013 at 2:26 AM, Thilina Rathnayake < > [email protected]> wrote: > >> Hi Ondrej, >> >> Thank you very much for the quick response. >> >> I copied these two functions from my current diophantine.py and I forgot >> to copy the lines containing >> the imports. Sorry about that. >> >> I think both the outputs [2] and [3] are correct. In [3], since we called >> `descent(2, 1)` the equation is >> `w**2 = 2*x**2 + y**2`, and the output `(x, y, w) = (0, 1, 1)` satisfies >> the above equation. As I noticed, >> outputs generated by `descent()` when either `A` or `B` equals 1 are >> correct. Things go wrong when >> both A and B are different from 1. For example: >> >> >>>descent(2, 7) # equation w**2 = 2*x**2 + 7*y**2 >> (1, 7, 3) >> >> Here,` (x, y, w) = (1, 7, 3)` and it is wrong because `3**2 != 2 * (1**2) >> + 7 * (7**2)`. >> >> I'll send you a notebook as soon as possible. >> >> Thank you and best regards, >> Thilina. >> >> >> >> On Tue, Jul 16, 2013 at 2:06 AM, Ondřej Čertík >> <[email protected]>wrote: >> >>> On Mon, Jul 15, 2013 at 2:31 PM, Ondřej Čertík <[email protected]> >>> wrote: >>> > Thilina, >>> > >>> > Thanks. You can use ipython notebook for exactly this exploration (let >>> > me know if you have any problems setting it up). Here it is: >>> > >>> > http://nbviewer.ipython.org/6003151 >>> >>> P.S. One has to put >>> >>> from math import floor >>> >>> on top, and also the "divisors" function is not defined. Thilina, in >>> the meantime, you can send >>> a notebook that produces the wrong results and is executable and we >>> can then debug it. >>> >>> IPython notebook is very useful for such exploratory debugging. >>> >>> Ondrej >>> >>> > >>> > So the output [2] is correct, isn't it? But the output [3] seems >>> > incorrect to me, isn't it? Let me know which other output you found >>> > incorrect. >>> > >>> > I'll have a look at the book later. >>> > >>> > Ondrej >>> > >>> > On Mon, Jul 15, 2013 at 1:22 PM, Thilina Rathnayake >>> > <[email protected]> wrote: >>> >> I attached the code herewith as a .py file so you can easily run it. >>> >> >>> >> >>> >> On Tue, Jul 16, 2013 at 12:47 AM, Thilina Rathnayake >>> >> <[email protected]> wrote: >>> >>> >>> >>> Hi Ondrej, >>> >>> >>> >>> Thank you very much for the words of encouragement. I am >>> implementing the >>> >>> descent method. >>> >>> I have used the solver you pointed out to validate answers for >>> quadratic >>> >>> Diophantine equations. >>> >>> We need something which solve the the ternary quadratic equations. >>> >>> >>> >>> I will continue implementing the algorithm as you suggested. I think >>> >>> that's the best way. >>> >>> I am confident that we will find some way to validate the results. I >>> >>> couldn't completely check PARI/GP >>> >>> and sage. They might contain a solver for ternary quadratic forms. >>> >>> >>> >>> I implemented the descent algorithm but it doesn't return the correct >>> >>> results. Here is the code. >>> >>> >>> >>> 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 a solution to the above >>> >>> equation. >>> >>> This solution is used as a base to calculate all the solutions. >>> >>> """ >>> >>> if abs(A) > abs(B): >>> >>> y, x, w = descent(B, A) >>> >>> return (x, y, w) >>> >>> >>> >>> if A == 1: >>> >>> return (1, 0, 1) >>> >>> >>> >>> if B == 1: >>> >>> return (0, 1, 1) >>> >>> >>> >>> if quadratic_congruence(A, B, 0.5) != None: >>> >>> r = quadratic_congruence(A, B, 0.5) >>> >>> Q = (r**2 - A)//B >>> >>> div = divisors(Q) >>> >>> >>> >>> B_0 = 0 >>> >>> >>> >>> for i in div: >>> >>> if isinstance(sqrt(Q // i), Integer): >>> >>> B_0, d = i, sqrt(Q // i) >>> >>> break >>> >>> >>> >>> if B_0 > 0: >>> >>> X, Y, W = descent(A, B_0) >>> >>> return (r*X + W, B*d*Y, A*X + r*W) >>> >>> >>> >>> def quadratic_congruence(a, m, k=1.0): >>> >>> """ >>> >>> 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. >>> >>> """ >>> >>> upper_bound = floor(abs(k*m)) + 1 >>> >>> >>> >>> for i in range(upper_bound): >>> >>> if (i**2 - a)% abs(m) == 0: >>> >>> return i >>> >>> >>> >>> return None >>> >>> >>> >>> Perhaps you could take a look at it and tell me where I have done >>> wrong. I >>> >>> tried to follow >>> >>> almost the same symbols as in the algorithm. Here B_0 corresponds to >>> B' in >>> >>> the algorithm. >>> >>> Please take a look at it when you have some free time. >>> >>> >>> >>> Thank you and best regards, >>> >>> Thilina >>> >>> >>> >>> >>> >>> On Mon, Jul 15, 2013 at 9:01 PM, Ondřej Čertík < >>> [email protected]> >>> >>> wrote: >>> >>>> >>> >>>> I see, it's the Chapter IV. >>> >>>> Are you going to implement the "Sieving method" (IV.3.1)? Or the >>> >>>> "descent method"? >>> >>>> >>> >>>> So, if I take the Pythagorean formula a^2 + b^2 = c^2, your new >>> solver >>> >>>> will be able to give me all integer solutions to this? That would be >>> >>>> really cool. >>> >>>> >>> >>>> As far as a solver, I found this one >>> >>>> (http://www.alpertron.com.ar/QUAD.HTM). >>> >>>> But it can only do equations for "x" an "y". You need equations for >>> >>>> "x", "y" and "z". >>> >>>> >>> >>>> I would suggest that you implement the solver, and then we just >>> check >>> >>>> the solutions numerically, that should be quite easy. We will not >>> know >>> >>>> for sure that you got all the solutions, but just cite the algorithm >>> >>>> reference in the docstring, so that users know what it does. That >>> >>>> would be good start. If somebody later finds some other software >>> that >>> >>>> can do this, or some other reference with examples and results, we >>> can >>> >>>> make sure that it returns all the solutions. >>> >>>> >>> >>>> Ondrej >>> >>>> >>> >>>> On Sun, Jul 14, 2013 at 11:21 PM, Thilina Rathnayake >>> >>>> <[email protected]> wrote: >>> >>>> > This is what I use. >>> >>>> > >>> >>>> > “The Algorithmic resolution of Diophantine Equations”, Nigel P. >>> Smart, >>> >>>> > Chapter IV, Quadratic ternary forms. >>> >>>> > >>> >>>> > Here is the book: >>> >>>> > >>> >>>> > >>> >>>> > >>> >>>> > On Mon, Jul 15, 2013 at 2:20 AM, Aaron Meurer <[email protected] >>> > >>> >>>> > wrote: >>> >>>> >> >>> >>>> >> On Sun, Jul 14, 2013 at 12:09 AM, Thilina Rathnayake >>> >>>> >> <[email protected]> wrote: >>> >>>> >> > Hi All, >>> >>>> >> > >>> >>>> >> > I am planning to implement ternary quadratic forms, i.e >>> equations of >>> >>>> >> > the >>> >>>> >> > form, >>> >>>> >> > a*x**2 + by**2 + cz**2 + fxy + gyz + hzx = 0. It would be >>> better If >>> >>>> >> > I >>> >>>> >> > can >>> >>>> >> > find >>> >>>> >> > a system which currently implement this so I can validate my >>> >>>> >> > results. If >>> >>>> >> > you >>> >>>> >> > know >>> >>>> >> > of any system which solves this or a source which have good >>> >>>> >> > literature >>> >>>> >> > on >>> >>>> >> > the problem, >>> >>>> >> > please let me know. >>> >>>> >> >>> >>>> >> What is the literature that you plan to use to get the algorithm >>> in >>> >>>> >> the first place? >>> >>>> >> >>> >>>> >> Aaron Meurer >>> >>>> >> >>> >>>> >> > >>> >>>> >> > Also I have to solve the quadratic congruence x**2 = D (mod m) >>> as a >>> >>>> >> > sub >>> >>>> >> > problem to implement the algorithm I found. I found a few >>> algorithms >>> >>>> >> > on >>> >>>> >> > this >>> >>>> >> > but >>> >>>> >> > none of them explain precisely how to solve the general >>> equation, >>> >>>> >> > all >>> >>>> >> > they >>> >>>> >> > do is >>> >>>> >> > solve the equation for m prime and gcd(D, m) = 1 and just ask >>> to use >>> >>>> >> > Chinese >>> >>>> >> > remainder theorem to combine the results to solve for the >>> general >>> >>>> >> > case >>> >>>> >> > where >>> >>>> >> > m >>> >>>> >> > is not prime and gcd(D, m) not necessarily equal to 1. Any >>> help on >>> >>>> >> > this >>> >>>> >> > would >>> >>>> >> > be highly appreciated. >>> >>>> >> > >>> >>>> >> > 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. >>> >>>> >> > >>> >>>> >> > >>> >>>> >> >>> >>>> >> -- >>> >>>> >> 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 >>> . >>> >>>> >> >>> >>>> >> >>> >>>> > >>> >>>> > -- >>> >>>> > 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. >>> >>>> > >>> >>>> > >>> >>>> >>> >>>> -- >>> >>>> 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. >>> >>>> >>> >>>> >>> >>> >>> >> >>> >> -- >>> >> 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. >>> >> >>> >> >>> >>> -- >>> 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. >>> >>> >>> >> > -- 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.
