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”).

A133388
Largest integer m such that n-m^2 is a square, or 0 if no such m exists.
18
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, 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, 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
OFFSET
1,4
COMMENTS
The sequence could be extended to a(0) = 0.
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.
LINKS
FORMULA
a(n) = max( sup { max(a,b) | a^2+b^2 = n ; a,b in Z }, 0 )
a(A022544(j))=0, j>0. - R. J. Mathar, Jun 17 2009
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
EXAMPLE
a(3) = 0 since 3 cannot be written as sum of 2 perfect squares;
a(5) = 2 since 5 = 2^2 + 1^2.
MAPLE
a:= proc(n) local t, d;
for t from 0 do d:= n-t^2;
if d<0 then break elif issqr(d) then return isqrt(d) fi
od; 0
end:
seq(a(n), n=1..100); # Alois P. Heinz, May 14 2015
MATHEMATICA
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];
Table[a[n], {n, 1, 100}] (* Jean-François Alcover, May 18 2018, after Alois P. Heinz *)
PROG
(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)
apply( A133388=n->if(n=sum2sqr(n), vecmax(Mat(n~))), [1..50]) \\ This sequence: maximum
(Python)
from sympy.solvers.diophantine.diophantine import diop_DN
def A133388(n): return max((a for a, b in diop_DN(-1, n)), default=0) # Chai Wah Wu, Sep 08 2022
CROSSREFS
Sequence in context: A031124 A063695 A081417 * A354643 A282516 A158092
KEYWORD
nonn
AUTHOR
M. F. Hasler, Nov 23 2007
STATUS
approved