A004431 - OEIS
5, 10, 13, 17, 20, 25, 26, 29, 34, 37, 40, 41, 45, 50, 52, 53, 58, 61, 65, 68, 73, 74, 80, 82, 85, 89, 90, 97, 100, 101, 104, 106, 109, 113, 116, 117, 122, 125, 130, 136, 137, 145, 146, 148, 149, 153, 157, 160, 164, 169, 170, 173, 178, 180, 181, 185, 193, 194, 197
COMMENTS
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
Reordering of A055096 by increasing values and without repetition. - Paul Curtz, Sep 08 2008
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
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
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.
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.
Numbers of the form j^2 + 2*j*k + 2*k^2 {j,k >= 1}. (End)
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
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
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.
The subsequence of primes is A002144. (End)
EXAMPLE
53 = 7^2 + 2^2, so 53 is in the sequence.
MAPLE
isA004431 := proc(n)
local a, b ;
for a from 2 do
if a^2>= n then
return false;
end if;
b := n -a^2 ;
if b < 1 then
return false ;
end if;
if issqr(b) then
if ( sqrt(b) <> a ) then
return true;
end if;
end if;
end do:
return false;
end proc:
option remember ;
local a;
if n = 1 then
5;
else
for a from procname(n-1)+1 do
if isA004431(a) then
return a;
end if;
end do:
end if;
MATHEMATICA
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 *)
Select[Range@ 200, Length[PowersRepresentations[#, 2, 2] /. {{0, _} -> Nothing, {a_, b_} /; a == b -> Nothing}] > 0 &] (* Michael De Vlieger, Mar 24 2016 *)
PROG
(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
(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
for(n=1, 1e3, if(is(n), print1(n, ", "))) \\ Altug Alkan, Dec 06 2015
(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
(Haskell)
import Data.List (findIndices)
a004431 n = a004431_list !! (n-1)
a004431_list = findIndices (> 1) a063725_list
(Python)
def aupto(limit):
s = [i*i for i in range(1, int(limit**.5)+2) if i*i < limit]
s2 = set(a+b for i, a in enumerate(s) for b in s[i+1:] if a+b <= limit)
return sorted(s2)