From Wikipedia, the free encyclopedia

return (function()

local builders = {}

local function register(name, f)

  buildersname = f

end

register('llpeg', function() return require [[Module:User:Cscott/llpeg]] end)



register('util', function(myrequire)

local function read_wiki_input(func)

    return function (frame, ...)

      if type(frame)=='string' then

        frame = { args = { frame, ... } }

      end

      local title = mw.title.new(frame.args1])

      local source = title:getContent()

      if source == nil then

        error("Can't find title " .. tostring(title))

      end

      source = source:gsub("^%s*<syntaxhighlight[^>]*>\n?", "", 1)

      source = source:gsub("</syntaxhighlight[^>]*>%s*$", "", 1)

      return func(source, frame.args2], frame.args3])

    end

end



return {

  read_wiki_input = read_wiki_input,

}



end)



register('advent.compat', function() return require [[Module:User:Cscott/compat]] end)



register('eq', function(myrequire)

-- "fix" lua's eq metamethod to be consistent with __add etc, that is:

-- try the metamethod if any operand is not a number

local function eq(a, b)

  local tya, tyb = type(a), type(b)

  if tya ~= 'number' or tyb ~= 'number' then

    local op = nil

    local mt = getmetatable(a)

    if mt ~= nil then op = mt.__eq end

    if op == nil then

      mt = getmetatable(b)

      if mt ~= nil then op = mt.__eq end

    end

    if op ~= nil then

      return op(a, b)

    end

  end

  return a == b

end



return eq



end)



register('lt', function(myrequire)

-- "fix" lua's lt metamethod to be consistent with __add etc, that is:

-- try the metamethod if any operand is not a number

local function lt(a, b)

  local tya, tyb = type(a), type(b)

  if tya ~= 'number' or tyb ~= 'number' then

    local op = nil

    local mt = getmetatable(a)

    if mt ~= nil then op = mt.__lt end

    if op == nil then

      mt = getmetatable(b)

      if mt ~= nil then op = mt.__lt end

    end

    if op ~= nil then

      return op(a, b)

    end

  end

  return a < b

end



return lt



end)



register('bignum', function(myrequire)

local compat = myrequire('advent.compat')

local eq = myrequire('eq')

local lt = myrequire('lt')



-- poor man's bignum library

local RADIX = 1000 -- power of 10 to make printout easier



local function maxdigits(a, b)

  if #a > #b then return #a end

  return #b

end



local function ltz(a)

  if type(a) == 'number' then

    return a < 0

  end

  return a.negative or false

end



local BigNum = {}

BigNum.__index = BigNum

function BigNum:new(n)

  if n == nil then n = 0 end

  assert(type(n)=='number')

  if n < 0 then

    return setmetatable( {-n, negative=true}, self):normalize()

  else

    return setmetatable( {n}, self):normalize()

  end

end

function BigNum:__tostring()

   local result = {}

   if self.negative then

      table.insert(result, "-")

   end

  local first = true

  for i=#self,1,-1 do

    local n = selfi

    if n ~= 0 or not first then

      local s = tostring(n)

      if first then

        first = false

      else

        while #s < 3 do s = '0' .. s end

      end

      table.insert(result, s)

    end

  end

  if #result == 0 then return "0" end

  return table.concat(result)

end

function BigNum:toNumber()

  local m = 1

  local sum = 0

  for i=1,#self do

    sum = sum + (selfi * m)

    m = m * RADIX

  end

  return sum

end

function BigNum:normalize()

  local i = 1

  local sawNonZero = false

  while selfi ~= nil do

    assert(selfi >= 0)

    if selfi > 0 then

      sawNonZero = true

    end

    if selfi >= 1000 then

      local carry = math.floor(selfi / 1000)

      selfi = selfi % 1000

      selfi+1 = (selfi+1 or 0) + carry

    end

    i = i + 1

  end

  if not sawNonZero then

    self.negative = nil -- -0 is 0

  end

  return self

end

function BigNum:copy()

  local r = BigNum:new()

  for i=1,#self do

    ri = selfi

  end

  r.negative = self.negative

  return r

end

function BigNum.__unm(a)

  if eq(a, 0) then return a end -- -0 is 0

  local r = a:copy()

  r.negative = not r.negative

  return r

end

function BigNum.__add(a, b)

  if ltz(b) then

    if ltz(a) then

      return -((-a) + (-b))

    end

    return a - (-b)

  end

  if ltz(a) then

    return b - (-a)

  end

  assert(not ltz(a))

  assert(not ltz(b))

  if type(a) == 'number' then

    a,b = b,a

  end

  assert(not a.negative)

  local r = a:copy()

  if type(b) == 'number' then

    assert(b >= 0)

    r1 = (r1 or 0) + b

  else

    assert(not b.negative)

    for i=1,#b do

      ri = (ri or 0) + bi

    end

  end

  return r:normalize()

end

function BigNum.__lt(a, b)

  if ltz(a) then

    if ltz(b) then

      return (-a) > (-b)

    end

    return true

  elseif ltz(b) then

    return false

  end

  if type(a) == 'number' then a = BigNum:new(a) end

  if type(b) == 'number' then b = BigNum:new(b) end

  for i=maxdigits(a,b),1,-1 do

    if (ai or 0) < (bi or 0) then return true end

    if (ai or 0) > (bi or 0) then return false end

  end

  return false -- they are equal

end

function BigNum.__le(a, b)

  return not (a > b)

end

function BigNum.__eq(a, b)

  if ltz(a) ~= ltz(b) then return false end

  if type(a) == 'number' then a = BigNum:new(a) end

  if type(b) == 'number' then b = BigNum:new(b) end

  for i=1,maxdigits(a,b) do

    if (ai or 0) ~= (bi or 0) then return false end

  end

  return true

end

function BigNum.__sub(a, b)

  if ltz(b) then

    return a + (-b)

  end

  if ltz(a) then

    return -((-a) + b)

  end

  if type(a) == 'number' then a = BigNum:new(a) end

  if type(b) == 'number' then b = BigNum:new(b) end

  if b > a then

    return -(b - a)

  end

  local r = a:copy()

  local borrow = 0

  for i=1,maxdigits(a,b) do

    ri = (ri or 0) - (bi or 0) - borrow

    borrow = 0

    while ri < 0 do

      ri = ri + RADIX

      borrow = borrow + 1

    end

    assert(ri >= 0)

  end

  assert(borrow == 0)

  return r:normalize() -- shouldn't actually be necessary?

end



function BigNum.__mul(a, b)

  if type(a) == 'number' then

    a,b = b,a

  end

  local r = BigNum:new()

  if type(b) == 'number' then

    if b < 0 then r.negative = true ; b = -b ; end

    assert(b>=0)

    for i=1,#a do

      ri = ai * b

    end

    if a.negative then r.negative = not r.negative end

    return r:normalize()

  end

  for i=1,#a do

    for j=1,#b do

      assert(ai >= 0)

      assert(bj >= 0)

      local prod = ai * bj

      ri+j-1 = (ri+j-1 or 0) + prod

    end

  end

  r.negative = a.negative

  if b.negative then r.negative = not r.negative end

  return r:normalize()

end



function toBinary(a)

  assert(not a.negative)

  local bits = {}

  local RADIX_DIV_2 = compat.idiv(RADIX, 2)

  while a1 ~= nil do

    table.insert(bits, a1 % 2)

    for i=1,#a do

      ai = compat.idiv(ai], 2) + ((ai+1 or 0) % 2) * RADIX_DIV_2

    end

    if a#a == 0 then a#a = nil end

  end

  return bits

