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 <thilina.r...@gmail.com>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 <ondrej.cer...@gmail.com>wrote: > >> On Mon, Jul 15, 2013 at 2:31 PM, Ondřej Čertík <ondrej.cer...@gmail.com> >> 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 >> > <thilina.r...@gmail.com> 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 >> >> <thilina.r...@gmail.com> 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 < >> ondrej.cer...@gmail.com> >> >>> 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 >> >>>> <thilina.r...@gmail.com> 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 <asmeu...@gmail.com> >> >>>> > wrote: >> >>>> >> >> >>>> >> On Sun, Jul 14, 2013 at 12:09 AM, Thilina Rathnayake >> >>>> >> <thilina.r...@gmail.com> 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 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. >> >>>> >> > >> >>>> >> > >> >>>> >> >> >>>> >> -- >> >>>> >> 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. >> >>>> >> >> >>>> >> >> >>>> > >> >>>> > -- >> >>>> > 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. >> >>>> > >> >>>> > >> >>>> >> >>>> -- >> >>>> 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. >> >>>> >> >>>> >> >>> >> >> >> >> -- >> >> 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. >> >> >> >> >> >> -- >> 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. >> >> >> > -- 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.
