I am trying to find an algorithm, similar to the Pell equation, that would solve for matricies of this form, that have a determinate of 1.
x nz ny
y x nz
z y x
Such a matrix would rapidly converge on the cube root of n. In the shown matrix, z is allowed to be fractional, divisible by some m, where m³ divides n². For example, the case of
n=12, x=55, y=24, z=10.5 works.
Some are quite large: n=1260, x=172226251, y=15945660, z=146337.5
The civilized approach to your problem is to use the built-in Pari/gp functions to calculate powers of fundamental units as mentioned in the comments. A more playful approach is to run the following script which finds a solution to the cubic pell equation $x^3+Ny^3+N^2z^3-3Nxyz=1$ for N=1260.
There are several observations about this algorithm that are perhaps worth noting here. First, rational matrices having your desired form can often be found quite easily. To see this, proceed as follows: on the way to finding a solution of determinant one, record two $(x,y,z)$ triples having equal norm (if they exist). Now put these triples into your special matrix form and multiply the larger by the inverse of the smaller. The result is a matrix of determinant one having all rational entries. For example, when $n=25$ at iteration 3 and 13, of 28, we get the pair of triples $(8,3,1)$ and $(38483,13161,4501)$ having norm 12. The resulting rational solution of norm one is $(6079/4,2079/4,711/4)$.
Second, while the algorithm is a two dimensional non-markovian continued fraction algorithm the "memory" requirement is not nearly as substantial as one might infer from looking at the code below. Experiments suggest that initializing the matrix B to B=1.0*A with a default real precision of, say, ten digits yields exactly the same results. So both the memory requirement to store B and the computational effort required to update B is quite a bit smaller than what is used below.
Third, precomputation of $n^{1/3}$ is not necessary. In the script below replace the vector D with powers of a matrix D' having the nice convergence properties you note above. Add more powers (by left multiplication) as needed. The matrices with rational entries found above are useful here as are others.
Finally, experiments suggest the norms of the intermediate approximations are always less than something near $n^6$. This fact can speed the computation significantly at times.
A_tmp=[[0, 0, 0; 1, 0, 0; 0, 1, 0], [1, 0, 0; 0, 0, 0; 0, 1, 0]]; B_tmp=A_tmp;
A=matid(3); B_inv=A; B=A;
alpha=vector(2);
default(realprecision,1000);
N=1260;
alpha[1]=N^(1/3); alpha[2]=N^(2/3);
D=vector(3);
D[1]=10^990;
for(i=2,3,D[i]=round(D[1]*alpha[i-1]));
done=0;
while(!done,
for(i=1,2,
for(j=1,3,A_tmp[i][j,3]=D[j]\D[i]);
B_tmp[i]=A_tmp[i]^(-1);
);
B_szs=vector(2,i,norml2(B_tmp[i]*B));
B_szs_prm=vecsort(B_szs,,1);
idx=B_szs_prm[1];
A=A*(A_tmp[idx]);
B=(B_tmp[idx]*B);
D=D*(B_tmp[idx])~;
done=(1==A[3,3]^3+N*A[2,3]^3+N^2*A[1,3]^3-3*N*A[1,3]*A[2,3]*A[3,3]);
);
print(A[,3]~);
Added 1/16:
I did find the following reference which may be helpful to you:
C. L. E. Wolfe, On the indeterminate cubic equation x 3 + Dy 3 + D 2 z 3 − 3Dxyz = 1, Univ. California pub. Math. 1 (1923) no. 16, 359–369. Available at: University of Michigan Historical Math Collection http://quod.lib.umich.edu/u/umhistmath/ACH6090.0001.001/381
As a check of the algorithm given above I computed the first 100 solutions (in 672 ms) and compared the results to the solutions given in the above paper. It seems that the solutions given for 72 and 82 are wrong. For D=85, I find a smaller solution.
72: 1263601, 303738, 73011
82: 19481839755856966061390701, 4484272707726892450472280, 1032176733268676625192495
85: 658895013725266441, 149856842965183254, 34082931143344968