end



local function divmod(a, b)

  if eq(b, 0) then error("division by zero") end

  if ltz(b) then

    if ltz(a) then return divmod(-a, -b) end

    local q,r = divmod(a, -b)

    if lt(0, r) then q = q + 1 end

    return -q, -r

  elseif ltz(a) then

    local q,r = divmod(-a, b)

    if lt(0, r) then q = q + 1 end

    return -q, r

  end

  -- ok, a and b are both positive now

  assert(not ltz(a))

  assert(not ltz(b))

  if type(a) == 'number' then a = BigNum:new(a) end

  if type(b) == 'number' then b = BigNum:new(b) end

  local N,D = a,b

  local Q,R = BigNum:new(0), BigNum:new(0)

  local Nbits = toBinary(N:copy())

  for i=#Nbits,1,-1 do

    --print("R=",R,"Q=",Q)

    for i=1,#R do Ri = Ri * 2 end

    if Nbitsi > 0 then R1 = R1 + 1 end

    R:normalize()

    for i=1,#Q do Qi = Qi * 2 end

    if R >= D then

      R = R - D

      Q1 = Q1 + 1

    end

    Q:normalize()

  end

  return Q,R

end



function BigNum.__idiv(a, b)

  local q,r = divmod(a, b)

  return q

end



function BigNum.__mod(a, b)

  local q,r = divmod(a, b)

  return r

end



--[[

print(BigNum:new(4) >= BigNum:new(2))

print(BigNum:new(4) > BigNum:new(2))

print(BigNum:new(2) >= BigNum:new(2))

print(BigNum:new(2) > BigNum:new(2))

print(BigNum:new(4653) // BigNum:new(210))

]]--



return BigNum



end)



register('gcd', function(myrequire)

local eq = myrequire('eq')



local function gcd(a, b)

  if eq(b, 0) then return a end

  return gcd(b, a % b) -- tail call

end



return gcd



end)



register('rational', function(myrequire)

-- poor man's rational number library

local eq = myrequire('eq')

local gcd = myrequire('gcd')

local compat = myrequire('advent.compat')



local Rational = {}

Rational.__index = Rational

function Rational:new(n, d, reduced)

  if d == nil then d = 1 end

  if d < 0 then

    n,d = -n,-d

  elseif d > 0 then

    -- no problem

  else

    error("divide by zero")

  end

  local r = nil

  if reduced then r = true end

  return setmetatable({n=n, d=d, reduced=r}, self)

end

function Rational:reduce()

  if self.reduced then return self end

  -- find gcd of numerator and denominator

  if eq(self.n, 0) then

    self.d = 1

  elseif self.d > 1 then

    local div

    if self.n > 0 then

      div = gcd(self.n, self.d)

    else

      div = gcd(-self.n, self.d)

    end

    if div ~= 1 then

      self.n = compat.idiv(self.n, div)

      self.d = compat.idiv(self.d, div)

    end

  end

  self.reduced = true

  return self

end



local function ensureRational(r)

  if type(r) == 'number' then return Rational:new(r, 1, true) end

  assert(getmetatable(r) == Rational)

  return r

  --[[

  if getmetatable(r) == Rational then return r end

  return Rational:newUnreduced(r, 1)

  ]]--

end



function Rational:isWholeNumber()

  self:reduce()

  return eq(self.d, 1)

end



function Rational:toWholeNumber()

  return compat.idiv(self.n, self.d)

end



function Rational:__tostring()

  self:reduce()

  if self:isWholeNumber() then

    return tostring(self.n)

  end

  return tostring(self.n) .. "/" .. tostring(self.d)

end

function Rational:__unm()

  if eq(self.n, 0) then return self end

  return Rational:new(-self.n, self.d, self.reduced)

end

function Rational.__add(a, b)

  a,b = ensureRational(a), ensureRational(b)

  return Rational:new(a.n * b.d + b.n * a.d, a.d * b.d)

end

function Rational.__sub(a, b)

  a,b = ensureRational(a), ensureRational(b)

  return Rational:new(a.n * b.d - b.n * a.d, a.d * b.d)

end

function Rational.__mul(a, b)

  a,b = ensureRational(a), ensureRational(b)

  return Rational:new(a.n*b.n, a.d*b.d)

end

function Rational.__div(a, b)

  a,b = ensureRational(a), ensureRational(b)

  if type(a.n) ~= 'number' then a.n:normalize() end

  if type(a.d) ~= 'number' then a.d:normalize() end

  if type(b.n) ~= 'number' then b.n:normalize() end

  if type(b.d) ~= 'number' then b.d:normalize() end

  return Rational:new(a.n*b.d, a.d*b.n)

end

function Rational.__lt(a, b)

  a,b = ensureRational(a), ensureRational(b)

  return (a.n * b.d) < (b.n * a.d)

end

function Rational.__le(a, b)

  return not (a > b)

end

function Rational.__eq(a, b)

  a,b = ensureRational(a), ensureRational(b)

  return eq(a.n * b.d, b.n * a.d)

end



return Rational



end)



register('matrix', function(myrequire)

local eq = myrequire('eq')



local Matrix = {}

Matrix.__index = Matrix



function Matrix:new(p)

  return setmetatable(p or {}, self)

end



function Matrix:newNxM(n, m)

  local m = {}

  for i=1,n do

    local row = {}

    for j=1,m do

      table.insert(row, 0)

    end

    table.insert(m, row)

  end

  return self:new(m)

end



-- destructive update

function Matrix:apply(f)

  for i=1,#self do

    for j=1,#selfi do

      selfi][j = f(selfi][j])

    end

  end

  return self

end



-- destructive update to self

