Denis Ivanov - cp's OEIS frontend

A368774 a(n) is the number of distinct real roots of the polynomial whose coefficients are the digits of the binary expansion of n.

0, 1, 1, 1, 0, 2, 0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 0, 2, 0, 1, 0, 2, 0, 2, 0, 2, 1, 1, 0, 2, 0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 0, 2, 0, 1, 0, 2, 0, 2, 0, 2, 2, 1, 0, 2, 0, 1, 0, 2, 0, 1, 0, 2, 0, 2, 0, 2, 0, 1, 0, 2, 0, 2, 0, 2, 1, 1

Offset: 1

Denis Ivanov, Jan 05 2024

In the full-form binary representation of n = c_m*2^m + c_{m-1}*2^(m-1) + ... + c_0*2^0 we replace 2 with x and count the distinct real roots of the polynomial f_n(x) = c_m*x^m + c_{m-1}*x^(m-1) + ... + c_0.
Multiple roots count only once.
The first occurrences a(n) = 3, 4, 5, 6 are at n = 150, 1686, 854646 and 437545878, respectively.
a(n) >= 1 if n is even or has an even number of binary digits. - Robert Israel, Jan 08 2024

			n = 1 = 1*2^0 -> f_n(x) = 1*x^0 = 1 -> no roots, a(1) = 0.
n = 6 = 1*2^2 + 1*2^1 + 0*2^0 -> f_n(x) = x^2 + x = x*(x + 1) -> roots {0,-1}, a(6) = 2.
n = 150 = 1*2^7 + 0*2^6 + 0*2^5 + 1*2^4 + 0*2^3 + 1*2^2 + 1*2^1 + 0*2^0 -> f_n(x) = x^7 + x^4 + x^2 + x -> roots {-1, -0.755..., 0}, a(150) = 3.

Robin Visser, Table of n, a(n) for n = 1..10000
M. Filaseta, C. Finch and C. Nicol, On three questions concerning 0,1-polynomials, Journal de théorie des nombres de Bordeaux, Tome 18 (2006) no. 2, pp. 357-370.
D. Ivanov et al., Number of real roots of 0,1 polynomial, MathOverflow.

Cf. A362344, A368824.

f:= proc(n) local L,p,i,z;
  L:= convert(n,base,2);
  p:= add(L[i]*z^(i-1),i=1..nops(L));
  sturm(sturmseq(p,z),z,-infinity,infinity);
end proc:
map(f, [$1..100]); # Robert Israel, Jan 08 2024

a[n_]:=Length@{ToRules@Reduce[FromDigits[IntegerDigits[n, 2], x] == 0, x, Reals]};

a(n) = #Set(polrootsreal(Pol(binary(n)))); \\ Michel Marcus, Jan 05 2024

from sympy.abc import x
from sympy import sturm, oo, sign, nan, LT
def A368774(n):
    if n == 1: return 0
    l = len(s:=bin(n)[2:])
    q = sturm(sum(int(s[i])*x**(l-i-1) for i in range(l)))
    a = [1 if (k:=LT(p).subs(x,-oo))==nan else sign(k) for p in q[:-1]]+[sign(q[-1])]
    b = [1 if (k:=LT(p).subs(x,oo))==nan else sign(k) for p in q[:-1]]+[sign(q[-1])]
    return sum(1 for i in range(len(a)-1) if a[i]!=a[i+1])-sum(1 for i in range(len(b)-1) if b[i]!=b[i+1]) # Chai Wah Wu, Feb 15 2024

def a(n):
    R. = PolynomialRing(ZZ); poly = R(list(bin(n)[:1:-1]))
    return poly.number_of_real_roots()  # Robin Visser, Aug 03 2024

A366544 a(n) is a lower bound for the number of distinct stable centroidal Voronoi tessellations (CVTs) of a square with n generators (seeds).

Original entry on oeis.org

1, 1, 1, 1, 2, 3, 3, 3, 2, 2, 3, 5, 8, 6, 5, 3, 4, 7, 10, 21, 21

Offset: 0

Denis Ivanov, Oct 12 2023

Stable CVTs are local minimizers of the CVT function (see first paper).
There are other CVTs which are saddle points.
Lloyd's process converges only to stable CVTs.
An efficient two-step semi-manual algorithm is used to recognize identical patterns and a fast code for the Lloyd's process.

			As initialization, clustering centers for a large number of points in the square are used. For every set of centers, Lloyd's algorithm is iterated and all variants symmetric with respect to rotations and reflections are removed.

Lin Lu, F. Sun, and H. Pan, Global optimization Centroidal Voronoi Tessellation with Monte Carlo Approach, 2012 IEEECS Log Number TVCG-2011-03-0067.

