%I #45 Aug 16 2024 08:36:44
%S 4,8,12,20,24,28,40,44,52,56,60,68,76,84,88,92,104,116,120,124,132,
%T 136,140,148,152,156,164,168,172,184,188,204,212,220,228,232,236,244,
%U 248,260,264,268,276,280,284,292,296,308,312,316,328,332,340,344,348,356
%N Positive integers k for which there is exactly one integer i in {1,2,3,...,k-1} such that i*k is a square.
%C It appears that all terms of this sequence are exactly four times those of the squarefree integers (A005117).
%C The observed behavior is true for all n. All positive integers n are written uniquely as k*m^2 where k is squarefree, k >=1, m >= 1. The square multiples of n are j^2*k*n, j >= 1. We seek n with exactly 1 multiple that is square and less than n^2. If m = 1, there are no such multiples as we have k = n, so the least square multiple is n^2. If m >= 2, k*n is square and less than n^2. However, 4*k*n also qualifies as square and less than n^2 if m > 2. So the qualifying values of n are those with m=2. - _Peter Munn_, Nov 28 2019
%C The asymptotic density of this sequence is 3/(2*Pi^2). - _Amiram Eldar_, Mar 08 2021
%H Reinhard Zumkeller, <a href="/A133466/b133466.txt">Table of n, a(n) for n = 1..1000</a>
%F A057918(a(n)) = 1. - _Reinhard Zumkeller_, Mar 27 2012
%F From _Peter Munn_, Nov 28 2019: (Start)
%F a(n) = 4 * A005117(n).
%F {a(n)} = {A225546(A007283(n)) : n >= 0}, where {a(n)} denotes the set of integers in the sequence.
%F (End)
%F Sum_{n>=1} 1/a(n)^s = zeta(s)/(4^s*zeta(2*s)), s>1. - _Amiram Eldar_, Sep 26 2023
%e 4 is in the sequence because among the products 1*4,2*4,3*4 = 4,8,12 there is exactly one square.
%t eoiQ[n_]:=Count[n*Range[n-1],_?(IntegerQ[Sqrt[#]]&)]==1; Select[Range[ 400],eoiQ] (* _Harvey P. Dale_, Mar 14 2015 *)
%o (Haskell)
%o a133466 n = a133466_list !! (n-1)
%o a133466_list = map (+ 1) $ elemIndices 1 a057918_list
%o -- _Reinhard Zumkeller_, Mar 27 2012
%o (PARI) isok(n) = sum(k=1, n-1, issquare(k*n)) == 1; \\ _Michel Marcus_, Nov 29 2019
%o (Magma) [k:k in [1..350]|#[m:m in [1..k-1]| IsSquare(m*k)] eq 1]; // _Marius A. Burtea_, Dec 03 2019
%o (Python)
%o from math import isqrt
%o from sympy import mobius
%o def A133466(n):
%o def f(x): return n+x-sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))
%o m, k = n, f(n)
%o while m != k:
%o m, k = k, f(k)
%o return int(m)<<2 # _Chai Wah Wu_, Aug 15 2024
%Y Cf. A005117, A007283, A225546.
%Y Cf. A057918, A195085.
%K nonn
%O 1,1
%A _John W. Layman_, Nov 28 2007