function Matrix:LUPDecompose(N, nopivots)

  local P = {}

  for i=1,N do

    Pi = i -- unit permutation matrix

  end

  local S = 0 -- counting pivots

  for i=1,N do

    if nopivots then

      if eq(selfi][i], 0) then

        return nil -- matrix is degenerate

      end

    else

      local maxA = 0

      local imax = i

      for k=i,N do

        local absA = selfk][i

        if absA < 0 then absA = -absA end

        if absA > maxA then

          maxA = absA

          imax = k

        end

      end

      if not (maxA > 0) then

        --print("i", i, "maxA", maxA)

        --self:print()

        return nil

      end -- failure, matrix is degenerate

      if imax ~= i then

        -- pivoting P

        Pi],Pimax = Pimax],Pi

        -- pivoting rows of A

        selfi],selfimax = selfimax],selfi

        -- counting pivots (for determinant)

        S = S + 1

      end

    end



    for j=i+1,N do

      --print(self[i][i])

      assert(not eq(selfi][i], 0))

      selfj][i = selfj][i / selfi][i

      for k=i+1,N do

        selfj][k = selfj][k - selfj][i * selfi][k

      end

    end

    --print("after pivot of",i)

    --self:print()

  end

  return P,S

end



-- destructive update to self

function Matrix:LUPSolve(b, N, nopivots)

  local P,S = self:LUPDecompose(N, nopivots)

  if P == nil then return nil end

  print("Decomposed!", P, S)

  local x = {}

  for i=1,N do

    xi = bPi]]

    for k=1,i-1 do

      xi = xi - selfi][k * xk

    end

  end

  for i=N,1,-1 do

    for k=i+1,N do

      xi = xi - selfi][k * xk

    end

    xi = xi / selfi][i

  end

  return x

end



function Matrix:print()

  local buf = {}

  for i=1,#self do

    for j=1,#selfi do

      table.insert(buf, tostring(selfi][j]))

      table.insert(buf, ' ')

    end

    table.insert(buf, "\n")

  end

  print(table.concat(buf))

end



--[[

local A = Matrix:new{

  { 1, 1, 1 },

  { 4, 3, -1 },

  { 3, 5, 3 }

}

local Rational = require 'rational'

A:apply(function(n) return Rational:new(n) end)

--local P,S = A:LUPDecompose(3)

--A:print()

--print(P,S)

local C = { 1, 6, 4 }

local X = A:LUPSolve(C, 3)

print(X[1], X[2], X[3])

]]--



return Matrix



end)



register('ring', function(myrequire)

-- modular arithmetic

local eq = myrequire('eq')

local compat = myrequire('advent.compat')



local Ring = {}

Ring.__index = Ring

function Ring:new(n, m)

  assert(n >=0 and m > 0)

  assert(n < m)

  return setmetatable({n=n, m=m}, self)

end

function Ring:__tostring()

  return string.format("%s mod %s", tostring(self.n), tostring(self.m))

end

function Ring:__unm()

  if self.n == 0 then return self end

  return Ring:new(self.m-self.n, self.m)

end

function Ring.__add(a, b)

  local r, m

  if type(a) == 'number' then

    m = b.m

    r = (a % m) + b.n

  elseif type(b) == 'number' then

    m = a.m

    r = a.n + (b % m)

  else

    assert(eq(a.m, b.m))

    m = a.m

    r = a.n + b.n

  end

  if r >= m then r = r - m end

  return Ring:new(r, m)

end

function Ring.__sub(a, b)

  return a + (-b)

end

function Ring.__eq(a, b)

  if type(a) == 'number' then

    return eq(a % b.m, b.n)

  elseif type(b) == 'number' then

    return eq(a.n, b % a.m)

  else

    assert(eq(a.m, b.m))

    return eq(a.n, b.n)

  end

end

function Ring.__mul(a, b)

  -- there are better algorithms, but this will do for now

  local r, m

  if type(a) == 'number' then

    m = b.m

    r = ((a % m) * b.n)

  elseif type(b) == 'number' then

    m = a.m

    r = (a.n * (b % m))

  else

    assert(eq(a.m, b.m))

    m = a.m

    r = (a.n * b.n)

  end

  return Ring:new(r % m, m)

end





function Ring.__div(a, b)

  -- compute modular inverse of b; this is only valid if modulus is prime

  local t, newt = 0, 1

  local r, newr = b.m, b.n

  while not eq(newr, 0) do

    local quotient = compat.idiv(r, newr)

    t, newt = newt, t - quotient * newt

    r, newr = newr, r - quotient * newr

  end

  if r > 1 then

    error("divisor is not invertible")

  elseif t < 0 then

    t = t + b.m

  end

  local inverse_b = Ring:new(t, b.m)

  if eq(a.n, 1) then return inverse_b end

  return a * inverse_b

end



function Ring.__idiv(a, b)

  return Ring.__div(a, b)

end



return Ring



end)



register('day24', function(myrequire)

--[[ DAY 24 ]]--

local l = myrequire('llpeg')

local read_wiki_input = myrequire('util').read_wiki_input

local BigNum = myrequire('bignum')

local Rational = myrequire('rational')

local Matrix = myrequire('matrix')

local Ring = myrequire('ring')

local eq = myrequire('eq')

local lt = myrequire('lt')

local gcd = myrequire('gcd')

local compat = myrequire('advent.compat')



local inspect = _G'print' or function() end



local function abs(n)

  if lt(n, 0) then return -n else return n end

end



local function primes(n)

  local result = {}

  local notA = {}

  local i = 2

  while i*i <= n do

    if notAi == nil then

      table.insert(result, i)

      local j = i*i

      while j <= n do

        notAj = true

        j = j + i

      end

    end

    i = i + 1

  end

  while i<=n do

    if notAi == nil then

      table.insert(result, i)

    end

    i = i + 1

  end

  return result

end



local function sortedKeys(tbl)

   local sorted = {}

   for k,_ in pairs(tbl) do

      table.insert(sorted, k)

   end

   table.sort(sorted)

   return sorted

end



--[[ PARSING ]]--



local function tobignum(n) return BigNum:new(tonumber(n)) end



local Hail = {}

Hail.__index = Hail

function make_hail(tbl)

  return setmetatable(tbl, Hail)

end

function Hail:__tostring()

  return string.format("%s,%s,%s @ %s,%s,%s",

                       self.position.x, self.position.y, self.position.z,

                       self.velocity.x, self.velocity.y, self.velocity.z)

end



local nl = l.P"\n"

local SKIP = l.S" \t"^0



local patt = l.P{

   "Hail",

   Hail = l.Ct( l.V"Hailstone" * (nl * l.V"Hailstone")^0 * nl^0) * -1,

   Hailstone = l.Ct( l.Cg(l.V"Coord", "position") * l.P"@" * SKIP * l.Cg(l.V"Coord", "velocity") ) / make_hail,

  Coord = l.Ct( l.Cg(l.V"Number", "x") * SKIP * l.P"," * SKIP *

                l.Cg(l.V"Number", "y") * SKIP * l.P"," * SKIP *

                l.Cg(l.V"Number", "z") * SKIP ),

  Number =  ((l.P"-" ^ -1) * (l.R"09"^1)) / tobignum,

}



local function parse(source)

   local ast, errlabel, pos = patt:match(source)

   if not ast then

      error(string.format("Error at pos %s label '%s'", pos, errlabel))

   end

   return ast

end





--[[ PART 1 ]]--



function determinant(a, b, c, d)

   return a*d - b*c

end



function sign(x)

   if lt(x, 0) then return -1 end

   if lt(0, x) then return  1 end

   return 0

end



function hailstone_intersection(a, b, min, max)

   local x1, x2 = a.position.x, a.position.x + a.velocity.x

   local y1, y2 = a.position.y, a.position.y + a.velocity.y

   local x3, x4 = b.position.x, b.position.x + b.velocity.x

   local y3, y4 = b.position.y, b.position.y + b.velocity.y

   local d = determinant(x1-x2, x3-x4, y1-y2, y3-y4)

   if d == 0 then return false end -- no intersection!

   local td = determinant(x1-x3, x3-x4, y1-y3, y3-y4)

   if sign(td) ~= sign(d) then return false end -- intersection in past

   local ud = determinant(x1-x3, x1-x2, y1-y3, y1-y2)

   if sign(ud) ~= sign(d) then return false end -- intersection in past

   local Px = a.position.x + compat.idiv(a.velocity.x * td, d)

   --print(Px)

   if lt(Px, min) or lt(max, Px) then return false end -- intersection not in range

   if eq(Px, max) then print("warning x") end

   local Py = a.position.y + compat.idiv(a.velocity.y * td, d)

   --print(Py)

   if lt(Py, min) or lt(max, Py) then return false end -- intersection not in race

   if eq(Py, max) then print("warning y") end

   return true

end



local function part1(source, min, max)

   min,max = tonumber(min),tonumber(max)

   local hail = parse(source)

   local count = 0

   for i=1,#hail do

      for j = i+1, #hail do

         if hailstone_intersection(haili], hailj], min, max) then

            --print("Intersection:", inspect(hail[i]), inspect(hail[j]))

            count = count + 1

         end

      end

   end

   return count

