1,2

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)

Considering representations where min(w,x)=0 shows that a(n) >= 2*A066839(n) - A038548(n).

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)

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.

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

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

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 *)

(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

Cf. A000203, A071561, A368276 (monotone), A368341 (fixed points), A368457, A368458 (semiprimes), A368580 (degenerated).

Sequence in context: A308160 A151729 A088652 * A266981 A045893 A071939

Adjacent sequences: A368204 A368205 A368206 * A368208 A368209 A368210

Don Knuth, Dec 16 2023

approved