OFFSET
1,2
COMMENTS
When n = p is an odd prime, Bacher proved that a(p) = (p+1)/2.
It appears that also a(k*p) = sigma(k)*(p+1)/2, for all prime p > 2k, where sigma(k) is the sum of the divisors of k (A000203).
It appears furthermore that a(p^2) = (p^2 + 3*p)/2; a(p^3) = (p^3 + p^2 + p + 1)/2; a(p^4) = (p^4 + p^3 + 3p^2 + p)/2, for all prime p.
Conjectures: (1) a(n) >= sigma(n)/2, with equality if and only if n has no middle divisors, i.e., if and only if n is in A071561. (2) a(n)/sigma(n) converges to 1/2. - Pontus von Brömssen, Dec 18 2023
From Chai Wah Wu, Dec 19 2023: (Start)
It appears that a(p^5) = (p^5 + p^4 + p^3 + p^2 + p + 1)/2 for all prime p > 2 and a(p^6) = (p^6 + p^5 + p^4 + 3p^3 + p^2 + p)/2 for all prime p.
Conjecture: a(p^m) = sigma(p^m)/2 for odd m and all prime p > 2. a(p^m) = (sigma(p^m)-1)/2 + p^(m/2) for even m and all prime p. a(2^m) = sigma(2^m)/2 + 1/2 for m odd. (End)
LINKS
Don Knuth, Table of n, a(n) for n = 1..10000
Roland Bacher, A quixotic proof of Fermat's two squares theorem for prime numbers, American Mathematical Monthly, Vol. 130, No. 9 (November 2023), 824-836; arXiv version, arXiv:2210.07657 [math.NT], 2022.
Michael S. Branicky, Python translation of Knuth's CWEB program
Pontus von Brömssen, Plot of (n, a(n) - sigma(n)/2) for n = 1..100000.
Pontus von Brömssen, Plot of (n, b(n)) for n = 1..100000, where a(n) = sigma(n)/2 + b(n)*sqrt(n).
Don Knuth, CWEB program.
Peter Luschny, A Formula for the Bacher Numbers.
FORMULA
Let t(n) = {d|n and d*d <= n}, and s(d, n) = 2*d - 1 if d*d = n, otherwise 4*d - 2. Then a(n) = (Sum_{d in t(n)} s(d, n)) + (Sum_{k=1..floor(n/2)} Sum_{w in t(k)} Sum_{y in t(n-k) and k < y*w} max(1, 2*([w*w < k] + [y*y < n - k]))), where [] denote the Iverson brackets. (See the 'Julia' implementation below.) - Peter Luschny, Dec 21 2023
EXAMPLE
For n = 13, the a(13) = 7 solutions are (w,x,y,z) = (0,0,1,13), (0,0,13,1), (1,1,2,6), (1,1,3,4), (1,1,4,3), (1,1,6,2), (2,2,3,3).
MATHEMATICA
t[n_] := t[n] = Select[Divisors[n], #^2 <= n&];
A368207[n_] := Sum[(1 + Boole[d^2 < n])(2d - 1), {d, t[n]}] + Sum[If[wx < y*w, Max[1, 2(Boole[w^2 < wx] + Boole[y^2 < n-wx])], 0], {wx, Floor[n/2]}, {w, t[wx]}, {y, t[n - wx]}];
Array[A368207, 100] (* Paolo Xausa, Jan 02 2024 *)
PROG
(CWEB) @ See Knuth link.
(Python) # See Branicky link for translation of Knuth's CWEB program.
(Python)
from math import isqrt
def A368207(n):
c, r = 0, isqrt(n)
for w in range(r+1):
for x in range(w, r+1):
wx = w*x
if wx>n:
break
for y in range(x+1, r+1):
for z in range(y, n+1):
yz = wx+y*z
if yz>n:
break
if yz==n:
m = 1
if w!=x:
m<<=1
if y!=z:
m<<=1
c+=m
return c # Chai Wah Wu, Dec 19 2023
(Python)
from sympy import divisors
# faster program
def A368207(n):
c = 0
for d2 in divisors(n):
if d2**2 > n:
break
c += (d2<<2)-2 if d2**2<n else (d2<<1)-1
for wx in range(1, (n>>1)+1):
for d1 in divisors(wx):
if d1**2 > wx:
break
for d2 in divisors(m:=n-wx):
if d2**2 > m:
break
if wx < d1*d2:
k = 1
if d1**2 != wx:
k <<=1
if d2**2 != m:
k <<=1
c+=k
return c # Chai Wah Wu, Dec 19 2023
(Julia)
function A368207(n)
t(n) = (d for d in divisors(n) if d * d <= n)
s(d) = d * d == n ? d * 2 - 1 : d * 4 - 2
c(y, w, wx) = max(1, 2*(Int(w*w < wx) + Int(y*y < n - wx)))
sum(sum(sum(c(y, w, wx) for y in t(n - wx) if wx < y * w; init=0)
for w in t(wx)) for wx in 1:div(n, 2); init=sum(s(d) for d in t(n)))
end
println([A368207(n) for n in 1:69]) # Peter Luschny, Dec 21 2023
CROSSREFS
KEYWORD
nonn,changed
AUTHOR
Don Knuth, Dec 16 2023
STATUS
approved