end



--[[ PART 2 ]]--

local Range = {}

Range.__index = Range

function Range:new()

  return setmetatable({}, self)

end

function Range:le(val)

  if self.max == nil or self.max > val then self.max = val end

end

function Range:ge(val)

  if self.min == nil or self.min < val then self.min = val end

end

function Range:__tostring(val)

  local buf = { '[' }

  if self.min == nil then

    table.insert(buf, "inf")

  else

    table.insert(buf, tostring(self.min))

  end

  table.insert(buf, ",")

  if self.max == nil then

    table.insert(buf, "inf")

  else

    table.insert(buf, tostring(self.max))

  end

  table.insert(buf, ']')

  return table.concat(buf)

end



local function check_bounds(hail, vPositive)

  local unk = { x=Range:new(), y=Range:new(), z=Range:new() }

  local coords = { "x", "y", "z" }

  for _,h in ipairs(hail) do

    --print(h.velocity.x, h.velocity.y, h.velocity.z)

    -- h.position + t*h.velocity = unk.position + u*unk.velocity

    -- since we know t and u have to be positive, if we know the sign of

    -- u (which we'll brute force) we can constrain unk.position by h.position

    for _,c in ipairs(coords) do

      if vPositivec then -- unk.position >= h.position

        if h.velocityc >= 0 then

          -- both velocities are positive, nothing more can be said

        else

          -- hail is moving - while unk is moving +

          -- thus h.position >= unk.position

          unkc]:le(h.positionc])

        end

      else

        if h.velocityc >= 0 then

          -- hail is moving + while unk is moving -

          -- thus h.position <= unk.position

          unkc]:ge(h.positionc])

        else

          -- both velocities are negative, nothing more can be said

        end

      end

    end

  end

  local buf = {"For "}

  for _,c in ipairs(coords) do

    table.insert(buf, c)

    if vPositivec then

      table.insert(buf, "+")

    else

      table.insert(buf, "-")

    end

    if c ~= 'z' then table.insert(buf, ", ") end

  end

  table.insert(buf, " limits are ")

  for _,c in ipairs(coords) do

    table.insert(buf, c)

    table.insert(buf, "=")

    table.insert(buf, tostring(unkc]))

    if c ~= 'z' then table.insert(buf, ", ") end

  end

  print(table.concat(buf))

  return unk

end



local function apply(tbl, f)

  for i,v in ipairs(tbl) do

    tbli = f(v)

  end

  return tbl

end



local function mkrat(n) return Rational:new(n, tobignum(1)) end

local function mkrat3(v)

  return {x=mkrat(v.x), y=mkrat(v.y), z=mkrat(v.z)}

end

local function bigrat(n,m) return Rational:new(tobignum(n),tobignum(m or 1)) end



function check8(hail)

  -- use x=az+b / y = cz + d form of each line

  -- z = pos.z + t*vel.z => t = z/vel.z - pos.z/vel.z

  -- x = pos.x + t*vel.x => x = pos.x + z(vel.x/vel.z) - (vel.x/vel.z)*pos.z

  --                        x = (vel.x/vel.z)*z + (pos.x - (vel.x/vel.z)*pos.z)

  local function abcd(pos, vel)

    local a = mkrat(vel.x) / mkrat(vel.z)

    local b = mkrat(pos.x) - (mkrat(pos.z) * mkrat(vel.x) / mkrat(vel.z))

    local c = mkrat(vel.y) / mkrat(vel.z)

    local d = mkrat(pos.y) - (mkrat(pos.z) * mkrat(vel.y) / mkrat(vel.z))

    return a,b,c,d

  end

  local a1,b1,c1,d1 = abcd(hail1].position, hail1].velocity)

  local a2,b2,c2,d2 = abcd(hail2].position, hail2].velocity)

  local a3,b3,c3,d3 = abcd(hail3].position, hail3].velocity)

  local a4,b4,c4,d4 = abcd(hail4].position, hail4].velocity)

  --[[

    a1*z1 + b1 = unk.a * z1 + unk.b

    c1*z1 + d1 = unk.c * z1 + unk.d



    unk.a*z1 - a1*z1 + unk.b - b1 = 0

    unk.c*z1 - c1*z1 + unk.d - d1 = 0

  ]]--

end