Denis Ivanov, Code, explanations and results (github).
J. C. Hateley, H. Wei, and L. Chen, Fast Methods for Computing Centroidal Voronoi Tessellations, J. Sci. Comput., 63, pp. 185-212, 2015.
Yang Liu, Wenping Wang, Bruno Lévy, Feng Sun, Dong-Ming Yan, Lin Lu, and Chenglei Yang, On centroidal Voronoi tessellation—Energy smoothness and fast computation, ACM Transactions on Graphics, Volume 28, Issue 4, Article No. 101, pp. 1-17, 2009.
Wikipedia, Centroidal Voronoi tessellation (unfortunately, article is a stub and contains inaccuracies).
Wikipedia, Lloyd's algorithm.

Cf. A363822 (disk).

A368824 a(n) is the smallest degree of (0,1)-polynomial with exactly n real distinct roots.

Original entry on oeis.org

1, 2, 7, 10, 19, 28

Offset: 1

Denis Ivanov, Jan 07 2024

(0,1) (or Newman) polynomials have coefficients from {0,1}.

P. Borwein, T. Erdélyi, and G. Kós, Littlewood-type problems on [0,1], Proc. London Math. Soc. 79 (1999), 22-46.
MathOverflow, Number of real roots of 0,1 polynomial.
A. Odlyzko and B. Poonen, Zeros of polynomials with 0,1 coefficients, L’Enseignement Mathématique 39 (1993), 317-348.

Cf. A368774, A362344.

(* Suitable only for 1 <= n <= 4;
for larger n, special Julia and Python packages are needed *)
Table[Exponent[Monitor[Catch[Do[
   poly = FromDigits[IntegerDigits[k, 2], x];
   res = Length@{ToRules@Reduce[poly == 0, x, Reals]};
   If[res == n, Throw@{res, Expand@poly}]
   , {k, 2000}]], k][[2]], x], {n, 1, 4}]

from itertools import count
from sympy.abc import x
from sympy import sturm, oo, sign, nan, LT
def A368824(n):
    for m in count(2):
        l = len(s:=bin(m)[2:])
        q = sturm(sum(int(s[i])*x**(l-i-1) for i in range(l)))
        a = [1 if (k:=LT(p).subs(x,-oo))==nan else sign(k) for p in q[:-1]]+[sign(q[-1])]
        b = [1 if (k:=LT(p).subs(x,oo))==nan else sign(k) for p in q[:-1]]+[sign(q[-1])]
        if n==sum(1 for i in range(len(a)-1) if a[i]!=a[i+1])-sum(1 for i in range(len(b)-1) if b[i]!=b[i+1]):
            return l-1 # Chai Wah Wu, Feb 15 2024

a(n) ~ C*n^2 (see Odlyzko and Poonen) with numerical estimate 0.7 < C < 0.9.

A363822 a(n) is the conjectured number of stable distinct centroidal Voronoi tessellations (CVTs) of a unit disk with n generators (seeds).

Original entry on oeis.org

1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 5, 6, 6, 5, 5, 5, 6, 9, 10

Offset: 0

Denis Ivanov, Oct 18 2023

Stable CVTs are local minimizers of the CVT function (see Hateley, Wei,and Chen article).
There are other CVTs which are saddle points.
Lloyd's process converges only to stable CVTs from which different with respect to rotation symmetry are selected.
An efficient two-step semi-manual algorithm is used to recognize identical patterns and a fast code for the Lloyd's process.
Code in Mathematica and details published on Github.

			As initialization, clustering centers for a large number of points in the unit disk are used. For every set of centers, Lloyd's algorithm is iterated and all variants symmetric with respect to rotations are removed.

J. C. Hateley, H. Wei, and L. Chen, Fast Methods for Computing Centroidal Voronoi Tessellations, 2014 J Sci Comput DOI 10.1007/10915-014-9894-1
Yang Liu, Wenping Wang, Bruno Lévy, Feng Sun, Dong-Ming Yan, Lin Lu, and Chenglei Yang, On centroidal Voronoi tessellation—Energy smoothness and fast computation, ACM Transactions on Graphics, Volume 28, Issue 4, Article No. 101, pp. 1-17, 2009, DOI 10.1145/1559755.1559758
Lin Lu, F. Sun, and H. Pan, Global optimization Centroidal Voronoi Tessellation with Monte Carlo Approach, 2012 IEEECS Log Number TVCG-2011-03-0067.

