160412

Ruby


Zagier's problems(7)

9 以下のm に対し、
U(m, n) = 0 mod ((m + 1)mn + m) となる10^6 以下のn を探してみた。

m = 1 のとき、解なし。

(http://math.stackexchange.com/questions/1736156/are-there-any-number-n-such-that-a-n-0-mod-2n-1-where-a-0-1-a-1)

m = 2 のとき、2755452, 4570452

m = 3 のとき、3006, 401781, 528531

m = 4 のとき、見つからず。

m = 5 のとき、196
m = 6 のとき、見つからず。
m = 7 のとき、26
m = 8 のとき、見つからず。
m = 9 のとき、見つからず。

その際使ったコードは以下のとおりです。
（一部haruya さんのコード(http://d.hatena.ne.jp/haruya12/20160131) を使用しています。）

require 'prime'
require 'matrix'

def f(n)
  return 1 if n <= 1
  (1..n).inject(:*)
end

def u(m, n)
  i = (0..n).inject(0){|s, k| s + f(m * n + k).to_r / (f(m * n - m * k) * f((m + 1) * k + 1))}
  j = i.denominator
  return i.to_i if j == 1
  i
end

def U(m, n)
  i = m * ((m + 1) * n + 1) * u(m, n)
  j = i.denominator
  return i.to_i if j == 1
  i
end

def berlekamp_massey(s, q)
  b, c = [1], [1] + [0] * (s.size - 1)
  l, m, a = 0, -1, 1
  s.size.times do |n|
    d = (0..l).inject(0) {|sum, i| (sum + c[i] * s[n - i]) % q}
    next if d == 0
    t = c[0..l]
    (0...[s.size - n + m, b.size].min).each do |j|
      c[n - m + j] = (c[n - m + j] - d * a * b[j]) % q
    end
    b, l, m, a = t, n + 1 - l, n, mod_inv(d, q) if 2 * l <= n
  end
  c[0..l]
end

def euclid(a, b)
  return [0, 1] if a == 0
  q, r = b.divmod(a)
  x, y = euclid(r, a)
  [y - q * x, x]
end

# x^(-1) (mod n)
def mod_inv(x, n)
  euclid(x, n)[0]
end

# x % n1 = r1, x % n2 = r2, |x| <= n1 * n2 / 2 となる x
def chinese(n1, r1, n2, r2)
  x = (n1 * (r2 - r1) * mod_inv(n1, n2) + r1) % (n = n1 * n2)
  2 * x > n ? x - n : x
end

# f を多項式として f=0 が数列 s を生成する漸化式の特性方程式となっているか
def test(f, s)
  (0..s.size - f.size).all? do |i|
    f.each_with_index.inject(0) {|sum, (fj, j)|
      sum + fj * s[f.size + i - j - 1]
    } == 0
  end
end

# 数列 s を生成する漸化式の特性方程式を返す
def polynomial(s)
  f, n = [], 1
  Prime.each do |q|
    c = berlekamp_massey(s, q)
    if c.size != f.size then
      f, n = c, q if c.size > f.size
      next
    end
    f = (0...f.size).map {|i| chinese(n, f[i], q, c[i])}
    return f if test(f, s)
    n *= q
  end
end

def power(a, n, mod)
  return Matrix.I(a.row_size) if n == 0
  m = power(a, n >> 1, mod)
  m = (m * m).map{|i| i % mod}
  return m if n & 1 == 0
  (m * a).map{|i| i % mod}
end

(1..9).each{|i|
  a = (0..24).map{|j| U(i, j)}
  ary0 = a[0..i].reverse
  v = Vector.elements(ary0)
  b = polynomial(a)[1..-1].map{|i| -i}
  ary1 = [b]
  (1..i).each{|j|
    c = Array.new(i + 1, 0)
    c[j - 1] = 1
    ary1 << c
  }
  a = Matrix[*ary1]
  p [i, v, a]
  (1..10 ** 6).each{|n|
    mod = (i + 1) * i * n + i
    p n if (power(a, n, mod) * v)[i] % mod == 0
  }
}

出力結果
[1, Vector[4, 1], Matrix[[3, -1], [1, 0]]]
[2, Vector[51, 10, 2], Matrix[[5, 2, 1], [1, 0, 0], [0, 1, 0]]]
[3, Vector[2395, 219, 18, 3], Matrix[[7, 39, 47, -1], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]]]
3006
401781
528531
[4, Vector[441221, 19620, 972, 28, 4], Matrix[[9, 244, 974, 44, 1], [1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0]]]
[5, Vector[320160441, 7408255, 154000, 4360, 40, 5], Matrix[[11, 1205, 12240, -6915, 1862, -1], [1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 1, 0]]]
196
[6, Vector[947687003440, 10353349563, 124825895, 1172628, 19401, 54, 6], Matrix[[13, 5466, 120994, -573455, 736546, -1401, 1], [1, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0]]]
[7, Vector[10849848134049807, 62681716148597, 329478589895, 2110500091, 8741607, 85211, 70, 7], Matrix[[15, 23919, 1050959, -22723953, 110533535, 1824095, 118911, -1], [1, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0]]]
26
[8, Vector[515038338884349001737, 1398347293343723656, 4149291508997544, 10350461259304, 35725792824, 64181272, 369672, 88, 8], Matrix[[17, 102824, 8460428, -660176392, 9869354006, 6085682408, 1508966588, -39688, 1], [1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0]]]
[9, Vector[93786705206867081003177917, 135131992303972351151403, 178716120071486628960, 275056422056846730, 321805755708813, 604220806080, 465906573, 1587420, 108, 9], Matrix[[19, 437409, 64949424, -16192124931, 648337672368, 2952342946923, 3894121358913, 4763358099, 11078999, -1], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1, 0]]]