function gradient_descent(hail)

  -- use x=az+b / y=cz+d form

  -- vel.x = a, vel.y = c, vel.z = 1

  -- pos.x = b, pos.y = d. pos.z = 0



  local function crossprod(a, b)

    return {

      x=a.y*b.z - a.z*b.y,

      y=a.z*b.x - a.x*b.z,

      z=a.x*b.y - a.y*b.x,

    }

  end

  local function dist2(hailstone, a, b, c, d)

    local unk_pos = { x=b, y=d, z=mkrat(0) }

    local unk_vel = { x=a, y=c, z=mkrat(1) }

    local hail_pos = mkrat3(hailstone.position)

    local hail_vel = mkrat3(hailstone.velocity)

    local n = crossprod(hail_vel, unk_vel)

    local r1mr2 = {

      x=hail_pos.x - unk_pos.x,

      y=hail_pos.y - unk_pos.y,

      z=hail_pos.z - unk_pos.z,

    }

    local dd = n.x * r1mr2.x + n.y * r1mr2.y + n.z * r1mr2.z

    local n_mag2 = n.x * n.x + n.y * n.y + n.z * n.z

    return (dd * dd) / n_mag2

  end



  -- arbitrarily selected starting position for search

  -- as hail[1] presumably doesn't have z velocity 1 or pos.z = 0

  -- this is *near* but not *actually* the same as hail[1]'s vector

  local start = {

    --[[

    a=hail[1].velocity.x, b=hail[1].position.x,

    c=hail[1].velocity.y, d=hail[1].position.y,

    ]]--

    --a=mkrat(-18), b=mkrat(206273907287337), c=mkrat(-17), d=mkrat(404536114336383),

--    a=bigrat(-1,10), b=bigrat(6893814583987912, 25), c=bigrat(-1,5), d=bigrat(15882393318613117,50)

    --    a=bigrat(-107,1000), b=bigrat(27575205905929417,100), c=bigrat(-37,250), d=bigrat(79411872811990497,250)

    -- a=bigrat(-107,1000), b=bigrat(27616905905929417,100), c=bigrat(-147,1000), d=bigrat(79307872811990497,250),

    -- a=bigrat(-541,5000), b=bigrat(27657105905929417,100), c=bigrat(-361,2500), d=bigrat(79207622811990497,250),

    --a=bigrat(-129,1000), b=bigrat(279241059059294,1), c=bigrat(-61,625), d=bigrat(309020491247961, 1)

--    a=bigrat(-1593,10000), b=bigrat(283241059059294,1), c=bigrat(-333,10000),d=bigrat(286020491247961,1),

    --a=bigrat(-1747,10000),b=bigrat(284841059059294,1),c=bigrat(73,10000),d=bigrat(278720491247961,1),

    --a=bigrat(-1857,10000),b=bigrat(286421059059294,1),c=bigrat(383,10000),d=bigrat(273340491247961,1),

    --a=bigrat(-1857,10000),b=bigrat(286419759059294,1),c=bigrat(383,10000),d=bigrat(273344691247961,1),

    --a=bigrat(-1857,10000),b=bigrat(286419759059294,1),c=bigrat(383,10000),d=bigrat(273344698247961,1),

    --a=bigrat(-1857,10000),b=bigrat(286419758729294,1),c=bigrat(383,10000),d=bigrat(273344698467961,1),

--a=bigrat(-1857,10000),b=bigrat(286419758728794,1),c=bigrat(383,10000),d=bigrat(273344698470961,1),

--a=bigrat(-1857,10000),b=bigrat(286419758728797,1),c=bigrat(383,10000),d=bigrat(273344698470938,1),

--a=bigrat(-9289,50000),b=bigrat(286419758728717,1),c=bigrat(1919,50000),d=bigrat(273344698470858,1),

--a=bigrat(-18763,100000),b=bigrat(286419718728717,1),c=bigrat(4237,100000),d=bigrat(273344658570858,1),

--a=bigrat(-241,1250),b=bigrat(286419463728717,1),c=bigrat(143,2500),d=bigrat(273344403570858,1),

--a=bigrat(-9693,50000),b=bigrat(286416553728717,1),c=bigrat(601,10000),d=bigrat(273341503570858,1),

--a=bigrat(-19439,100000),b=bigrat(286265553728717,1),c=bigrat(77,1250),d=bigrat(273191503570858,1),

--a=bigrat(-19439,100000),b=bigrat(285835553728717,1),c=bigrat(77,1250),d=bigrat(272001503570858,1),

--a=bigrat(-1949,10000),b=bigrat(285945553728717,1),c=bigrat(3149,50000),d=bigrat(268791503570858,1),

--a=bigrat(-39029,200000),b=bigrat(285795553728717,1),c=bigrat(63691,1000000),d=bigrat(268652503570858,1),



--a=bigrat(-39029,200000),b=bigrat(280000000000000,1),c=bigrat(63691,1000000),d=bigrat(260000000000000,1),

--a=bigrat(-122,625),b=bigrat(285200000000000,1),c=bigrat(319,5000),d=bigrat(268500000000000,1),

--a=bigrat(-19519,100000),b=bigrat(285830000000000,1),c=bigrat(6383,100000),d=bigrat(268630000000000,1),

--a=bigrat(-48823,250000),b=bigrat(285795000000000,1),c=bigrat(32067,500000),d=bigrat(268575000000000,1),

--a=bigrat(-97701,500000),b=bigrat(285794400000000,1),c=bigrat(64463,1000000),d=bigrat(268545000000000,1),

--a=bigrat(-3054,15625),b=bigrat(285846500000000,1),c=bigrat(16157,250000),d=bigrat(268488400000000,1),



a=bigrat(3054,15625),b=bigrat(285846500000000,1),c=bigrat(16157,250000),d=bigrat(268488400000000,1),

--[[

  206273907288897 - 18t = x - 97701/500000u

  404536114337943 + 6t  = y + 64463/1000000u

  197510451330134 + 92t = z + u



]]--

  }

  local vars = { "a", "b", "c", "d" }

  local function score(guess)

    local sum = mkrat(0)

    for i=1,4 do

      sum = sum + dist2(haili], guess.a, guess.b, guess.c, guess.d)

    end

    return sum

  end

  local function quantize(n, amt)

    local m = (n*amt):toWholeNumber()

    return Rational:new(m, amt)

  end



  local guess = start

  local dim = 1

  local epsilon = bigrat(1, 1000000)

  local beta = bigrat(1,5)

  while true do

    local s = score(guess)

    --print(s:toWholeNumber(), guess.a, guess.b, guess.c, guess.d)

    local buf = {}

    for _,v in ipairs(vars) do

      table.insert(buf, v) ; table.insert(buf, "=bigrat(")

      guessv]:reduce()

      table.insert(buf, tostring(guessv].n))

      table.insert(buf, ",")

      table.insert(buf, tostring(guessv].d))

      table.insert(buf, "),")

    end

    print(s:toWholeNumber(), table.concat(buf))

    if not (s > 0 or s < 0) then break end

    local orig = guessvarsdim]]

    guessvarsdim]] = orig + epsilon

    local s2 = score(guess)

    guessvarsdim]] = orig - epsilon

    local s3 = score(guess)

    guessvarsdim]] = orig

    --print("Score 2", vars[dim], s2:toWholeNumber())

    --print("Score 3", vars[dim], s3:toWholeNumber())

    -- compute the derivative in this dimension

    local adjustment = 0

    local deriv = 0

    if s2 < s3 then

      if s2 < s then adjustment = mkrat(1) ; deriv = s - s2; end

    else

      if s3 < s then adjustment = mkrat(-1) ; deriv = s - s3; end

    end

    local q

    if dim == 2 or dim == 4 then

      --adjustment = adjustment*guess[vars[dim]]/1000

      --adjustment = adjustment * guess[vars[dim]]/100000

      --adjustment = adjustment * bigrat(1000000, 1)

      --adjustment = adjustment * s / 100*(deriv / epsilon)

      adjustment = adjustment * 100000000

      q = 1

    else

      q = 1000000

      adjustment = adjustment / q

    end

      --[[

    local deriv = (s2 - s)

    local adjustment = (-s * beta / deriv):toWholeNumber()

    --print("Score 1", s:toWholeNumber())

    --print("Score 2", s2:toWholeNumber())

    --print("Derivative", vars[dim], deriv:toWholeNumber())

    if not (adjustment<0 or adjustment>0) then

      if deriv > 0 then adjustment = -1 else adjustment = 1 end

    end

        --print("Adjusting",vars[dim],"by",adjustment)

      ]]--

    guessvarsdim]] = quantize(guessvarsdim]] + adjustment, q)

    local s3 = score(guess)

    --print("Score 3", s3:toWholeNumber())

    if s3 > s then -- that didn't work, undo

      guessvarsdim]] = orig

    end

    dim = dim + 1

    if dim > #vars then dim = 1 end

  end

  print("Found it!", guess.a, guess.b, guess.c, guess.d)

