The OEIS is supported by the many generous donors to the OEIS Foundation.
%I #130 Apr 14 2022 09:42:15
%S 5,10,13,17,20,25,26,29,34,37,40,41,45,50,52,53,58,61,65,68,73,74,80,
%T 82,85,89,90,97,100,101,104,106,109,113,116,117,122,125,130,136,137,
%U 145,146,148,149,153,157,160,164,169,170,173,178,180,181,185,193,194,197
%N Numbers that are the sum of 2 distinct nonzero squares.
%C Numbers whose prime factorization includes at least one prime congruent to 1 mod 4 and any prime factor congruent to 3 mod 4 has even multiplicity. - _Franklin T. Adams-Watters_, May 03 2006
%C Reordering of A055096 by increasing values and without repetition. - _Paul Curtz_, Sep 08 2008
%C A063725(a(n)) > 1. - _Reinhard Zumkeller_, Aug 16 2011
%C The square of these numbers is also the sum of two nonzero squares, so this sequence is a subsequence of A009003. - _Jean-Christophe Hervé_, Nov 10 2013
%C Closed under multiplication. Primitive elements are those with exactly one prime factor congruent to 1 mod 4 with multiplicity one (A230779). - _Jean-Christophe Hervé_, Nov 10 2013
%C From _Bob Selcoe_, Mar 23 2016: (Start)
%C Numbers c such that there is d < c, d >= 1 where c + d and c - d are square. For example, 53 + 28 = 81, 53 - 28 = 25.
%C Given a prime p == 1 mod 4, a term appears if and only if it is of the form p^i, p*2^j or p*k^2 {i,j,k >= 1}, or a product of any combination of these forms. Therefore, the products of any terms to any powers also are terms. For example, p(1) = 5 and p(2) = 13 so term 45 appears because 5*3^2 = 45 and term 416 appears because 13*2^5 = 416; therefore 45 * 416 = 18720 appears, as does 45^3 * 416^11 = 18720^3 * 416^8.
%C Numbers of the form j^2 + 2*j*k + 2*k^2 {j,k >= 1}. (End)
%C Suppose we have a term t = x^2 + y^2. Then s^2*t = (s*x)^2 + (s*y)^2 is a term for any s > 0. Also 2*t = (y+x)^2 + (x-y)^2 is a term. It follows that q*s^2*t is a term for any s > 0 and q=1 or 2. Examples: 2*7^2*26 = 28^2 + 42^2; 6^2*17 = 6^2 + 24^2. - _Jerzy R Borysowicz_, Aug 11 2017
%C To find terms up to some upper bound u, we can search for x^2 + y^2 = t where x is odd and y is even. Then we add all numbers of the form 2^m * t <= u and then remove duplicates. - _David A. Corneth_, Oct 04 2017
%C From _Bernard Schott_, Apr 13 2022: (Start)
%C The 5th comment "Closed under multiplication" can be proved with Brahmagupta's identity: (a^2+b^2) * (c^2+d^2) = (ac + bd)^2 + (ad - bc)^2.
%C The subsequence of primes is A002144. (End)
%H T. D. Noe, <a href="/A004431/b004431.txt">Table of n, a(n) for n = 1..10000</a>.
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Brahmagupta's_identity">Brahmagupta's identity</a>.
%H <a href="/index/Su#ssq">Index entries for sequences related to sums of squares</a>.
%e 53 = 7^2 + 2^2, so 53 is in the sequence.
%p isA004431 := proc(n)
%p local a,b ;
%p for a from 2 do
%p if a^2>= n then
%p return false;
%p end if;
%p b := n -a^2 ;
%p if b < 1 then
%p return false ;
%p end if;
%p if issqr(b) then
%p if ( sqrt(b) <> a ) then
%p return true;
%p end if;
%p end if;
%p end do:
%p return false;
%p end proc:
%p A004431 := proc(n)
%p option remember ;
%p local a;
%p if n = 1 then
%p 5;
%p else
%p for a from procname(n-1)+1 do
%p if isA004431(a) then
%p return a;
%p end if;
%p end do:
%p end if;
%p end proc: # _R. J. Mathar_, Jan 29 2013
%t A004431 = {}; Do[a = 2 m * n; b = m^2 - n^2; c = m^2 + n^2; AppendTo[A004431, c], {m, 100}, {n, m - 1}]; Take[Union@A004431, 63] (* _Robert G. Wilson v_, May 02 2009 *)
%t Select[Range@ 200, Length[PowersRepresentations[#, 2, 2] /. {{0, _} -> Nothing, {a_, b_} /; a == b -> Nothing}] > 0 &] (* _Michael De Vlieger_, Mar 24 2016 *)
%o (PARI) select( isA004431(n)={n>1 && vecmin((n=factor(n)%4)[,1])==1 && ![f[1]>2 && f[2]%2 | f<-n~]}, [1..199]) \\ _M. F. Hasler_, Feb 06 2009, updated Nov 24 2019
%o (PARI) is(n)=if(n<5, return(0)); my(f=factor(n)%4); if(vecmin(f[, 1])>1, return(0)); for(i=1, #f[, 1], if(f[i, 1]==3 && f[i, 2]%2, return(0))); 1
%o for(n=1, 1e3, if(is(n), print1(n, ", "))) \\ _Altug Alkan_, Dec 06 2015
%o (PARI) upto(n) = {my(res = List(), s); forstep(i=1, sqrtint(n), 2, forstep(j = 2, sqrtint(n - i^2), 2, listput(res, i^2 + j^2))); s = #res; for(i = 1, s, t = res[i]; for(e = 1, logint(n \ res[i], 2), listput(res, t<<=1))); listsort(res, 1); res} \\ _David A. Corneth_, Oct 04 2017
%o (Haskell)
%o import Data.List (findIndices)
%o a004431 n = a004431_list !! (n-1)
%o a004431_list = findIndices (> 1) a063725_list
%o -- _Reinhard Zumkeller_, Aug 16 2011
%o (Python)
%o def aupto(limit):
%o s = [i*i for i in range(1, int(limit**.5)+2) if i*i < limit]
%o s2 = set(a+b for i, a in enumerate(s) for b in s[i+1:] if a+b <= limit)
%o return sorted(s2)
%o print(aupto(197)) # _Michael S. Branicky_, May 10 2021
%Y Complement of A004439.
%Y Cf. A000404, A002144, A007692, A009000, A009003, A009177, A024507, A081324, A081325, A118882, A230779.
%K nonn
%O 1,1
%A _N. J. A. Sloane_