linear-ineq.rb

require "mathn"

SHIFT = 16
SHIFT2 = SHIFT*2
A = 0x41c64e6d % 2**SHIFT2
B = 0x6073 % 2**SHIFT2

def f(x)
  (A * x + B) % 2**SHIFT2
end

def g(x)
  f(x) >> SHIFT
end

def extgcd(a, b)
  if b == 0 then
    return [a > 0 ? 1 : -1, 0]
  end
  (q, r) = a.divmod(b)
  (x, y) = extgcd(b, r)
  [y, x - q * y]
end

def solve_ineq(a, b, t0, t1, x0, x1, y0, y1)
  # 不等式 t0 <= ax + by <= t1
  #       かつ x0 <= x <= x1
  #       かつ y0 <= y <= y1 を解く
  if b == 0 then
    #    t0 <= ax <= t1
    # ~> t0 / a <= x <= t1 / a
    xx0 = (a / t0).ceil
    xx1 = (a / t1).floor
    return [[xx0, x0].max, [xx1, x1].min, y0, y1]
  end
  (q, r) = a.divmod(b)
  if r == 0
    #    t0 <= ax + qay <= t1
    # ~> t0 / a <= x + qy <= t1 / a
    return solve_ineq(1, q, (t0 / a).ceil, (t1 / a).floor, x0, x1)
  end
  #    t0 <= ax + (qa + r)y <= t1
  # ~> t0 <= a(x + qy) + ry <= t1
  # ~> t0 <= ry + az <= t1
  z1 = [x0.abs, x1.abs].max + q.abs * [y0.abs, y1.abs].max
  z0 = -z1
  sols = solve_ineq(r, a, t0, t1, y0, y1, z0, z1)
  sols.flat_map {|ymin, ymax, zmin, zmax|
    # あ、だめだ
  }
end

x0 = extgcd(A, -2**SHIFT)[0]
z0 = extgcd(A, -2**SHIFT)[1]
c = 4
d = 2**SHIFT * c - B
p x0
p z0
p x0 * A + z0 * -(2**SHIFT)
p -x0 * d
p (65535 - x0 * d)

table = []
(2**SHIFT).times do |i|
  table[g(i)] = [] if !table[g(i)]
  table[g(i)] << i
end
p table[0..10]