OFFSET
1,2
COMMENTS
Number of integer-sided triangles with at least two sides <= n and sides relatively prime. - Henry Bottomley, Sep 29 2006
LINKS
R. J. Mathar, Table of n, a(n) for n = 1..10000
FORMULA
Partial sums of the Moebius transform of the triangular numbers (A007438). - Steve Butler, Apr 18 2006
Row sums of triangle A134543. - Gary W. Adamson, Oct 31 2007
a(n) ~ n^3 / (6*Zeta(3)). - Vaclav Kotesovec, Jan 31 2019
G.f.: (1/(1 - x)) * Sum_{k>=1} mu(k) * x^k / (1 - x^k)^3. - Ilya Gutkovskiy, Feb 14 2020
a(n) = n*(n+1)*(n+2)/6 - Sum_{j=2..n} a(floor(n/j)) = A000292(n) - Sum_{j=2..n} a(floor(n/j)). - Chai Wah Wu, Mar 30 2021
EXAMPLE
a(4) = 15 because the 15 triples in question are in lexicographic order: [1,1,1], [1,1,2], [1,1,3], [1,1,4], [1,2,2], [1,2,3], [1,2,4], [1,3,3], [1,3,4], [1,4,4], [2,2,3], [2,3,3], [2,3,4], [3,3,4] and [3,4,4]. - Wolfdieter Lang, Apr 04 2013
The a(4) = 15 triangles with at least two sides <= 4 and sides relatively prime (see Henry Bottomley's comment above) are: [1,1,1], [1,2,2], [2,2,3], [1,3,3], [2,3,3], [2,3,4], [3,3,4], [3,3,5], [1,4,4], [2,4,5], [3,4,4], [3,4,5], [3,4,6], [4,4,5], [4,4,7]. - Alois P. Heinz, Feb 14 2020
MAPLE
with(numtheory):
b:= proc(n) option remember;
add(mobius(n/d)*d*(d+1)/2, d=divisors(n))
end:
a:= proc(n) option remember;
b(n) + `if`(n=1, 0, a(n-1))
end:
seq(a(n), n=1..60); # Alois P. Heinz, Feb 09 2011
MATHEMATICA
a[1] = 1; a[n_] := a[n] = Sum[MoebiusMu[n/d]*d*(d+1)/2, {d, Divisors[n]}] + a[n-1]; Table[a[n], {n, 1, 60}] (* Jean-François Alcover, Jan 20 2014, after Maple *)
Accumulate[Table[Sum[MoebiusMu[n/d]*d*(d + 1)/2, {d, Divisors[n]}], {n, 1, 50}]] (* Vaclav Kotesovec, Jan 31 2019 *)
PROG
(Magma) [n eq 1 select 1 else Self(n-1)+ &+[MoebiusMu(n div d) *d*(d+1)/2:d in Divisors(n)]:n in [1..50]]; // Marius A. Burtea, Feb 14 2020
(Python)
from functools import lru_cache
@lru_cache(maxsize=None)
def A015631(n):
if n == 0:
return 0
c, j = 1, 2
k1 = n//j
while k1 > 1:
j2 = n//k1 + 1
c += (j2-j)*A015631(k1)
j, k1 = j2, n//j2
return n*(n-1)*(n+4)//6-c+j # Chai Wah Wu, Mar 30 2021
(PARI) a(n) = sum(k=1, n, sumdiv(k, d, moebius(k/d)*binomial(d+1, 2))); \\ Seiichi Manyama, Jun 12 2021
(PARI) a(n) = binomial(n+2, 3)-sum(k=2, n, a(n\k)); \\ Seiichi Manyama, Jun 12 2021
(PARI) my(N=66, x='x+O('x^N)); Vec(sum(k=1, N, moebius(k)*x^k/(1-x^k)^3)/(1-x)) \\ Seiichi Manyama, Jun 12 2021
CROSSREFS
KEYWORD
nonn
AUTHOR
STATUS
approved