%I #77 Jan 03 2023 02:10:08
%S 1,2,3,16,17,18,19,56,105,106,107,144,145,146,147,208,209,306,307,320,
%T 321,322,323,432,529,530,723,736,737,738,739,968,969,970,971,1176,
%U 1177,1178,1179,1288,1289,1290,1291,1304,1401,1402,1403,1608,1777,2018,2019,2032
%N Number of length four 1..n vectors that contain their geometric mean.
%C From _David A. Corneth_, Aug 27 2020: (Start)
%C Let (a, b, g, n) be a tuple where g is the geometric mean of a*b*g*n and a, b, g <= n. Then there are two cases: g < n and g = n.
%C If g = n then (a, b, g, n) = (n, n, n, n) adding 1 to the number of permutations.
%C If g < n then a*b = g^4 / (g*n) = g^3 / n. Furthermore let s be the squarefree part of n (Cf. A007947). Then s | g and so candidates for g are (s*i) where 1 <= i < floor(n/s), depending on whether g^3 / n = (s*i)^3 / n is an integer.
%C It follows that for suitable values of i, (a, b) are a pair of divisors of (s*i)^3 / n where a*b = (s*i)^3 / n and max(a, b) <= n. (End)
%C Bounds: n + 12*(floor(n/4) + 2*floor(n/8) + 4*floor(n/9) + 2*floor(n/12) + 2*floor(n/16) + 4*floor(n/18)) <= a(n) <= n^4 - 14*(n-1). - _Hywel Normington_, Jan 25 2021
%H David A. Corneth, <a href="/A337111/b337111.txt">Table of n, a(n) for n = 1..10000</a>
%H Hywel Normington, <a href="https://github.com/Horep/Number-of-vectors-that-contain-their-average/blob/master/A337111.py">Python code</a>, 2020.
%H Hywel Normington, <a href="https://github.com/Horep/Number-of-vectors-that-contain-their-average/blob/master/Julia_Edition/A337111.jl">Julia code</a>, 2023.
%F Empirical: if A000189(n) = 1 then a(n) = a(n-1) + 1.
%F From _David A. Corneth_, Aug 25 2020: (Start)
%F The above holds. That is: if x^3 == 0 (mod n) has only one solution then a(n) = a(n-1) + 1. Proof:
%F Let (a, b, c, n) be such a tuple. Let without loss of generality c be the geometric mean of the tuple. Then a*b*c*n = c^4 and as c is not 0 we have c^3 = a*b*n. So then c^3 == 0 (mod n). If c^3 == 0 (mod n) has only 1 solution then c = n. This gives the tuple (n, n, n, n) which has 1 permutation. So giving a(n) = a(n-1) + 1. (End)
%F a(n) - a(n-1) == 1 (mod 12). - _Hywel Normington_, Sep 28 2020
%e For n = 2 the a(2) = 2 solutions are: (1,1,1,1) and (2,2,2,2).
%e For n = 4 the a(4) = 16 solutions are:
%e (1, 1, 1, 1), (2, 2, 2, 2), (3, 3, 3, 3), (4, 4, 4, 4),
%e and the 12 permutations of (1, 2, 2, 4).
%e For n = 40, the a(40)-a(39) = 109 new solutions are:
%e (40,40,40,40),
%e the 24 permutations of (1, 10, 25, 40),
%e the 12 permutations of (5, 5, 10, 40),
%e the 12 permutations of (5, 20, 40, 40),
%e the 24 permutations of (8, 20, 25, 40),
%e the 12 permutations of (10, 20, 20, 40),
%e and the 24 permutations of (25, 27, 30, 40).
%o (PARI) first(n) = { my(res = vector(n)); s = 0; for(i = 1, n, s += b(i); res[i] = s; ); res }
%o b(n) = {my(resa = 1); my(s = factorback(factor(n)[, 1])); for(i = 1, n \ s - 1, s4 = (s*i)^3; if(s4 % n == 0, c = tuples((s*i)^3/n, s*i, n); for(i = 1, #c, resa+=qperms(c[i]) ) ) ); resa }
%o qperms(v) = {my(r=1,t); v = vecsort(v); for(i=1,#v-1,if(v[i]==v[i+1],t++,r*=binomial(i,t+1);t=0));r*=binomial(#v,t+1)}
%o tuples(n, s, u) = {my(res = List(), u4n, d, i); d = divisors(n); i = (#d + 1) \ 2; while(i > 0 && d[#d - i + 1] <= u, c = vecsort([d[i], d[#d - i + 1], s, u]); listput(res, c); i--); res} \\ _David A. Corneth_, Aug 28 2020
%Y Cf. A007947, A000189, A248435, A337110.
%K nonn,easy
%O 1,2
%A _Hywel Normington_, Aug 16 2020