Denis Ivanov, Github with code, explanations and results.
Wikipedia, Centroidal Voronoi tessellation (unfortunately, article is a stub and contains inaccuracies).
Wikipedia, Lloyd's algorithm.

Cf. A366544 (square).

A364200 Minimal number of terms of mixed-sign Egyptian fraction f such that H(n) + f is an integer, where H(n) is the n-th harmonic number.

Original entry on oeis.org

0, 1, 1, 1, 2, 2, 3, 3, 3, 2, 2, 3, 4, 3, 4, 4, 4, 4, 4, 5, 5, 4, 5, 5, 5, 4, 5, 5, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6

Offset: 1

Denis Ivanov, Jul 13 2023

For H(n) - floor(H(n)) and ceiling(H(n)) - H(n), the shortest mixed-sign Egyptian fractions are calculated, and the smaller length of fractions is selected.
Similar to A106394 and A224820. But those sequences use the greedy algorithm, which does not guarantee the shortest length of expansion.
For 1 < n < 41, a(n) < A363937(n) only for n = 10 and n = 22.

			For n=10: H(10) = 7381/2520 = 2.928...; H(10) - floor(H(10)) = 7381/2520 - 2 = 2341/2520 = 1/2 + 1/7 + 1/8 + 1/9 + 1/20, which cannot be expressed as the sum of fewer than 5 reciprocals, and ceiling(H(10)) - H(10) = 3 - 7381/2520 = 179/2520 = 1/30 + 1/42 + 1/72, which cannot be expressed as the sum of fewer than 3 reciprocals, so A363937(10) = 3.
But 179/2520 = 1/14 - 1/2520 (a "mixed-sign Egyptian fraction"), so a(10) = 2.

Ron Knott, Egyptian Fractions

Cf. A363937.

check[f_, k_] := (If[Numerator@f == 1, Return@True];
   If[k == 1, Return@False];
   Catch[Do[If[check[f - 1/i, k - 1], Throw@True],
     {i, Range[Ceiling[1/f], Floor[k/f]]}];
    Throw@False]);
checkMixed[f_, k_, m_] := If[m == 1,
   Catch[Do[If[check[1/i - f, k], Throw@True],
     {i, Range[2, Floor[1/f]]}];
    Throw@False],
   checkMixed[f, k, m - 1]];
a[n_] := (h = HarmonicNumber[n];
  d = Min[h - Floor@h, Ceiling@h - h];
  j = 1;
  While[Not@check[d, j], j++];
  res = j;
  Do[
   If[checkMixed[d, i - m, m], res = i],
   {i, 2, j - 1}, {m, 1, i - 1}];
  res);

a(n) <= A363937(n).

A363937 Minimal number of terms of an Egyptian fraction to be added to, or subtracted from, harmonic number H(n) to get an integer.

Original entry on oeis.org

0, 1, 1, 1, 2, 2, 3, 3, 3, 3, 2, 3, 4, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 4, 5, 5, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6

Offset: 1

Denis Ivanov, Jun 29 2023

The shortest Egyptian fractions for H(n) - floor(H(n)) and ceiling(H(n)) - H(n) are calculated and the smaller length of those fractions is a(n).

			For n = 2: H(2) = 3/2 which is between 1 and 2 and they are reached by the same H(2) - 1 = 2 - H(2) = 1/2 which is 1 term, so a(2) = 1.
For n = 5: H(5) = 137/60 is between 2 and 3; going up 3 - H(5) = 1/2 + 1/6 + 1/20 is 3 terms but going down H(5) - 2 = 1/5 + 1/12 is 2 terms, so the latter is shorter and a(5) = 2 terms.

Ron Knott, Egyptian Fractions

Cf. A281530

(* Thanks to Ron Knott for the algorithm. Slow for n>15. *)
check[f_, k_]:= (If[Numerator@f == 1, Return@True];
  If[k == 1, Return@False];
  Catch[Do[
    If[check[f - 1/i, k - 1], Throw@True],
    {i, Range[Ceiling[1/f],Floor[k/f]]}];
   Throw@False]
  );
