cp's OEIS Frontend

This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.

Showing 1-2 of 2 results.

A018806 Sum of gcd(x, y) for 1 <= x, y <= n.

Original entry on oeis.org

1, 5, 12, 24, 37, 61, 80, 112, 145, 189, 220, 288, 325, 389, 464, 544, 593, 701, 756, 880, 989, 1093, 1160, 1336, 1441, 1565, 1700, 1880, 1965, 2205, 2296, 2488, 2665, 2829, 3028, 3328, 3437, 3621, 3832, 4152, 4273, 4621, 4748, 5040, 5373, 5597, 5736, 6168
Offset: 1

Views

Author

David W. Wilson

Keywords

Comments

a(n) is also the entrywise 1-norm of the n X n GCD matrix.

Links

Crossrefs

Cf. A000010, A018805, A064951, A344479.

Programs

  • Maple
    N:= 1000 # to get a(1) to a(N)
g:= add(numtheory:-phi(k)*x^k*(1+x^k)/((1-x^k)^2*(1-x)),k=1..N):
S:= series(g, x, N+1):
seq(coeff(S,x,j), j=1..N); # Robert Israel, Jan 14 2015
  • Mathematica
    Table[nn = n;Total[Level[Table[Table[GCD[i, j], {i, 1, nn}], {j, 1, nn}], {2}]], {n, 1, 48}] (* Geoffrey Critzer, Jan 14 2015 *)
  • PARI
    a(n)=2*sum(i=1,n,sum(j=1,i-1,gcd(i,j)))+n*(n+1)/2 \\ Charles R Greathouse IV, Jun 21 2013
  • PARI
    a(n)=sum(k=1,n,eulerphi(k)*(n\k)^2) \\ Charles R Greathouse IV, Jun 21 2013
  • Python
    from sympy import totient
def A018806(n): return sum(totient(k)*(n//k)**2 for k in range(1,n+1)) # Chai Wah Wu, Aug 05 2024

Formula

Sum_{k=1..n} phi(k)*(floor(n/k))^2. - Vladeta Jovovic, Nov 10 2002
a(n) ~ kn^2 log n, with k = 6/Pi^2. - Charles R Greathouse IV, Jun 21 2013
G.f.: Sum_{k >= 1} phi(k)*x^k*(1+x^k)/((1-x^k)^2*(1-x)). - Robert Israel, Jan 14 2015

A367690 Total number of steps of Euclid's GCD algorithm to calculate gcd(x,y) for all pairs x,y in the range 1 <= x,y <= n.

Original entry on oeis.org

1, 5, 14, 26, 47, 67, 100, 136, 177, 221, 286, 338, 415, 491, 568, 652, 761, 861, 990, 1098, 1219, 1351, 1520, 1652, 1813, 1985, 2162, 2334, 2555, 2727, 2960, 3172, 3397, 3641, 3878, 4098, 4383, 4659, 4944, 5204, 5537, 5805, 6158, 6466, 6775, 7123, 7524, 7848
Offset: 1

Views

Author

Darío Clavijo, Nov 26 2023

Keywords

Comments

The number of steps to calculate gcd(x,y) is A107435(x,y) for x<=y, and is the length of the continued fraction expansion of x/y (including 0 for the integer part).
A018806(n) <= a(n) <= A064951(n).

Links

Crossrefs

Cf. A018806, A049826, A064951, A107435.

Programs

  • Maple
    g:= proc(x, y) option remember;
      `if`(y=0, 0, 1+g(y, irem(x, y)))
    end:
a:= proc(n) option remember; `if`(n=0, 0,
      a(n-1)+n+2*add(g(n, j), j=1..n-1))
    end:
seq(a(n), n=1..100);  # Alois P. Heinz, Nov 27 2023
  • Mathematica
    g[x_, y_] := g[x, y] = If[y == 0, 0, 1 + g[y, Mod[x, y]]];
a[n_] := a[n] = If[n == 0, 0, a[n-1] + n + 2*Sum[g[n, j], {j, 1, n-1}]];
Table[a[n], {n, 1, 100}] (* Jean-François Alcover, Apr 19 2025, after Alois P. Heinz *)
  • PARI
    A107435(n, k) = if (k == 0, 0, 1 + A107435(k, n % k));
a(n) = sum(x=1, n, sum(y=1, n, A107435(x,y)));
print(vector(49,n,a(n)));
  • Python
    from functools import cache
A107435 = lambda x,y: 0 if y == 0 else 1 + A107435(y, x % y)
A049826 = lambda n: sum(A107435(n,j) for j in range(1, n))
@cache
def a(n):
  # Code after Alois P. Heinz
  if n == 0: return 0
  return a(n-1) + n + A049826(n) * 2
print([a(n) for n in range(1,49)])

Formula

a(n) = a(n-1) + n + A049826(n) * 2 and a(1) = 1.
a(n) = a(n-1) + n + (Sum_{j=1..n-1} A107435(n,j)) * 2 and a(1) = 1.
a(n) = Sum_{x=1..n} Sum_{y=1..n} A107435(x,y). - Michel Marcus, Nov 28 2023
Showing 1-2 of 2 results.

Started in 1964 by Neil J. A. Sloane | Maintained by The OEIS Foundation Inc.

Content is available under The OEIS End-User License Agreement.