end



function check(hail, vx, vy, vz)

  -- use first 3 hailstones to constrain the three+(2*n) unknowns

  --[[

    [1 0 0 -h1.vx vx 0 0 0 0][unk.x ]   [h1.x]

    [0 1 0 -h1.vy vy 0 0 0 0][unk.y ]   [h1.y]

    [0 0 1 -h1.vz vz 0 0 0 0][unk.z ]   [h1.x]

    [1 0 0 0 0 -h2.vx vx 0 0][t1    ]   [h2.x]

    [0 1 0 0 0 -h2.vy vy 0 0][u1    ]   [h2.y]

    [0 0 1 0 0 -h2.vz vz 0 0][t2    ]   [h2.z]

    [1 0 0 0 0 0 0 -h3.vx vx][u2    ] = [h3.x]

    [0 1 0 0 0 0 0 -h3.vy vy][t3    ] = [h3.y]

    [0 0 1 0 0 0 0 -h3.vz vz][u3    ] = [h3.z]

  ]]--

  local M = Matrix:new{

    {1, 0, 0, -hail1].velocity.x, vx, 0, 0, 0, 0, },

    {0, 1, 0, -hail1].velocity.y, vy, 0, 0, 0, 0, },

    {0, 0, 1, -hail1].velocity.z, vz, 0, 0, 0, 0, },

    {1, 0, 0, 0, 0, -hail2].velocity.x, vx, 0, 0, },

    {0, 1, 0, 0, 0, -hail2].velocity.y, vy, 0, 0, },

    {0, 0, 1, 0, 0, -hail2].velocity.z, vz, 0, 0, },

    {1, 0, 0, 0, 0, 0, 0, -hail3].velocity.x, vx, },

    {0, 1, 0, 0, 0, 0, 0, -hail3].velocity.y, vy, },

    {0, 0, 1, 0, 0, 0, 0, -hail3].velocity.z, vz, },

  }

  --print("originally")

  --M:print()

  M:apply(function(n)

      if type(n)=='number' then return bigrat(n,1) end

      if getmetatable(n)==BigNum then return Rational:new(n, tobignum(1)) end

      return n

  end)

  local b = {

    hail1].position.x,

    hail1].position.y,

    hail1].position.z,

    hail2].position.x,

    hail2].position.y,

    hail2].position.z,

    hail3].position.x,

    hail3].position.y,

    hail3].position.z,

  }

  apply(b, mkrat)

  local x = M:LUPSolve(b, 9)

  print(x)

  if x == nil then return false end

  print("Solved!", inspect(x))

  -- check for non-integer values in solution array

  for i=1,9 do

    if not xi]:isWholeNumber() then return false end

  end

  -- check for negative time values

  for i=4,9 do

    if xi < 0 then return false end

  end

  -- okay, now verify solutions for all other hailstones

  local unk = {

    x=x1]:toWholeNumber(),

    y=x2]:toWholeNumber(),

    z=x3]:toWholeNumber(),

  }

  for i=4, #hail do

    -- 2 unknowns, 2 equations

    -- -h.vx*t + unk.vx*u = h.x - unk.x

    local M2 = Matrix:new{

      { -haili].velocity.x, vx },

      { -haili].velocity.y, vy },

    }

    M2:apply(mkrat)

    local b = {

      haili].position.x - unk.x,

      haili].position.y - unk.y,

    }

    apply(b, mkrat)

    local x = M2:LUPSolve(b, 2)

    if x == nil then return false end

    -- check for non-integer or non-positive values

    for i=1,2 do

      if xi < 0 or not xi]:isWholeNumber() then return false end

    end

    -- check that t/u values also work for z axis

    local t,u = x1]:toWholeNumber(), x2]:toWholeNumber()

    local hz = haili].position.z + t*haili].velocity.z

    local uz = unk.z + u*vz

    if hz < uz then return false end

    if hz > uz then return false end

  end

  -- found a solution!

  print(vz,vy,vz, "works!")

  return true

end



local function check_all_bounds(source)

  local hail = parse(source)

  for xsign=0,1 do

    for ysign=0,1 do

      for zsign=0,1 do

        check_bounds(hail, {x=(xsign==1), y=(ysign==1), z=(zsign==1)})

      end

    end

  end

end



local function part2_brute(source)

  local hail = parse(source)

  assert(not check(hail, 0, 0, 0))

  -- brute force search through small vx/vy/vz

  for dist=72,1000 do

    print("dist=",dist)

    local seen = {} -- hack to avoid checking edges twice

    local function mycheck(x,y,z)

      if x==dist or x==-dist or y==dist or y==-dist or z==dist or z==-dist then

        local key = string.format("%d,%d,%d", x, y, z)

        if seenkey then return false end -- already checked

        seenkey = true

      end

      return check(hail, x, y, z)

    end

    for i=-dist, dist do

      for j=-dist, dist do

        if mycheck( dist,i,j) then return end

        if mycheck(-dist,i,j) then return end

        if mycheck(i, dist,j) then return end

        if mycheck(i,-dist,j) then return end

        if mycheck(i,j, dist) then return end

        if mycheck(i,j,-dist) then return end

      end

    end

  end

end



