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.
a(3) = 66 because 66 = 11*12/2 is the 11th triangular number and is the product of 3 distinct primes 2*3*11. a(4) = 0 because a 4-gonal number is a square, and thus not the product of distinct primes.
f:= proc(s) local n,p,F;
for n from 1 do
p:= (s-2)*n*(n-1)/2 + n;
F:= ifactors(p)[2];
if nops(F) = s and andmap(t -> t[2]=1, F) then return p fi
od
end proc:
f(2):= 0:
map(f, [$2..11]);
f[s_] := f[s] = Module[{n, p, F}, For[n = 1, True, n++, p = (s - 2)*n*(n-1)/2 + n; F = FactorInteger[p]; If[Length[F] == s && AllTrue[F, #[[2]] == 1&], Return[ p]]]];
f[4] = 0;
Table[Print[n, " ", f[n]]; f[n], {n, 2, 11}] (* Jean-François Alcover, Jan 24 2023, after Maple program *)
squarefree_omega_polygonals(A, B, n, k) = A=max(A, vecprod(primes(n))); (f(m, p, j) = my(list=List()); my(s=sqrtnint(B\m, j)); if(j==1, forprime(q=max(p, ceil(A/m)), s, if(ispolygonal(m*q, k), listput(list, m*q))), forprime(q=p, s, my(t=m*q); if(ceil(A/t) <= B\t, list=concat(list, f(t, q+1, j-1))))); list); vecsort(Vec(f(1, 2, n))); a(n, k=n) = if(n < 2, return()); if(n==2, return(6)); if(n==4, return(0)); my(x=vecprod(primes(n)), y=2*x); while(1, my(v=squarefree_omega_polygonals(x, y, n, k)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ Daniel Suteu, Jan 18 2023
a(3) = 37 = 1 + 15 + 21 = 3 + 6 + 28 = 6 + 10 + 21. a(4) = 161 = 1^2 + 4^2 + 12^2 = 2^2 + 6^2 + 11^2 = 4^2 + 8^2 + 9^2 = 5^2 + 6^2 + 10^2.
a(3) = 4 subsets: {1}, {3}, {6}, {1, 3, 6}.
a(4) = 5 subsets: {1}, {4}, {9}, {16}, {9, 16}.
a(5) = 7 subsets: {1}, {5}, {12}, {22}, {35}, {1, 12, 22}, {1, 12, 22, 35}.
from functools import cache
from itertools import count, takewhile
def ngonal(n, k): return k*((n-2)*k - (n-4))//2
def a(n):
@cache
def b(i, s):
if i == 0: return 1 if s > 0 and s in ISNGONAL else 0
return b(i-1, s) + b(i-1, s+NGONAL[i-1])
NGONAL = [ngonal(n, i) for i in range(1, n+1)]
BOUND = sum(NGONAL)
ISNGONAL = set(takewhile(lambda x: x<=BOUND, (ngonal(n, i) for i in count(1))))
b.cache_clear()
return b(n, 0)
print([a(n) for n in range(2, 23)]) # Michael S. Branicky, Dec 21 2024
a(0) = 3; a(1) = 6 and 6 is the 3rd triangular number; a(2) = 66 and 66 is the 6th hexagonal number; a(3) = 137346 and 137346 is the 66th 66-gonal number, etc.
RecurrenceTable[{a[0] == 3, a[n] == a[n - 1] (a[n - 1]^2 - 3 a[n - 1] + 4)/2}, a[n], {n, 5}]
For n = 5 we have: ---------------------------- 0 1 2 3 4 [5] ---------------------------- 0, 1, 5, 12, 22, 35, ... A000326 (pentagonal numbers) 0, 1, 6, 18, 40, 75, ... A002411 (pentagonal pyramidal numbers) 0, 1, 7, 25, 65, 140, ... A001296 (4-dimensional pyramidal numbers) 0, 1, 8, 33, 98, 238, ... A051836 (partial sums of A001296) 0, 1, 9, 42, 140, 378, ... A051923 (partial sums of A051836) 0, 1, 10, 52, 192, [570], ... A050494 (partial sums of A051923) ---------------------------- therefore a(5) = 570.
Table[n (n^2 - 2 n + 4) Binomial[2 n, n]/((n + 1) (n + 2)), {n, 0, 27}]
nmax = 27; CoefficientList[Series[(-4 + 31 x - 66 x^2 + 28 x^3 + (4 - 7 x) (1 - 4 x)^(3/2))/(2 x^2 (1 - 4 x)^(3/2)), {x, 0, nmax}], x]
nmax = 27; CoefficientList[Series[Exp[2 x] (4 - x + 2 x^2) BesselI[1, 2 x]/x - 2 Exp[2 x] (2 - x) BesselI[0, 2 x], {x, 0, nmax}], x] Range[0, nmax]!
Table[SeriesCoefficient[x (1 - 3 x + n x)/(1 - x)^(n + 3), {x, 0, n}], {n, 0, 27}]
A308488_vec(lim,J=10^6)={my( pyramid(s,n)=(3*n^2 + n^3*(s-2)-n*(s-5))/6, check(s)=j=if(lift(Mod(s,3))==2,((s-2)^2)/3-2,J);m=3;while(m<=j,if(ispolygonal(pyramid(s,m),s),return(pyramid(s,m)),m++));0); vector(lim,s,check(s+2))}
For n = 3: 1 + 6 + 21 = 3 + 10 + 15 = 28.
a(3) = 19 = 1 + 3 + 15 = 3 + 6 + 10. a(4) = 90 = 1^2 + 2^2 + 6^2 + 7^2 = 1^2 + 3^2 + 4^2 + 8^2.
Comments