Squfof function in C or C++

[Rosario19]

i would like to write the Squfof function in C or C++, i see in

http://homes.cerias.purdue.edu/~ssw/squfof.pdf

http://guan.cse.nsysu.edu.tw/note/squfof.pdf

http://pari.math.u-bordeaux.fr/dochtml/scripts/squfof.gp.html

all these paper propose something it seems different in its
implementation, at my first seen

than i think one could use one 64 bit type, signed
for build the "form" f=(a,b,c) type [with a,b,c i64] (or struct) and
doing all

for sqrt() of one 64 bit signed integer, can i use the
long double sqrtl() ?

something as

#define sqrti64(x) (i64)sqrtl((long double)x)

i64 a,b;
...
b=sqrti64(a);

[where i64 is a define of one interger type of 64 bit]i think yes...
what would be the link it is better to follow?

thank you

[Rosario19]

this is the first try but seems something goes wrong...

void print(i64& n)
{char p[128];
*p=0;
_i64toa(n, p, 10);
P("%s", p);
}

#include <math.h>

#define sqrti64(x) ((i64) sqrtl((long double) x) )
#define absi64(x) (x>=0?x:-x)


class form{
i64 a;
i64 b;
i64 c;
i64 dD;

i32 IsReduced(void){i64 t;t=dD-2*absi64(a);R absi64(t)<b && b<dD;}

form(){a=b=c=dD=0;}
form(i64 fa, i64 fb, i64 fc){a=fa; b=fb; c=fc; dD=b*b-4*a*c;}

void assign(i64 fa, i64 fb, i64 fc){a=fa; b=fb; c=fc;
dD=b*b-4*a*c;}

friend ostream& operator<<(ostream& o,const form& x)
{o<<"("<<(long double)x.a<<", "<<(long double)x.b<<", "<<(long
double)x.c<<", "<<(long double)x.dD<<")";
R o;
}

friend i32 fr(i64* res, i64& xb, i64& xa, i64& xdD)
{i64 a2,absa,r,k;

if(res==0) R -1;
*res=0;
if(xa ==0) R -1;
r=xb; a2=2*xa; absa=absi64(xa);
if(absa==xdD)
{P("errore0"); R -1;}
if(absa>xdD)
{if(xa>0) k=((-xa-xb)/a2)+1;
else k=(( xa-xb)/a2)-1;
r=xb+k*a2;
if(r> absa)
r-=absi64(a2);
if(r> absa)
{ P("errore1 r="); print(r); P(" |a|="); print(absa);
P("\n");
P("k="); print(k); P(" b="); print(xb); P(" a2=");
print(a2); P("\n");
R -1;
}
if(r<=-absa)
{P("errore2");R -1;}
}
else
{if(xa>0) k=((xdD-a2-xb)/a2)+1;
else k=((xdD+a2-xb)/a2)-1;
r=xb+k*a2;
if(r> xdD)
r-=absi64(a2);
if(r> xdD)
{P("errore3"); R -1;}
if(r<=xdD-2*absa)
{P("errore4"); R -1;}
}
*res=r;
R 0;
}


friend i32 rhoF(form* res, form& f)
{i64 r, r2, fc4;
i32 v;

if(res== 0)R -1;
if(f.c== 0)R -1;
res->a=res->b=res->c=res->dD=0;
r2=-f.b;
v=fr(&r, r2, f.c, f.dD);
if(v ==-1)R -1;
res->a =f.c; res->b=r;
r2=r*r; fc4=4*f.c; res->c =(r2-f.dD)/fc4;
res->dD=r2-fc4*res->c;
R 0;
}


// suppone che i numeri a esso passato nn siano primi
friend i32 Squfof(i64* res, i64& n)
{i64 xdD,xb,dd,t,x, m, bound, l, gx, b1, ve[512];
i32 r;
u32 i;
form f,g, h;
Mset<i64> buf;

if(n<=1||res==0)R -1;
*res=0; dd=sqrti64(n);
if(dd*dd==n){*res=dd; R 0;}

if(n%4==1){xdD= n; m=1; xb=(dd-1)/2; xb=2*xb+1;}
else {xdD=4*n; m=2;dd=sqrti64(xdD); xb= dd /2; xb=2*xb ;}
f.a=1; f.b=xb; f.c=(xb*xb-xdD)/4; l=sqrti64(dd); bound=4*l;
f.dD=xb*xb-4*f.c;
t=absi64( f.c );
if(t==0)R -1;
gx=gcdi64(t,m); gx=t/gx;
buf.assegnaMem(ve, 512);

if(gx<=l)
{r=buf.add1(gx);
if(r==-1)R -1;
}
bound=158;
F(i=2;i<bound;++i)
{
r=rhoF(&g,f);
if(r==-1)
R -1;
f=g;
ooo<<"i="<<i<<" f="<<f<<"\n";
if(i%2==0) continue;
if(f.c>=0)
{x=sqrti64(f.c);
if(x<=1) R -1;
if(x*x==abs(f.c) && !buf.is(x))
break;
}
if(f.c==0)
R -1;
gx=absi64(f.c); x=gcdi64(gx,m); gx=gx/x;
if(gx<=l){r=buf.add1(gx);
if(r==-1)
R -1;
}
}
if(i>=bound)
R -1;
h.a=x*f.a; h.b=-f.b; h.c=-f.c; h.dD=h.b*h.b-4*h.a*h.c;
r=rhoF(&g, h);
if(r==-1)
R -1;
la: b1=g.b; r=rhoF(&h, g);
if(r==-1)
R -1;
g=h;if(b1!=g.b)G la;
if(g.c%2==0) g.c=g.c/2;
*res=absi64(g.c); R 0;
}

friend i32 squfofM(i64* res, i64 n)
{i64 xdD,xb,dd,t,x, m, bound, l, gx, b1, ve[512];
i32 r;
u32 j, i;
form f,g, h;
Mset<i64> buf;

if(n<=1||res==0)R -1;
*res=0; dd=sqrti64(n);
if(dd*dd==n){*res=dd; R 0;}
if(n%4==1){xdD= n; m=1; xb=(dd-1)/2; xb=2*xb+1;}
else {xdD=4*n; m=2;dd=sqrti64(xdD); xb= dd /2; xb=2*xb ;}
l=sqrti64(dd);
f.a=1; f.b=xb; f.c=(xb*xb-xdD)/4; f.dD=xb*xb-4*f.c;
buf.assegnaMem(ve, 512);
F(i=1;;++i)
{r=rhoF(&g,f);if(r==-1)R -1;f=g;
if(i%2==0&&f.a>=0)
{x=sqrti64(f.a);
if(x*x==f.a)
{if(!buf.is(x)) break; // nn trovato
}
}
x=absi64(f.a);
if(x<=l)
buf.add1(x);
if(i>300)R -1;
}
m=gcdi64(x, f.b); m=gcdi64(m, xdD);
if(m>1){*res=m; R 0;}
g.a=x; g.b=-f.b; g.c=x*f.c; g.dD=g.b*g.b-4*x*g.c;
r=rhoF(&h,g);if(r==-1)R -1;g=h;
la: b1=g.b; r=rhoF(&h, g);if(r==-1)R -1;g=h;if(b1!=g.b)G la;
if(g.a%2==0) g.a=g.a/2;
*res=absi64(g.a);
R 0;
}

};

