123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149 |
- #!/usr/bin/ruby
- # a(n) is the least integer m such that M(m) is divisible by prime(n)^2 or -1 if no such m exists.
- # https://oeis.org/A354293
- # See also:
- # https://oeis.org/A354292
- # Known terms:
- # 3, 4, -1, 23, 21, -1, 188, 65, 1010, 2231, -1, -1, 1326, 389, 1092, 13196, 1450, -1, 40466, 85553, 665, -1, 5139193, 333, -1, 408241, -1, 3072, 6702, 1393, 5832, 935, 1071, 77421, 292187, 775383, 493135, 4185, 1784560, 10632, 7935, 743003, 13418, 64499
- var known_terms = Hash(
- 1 => 3
- 2=> 4
- 5=> 21
- 4=> 23
- 8=> 65
- 7=> 188
- 24=> 333
- 14=> 389
- 21=> 665
- 9=> 1010
- 15=> 1092
- 13=> 1326
- 17=> 1450
- 10=> 2231
- 16=> 13196
- 19=> 40466
- 20=> 85553
- 26=> 408241
- 31 => 5832
- 32 => 935
- 33 => 1071
- 34 => 77421
- 35 => 292187
- 30 => 1393,
- 36 => 775383,
- 37 => 493135,
- 28 => 3072,
- 29 => 6702,
- 38 => 4185,
- 41 => 7935,
- 40 => 10632,
- 46 => 12176,
- 43 => 13418,
- 44 => 64499,
- 1 => 3
- 2 => 4
- 3 => -1
- 4 => 23
- 5 => 21
- 6 => -1
- 7 => 188
- 8 => 65
- 9 => 1010
- 10 => 2231
- 11 => -1
- 12 => -1
- 13 => 1326
- 14 => 389
- 15 => 1092
- 39 => 1784560,
- 42 => 743003,
- # Term for p = 83 is given in the paper:
- # https://arxiv.org/pdf/2205.09902.pdf
- 83.pi => 5139193,
- 5.pi => -1,
- 13.pi => -1,
- 31.pi => -1,
- 37.pi => -1,
- 61.pi => -1,
- 79.pi => -1,
- 97.pi => -1,
- 103.pi => -1,
- )
- say "Known terms:\n"
- var prev_k = 0
- known_terms.keys.map{.to_i}.sort.each {|k|
- k - prev_k == 1 || break
- if (known_terms{k} != -1) {
- if (known_terms{k} < 1e6) {
- assert_eq(motzkin(known_terms{k}) % k.prime**2, 0)
- }
- }
- print(known_terms{k}, ", ")
- prev_k = k
- }
- say "\n"
- #primes -= [ 5, 13, 31, 37, 61, 79, 97, 103 ]
- #var primes = [39.prime]
- var primes = 200.primes
- primes = primes.grep {|p|
- !known_terms.has(p.pi)
- }
- say primes
- func catalanmod(n,m) {
- if (is_coprime(n+1, m)) {
- divmod(binomialmod(2*n, n, m), n+1, m)
- }
- else {
- idiv(binomial(2*n, n), (n+1)) % m
- }
- }
- func isok(n, m) {
- sum(0..n, {|k|
- (-1)**(n-k) * mulmod(binomialmod(n,k,m), catalanmod(k+1, m), m)
- }) % m == 0
- }
- #say isok(2231, 10.prime**2)
- #say isok(5139193, 83**2)
- var (a, b, n) = (0, 1, 1)
- Inf.times {
- var M = idiv(b,n)
- n += 1
- (a, b) = (b, idiv((3*(n-1)*n*a + (2*n - 1)*n*b), ((n+1)*(n-1))))
- var found = []
- for p in (primes) {
- if (p.sqr `divides` M) {
- say "#{p.pi} => #{n-2}"
- found << p
- }
- }
- if (found) {
- primes -= found
- primes || break
- }
- }
|