a[n_]:= (h = HarmonicNumber[n];
  d = {h - Floor[h], Ceiling[h] - h};
  j = 1;
  While[Not[Or @@ (check[#, j] & /@ d)], j++]; j);

a(31)-a(45) from Dmitry Petukhov, Jul 24 2023

A363501 a(n) = smallest product > n of some subset of the divisors of n, or if no product > n exists then a(n) = n.

Original entry on oeis.org

1, 2, 3, 8, 5, 12, 7, 16, 27, 20, 11, 18, 13, 28, 45, 32, 17, 27, 19, 40, 63, 44, 23, 32, 125, 52, 81, 56, 29, 36, 31, 64, 99, 68, 175, 48, 37, 76, 117, 50, 41, 63, 43, 88, 75, 92, 47, 64, 343, 100, 153, 104, 53, 81, 275, 64, 171, 116, 59, 72, 61, 124, 147

Offset: 1

Denis Ivanov, Jun 06 2023

nonn

a(n) = n iff n=1 or n is prime.

For composite n, A363627(n) < n < a(n) and where both bounds are products of divisors of n and as tight as possible.

			n = 4; divisors: [1,2,4]; subsets: [[], [1], [2], [4], [1, 2], [1, 4], [2, 4], [1, 2, 4]]; products: [1, 1, 2, 4, 2, 4, 8, 8]; minimal product that greater than 4 is 8, so a(4) = 8.
n = 5; divisors: [1,5]; subsets: [[], [1], [5], [1, 5]]; products: [1, 1, 5, 5]; no products greater than 5, so a(5) = 5.

Cf. A027750, A363627.

a[n_] := If[PrimeQ@n || n == 1, n,
  First@Select[Union[Times @@@ Subsets[Divisors@n]], # > n &]];

a(n) = my(d=divisors(n), nb = #d, m=oo); forsubset(nb, s, my(p=vecprod(vector(#s, k, d[s[k]]))); if (p>n, m=min(m, p))); if (mMichel Marcus, Jun 17 2023

A363627 a(n) = greatest product < n of some subset of the divisors of n, or if n is in A008578 then a(n) = n.

Original entry on oeis.org

1, 2, 3, 2, 5, 3, 7, 4, 3, 5, 11, 8, 13, 7, 5, 8, 17, 12, 19, 10, 7, 11, 23, 18, 5, 13, 9, 14, 29, 20, 31, 16, 11, 17, 7, 27, 37, 19, 13, 32, 41, 36, 43, 22, 27, 23, 47, 36, 7, 25, 17, 26, 53, 36, 11, 32

Offset: 1

Denis Ivanov, Jun 12 2023

nonn

a(n) = n <=> n in A008578.

For composite n, a(n) < n < A363501(n) and where both bounds are products of divisors of n and as tight as possible.

			n = 4; divisors: [1,2,4]; subsets: [[], [1], [2], [4], [1, 2], [1, 4], [2, 4], [1, 2, 4]]; products: [1, 1, 2, 4, 2, 4, 8, 8]; the maximal product that is lesser than 4 is 2, so a(4) = 2.

Cf. A008578, A363501.

If[PrimeQ@n || n == 1, n,
  Last@Select[Union[Times @@@ Subsets[Divisors@n]], # < n &]];

a(n) = my(d=divisors(n), nb = #d, m=1); forsubset(nb, s, my(p=vecprod(vector(#s, k, d[s[k]]))); if (p1, m, n); \\ Michel Marcus, Jun 17 2023

cp's OEIS Frontend

User: Denis Ivanov

A368774 a(n) is the number of distinct real roots of the polynomial whose coefficients are the digits of the binary expansion of n.

Views

Author

Keywords

Comments

Examples

Links

Crossrefs

Programs

Maple

Mathematica

PARI

Python

Sage

A366544 a(n) is a lower bound for the number of distinct stable centroidal Voronoi tessellations (CVTs) of a square with n generators (seeds).

Views

Author

Keywords

Comments

Examples

References

Links

Crossrefs

A368824 a(n) is the smallest degree of (0,1)-polynomial with exactly n real distinct roots.

Views

Author

Keywords

Comments

Links

Crossrefs

Programs

Mathematica

Python

Formula

A363822 a(n) is the conjectured number of stable distinct centroidal Voronoi tessellations (CVTs) of a unit disk with n generators (seeds).

Views

Author

Keywords

Comments

Examples

References

Links

Crossrefs

A364200 Minimal number of terms of mixed-sign Egyptian fraction f such that H(n) + f is an integer, where H(n) is the n-th harmonic number.

Views

Author

Keywords

Comments

Examples

Links

Crossrefs

Programs

Mathematica

Formula

A363937 Minimal number of terms of an Egyptian fraction to be added to, or subtracted from, harmonic number H(n) to get an integer.

Views

Author

Keywords

Comments

Examples

Links

Crossrefs

Programs

Mathematica

Extensions

A363501 a(n) = smallest product > n of some subset of the divisors of n, or if no product > n exists then a(n) = n.

Views

Author

Keywords

Comments

Examples

Crossrefs

Programs

Mathematica

PARI

A363627 a(n) = greatest product < n of some subset of the divisors of n, or if n is in A008578 then a(n) = n.

Views

Author