#undef sqrti64
#undef absi64

i32 main96()
{i64 m1, m2, mm;
u32 i,j;
i32 r;
u8 *pc;

mm=633003781;
r=Squfof(&m1, mm);
ooo<<"r="<<r<<" res=";print(m1);P("\n");
R 0;
}

but print it give always 2 results

i=2 f=(-7125.000000, 633002591.000000, -14059378230867.000,
633003781.000000)
i=3 f=(-14059378230867.000, -633002591.000000, -7125.000000,
633003781.000000)
i=4 f=(-7125.000000, 633002591.000000, -14059378230867.000,
633003781.000000)
i=5 f=(-14059378230867.000, -633002591.000000, -7125.000000,
633003781.000000)
i=6 f=(-7125.000000, 633002591.000000, -14059378230867.000,
633003781.000000)
i=7 f=(-14059378230867.000, -633002591.000000, -7125.000000,
633003781.000000)
i=8 f=(-7125.000000, 633002591.000000, -14059378230867.000,
633003781.000000)
i=9 f=(-14059378230867.000, -633002591.000000, -7125.000000,
633003781.000000)

and continue always in this way...
someone know where is the problem?

[Rosario19]

one error was here because the above should be something as
"v=fr(&r, r2, f.c, sqrti64(f.dD))"

