12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667 |
- #!/usr/bin/ruby
- # a(n) is the smallest centered n-gonal number with exactly n prime factors (counted with multiplicity).
- # https://oeis.org/A358926
- # Known terms:
- # 316, 1625, 456, 3964051, 21568, 6561, 346528
- # PARI/GP program:
- # a(n) = if(n<3, return()); for(k=1, oo, my(t=((n*k*(k+1))/2+1)); if(bigomega(t) == n, return(t))); \\ ~~~~
- # PARI/GP generation program:
- #`(
- bigomega_polygonals(A, B, n, k) = A=max(A, 2^n); (f(m, p, n) = my(list=List()); if(n==1, forprime(q=max(p,ceil(A/m)), B\m, my(t=m*q); if(issquare((8*(t-1))/k + 1) && ((sqrtint((8*(t-1))/k + 1)-1)%2 == 0), listput(list, t))), forprime(q = p, sqrtnint(B\m, n), my(t=m*q); if(ceil(A/t) <= B\t, list=concat(list, f(t, q, n-1))))); list); vecsort(Vec(f(1, 2, n)));
- a(n, k=n) = if(k < 3, return()); my(x=2^n, y=2*x); while(1, my(v=bigomega_polygonals(x, y, n, k)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ ~~~~
- # Version with primes in certain residue classes
- bigomega_polygonals(A, B, n, k) = A=max(A, 2^n); (f(m, p, n) = my(list=List()); if(n==1, forprime(q=max(p,ceil(A/m)), B\m, my(s=q%k); if(s==1 || s==7 || s==9 || s==15, my(t=m*q); if(issquare((8*(t-1))/k + 1) && ((sqrtint((8*(t-1))/k + 1)-1)%2 == 0), listput(list, t)))), forprime(q = p, sqrtnint(B\m, n), my(s=q%k); if(s==1 || s==7 || s==9 || s==15, my(t=m*q); if(ceil(A/t) <= B\t, list=concat(list, f(t, q, n-1)))))); list); vecsort(Vec(f(1, 2, n)));
- a(n, k=n) = if(k < 3, return()); my(x=2^n, y=2*x); while(1, my(v=bigomega_polygonals(x, y, n, k)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ ~~~~
- )
- func upper_bound(n, from = 2, upto = 2*from) {
- say "\n:: Searching an upper-bound for a(#{n})\n"
- loop {
- var count = n.almost_prime_count(from, upto)
- if (count > 0) {
- say "Sieving range: [#{from}, #{upto}]"
- say "This range contains: #{count.commify} elements\n"
- n.almost_primes_each(from, upto, {|v|
- if (v.is_centered_polygonal(n)) {
- say "a(#{n}) = #{v}"
- return v
- }
- })
- }
- from = upto+1
- upto *= 2
- }
- }
- upper_bound(16)
- __END__
- a(3) = 316
- a(4) = 1625
- a(5) = 456
- a(6) = 3964051
- a(7) = 21568
- a(8) = 6561
- a(9) = 346528
- a(10) = 3588955448828761
- a(11) = 1684992
- a(12) = 210804461608463437
- a(13) = 36865024
- a(14) = 835904150390625
- a(15) = 2052407296
|