login

Reminder: The OEIS is hiring a new managing editor, and the application deadline is January 26.

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

%I #135 Jan 06 2025 22:25:11

%S 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,

%T 30,15,44,16,32,24,27,30,51,19,30,28,49,21,58,22,42,41,36,24,70,35,47,

%U 36,49,27,66,36,72,40,45,30,88,31,48,62,71,42,74,34,63,48

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

%C When n = p is an odd prime, Bacher proved that a(p) = (p+1)/2.

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

%C 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.

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

%C From _Chai Wah Wu_, Dec 19 2023: (Start)

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

%C 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.

%C 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)

%H Don Knuth, <a href="/A368207/b368207.txt">Table of n, a(n) for n = 1..10000</a>

%H Roland Bacher, <a href="https://doi.org/10.1080/00029890.2023.2242034">A quixotic proof of Fermat's two squares theorem for prime numbers</a>, American Mathematical Monthly, Vol. 130, No. 9 (November 2023), 824-836; <a href="https://arxiv.org/abs/2210.07657">arXiv version</a>, arXiv:2210.07657 [math.NT], 2022.

%H Michael S. Branicky, <a href="/A368207/a368207_2.txt">Python translation of Knuth's CWEB program</a>

%H Pontus von Brömssen, <a href="/A368207/a368207.png">Plot of (n, a(n) - sigma(n)/2) for n = 1..100000</a>.

%H Pontus von Brömssen, <a href="/A368207/a368207_1.png">Plot of (n, b(n)) for n = 1..100000, where a(n) = sigma(n)/2 + b(n)*sqrt(n)</a>.

%H Don Knuth, <a href="/A368207/a368207.txt">CWEB program</a>.

%H Peter Luschny, <a href="https://github.com/PeterLuschny/Gists/blob/main/BacherNumbers.pdf">A Formula for the Bacher Numbers</a>.

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

%e 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 t[n_] := t[n] = Select[Divisors[n], #^2 <= n&];

%t 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]}];

%t Array[A368207, 100] (* _Paolo Xausa_, Jan 02 2024 *)

%o (CWEB) @ See Knuth link.

%o (Python) # See Branicky link for translation of Knuth's CWEB program.

%o (Python)

%o from math import isqrt

%o def A368207(n):

%o c, r = 0, isqrt(n)

%o for w in range(r+1):

%o for x in range(w,r+1):

%o wx = w*x

%o if wx>n:

%o break

%o for y in range(x+1,r+1):

%o for z in range(y,n+1):

%o yz = wx+y*z

%o if yz>n:

%o break

%o if yz==n:

%o m = 1

%o if w!=x:

%o m<<=1

%o if y!=z:

%o m<<=1

%o c+=m

%o return c # _Chai Wah Wu_, Dec 19 2023

%o (Python)

%o from sympy import divisors

%o # faster program

%o def A368207(n):

%o c = 0

%o for d2 in divisors(n):

%o if d2**2 > n:

%o break

%o c += (d2<<2)-2 if d2**2<n else (d2<<1)-1

%o for wx in range(1,(n>>1)+1):

%o for d1 in divisors(wx):

%o if d1**2 > wx:

%o break

%o for d2 in divisors(m:=n-wx):

%o if d2**2 > m:

%o break

%o if wx < d1*d2:

%o k = 1

%o if d1**2 != wx:

%o k <<=1

%o if d2**2 != m:

%o k <<=1

%o c+=k

%o return c # _Chai Wah Wu_, Dec 19 2023

%o (Julia)

%o function A368207(n)

%o t(n) = (d for d in divisors(n) if d * d <= n)

%o s(d) = d * d == n ? d * 2 - 1 : d * 4 - 2

%o c(y, w, wx) = max(1, 2*(Int(w*w < wx) + Int(y*y < n - wx)))

%o sum(sum(sum(c(y, w, wx) for y in t(n - wx) if wx < y * w; init=0)

%o for w in t(wx)) for wx in 1:div(n, 2); init=sum(s(d) for d in t(n)))

%o end

%o println([A368207(n) for n in 1:69]) # _Peter Luschny_, Dec 21 2023

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

%K nonn

%O 1,2

%A _Don Knuth_, Dec 16 2023