local function part2_grad(source)

  local hail = parse(source)

  gradient_descent(hail)

end



local function part2_another_attempt(source)

  local hail = parse(source)

  --check(hail, bigrat(-97701,500000), bigrat(64463,1000000), bigrat(1,1))

  local function abcd(pos, vel)

    local a = mkrat(vel.x) / mkrat(vel.z)

    local b = mkrat(pos.x) - (mkrat(pos.z) * mkrat(vel.x) / mkrat(vel.z))

    local c = mkrat(vel.y) / mkrat(vel.z)

    local d = mkrat(pos.y) - (mkrat(pos.z) * mkrat(vel.y) / mkrat(vel.z))

    return a,b,c,d

  end

  local unk_a = bigrat(-3054,15625)

  local unk_c = bigrat(16157,250000)

  local a1,b1,c1,d1 = abcd(hail1].position, hail1].velocity)

  local a2,b2,c2,d2 = abcd(hail2].position, hail2].velocity)

  --[[

    a1*z1 + b1 = unk.a * z1 + unk.b

    c1*z1 + d1 = unk.c * z1 + unk.d



    unk.a*z1 - a1*z1 + unk.b - b1 = 0

    unk.c*z1 - c1*z1 + unk.d - d1 = 0



    [unk.a-a1 0 1 0][z1]      [ b1 ]

    [unk.c-c1 0 0 1][z2]      [ d1 ]

    [0 unk.a-a2 1 0][unk.b] = [b2]

    [0 unk.c-c2 0 1][unk.d]   [d2]



                             |        0 0 1 |              |       0 1 0 |

    determinant = (unk.a-a1)*| unk.a-a2 1 0 | - (unk.c-c1)*|unk.a-a2 1 0 |

                             | unk.c-c2 0 1 |              |unk.c-c2 0 1 |



    = (unk.a-a1)*(c2-unk.c) + (c1-unk.c)*(a2-unk.a)

  ]]--

  print("unk_a",unk_a) print("unk_c",unk_c)

  print("a1",a1) print("c1", c1)

  print("a2",a2) print("c2", c2)

  print((unk_a-a1)*(c2-unk_c))

  print((c1-unk_c)*(a2-unk_a))

  local M2 = Matrix:new{

    { (unk_a - a1), bigrat(0), bigrat(1), bigrat(0) },

    { (unk_c - c1), bigrat(0), bigrat(0), bigrat(1) },

    { bigrat(0), (unk_a - a1), bigrat(1), bigrat(0) },

    { bigrat(0), (unk_c - c1), bigrat(0), bigrat(1) },

  }

  local b = { b1, d1, b2, d2 }

  local x = M2:LUPSolve(b, 4)

  print(x)

end



local function part2_residue(source)

  local coords = { "x", "y", "z" }

  local max = 0

  local p = primes(689) -- maximum common factor is 689



  local vResults = {}

  local pResults = {}

  local p2Results = {}

  local hail = parse(source)

  local HACK = 2

  -- for all pairs of hailstones:

  for i=1,#hail do

    for j=i+1,#hail do

      local hail1,hail2 = haili], hailj

      -- for all pairs of coords:

      for k=1,#coords do

        for l=k+1,#coords do

          local d1,d2 = coordsk], coordsl

          -- for all prime factors in common between the velocity vectors

          local common1 = gcd(abs(hail1.velocityd1]), abs(hail1.velocityd2]))

          local common2 = gcd(abs(hail2.velocityd1]), abs(hail2.velocityd2]))

          local common = gcd(common1, common2)

          if common > max then max = common end -- track max common factor

          if not eq(common, 1) then

            for _,pp in ipairs(p) do

              --pp = pp * pp * pp * pp * pp -- HACK

              if pp*pp > common then break end

              if eq(common % pp, 0) then

                -- check for usable residue!

                local pos1d1 = hail1.positiond1 % pp

                local pos2d1 = hail2.positiond1 % pp

                local pos1d2 = hail1.positiond2 % pp

                local pos2d2 = hail2.positiond2 % pp

                if eq(pos1d1, pos2d1) and not eq(pos1d2, pos2d2) then

                  --print("vel", d1,"=0 mod", pp)

                  local key = string.format("v%s=0mod%s", d1, pp)

                  vResultskey = { coord=d1, mod=pp }

                  --print("pos", d1,"=", pos1d1, "mod", pp)

                  key = string.format("p%s=%smod%s", d1, pos1d1, pp)

                  pResultskey = { coord=d1, rem=pos1d1, mod=pp }

                  if d1=='x' and eq(pp, 5) then

                    print(string.format("Looking at %s and %s", d1, d2))

                    print(hail1)

                    print(hail2)

                    print("pos", d1,"=", pos1d1, "mod", pp)

                  end

                end

                if eq(pos1d2, pos2d2) and not eq(pos1d1, pos2d1) then

                  --print("vel", d2,"=0 mod", pp)

                  local key = string.format("v%s=0mod%s", d2, pp)

                  vResultskey = { coord=d2, mod=pp }

                  --print("pos", d2,"=", pos1d2, "mod", pp)

                  key = string.format("p%s=%smod%s", d2, pos1d2, pp)

                  pResultskey = { coord=d2, rem=pos1d2, mod=pp }

                end

              end

            end

          end

        end

      end

    end

  end

  print("Largest common factor is ", max)

  local v = { x=1, y=1, z=1 }

  for _,k in ipairs(sortedKeys(vResults)) do

    local r = vResultsk

    print(string.format("Velocity %s = 0 mod %s", r.coord, r.mod))

    vr.coord = vr.coord * r.mod

  end

  print("So:")

  for _,c in ipairs(coords) do

    print(string.format("Velocity %s = <some constant> * %s", c, vc]))

  end

  for _,k in ipairs(sortedKeys(pResults)) do

    local r = pResultsk

    print(string.format("Position %s = %s mod %s", r.coord, r.rem, r.mod))

  end

  return v.x, v.y, v.z

end



--[[

Largest common factor is 	44

Velocity x = 0 mod 3

Velocity y = 0 mod 2

Velocity y = 0 mod 3

Velocity z = 0 mod 4

So:

Velocity x = <some constant> * 3

Velocity y = <some constant> * 6

Velocity z = <some constant> * 2

Position x = 0 mod 3

Position y = 0 mod 3

Position y = 1 mod 2

Position z = 2 mod 4

]]--



