|
|
Re: cube root of a given number
Posted:
Aug 17, 2007 12:49 AM
|
|
On Aug 12, 2:33?pm, "sttscitr...@tesco.net" <sttscitr...@tesco.net>
wrote:
> On 12 Aug, 00:56, rich burge <r3...@aol.com> wrote:
>
>
>
>
>
> > On Aug 11, 4:11?pm, "sttscitr...@tesco.net" <sttscitr...@tesco.net>
>
> > \\
> > \\
> > \\ This script finds a solution to
> > P_3(k;x,y,z)=x^3+ky^3+k^2z^3-3kxyz=1 for k=1000700. Here k,x,y,z are
> > \\ positive integers.
> > \\
>
> > nm1(k,x)=
> > {
> > \\
> > \\ k is k. x is an integer triple, where x[3]/x[1] and x[2]/x[1] are
> > the reduced approximations
> > \\ to k^(2/3) and k^(1/3) found below.
> > \\
> > return((x[3]^3+k*x[2]^3+k^2*x[1]^3-3*k*x[1]*x[2]*x[3])%modk); \\
> > return the norm mod k^6
> > \\
>
> > }
>
> > lessthan(A,B)=
> > \\
> > \\ This routine determines which of the two sets of column norms
> > corresponding to the reduced approximation
> > \\ matricies A and B is 'smallest'. Returns true iff A <* B.
> > \\
> > \\ The idea here is this: given two matricies, A and B, the columns of
> > which are reduced approximation
> > \\ vectors, determine which one has the smallest maximum column norm.
> > In the case where the largest column
> > \\ norms are equal, compare the second largest, ect.
> > \\
> > {
> > local(i);
>
> > \\(re)compute the norms of the columns of A and B
> > for(i=1,3,nms_A[i]=nm1(k,A[,i]~);nms_B[i]=nm1(k,B[,i]~)); \\ only 2
> > of these are new!!!
>
> > \\sortem,
> > nms_A_srt=vecsort(nms_A);
> > nms_B_srt=vecsort(nms_B);
>
> > \\ sort through the ties,
> > i=3;
> > while((i>1)&&nms_A_srt[i]==nms_B_srt[i],i=i-1);
>
> > \\and return the smallest one.
> > return(nms_A_srt[i]<nms_B_srt[i])
>
> > }
>
> > iter()=
> > \\
> > \\ This routine acts on D and A. D[1,] is initially an integer
> > approximation vector to (1,k^(1/3),k^(2/3)). A is
> > \\ the identity matrix. Given D, two reduction matricies B_s and B
> > are computed, along with the two possible
> > \\ new values of A. lessthan is used to select the correct new A and
> > B then D is reduced.
> > \\
> > \\ The entries of A, if not controlled, grow rapidly and result in an
> > increasing execution time. I.e if a large
> > \\ number of iterations are needed to find a solution the process of
> > determining what swap to make is painfully slow.
> > \\ Note, however, that lessthan(A,B)=lessthan(A mod k^6,B mod k^6),
> > given each column norm is bounded by k^6.
> > \\ So the computational effort required to make the swap decision is
> > bounded.
> > \\
>
> > {
> > B_s=[1,-(D[1,1]\D[1,2]),0;0,-(D[1,3]\D[1,2]),1;0,1,0];
> > A_s=(A*(B_s)^(-1))%modk;
> > \\
> > B=[-(D[1,2]\D[1,1]),1,0;-(D[1,3]\D[1,1]),0,1;1,0,0];
> > A=(A*B^(-1))%modk;
> > \\
> > \\ this is the key: directed swaps
> > \\
> > if(lessthan(A_s,A),A=A_s;B=B_s);
> > \\
> > \\ now reduce D accordingly
> > \\
> > D=D*B~;
> > \\
> > \\ Every fifteenth iteration, add some more precision. 15 is a
> > guess!!
> > \\
> > if(itercnt%15==0,D=C10*D);
> > \\
> > \\ Keep the unreduced values in A_act
> > \\
> > A_act=A_act*B^(-1);
> > \\
> > \\ (Re)compute the norm and stop if = 1, printing z where
> > x=round(z*k^(2/3)) and y=round(z*k^(1/3)).
> > \\
> > tmp=nm1(k,A[,3]~);
> > if((tmp==1)&&(A_act[1,3]>10),print(vecmin(A_act[,
> > 3]));return(0),return(1));
>
> > }
>
> > find_p3()=
> > {
> > default(timer,1);
> > k=1000700;
> > modk=k^6;
> > A=matid(3);
> > A_act=A;
> > nms_A=vector(3);
> > nms_B=nms_A;
> > \\ C is the matrix used to approximate [1,k^(1/3),k^(2/3)]
> > \\ The entries are derived from the small rational solution
> > [900630049/49,9004200/49,90021/49] to
> > \\ P_3(k,x,y,z)=1. These small rational solutions are fairly easy to
> > find.
> > \\
> > C=[900630049,9004200,90021;
> > k*90021,900630049,9004200;
> > k*9004200,k*90021,900630049]~;
> > C10=C;
> > D=C;
> > itercnt=0;
> > while(iter(),itercnt=itercnt+1);
> > print();
>
> > }
>
> > find_p3;
>
> > \\ end script
>
> Very impressive. But I can't say I understand
> how it all works at the moment.
> What's the largest k you have found solutions for ?
The largest solution I have found so far is for k=1000007000 . It
takes about 11 million iterations and can be done in about 90 minutes
(w/o updating A_act). This requires a real good choice for C. I was
hoping to find enough particular solutions to the family 10^(3*e)
+7*10^e to find a parametric solution, but this does not seem feasible
at this point.
Rich
|
|
