login

Year-end appeal: Please make a donation to the OEIS Foundation to support ongoing development and maintenance of the OEIS. We are now in our 61st year, we have over 378,000 sequences, and we’ve reached 11,000 citations (which often say “discovered thanks to the OEIS”).

a(n) is the greatest prime q such that A367798(n)^2 is the sum of q and its reversal.
2

%I #10 Dec 30 2023 23:12:45

%S 2,83,81131,894500063,88990607813,8499228209501,8597793891803,

%T 800072140300001,859981720058603,899969843983163,82943190509220401,

%U 86999838571212401,88290616680100001,89991996902408171,83909667566050103,89690298128004023,89919974791600043,8069990701280128001,8299959944574088001

%N a(n) is the greatest prime q such that A367798(n)^2 is the sum of q and its reversal.

%C a(n) is the last term q of A367796 such that A056964(q) = A367798(n)^2.

%F A056964(a(n)) = A367798(n)^2.

%e a(4) = 894500063 because A367798(3) = 35419 and 35419^2 = 1254505561 = 894500063 + 360005498 and 894500063 is the greatest prime that works.

%p f:= proc(n) local y, c, d, dp, i, delta, m;

%p y:= convert(n^2, base, 10);

%p d:= nops(y);

%p if d::even then

%p if y[-1] <> 1 then return false fi;

%p dp:= d-1;

%p y:= y[1..-2];

%p c[dp]:= 1;

%p else

%p dp:= d;

%p c[dp]:= 0;

%p fi;

%p c[0]:= 0;

%p for i from 1 to floor(dp/2) do

%p delta:= y[i] - y[dp+1-i] - c[i-1] - 10*c[dp+1-i];

%p if delta = 0 then c[dp-i]:= 0; c[i]:= 0;

%p elif delta = -1 then c[dp-i]:= 1; c[i]:= 0;

%p elif delta = -10 then c[dp-i]:= 0 ; c[i]:= 1;

%p elif delta = -11 then c[dp-i]:= 1; c[i]:= 1;

%p else return false

%p fi;

%p if y[i] + 10*c[i] - c[i-1] < 0 or (i=1 and y[i]+10*c[i]-c[i-1]=1) then return false fi;

%p od;

%p m:= (dp+1)/2;

%p delta:= y[m] + 10*c[m] - c[m-1];

%p if not member(delta, [seq(i, i=0..18, 2)]) then return false fi;

%p [seq(y[i]+ 10*c[i]-c[i-1], i=1..m)]

%p end proc:

%p g:= proc(L) local T, d, t, p, x, i; uses combinat;

%p d:= nops(L);

%p T:= cartprod([select(t -> t[1]::odd, [seq([L[1]-x, x], x=min(L[1], 9)..max(1, L[1]-9),-1)]),

%p seq([seq([L[i]-x, x], x=min(9, L[i])..max(0, L[i]-9),-1)], i=2..d-1)]);

%p while not T[finished] do

%p t:= T[nextvalue]();

%p p:= add(t[i][1]*10^(i-1), i=1..d-1) + L[-1]/2 * 10^(d-1) +

%p add(t[i][2]*10^(2*d-i-1), i=1..d-1);

%p if isprime(p) then return p fi;

%p od;

%p -1

%p end proc:

%p p:= 2, 11: Q:= 83:

%p while p < 10^10 do

%p p:= nextprime(p);

%p d:= 1+ilog10(p^2);

%p if d::even and p^2 >= 2*10^(d-1) then p:= nextprime(floor(10^(d/2))); fi;

%p v:= f(p);

%p if v = false then next fi;

%p q:= g(v);

%p if q = -1 then next fi;

%p Q:= Q, q;

%p od:

%p Q;

%Y Cf. A056964, A367796, A367798, A367871.

%K nonn,base

%O 1,1

%A _Robert Israel_, Dec 04 2023