Нашел таким поиском:
(код pari/gp)
Код:
 my(V,P,p1,c,r,m,x,y,z,x0,y0,x1,y1,x2,y2,x3,y3,z3);
 V= [];
 forprime(p=3,100000,
  P=polrootsmod(1+'x^2+'x^3, p)~;
  for(i=1,#P,
   V= concat(V,P[i]);
   parfor(j=1,#V,
    p1= V[j].mod;
    if(p1<p,
     c= chinese(P[i],V[j]);
     c= chinese(Mod(4,9),c);
     r= lift(c); m= c.mod;
     for(k=0, 16, forstep(ksign=-1,1,2,
      x= r+k*ksign*m;
      y= 9*p1*p;
      z= -(1+x^2+x^3+y^2)/(x*y);
      if(z==floor(z),
       print("("x","y","z")    ",factorint(y),"    ",k"\n");
       x0= x; y0= y;
       x1= (y0^2+1)/x0; y1= y0;
       x2= x1; y2= (x1^3+x1+1)/y1;
       x3= (y2^2+1)/x2; y3= y2;
       z3= -(1+x3^2+x3^3+y3^2)/(x3*y3);
       if(z3==floor(z3), if(z3%9==0,
        print("(*)    ("x3","y3","z3/9")\n");
       ));
       y= (1+x^2+x^3)/y;
       print("("x","y","z")    ",factorint(y),"\n");
       x0= x; y0= y;
       x1= (y0^2+1)/x0; y1= y0;
       x2= x1; y2= (x1^3+x1+1)/y1;
       x3= (y2^2+1)/x2; y3= y2;
       z3= -(1+x3^2+x3^3+y3^2)/(x3*y3);
       if(z3==floor(z3), if(z3%9==0,
        print("(*)    ("x3","y3","z3/9")\n");
       ))
      )
     ))
    )
   )
  )
 )
Было бы здорово, если кто-то сможет объяснить магию этого 
метода.
Происходит примерно следующее. Начать нужно с более простого уравнения 

, для которого сравнительно легко находить малые целые решения 

, которые будут стартовыми в итерациях алгоритма. Легко видеть, что в дивизорах 

 два целых значения 

, которые участвуют в решениях. А в дивизорах 

 лишь одно целое значение 

, два других комплексные/иррациональные сопряженные, т.к. по 

 уравнение кубическое. Волшебство случается, когда метод прыгает на "соседнюю" поверхность 

, делаются три шага 

. В результате третья целая точка 

 почему-то оказывается на исходной поверхности. Если при этом 

, то решение найдено.