Alasdair, Thanks! On 9 feb, 11:54, Alasdair <amc...@gmail.com> wrote: > Do a google search on "Cornacchia's algorithm". Shouldn't be too hard > to program in Sage (if it isn't there already). > > Alasdair > > On Feb 8, 4:07 pm, Rolandb <rola...@planet.nl> wrote: > > > Thanks William, > > > Actually I try to solve it for different x and n. A typical example of > > (x,n) is: > > > %time > > bruteforce(7^10*29^5,973) > > [(3899224, 2437015)] > > CPU time: 25.88 s, Wall time: 26.12 s > > > Roland > > > On 7 feb, 22:29, William Stein <wst...@gmail.com> wrote: > > > > On Sun, Feb 7, 2010 at 12:35 PM, Rolandb <rola...@planet.nl> wrote: > > > > Hi, > > > > Consider Euler’s famous expression x=a^2+n*b^2. I want to solve a and > > > > b, given x and n. Because solve_mod is totally broken, I tried to > > > > develop a clever method myself. Counterintuitively, I found that the > > > > most simple brute force method in SAGE is relatively fast. > > > > > def bruteforce(number,n): > > > > out=[] > > > > for b in xrange(isqrt(number/n)): > > > > (bool,a)=is_square(number-n*b^2,True) > > > > if bool: out.append((a,b)) > > > > return out > > > > > My questions are: > > > > - Why is the brute force method relatively fast? > > > > - Is there a clever approach (not dependent on solve_mod)? > > > > > Thanks in advance! Roland > > > > We have > > > > x=a^2+n*b^2 = (a-sqrt(n)*b)*(a+sqrt(n)*b) = Norm(a+sqrt(n)*b), > > > > where the norm is in the quadratic field Q(sqrt(n)). Thus a solution > > > corresponds to a factorization of the ideal generated by x in the ring > > > of integers of Q(sqrt(n)). So here is some code related to your > > > question above, which you may or may not find relevant, depending on > > > whether you're making n big or x: > > > > def eqn(x, n): > > > """ > > > Solve x = a^2+n*b^2 when x is prime, n is not a square, and there > > > is a solution. > > > > EXAMPLES: > > > > sage: a,b = eqn(13, -3); a,b > > > (4, 1) > > > sage: a^2 + -3*b^2 > > > 13 > > > sage: a,b = eqn(113, 97); a,b > > > (-4, 1) > > > sage: a^2 + 97*b^2 > > > 113 > > > """ > > > proof.number_field(False) > > > assert is_prime(x), "x must be prime" > > > assert not is_square(n), "n must not be a square" > > > R.<X> = QQ[] > > > K.<a> = NumberField(X^2 + n) > > > F = K.factor(x) > > > if len(F) == 1 and F[0][1] == 1: > > > raise ValueError, "no solution (x is inert)" > > > A = F[0][0].gens_reduced() > > > if len(A) > 1: > > > raise ValueError, "no solution (ideal isn't principal)" > > > return A[0][0], A[0][1]
-- To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org
