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.

A385133 The number of integers k from 1 to n such that gcd(n, k) has an odd maximum exponent in its prime factorization.

This page as a plain text file.
%I A385133 #12 Jun 25 2025 01:26:16
%S A385133 0,1,1,1,1,4,1,3,2,6,1,5,1,8,7,5,1,10,1,7,9,12,1,13,4,14,7,9,1,22,1,
%T A385133 11,13,18,11,12,1,20,15,19,1,30,1,13,16,24,1,23,6,28,19,15,1,32,15,25,
%U A385133 21,30,1,29,1,32,20,21,17,46,1,19,25,46,1,33,1,38,32
%N A385133 The number of integers k from 1 to n such that gcd(n, k) has an odd maximum exponent in its prime factorization.
%H A385133 Amiram Eldar, <a href="/A385133/b385133.txt">Table of n, a(n) for n = 1..10000</a>
%F A385133 a(n) = Sum_{k=1..n} (A051903(gcd(n, k)) mod 2).
%F A385133 a(n) = n - A385132(n).
%F A385133 a(n) = Sum_{k=1..kmax(n)} (-1)^k * Product_{i=1..r} f(p_i^e_i, k), for n >= 2; if n = Product_{i=1..r} p_i^e_i, r = omega(n) = A001221(n), then emax(n) = max(e_i) = A051903(n), kmax(n) = emax(n)+1 if emax(n) is odd, and emax(n) otherwise, f(p^e, k) = p^e - p^(e-k) if e >= k, and f(p^e, k) = p^e if e < k.
%F A385133 Sum_{k=1..n} a(k) ~ c * n^2 / 2, where c = Sum_{k>=1} (-1)^(k+1)*(1-1/zeta(2*k)) = 0.32979448728905469699... .
%t A385133 e[n_] := If[n == 1, 0, Max[FactorInteger[n][[;; , 2]]]]; a[n_] := Sum[Boole[OddQ[e[GCD[n, k]]]], {k, 1, n}]; Array[a, 100]
%t A385133 (* second program: *)
%t A385133 a[n_] := Module[{f = FactorInteger[n], p, e, emax, kmax}, p = f[[;;, 1]]; e = f[[;;, 2]]; emax = Max[e]; kmax = emax + Mod[emax, 2]; Sum[(-1)^k * Product[p[[i]]^e[[i]] - If[e[[i]] < k, 0, p[[i]]^(e[[i]]-k)], {i, 1, Length[p]}], {k, 1, kmax}]]; a[1] = 0; Array[a, 100]
%o A385133 (PARI) q(n) = if(n == 1, 0, vecmax(factor(n)[,2]) % 2);
%o A385133 a(n) = sum(k = 1, n, q(gcd(n, k)));
%o A385133 (PARI) a(n) = if(n == 1, 0, my(f = factor(n), p = f[,1], e = f[,2], emax = vecmax(e), kmax = emax + emax % 2); sum(k = 1, kmax, (-1)^k * prod(i = 1, #p, p[i]^e[i] - if(e[i] < k, 0, p[i]^(e[i]-k)))));
%Y A385133 Cf. A001221, A051903, A368714, A384655, A385132.
%K A385133 nonn,easy
%O A385133 1,6
%A A385133 _Amiram Eldar_, Jun 24 2025