> if(v ==-1)R -1;
> res->a =f.c; res->b=r;
> r2=r*r; fc4=4*f.c; res->c =(r2-f.dD)/fc4;
> res->dD=r2-fc4*res->c;
> R 0;
>}
>

i64 gcdi64(i64 u, i64 v)
{i64 t;
if(u==0||v==0) R 0;
W(u>0){if(u<v){t=u; u=v; v=t;}
u = u%v;
}
R v;

}

void print(i64& n)
{char p[128];
*p=0;
_i64toa(n, p, 10);
P("%s", p);
}

#include <math.h>

#define sqrti64(x) ((i64) sqrtl((long double) x) )
#define absi64(x) (x>=0?x:-x)


class form{
i64 a;
i64 b;
i64 c;

i32 IsReduced(void){i64 t, dD; dD=b*b-4*a*c;t=dD-2*absi64(a);R

absi64(t)<b && b<dD;}

form(){a=b=c=0;}
form(i64& fa,i64& fb,i64& fc){a=fa; b=fb; c=fc;}
void assign(i64& fa,i64& fb,i64& fc){a=fa; b=fb; c=fc;}

friend ostream& operator<<(ostream& o,const form& x)

{i64 dD=x.b*x.b-4*x.a*x.c;

o<<"("<<(long double)x.a<<", "<<(long double)x.b<<", "<<(long

double)x.c<<", "<<(long double)dD<<")";

R o;
}

friend i32 fr(i64* res, i64& xb, i64& xa, i64& xdD)
{i64 a2,absa,r,k;

if(res==0) R -1;
*res=0;
if(xa ==0) R -1;
r=xb; a2=2*xa; absa=absi64(xa);

if(absa==xdD)R -1;

if(absa>xdD)
{if(xa>0) k=((-xa-xb)/a2)+1;
else k=(( xa-xb)/a2)-1;
r=xb+k*a2;
if(r> absa)
r-=absi64(a2);

//if(r> absa)R -1;
//if(r<=-absa)R -1;

}
else
{if(xa>0) k=((xdD-a2-xb)/a2)+1;
else k=((xdD+a2-xb)/a2)-1;
r=xb+k*a2;
if(r> xdD) r-=absi64(a2);

//if(r> xdD) R -1;
//if(r<=xdD-2*absa) R -1;
}
*res=r;
R 0;
}


friend i32 rhoF(form* res, form& f, i64& delta, i64& sqrtdelta)

{i64 r, r2, fc4;
i32 v;

if(res== 0)R -1;
if(f.c== 0)R -1;

res->a=res->b=res->c=0;
r2=-f.b;
v=fr(&r, r2, f.c, sqrtdelta);

if(v ==-1)R -1;
res->a =f.c; res->b=r;
r2=r*r; fc4=4*f.c;

res->c =(r2-delta)/fc4;

R 0;
}


// suppone che i numeri a esso passato nn siano primi

// Ritorna -1 per errore, altrimenti ritorna il numero dei loop
effettuati e in res,
// un fattore di n

friend i32 Squfof(i64* res, i64& n)
{i64 xdD,xb,dd,t,x, m, bound, l, gx, b1, ve[512];
i32 r;

u32 i,v;

form f,g, h;
Mset<i64> buf;

if(n<=1||res==0)R -1;
*res=0; dd=sqrti64(n);
if(dd*dd==n){*res=dd; R 0;}

if(n%4==1){xdD= n; m=1; xb=(dd-1)/2; xb=2*xb+1;}
else {xdD=4*n; m=2;dd=sqrti64(xdD); xb= dd /2; xb=2*xb ;}
f.a=1; f.b=xb; f.c=(xb*xb-xdD)/4; l=sqrti64(dd); bound=4*l;

t=absi64( f.c );
if(t==0)R -1;
gx=gcdi64(t,m); gx=t/gx;
buf.assegnaMem(ve, 512);

if(gx<=l)
{r=buf.add1(gx);
if(r==-1)R -1;
}

F(i=2;i<bound;++i)
{
r=rhoF(&g,f, xdD, dd);if(r==-1)R -1;
f=g;
t=absi64(f.c);
if(i%2==0)
{ooo<<"i="<<i<<" f="<<f<<"\n";
x=sqrti64(t);
if(x<=1) R -1;
if(x*x==t && !buf.is(x))
break;
}
gx=gcdi64(t,m); gx=t/gx;

if(gx<=l){r=buf.add1(gx);
if(r==-1)
R -1;
}
}
if(i>=bound)
R -1;

v=i;
if(f.c>0)
{h.a= x*f.a; h.b=-f.b; h.c=x;}
else {h.a=-x*f.a; h.b=-f.b; h.c=x;}
// h.dD=h.b*h.b-4*h.a*h.c= f.b^2-4*x^2*f.a
r=rhoF(&g, h, xdD, dd);if(r==-1) R -1;
i=0;
la: b1=g.b; r=rhoF(&h, g, xdD, dd);
if(i%2==0)ooo<<"i="<<i<<" g="<<g<<"\n";
if(r==-1) R -1;
g=h;++i;if(b1!=g.b)G la;
ooo<<"i="<<i<<" g="<<g<<"\n";
t=absi64(g.a); v+=i;
if(t<=1) R -1;
if(t%2==0) t=t/2;
x=n/t; m=x>t?t:x;
if(x*t!=n)R -1;
*res=m; R v;

}


};

