%I #38 Sep 24 2024 02:11:56
%S 1,1,0,2,2,0,0,2,3,3,0,0,3,0,0,4,4,3,0,4,0,0,0,0,5,5,0,0,5,0,0,4,0,5,
%T 0,6,6,0,0,6,5,0,0,0,6,0,0,0,7,7,0,6,7,0,0,0,0,7,0,0,6,0,0,8,8,0,0,8,
%U 0,0,0,6,8,7,0,0,0,0,0,8,9,9,0,0,9,0,0,0,8,9,0,0,0,0,0,0,9,7,0,10
%N Largest integer m such that n-m^2 is a square, or 0 if no such m exists.
%C The sequence could be extended to a(0) = 0.
%C We could have defined a(n) = -1 instead of 0 if n is not sum of two squares, and then include unambiguously a(0) = 0. At present, a(n) = 0 <=> A000161(n) = 0.
%H Alois P. Heinz, <a href="/A133388/b133388.txt">Table of n, a(n) for n = 1..20000</a>
%F a(n) = max( sup { max(a,b) | a^2+b^2 = n ; a,b in Z }, 0 )
%F a(A022544(j))=0, j>0. - _R. J. Mathar_, Jun 17 2009
%F a(n^2) = a(n^2 + 1) = n, for all n. Conversely, whenever a(n) = a(n+1), then n = k^2. - _M. F. Hasler_, Sep 02 2018
%e a(3) = 0 since 3 cannot be written as sum of 2 perfect squares;
%e a(5) = 2 since 5 = 2^2 + 1^2.
%p a:= proc(n) local t, d;
%p for t from 0 do d:= n-t^2;
%p if d<0 then break elif issqr(d) then return isqrt(d) fi
%p od; 0
%p end:
%p seq(a(n), n=1..100); # _Alois P. Heinz_, May 14 2015
%t a[n_] := Module[{m, d, s}, For[m = 0, True, m++, d = n - m^2; If[d < 0, Break[], s = Sqrt[d]; If[IntegerQ[s], Return[s]]]]; 0];
%t Table[a[n], {n, 1, 100}] (* _Jean-François Alcover_, May 18 2018, after _Alois P. Heinz_ *)
%o (PARI) sum2sqr(n)={ if(n>1, my(L=List(), f, p=1); for(i=1, matsize(f=factor(n))[1], if(f[i,1]%4==1, listput(L, [qfbsolve(Qfb(1,0,1), f[i,1])*[1, I]~, f[i,2]] ),/*elseif*/ f[i,1]==2, p = (1+I)^f[i,2],/*elseif*/ bittest(f[i,2],0), return([]),/*else*/ p *= f[i,1]^(f[i,2]\2))); L=apply(s->vector(s[2]+1, j, s[1]^(s[2]+1-j)*conj(s[1])^(j-1)), L); my(S=List()); forvec(T=vector(#L, i, [1,#L[i]]), listput(S, prod( j=1, #T, L[j][T[j]] ))); Set(apply(f->vecsort(abs([real(f), imag(f)])), Set(S)*p)), if(n<0, [], [[0, n]]))} \\ updated by _M. F. Hasler_, May 12 2018. (If PARI version 2.12.x returns an error, append [1] to qfbsolve(...) above. - _M. F. Hasler_, Dec 12 2019)
%o apply( A133388=n->if(n=sum2sqr(n),vecmax(Mat(n~))), [1..50]) \\ This sequence: maximum
%o (Python)
%o from sympy.solvers.diophantine.diophantine import diop_DN
%o def A133388(n): return max((a for a, b in diop_DN(-1,n)),default=0) # _Chai Wah Wu_, Sep 08 2022
%Y Cf. A000161, A001481.
%K nonn
%O 1,4
%A _M. F. Hasler_, Nov 23 2007