A368207 Bacher numbers: number of nonnegative representations of n = w*x+y*z with max(w,x) < min(y,z).

1, 2, 2, 5, 3, 8, 4, 8, 9, 9, 6, 18, 7, 12, 14, 19, 9, 20, 10, 27, 16, 18, 12, 34, 20, 21, 20, 30, 15, 44, 16, 32, 24, 27, 30, 51, 19, 30, 28, 49, 21, 58, 22, 42, 41, 36, 24, 70, 35, 47, 36, 49, 27, 66, 36, 72, 40, 45, 30, 88, 31, 48, 62, 71, 42, 74, 34, 63, 48

Offset: 1

Author

Don Knuth, Dec 16 2023

Keywords

• nonn

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)

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)

Examples

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

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.

Crossrefs

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

Programs

CWEB

@ See Knuth link.

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

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

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

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