cp's OEIS Frontend

A380175 Greedy sums of distinct squares.

%I A380175 #37 Feb 15 2025 02:06:05
%S A380175 0,1,4,5,9,10,13,14,16,17,20,21,25,26,29,30,34,35,36,37,40,41,45,46,
%T A380175 49,50,53,54,58,59,62,63,64,65,68,69,73,74,77,78,80,81,82,85,86,90,91,
%U A380175 94,95,97,98,100,101,104,105,109,110,113,114,116,117,120,121,122,125,126,130
%N A380175 Greedy sums of distinct squares.
%C A380175 Let n be a positive integer. Suppose that s1^2 is the largest square not exceeding n, that s2^2 is the largest square not exceeding n-s1^2, and so on, so that n=s1^2+s2^2+...+sr^2 for some r. Clearly the si are weakly decreasing, but if they are strictly decreasing, s1>s2>...>sr, then we say that n is a greedy sum of distinct squares.
%C A380175 The set of integers for which the summands are distinct does not have a natural density, but the counting function oscillates in a predictable way (see Montgomery link).
%H A380175 Robert Israel, <a href="/A380175/b380175.txt">Table of n, a(n) for n = 1..10000</a>
%H A380175 Hugh Montgomery and Ulrike Vorhauer, <a href="https://doi.org/10.1090/S0025-5718-03-01513-8">Greedy sums of distinct squares</a>, Mathematics of Computation 73.245 (2004): 493-513.
%e A380175 21 is a term since greedily taking squares is 21 = 4^2 + 2^2 + 1^2 and all are distinct.
%e A380175 38 is not a term since greedily 38 = 6^2 + 1^2 + 1^2 and 1^2 has repeated (it can be distinct squares 38 = 5^2 + 3^2 + 2^2 but that's not the greedy sum).
%p A380175 filter:= proc(n) local s,x;
%p A380175   s:= n; x:= floor(sqrt(s));
%p A380175   do
%p A380175     s:= s - x^2;
%p A380175     if s >= x^2 then return false fi;
%p A380175     if s = 0 then return true fi;
%p A380175     x:= floor(sqrt(s))
%p A380175   od;
%p A380175 end proc:
%p A380175 filter(0):= true:
%p A380175 select(filter, [$0..1000]); # _Robert Israel_, Feb 14 2025
%t A380175 Select[Range[0,200], DuplicateFreeQ[#] &@ Differences[NestWhileList[# - Floor[Sqrt[#]]^2 &, #, # > 0 &]] &]
%o A380175 (Python)
%o A380175 from math import isqrt
%o A380175 def ok(n):
%o A380175     rprev = n+1
%o A380175     while 0 < (r:=isqrt(n)) < rprev: n, rprev = n-r**2, r
%o A380175     return r == 0
%o A380175 print([k for k in range(131) if ok(k)]) # _Michael S. Branicky_, Jan 15 2025
%Y A380175 Cf. A003995 (if not taken greedily), A380177.
%K A380175 nonn
%O A380175 1,3
%A A380175 _Mike Sheppard_, Jan 14 2025