#undef sqrti64
#undef absi64

i32 main96()
{i64 m1, m2, mm;
u32 i,j;
i32 r;
u8 *pc;

Mset<snum> bu;

mm=633003781;
r=Squfof(&m1, mm);
if(r==-1) ooo<<"Error\n";
else ooo<<"loop Conter="<<r<<" res=";print(m1);P("\n");
R 0;
}

i=2 f=(-7125.000000, 17591.000000, 11353.000000, 633003781.000000)
i=4 f=(-13363.000000, 21611.000000, 3105.000000, 633003781.000000)
i=6 f=(-12495.000000, 3131.000000, 12469.000000, 633003781.000000)
i=8 f=(-3157.000000, 22391.000000, 10425.000000, 633003781.000000)
i=10 f=(-6239.000000, 18125.000000, 12201.000000, 633003781.000000)
i=12 f=(-12163.000000, 18049.000000, 6315.000000, 633003781.000000)
i=14 f=(-9475.000000, 18059.000000, 8097.000000, 633003781.000000)
i=16 f=(-13205.000000, 12081.000000, 9221.000000, 633003781.000000)
i=18 f=(-483.000000, 24463.000000, 17891.000000, 633003781.000000)
i=20 f=(-7055.000000, 16901.000000, 12309.000000, 633003781.000000)
i=22 f=(-11647.000000, 15577.000000, 8379.000000, 633003781.000000)
i=24 f=(-9285.000000, 19201.000000, 7117.000000, 633003781.000000)
i=26 f=(-2835.000000, 21859.000000, 13685.000000, 633003781.000000)
i=28 f=(-11009.000000, 16507.000000, 8187.000000, 633003781.000000)
i=30 f=(-11275.000000, 6309.000000, 13153.000000, 633003781.000000)
i=32 f=(-4431.000000, 24313.000000, 2363.000000, 633003781.000000)
i=34 f=(-11261.000000, 22097.000000, 3213.000000, 633003781.000000)
i=36 f=(-8503.000000, 11127.000000, 14971.000000, 633003781.000000)
i=38 f=(-4659.000000, 18457.000000, 15687.000000, 633003781.000000)
i=40 f=(-7429.000000, 16799.000000, 11805.000000, 633003781.000000)
i=42 f=(-12423.000000, 18035.000000, 6193.000000, 633003781.000000)
i=44 f=(-10791.000000, 24041.000000, 1275.000000, 633003781.000000)
i=46 f=(-7295.000000, 19361.000000, 8847.000000, 633003781.000000)
i=48 f=(-10629.000000, 5231.000000, 14245.000000, 633003781.000000)
i=50 f=(-1615.000000, 21961.000000, 23331.000000, 633003781.000000)
i=52 f=(-245.000000, 24789.000000, 18887.000000, 633003781.000000)
i=54 f=(-6147.000000, 23897.000000, 2519.000000, 633003781.000000)
i=56 f=(-17181.000000, 12917.000000, 6783.000000, 633003781.000000)
i=58 f=(-15883.000000, 17551.000000, 5115.000000, 633003781.000000)
i=60 f=(-4247.000000, 19101.000000, 15785.000000, 633003781.000000)
i=62 f=(-7563.000000, 17783.000000, 10471.000000, 633003781.000000)
i=64 f=(-1245.000000, 23209.000000, 18945.000000, 633003781.000000)
i=66 f=(-5509.000000, 18373.000000, 13407.000000, 633003781.000000)
i=68 f=(-10475.000000, 12509.000000, 11373.000000, 633003781.000000)
i=70 f=(-11611.000000, 12985.000000, 9999.000000, 633003781.000000)
i=72 f=(-14597.000000, 22181.000000, 2415.000000, 633003781.000000)
i=74 f=(-18611.000000, 15933.000000, 5093.000000, 633003781.000000)
i=76 f=(-855.000000, 24779.000000, 5557.000000, 633003781.000000)
i=78 f=(-11059.000000, 24559.000000, 675.000000, 633003781.000000)
i=80 f=(-20383.000000, 16725.000000, 4333.000000, 633003781.000000)
i=82 f=(-17955.000000, 17971.000000, 4317.000000, 633003781.000000)
i=84 f=(-20767.000000, 24969.000000, 115.000000, 633003781.000000)
i=86 f=(-23805.000000, 22669.000000, 1251.000000, 633003781.000000)
i=88 f=(-2905.000000, 21611.000000, 14283.000000, 633003781.000000)
i=90 f=(-10233.000000, 13511.000000, 11005.000000, 633003781.000000)
i=92 f=(-12739.000000, 16979.000000, 6765.000000, 633003781.000000)
i=94 f=(-2791.000000, 21045.000000, 17029.000000, 633003781.000000)
i=96 f=(-6807.000000, 14215.000000, 15827.000000, 633003781.000000)
i=98 f=(-5195.000000, 24121.000000, 2463.000000, 633003781.000000)
i=100 f=(-105.000000, 25051.000000, 12979.000000, 633003781.000000)
i=102 f=(-12177.000000, 23447.000000, 1709.000000, 633003781.000000)
i=104 f=(-5471.000000, 19363.000000, 11793.000000, 633003781.000000)
i=106 f=(-13041.000000, 21859.000000, 2975.000000, 633003781.000000)
i=108 f=(-20279.000000, 20767.000000, 2487.000000, 633003781.000000)
i=110 f=(-5735.000000, 21881.000000, 6723.000000, 633003781.000000)
i=112 f=(-10871.000000, 25027.000000, 153.000000, 633003781.000000)
i=114 f=(-211.000000, 25061.000000, 5865.000000, 633003781.000000)
i=116 f=(-6615.000000, 17831.000000, 11907.000000, 633003781.000000)
i=118 f=(-12539.000000, 19095.000000, 5351.000000, 633003781.000000)
i=120 f=(-3303.000000, 22529.000000, 9495.000000, 633003781.000000)
i=122 f=(-10381.000000, 5311.000000, 14565.000000, 633003781.000000)
i=124 f=(-1127.000000, 23515.000000, 17757.000000, 633003781.000000)
i=126 f=(-6885.000000, 15541.000000, 14215.000000, 633003781.000000)
i=128 f=(-8211.000000, 19955.000000, 7149.000000, 633003781.000000)
i=130 f=(-3735.000000, 21881.000000, 10323.000000, 633003781.000000)
i=132 f=(-6205.000000, 17819.000000, 12711.000000, 633003781.000000)
i=134 f=(-11313.000000, 15023.000000, 9001.000000, 633003781.000000)
i=136 f=(-5355.000000, 21859.000000, 7245.000000, 633003781.000000)
i=138 f=(-5727.000000, 24205.000000, 2057.000000, 633003781.000000)
i=140 f=(-23085.000000, 25121.000000, 21.000000, 633003781.000000)
i=142 f=(-3933.000000, 22043.000000, 9351.000000, 633003781.000000)
i=144 f=(-10615.000000, 5869.000000, 14097.000000, 633003781.000000)
i=146 f=(-2387.000000, 20641.000000, 21675.000000, 633003781.000000)
i=148 f=(-1353.000000, 23293.000000, 16711.000000, 633003781.000000)
i=150 f=(-7935.000000, 21611.000000, 5229.000000, 633003781.000000)
i=152 f=(-10715.000000, 22639.000000, 2811.000000, 633003781.000000)
i=154 f=(-11923.000000, 1509.000000, 13225.000000, 633003781.000000)
i=156 f=(-207.000000, 25153.000000, 399.000000, 633003781.000000)
i=158 f=(-1215.000000, 23479.000000, 16819.000000, 633003781.000000)
i=160 f=(-7875.000000, 21341.000000, 5637.000000, 633003781.000000)
i=162 f=(-3047.000000, 24997.000000, 669.000000, 633003781.000000)
i=164 f=(-12075.000000, 23791.000000, 1387.000000, 633003781.000000)
i=166 f=(-15679.000000, 7991.000000, 9075.000000, 633003781.000000)
i=168 f=(-14595.000000, 19031.000000, 4639.000000, 633003781.000000)
i=170 f=(-16495.000000, 14909.000000, 6225.000000, 633003781.000000)
i=172 f=(-5197.000000, 19135.000000, 12837.000000, 633003781.000000)
i=174 f=(-11495.000000, 16451.000000, 7881.000000, 633003781.000000)
i=176 f=(-12873.000000, 10673.000000, 10081.000000, 633003781.000000)
i=178 f=(-13465.000000, 17441.000000, 6105.000000, 633003781.000000)
i=180 f=(-10843.000000, 24183.000000, 1111.000000, 633003781.000000)
i=182 f=(-5145.000000, 16459.000000, 17595.000000, 633003781.000000)
i=184 f=(-4009.000000, 21359.000000, 11025.000000, 633003781.000000)
i=186 f=(-2627.000000, 24545.000000, 2907.000000, 633003781.000000)
i=188 f=(-12939.000000, 3911.000000, 11935.000000, 633003781.000000)
i=190 f=(-4915.000000, 19361.000000, 13131.000000, 633003781.000000)
i=192 f=(-11145.000000, 15389.000000, 8887.000000, 633003781.000000)
i=194 f=(-6375.000000, 18091.000000, 11989.000000, 633003781.000000)
i=196 f=(-12477.000000, 19067.000000, 5399.000000, 633003781.000000)
i=198 f=(-2361.000000, 23095.000000, 10549.000000, 633003781.000000)
i=200 f=(-6355.000000, 19029.000000, 10657.000000, 633003781.000000)
i=202 f=(-1785.000000, 22811.000000, 15779.000000, 633003781.000000)
i=204 f=(-8817.000000, 8887.000000, 15709.000000, 633003781.000000)
i=206 f=(-1995.000000, 21359.000000, 22155.000000, 633003781.000000)
i=208 f=(-1199.000000, 25009.000000, 1575.000000, 633003781.000000)
i=210 f=(-21959.000000, 21677.000000, 1857.000000, 633003781.000000)
i=212 f=(-14675.000000, 6459.000000, 10073.000000, 633003781.000000)
i=214 f=(-11061.000000, 8435.000000, 12699.000000, 633003781.000000)
i=216 f=(-6797.000000, 23819.000000, 2415.000000, 633003781.000000)
i=218 f=(-3487.000000, 24337.000000, 2919.000000, 633003781.000000)
i=220 f=(-11367.000000, 23101.000000, 2185.000000, 633003781.000000)
i=222 f=(-1093.000000, 23123.000000, 22491.000000, 633003781.000000)
i=224 f=(-1725.000000, 22991.000000, 15133.000000, 633003781.000000)
i=226 f=(-9583.000000, 11891.000000, 12825.000000, 633003781.000000)
i=228 f=(-8649.000000, 20837.000000, 5747.000000, 633003781.000000)
i=230 f=(-45.000000, 25081.000000, 21929.000000, 633003781.000000)
i=232 f=(-3197.000000, 19587.000000, 19499.000000, 633003781.000000)
i=234 f=(-3285.000000, 20009.000000, 17705.000000, 633003781.000000)
i=236 f=(-5589.000000, 18133.000000, 13607.000000, 633003781.000000)
i=238 f=(-10115.000000, 11149.000000, 12573.000000, 633003781.000000)
i=240 f=(-8691.000000, 20767.000000, 5803.000000, 633003781.000000)
i=242 f=(-18765.000000, 23479.000000, 1089.000000, 633003781.000000)
i=0 g=(33.000000, 25129.000000, -11645.000000, 633003781.000000)
i=2 g=(3711.000000, 23081.000000, -6755.000000, 633003781.000000)
i=4 g=(12159.000000, 6869.000000, -12045.000000, 633003781.000000)
i=6 g=(6983.000000, 24677.000000, -861.000000, 633003781.000000)
i=8 g=(22915.000000, 22291.000000, -1485.000000, 633003781.000000)
i=10 g=(23155.000000, 24051.000000, -589.000000, 633003781.000000)
i=12 g=(19137.000000, 14027.000000, -5699.000000, 633003781.000000)
i=14 g=(9927.000000, 19541.000000, -6325.000000, 633003781.000000)
i=16 g=(11625.000000, 4841.000000, -13109.000000, 633003781.000000)
i=18 g=(3357.000000, 18907.000000, -20519.000000, 633003781.000000)
i=20 g=(1745.000000, 23239.000000, -13317.000000, 633003781.000000)
i=22 g=(11667.000000, 19939.000000, -5045.000000, 633003781.000000)
i=24 g=(10703.000000, 22391.000000, -3075.000000, 633003781.000000)
i=26 g=(16765.000000, 12871.000000, -6969.000000, 633003781.000000)
i=28 g=(14631.000000, 14257.000000, -7343.000000, 633003781.000000)
i=30 g=(13773.000000, 12431.000000, -8685.000000, 633003781.000000)
i=32 g=(3895.000000, 24431.000000, -2319.000000, 633003781.000000)
i=34 g=(16305.000000, 10661.000000, -7963.000000, 633003781.000000)
i=36 g=(5775.000000, 25009.000000, -327.000000, 633003781.000000)
i=38 g=(17707.000000, 10719.000000, -7315.000000, 633003781.000000)
i=40 g=(9885.000000, 20999.000000, -4857.000000, 633003781.000000)
i=42 g=(16169.000000, 14481.000000, -6545.000000, 633003781.000000)
i=44 g=(707.000000, 24701.000000, -8085.000000, 633003781.000000)
i=46 g=(2045.000000, 21181.000000, -22539.000000, 633003781.000000)
i=48 g=(687.000000, 24193.000000, -17359.000000, 633003781.000000)
i=50 g=(7521.000000, 19559.000000, -8325.000000, 633003781.000000)
i=52 g=(13339.000000, 12937.000000, -8727.000000, 633003781.000000)
i=54 g=(4305.000000, 21079.000000, -10957.000000, 633003781.000000)
i=56 g=(2635.000000, 24681.000000, -2263.000000, 633003781.000000)
i=58 g=(303.000000, 24587.000000, -23501.000000, 633003781.000000)
i=60 g=(1389.000000, 24811.000000, -3135.000000, 633003781.000000)
i=62 g=(21451.000000, 23823.000000, -763.000000, 633003781.000000)
i=64 g=(2475.000000, 24491.000000, -3353.000000, 633003781.000000)
i=66 g=(9615.000000, 16009.000000, -9795.000000, 633003781.000000)
i=68 g=(2453.000000, 20983.000000, -19641.000000, 633003781.000000)
i=70 g=(3795.000000, 19651.000000, -16261.000000, 633003781.000000)
i=72 g=(7185.000000, 15869.000000, -13263.000000, 633003781.000000)
i=74 g=(9791.000000, 8925.000000, -14129.000000, 633003781.000000)
i=76 g=(4587.000000, 17363.000000, -18069.000000, 633003781.000000)
i=78 g=(3881.000000, 20035.000000, -14919.000000, 633003781.000000)
i=80 g=(8997.000000, 8191.000000, -15725.000000, 633003781.000000)
i=82 g=(1463.000000, 23557.000000, -13341.000000, 633003781.000000)
i=84 g=(11679.000000, 20233.000000, -4787.000000, 633003781.000000)
i=86 g=(16019.000000, 13975.000000, -6831.000000, 633003781.000000)
i=88 g=(16645.000000, 19941.000000, -3535.000000, 633003781.000000)
i=90 g=(9031.000000, 13645.000000, -12369.000000, 633003781.000000)
i=92 g=(10307.000000, 9521.000000, -13155.000000, 633003781.000000)
i=94 g=(6673.000000, 23249.000000, -3465.000000, 633003781.000000)
i=96 g=(21427.000000, 24523.000000, -369.000000, 633003781.000000)
i=98 g=(8027.000000, 23239.000000, -2895.000000, 633003781.000000)
i=100 g=(8659.000000, 11555.000000, -14421.000000, 633003781.000000)
i=102 g=(5793.000000, 17471.000000, -14145.000000, 633003781.000000)
i=104 g=(9119.000000, 7419.000000, -15845.000000, 633003781.000000)
i=106 g=(693.000000, 24239.000000, -16405.000000, 633003781.000000)
i=108 g=(8527.000000, 8483.000000, -16449.000000, 633003781.000000)
i=110 g=(561.000000, 24953.000000, -4613.000000, 633003781.000000)
i=112 g=(10001.000000, 18827.000000, -6963.000000, 633003781.000000)
i=114 g=(3815.000000, 22829.000000, -7329.000000, 633003781.000000)
i=116 g=(6341.000000, 16901.000000, -13695.000000, 633003781.000000)
i=118 g=(9547.000000, 8605.000000, -14637.000000, 633003781.000000)
i=120 g=(3515.000000, 21511.000000, -12111.000000, 633003781.000000)
i=122 g=(12915.000000, 23119.000000, -1907.000000, 633003781.000000)
i=124 g=(15735.000000, 8821.000000, -8821.000000, 633003781.000000)
i=125 g=(-8821.000000, 8821.000000, 15735.000000, 633003781.000000)
loop Conter=367 res=8821

as in one paper say (changing something)

some other result

x=127951157 y=87366371 z=11178628252341247
x=127951157 y=87366371 z=11178628252341247
mm=11178628252341247
r=1343 res=87366371
All ok

x=120959249 y=114609553 z=13863085459105697
x=120959249 y=114609553 z=13863085459105697
mm=13863085459105697
r=3082 res=114609553
All ok

i see loop of 30000 cicles too, i make some other error?

than there are the case of fail for bound detection...

x=85568101 y=69332609 z=5932659689505509
x=85568101 y=69332609 z=5932659689505509
mm=5932659689505509
r=-1 res=0