|
|
PROG
| (PARI) admissibleMod(LIM=10000)={ local( t=[4], tt=1 ); while( LIM > tt*=10, t=concat([t, t+vector(#t, i, 3)*tt, t+vector(#t, i, 4)*tt])); /*print("t="t); */ t=Set(t); tt=[]; for(i=1, LIM, if( setsearch(t, i^2%LIM), tt=concat(tt, i))); concat(tt, LIM+tt[1])} A058429(Nmax=1e19, N=2, addMod=100000, debug=1)={local( a=[], Nnext, N2, numDigits, place, addNext=admissibleMod(addMod=round(addMod)), d=1, add=vector(addMod, i, if(i-1>addNext[d], d++); addNext[d]-i+1), nextOK=[0, 2, 1, 0, 0, 5, 4, 3, 2, 1], pow10 = vector( d=#Str((Nmax=round(Nmax))^2), i, 10^(i-1)) ); nextOK = vector( #nextOK, i, if( nextOK[i], nextOK[i]*pow10)); N=round(N); while( Nmax >= N, numDigits = #Str(N2=N^2); if( 3 > N2 \ pow10[numDigits], N = sqrtint( 3*pow10[numDigits]+3 )+1); Nnext = min( Nmax, sqrtint(pow10[numDigits]*10)\3*2); if( debug, print( "checking from "N" to "Nnext": <= ", 1+max(0, Nnext-N)*(#addNext-1) \ addMod, " candidates.")); N += add[1+ N%addMod]; place=1; while( Nnext >= N, dr = divrem( N2=N^2, pow10[ place=numDigits ] ); while( place-- & !d=nextOK[1+ (dr = divrem( dr[2], pow10[ place ] ))[1]], ); if( !place, break); N = sqrtint( N2 - dr[2] + d[ place ])+1; N+=add[1+N%addMod]; ); if( !place, if( debug, print( N, "^2 = ", N^2 )); a=concat(a, N); N=Nnext+1 ); N = floor(Nnext*5/2); ); a } - M. F. Hasler (Maximilian.Hasler(AT)gmail.com), May 14 2007
|
