|
%I #40 Sep 14 2021 10:23:41
%S 1,7,25,55,115,181,307,439,637,841,1171,1447,1915,2329,2881,3433,4249,
%T 4879,5905,6745,7861,8911,10429,11557,13297,14773,16663,18355,20791,
%U 22495,25285,27541,30361,32905,36289,38845,42841,46027,49987,53395
%N Number of ordered triples (a, b, c) with gcd(a, b, c) = 1 and 1 <= {a, b, c} <= n.
%H Charles R Greathouse IV, <a href="/A071778/b071778.txt">Table of n, a(n) for n = 1..10000</a>
%H IBM Ponder This, <a href="https://www.research.ibm.com/haifa/ponderthis/challenges/June2002.html">Coin-weighing problem</a>, Jun 01 2002
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/GreatestCommonDivisor.html">Greatest Common Divisor</a>
%F a(n) = Sum_{k=1..n} mu(k)*floor(n/k)^3. - _Benoit Cloitre_, May 11 2003
%F a(n) = n^3 - Sum_{j=2..n} a(floor(n/j)). - _Vladeta Jovovic_, Nov 30 2004
%F G.f.: (1/(1 - x)) * Sum_{k >= 1} mu(k) * x^k * (1 + 4*x^k + x^(2*k))/(1 - x^k)^3. - _Seiichi Manyama_, May 22 2021
%F a(n) ~ n^3/zeta(3). - _Vaclav Kotesovec_, Sep 14 2021
%p f:=proc(n) local i,j,k,t1,t2,t3; t1:=0; for i from 1 to n do for j from 1 to n do t2:=gcd(i,j); for k from 1 to n do t3:=gcd(t2,k); if t3 = 1 then t1:=t1+1; fi; od: od: od: t1; end;
%t a[n_] := Sum[MoebiusMu[k]*Quotient[n, k]^3, {k, 1, n}]; Array[a, 40] (* _Jean-François Alcover_, Apr 14 2014, after _Benoit Cloitre_ *)
%o (Java) public class Triples { public static void main(String[] argv) { int i, j, k, a, m, n, d; boolean cf; try {a = Integer.parseInt(argv[0]);} catch (Exception e) {a = 10;}
%o for (m = 1; m <= a; m++) { n = 0; for (i = 1; i <= m; i++) for (j = 1; j <= m; j++) for (k = 1; k <= m; k++) { cf = false; for (d = 2; d <= m; d++) cf = cf || ((i % d == 0) && (j % d == 0) && (k % d == 0)); if (!cf) n++; } System.out.println(m + ": " + n); } } }
%o (PARI) a(n)=sum(k=1,n,moebius(k)*(n\k)^3)
%o (PARI) a(n)=my(s); forsquarefree(k=1,n, s+=moebius(k)*(n\k[1])^3); s \\ _Charles R Greathouse IV_, Jan 08 2018
%o (PARI) my(N=66, x='x+O('x^N)); Vec(sum(k=1, N, moebius(k)*x^k*(1+4*x^k+x^(2*k))/(1-x^k)^3)/(1-x)) \\ _Seiichi Manyama_, May 22 2021
%o (Python)
%o from functools import lru_cache
%o @lru_cache(maxsize=None)
%o def A071778(n):
%o if n == 0:
%o return 0
%o c, j = 1, 2
%o k1 = n//j
%o while k1 > 1:
%o j2 = n//k1 + 1
%o c += (j2-j)*A071778(k1)
%o j, k1 = j2, n//j2
%o return n*(n**2-1)-c+j # _Chai Wah Wu_, Mar 29 2021
%Y Cf. A018805 (ordered pairs), A082540, A082544, A343978, A344522.
%K nonn
%O 1,2
%A Michael Malak (mmalak(AT)alum.mit.edu), Jun 04 2002
|