local function part2_search2(source)

  local hail = parse(source)

  --check(hail, 3, 6, 4)



  local MOD = 5

  for vx=0,4 do

    for vy=0,4 do

      for vz=0,4 do



  local hail1, hail2, hail3 = hail1], hail2], hail3



  local M = Matrix:new{

    {1, 0, 0, -hail1.velocity.x, vx, 0, 0, 0, 0, },

    {0, 1, 0, -hail1.velocity.y, vy, 0, 0, 0, 0, },

    {0, 0, 1, -hail1.velocity.z, vz, 0, 0, 0, 0, },

    {1, 0, 0, 0, 0, -hail2.velocity.x, vx, 0, 0, },

    {0, 1, 0, 0, 0, -hail2.velocity.y, vy, 0, 0, },

    {0, 0, 1, 0, 0, -hail2.velocity.z, vz, 0, 0, },

    {1, 0, 0, 0, 0, 0, 0, -hail3.velocity.x, vx, },

    {0, 1, 0, 0, 0, 0, 0, -hail3.velocity.y, vy, },

    {0, 0, 1, 0, 0, 0, 0, -hail3.velocity.z, vz, },

  }

  local function make_ring(n)

    if type(n)=='number' then

      n = n % MOD

      return Ring:new(n, MOD)

    end

    if getmetatable(n)==BigNum then

      n = n % MOD

      n = n:toNumber()

      return Ring:new(n, MOD)

    end

    error("what is this?")

  end

  M:apply(make_ring)

  --M:print()

  local b = {

    hail1.position.x,

    hail1.position.y,

    hail1.position.z,

    hail2.position.x,

    hail2.position.y,

    hail2.position.z,

    hail3.position.x,

    hail3.position.y,

    hail3.position.z,

  }

  apply(b, make_ring)

  local x = M:LUPSolve(b, 9, true)

  print(x)

  if x ~= nil then

    print("Solved!", inspect(x))

  end

      end

    end

  end

end



local function part2_solve3(source)

  local hail = parse(source)

  --[[

    h[i].pos.x + t[i]*h[i].vel.x = unk.pos.x + t[i]*unk.vel.x

    =>

    t[i]*(h[i].vel.x - unk.vel.x) = unk.pos.x - h[i].pos.x

    =>

    t[i] = (unk.pos.x - h[i].pos.x) / (h[i].vel.x - unk.vel.x)

    ...now equate x and y for the same hailstone (aka, same t[i])...

    (unk.pos.x - h[i].pos.x)/(h[i].vel.x - unk.vel.x) =

      (unk.pos.y - h[i].pos.y)/(h[i].vel.y - unk.vel.y)

    =>

    (unk.pos.x - h[i].pos.x) * (h[i].vel.y - unk.vel.y) =

      (unk.pos.y - h[i].pos.y) * (h[i].vel.x - unk.vel.x)

    =>

    u.p.x * h.v.y - u.p.x * u.v.y - h.p.x * h.v.y + h.p.x * u.v.y =

       u.p.y * h.v.x - u.p.y * u.v.x - h.p.y * h.v.x + h.p.y * u.v.x

    =>

    u.p.y * u.v.x - u.p.x * u.v.y =

      -u.p.x * h.v.y + h.p.x * h.v.y - h.p.x * u.v.y + u.p.y * h.v.x - h.p.y * h.v.x + h.p.y * u.v.x

    = -h.v.y * u.p.x + h.v.x * u.p.y + h.p.y * u.v.x - h.p.x * u.v.y  + (h.p.x * h.v.y - h.p.y * h.v.x)



    [ (h2.v.y-h1.v.y) (h1.v.x-h2.v.x) (h1.p.y-h2.p.y) (h2.p.x-h1.p.x) ][u.p.x]

    [ (h3.v.y-h1.v.y) (h1.v.x-h3.v.x) (h1.p.y-h3.p.y) (h3.p.x-h1.p.x) ][u.p.y]

    [ (h4.v.y-h1.v.y) (h1.v.x-h4.v.x) (h1.p.y-h4.p.y) (h4.p.x-h1.p.x) ][u.v.x]

    [ (h5.v.y-h1.v.y) (h1.v.x-h5.v.x) (h1.p.y-h5.p.y) (h5.p.x-h1.p.x) ][u.v.y]



    =[ h2.p.x * h2.v.y - h2.p.y * h2.v.x - h1.p.x * h1.v.y + h1.p.y * h1.v.x]

    =[ h3.p.x * h3.v.y - h3.p.y * h3.v.x - h1.p.x * h1.v.y + h1.p.y * h1.v.x]

    =[ h4.p.x * h4.v.y - h4.p.y * h4.v.x - h1.p.x * h1.v.y + h1.p.y * h1.v.x]

    =[ h5.p.x * h5.v.y - h5.p.y * h5.v.x - h1.p.x * h1.v.y + h1.p.y * h1.v.x]

  ]]--



  local ans = { position={}, velocity={} }

  for _,d in ipairs{ "y", "z" } do

    local h1 = hail1

    local rows = {}

    local rhs = {}

    for i=2,5 do

      local h2 = haili

      table.insert(rows, {

                     (h2.velocityd - h1.velocityd]),

                     (h1.velocity.x - h2.velocity.x),

                     (h1.positiond - h2.positiond]),

                     (h2.position.x - h1.position.x),

      })

      table.insert(rhs,

                   h2.position.x * h2.velocityd -

                   h2.positiond * h2.velocity.x -

                   h1.position.x * h1.velocityd +

                   h1.positiond * h1.velocity.x)

    end

    local M = Matrix:new(rows)

    M:apply(mkrat)

    apply(rhs, mkrat)

    local x = M:LUPSolve(rhs, 4)

    if x ~= nil then

      --print("Solved!")

      --print("pos=", x[1], d, x[2])

      --print("vel=", x[3], d, x[4])

      ans.position.x = x1

      ans.positiond = x2

      ans.velocity.x = x3

      ans.velocityd = x4

    end

  end

  print("Position", ans.position.x, ans.position.y, ans.position.z)

  print("Velocity", ans.velocity.x, ans.velocity.y, ans.velocity.z)

  return ans.position.x + ans.position.y + ans.position.z

end



local function part2_manual(input)

   return "Calculation done with paper and pencil"

end



local part2 = part2_manual



--[[ CLI ] ]--

local do_part1 = true

local use_example = true



local filename

if use_example then

  filename = "day24.example"

else

  filename = "day24.input"

end

local source = io.input(filename):read("*a")

if do_part1 then

  -- part 1

  if use_example then

    print('Intersecting hailstones:', part1(source, 7, 27))

  else

    print('Intersecting hailstones:', part1(source, 200000000000000, 400000000000000))

  end

else

  -- part 2

  print("Sum:", part2(source))

end

--[ [ END CLI ]]--



return {

  part1 = read_wiki_input(part1),

  part2 = read_wiki_input(part2),

}



end)



local modules = {}

modules'bit32' = require('bit32')

modules'string' = require('string')

modules'strict' = {}

modules'table' = require('table')

local function myrequire(name)

  if modulesname == nil then

    modulesname = true

    modulesname = (buildersname])(myrequire)

  end

  return modulesname

end

return myrequire('day24')

end)()