OFFSET
1,1
COMMENTS
No more terms up to 10^23. - Charles R Greathouse IV, Jul 27 2009
LINKS
P. De Geest, Index to related sequences
Hisanori Mishima, Sporadic tridigital solutions
PROG
(PARI) /* helper function: */
admissibleMod(M=10^5, t=[3, 9], debug=0)={ my(p=1, v); while(M > p *= 10,
t = concat([t, t + t[1]*v=vector(#t, i, p), t + t[2]*v]));
debug && print("t="t); t=Set(t); v=[];
for(k=1, M, setsearch(t, k^2 % M) && v=concat(v, k)); concat(v, M+v[1])}
/* optional arguments: Nmax = upper search limit, N = start / lower limit,
addMod = step/chunk size, debug: 0=silent, 1=verbose */
A058433(Nmax=1e10, N=1, addMod=10^5, debug=1)={ my(a=[], d=1,
addNext = admissibleMod(addMod=10^logint(addMod\/1, 10), [3, 9]),
add = vector(addMod, i, i-1 > addNext[d] && d++; addNext[d]-i+1),
pow10 = [10^k | k<-[0..logint((Nmax \/= 1)^2, 10)]],
nextOK = [if(n, n*pow10) | n<-[0, 2, 1, 0, 5, 4, 3, 2, 1, 0]]); N \/= 1;
while( Nmax >= N, my(N2 = N^2, numDigits = logint(N2, 10)+1,
place = nextOK[1 + d = N2 \ pow10[numDigits]]);
if( place, N = max(sqrtint(place[numDigits] + d*pow10[numDigits]), N+1);
next); place = 1;
my(Nnext = min(sqrtint((d+1)*pow10[numDigits]), Nmax));
debug && print("checking from "N" to "Nnext": <= ",
1+max(0, Nnext-N)*(#addNext-1)\ addMod, " candidates.");
while( Nnext >= N += add[1 + N % addMod],
my(dr = divrem( N2 = N^2, pow10[place = numDigits] ));
while( place-- && !d=nextOK[1+ (dr = divrem(dr[2], pow10[place]))[1]], );
place || break; N = sqrtint(N2 - dr[2] + d[place]) + 1;
); if( !place, debug && print(N "^2 = ", N^2); a=concat(a, N));
N = Nnext*3\2+1); a} \\ M. F. Hasler, May 14 2007
CROSSREFS
KEYWORD
nonn,base,hard,more,bref
AUTHOR
Patrick De Geest, Nov 15 2000
STATUS
approved