Module:User:Cscott/Advent Of Code 2023/Day 24

return (function local builders = {} local function register(name, f) builders[name] = 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.args[1])      local source = title:getContent      if source == nil then        error("Can't find title " .. tostring(title))     end      source = source:gsub("^%s*]*>\n?", "", 1)      source = source:gsub("]*>%s*$", "", 1)      return func(source, frame.args[2], frame.args[3])    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 = self[i] 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 + (self[i] * m)    m = m * RADIX end return sum end function BigNum:normalize local i = 1 local sawNonZero = false while self[i] ~= nil do   assert(self[i] >= 0) if self[i] > 0 then sawNonZero = true end if self[i] >= 1000 then local carry = math.floor(self[i] / 1000) self[i] = self[i] % 1000 self[i+1] = (self[i+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   r[i] = self[i] 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) r[1] = (r[1] or 0) + b else assert(not b.negative) for i=1,#b do     r[i] = (r[i] or 0) + b[i] 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 (a[i] or 0) < (b[i] or 0) then return true end if (a[i] or 0) > (b[i] 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 (a[i] or 0) ~= (b[i] 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   r[i] = (r[i] or 0) - (b[i] or 0) - borrow borrow = 0 while r[i] < 0 do     r[i] = r[i] + RADIX borrow = borrow + 1 end assert(r[i] >= 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     r[i] = a[i] * 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(a[i] >= 0) assert(b[j] >= 0) local prod = a[i] * b[j] r[i+j-1] = (r[i+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 a[1] ~= nil do   table.insert(bits, a[1] % 2) for i=1,#a do     a[i] = compat.idiv(a[i], 2) + ((a[i+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 R[i] = R[i] * 2 end if Nbits[i] > 0 then R[1] = R[1] + 1 end R:normalize for i=1,#Q do Q[i] = Q[i] * 2 end if R >= D then R = R - D     Q[1] = Q[1] + 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,#self[i] do      self[i][j] = f(self[i][j]) end end return self end

-- destructive update to self function Matrix:LUPDecompose(N, nopivots) local P = {} for i=1,N do   P[i] = i -- unit permutation matrix end local S = 0 -- counting pivots for i=1,N do   if nopivots then if eq(self[i][i], 0) then return nil -- matrix is degenerate end else local maxA = 0 local imax = i     for k=i,N do        local absA = self[k][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       P[i],P[imax] = P[imax],P[i] -- pivoting rows of A       self[i],self[imax] = self[imax],self[i] -- counting pivots (for determinant) S = S + 1 end end

for j=i+1,N do     --print(self[i][i]) assert(not eq(self[i][i], 0)) self[j][i] = self[j][i] / self[i][i] for k=i+1,N do       self[j][k] = self[j][k] - self[j][i] * self[i][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   x[i] = b[P[i]] for k=1,i-1 do     x[i] = x[i] - self[i][k] * x[k] end end for i=N,1,-1 do   for k=i+1,N do      x[i] = x[i] - self[i][k] * x[k] end x[i] = x[i] / self[i][i] end return x end

function Matrix:print local buf = {} for i=1,#self do   for j=1,#self[i] do      table.insert(buf, tostring(self[i][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 notA[i] == nil then table.insert(result, i)     local j = i*i while j <= n do       notA[j] = true j = j + i     end end i = i + 1 end while i<=n do   if notA[i] == 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(hail[i], hail[j], 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 vPositive[c] then -- unk.position >= h.position if h.velocity[c] >= 0 then -- both velocities are positive, nothing more can be said else -- hail is moving - while unk is moving + -- thus h.position >= unk.position unk[c]:le(h.position[c]) end else if h.velocity[c] >= 0 then -- hail is moving + while unk is moving - -- thus h.position <= unk.position unk[c]:ge(h.position[c]) 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 vPositive[c] 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(unk[c])) 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    tbl[i] = 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(hail[1].position, hail[1].velocity) local a2,b2,c2,d2 = abcd(hail[2].position, hail[2].velocity) local a3,b3,c3,d3 = abcd(hail[3].position, hail[3].velocity) local a4,b4,c4,d4 = abcd(hail[4].position, hail[4].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(hail[i], 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(")      guess[v]:reduce      table.insert(buf, tostring(guess[v].n))      table.insert(buf, ",")      table.insert(buf, tostring(guess[v].d))      table.insert(buf, "),") end print(s:toWholeNumber, table.concat(buf)) if not (s > 0 or s < 0) then break end local orig = guess[vars[dim]] guess[vars[dim]] = orig + epsilon local s2 = score(guess) guess[vars[dim]] = orig - epsilon local s3 = score(guess) guess[vars[dim]] = 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) then      if deriv > 0 then adjustment = -1 else adjustment = 1 end    end        --print("Adjusting",vars[dim],"by",adjustment)      -- guess[vars[dim]] = quantize(guess[vars[dim]] + adjustment, q)   local s3 = score(guess) --print("Score 3", s3:toWholeNumber) if s3 > s then -- that didn't work, undo guess[vars[dim]] = 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, -hail[1].velocity.x, vx, 0, 0, 0, 0, }, {0, 1, 0, -hail[1].velocity.y, vy, 0, 0, 0, 0, }, {0, 0, 1, -hail[1].velocity.z, vz, 0, 0, 0, 0, }, {1, 0, 0, 0, 0, -hail[2].velocity.x, vx, 0, 0, }, {0, 1, 0, 0, 0, -hail[2].velocity.y, vy, 0, 0, }, {0, 0, 1, 0, 0, -hail[2].velocity.z, vz, 0, 0, }, {1, 0, 0, 0, 0, 0, 0, -hail[3].velocity.x, vx, }, {0, 1, 0, 0, 0, 0, 0, -hail[3].velocity.y, vy, }, {0, 0, 1, 0, 0, 0, 0, -hail[3].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 = { hail[1].position.x,   hail[1].position.y,    hail[1].position.z,    hail[2].position.x,    hail[2].position.y,    hail[2].position.z,    hail[3].position.x,    hail[3].position.y,    hail[3].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 x[i]:isWholeNumber then return false end end -- check for negative time values for i=4,9 do   if x[i] < 0 then return false end end -- okay, now verify solutions for all other hailstones local unk = { x=x[1]:toWholeNumber, y=x[2]:toWholeNumber, z=x[3]:toWholeNumber, } for i=4, #hail do    -- 2 unknowns, 2 equations -- -h.vx*t + unk.vx*u = h.x - unk.x   local M2 = Matrix:new{ { -hail[i].velocity.x, vx }, { -hail[i].velocity.y, vy }, }   M2:apply(mkrat) local b = { hail[i].position.x - unk.x,     hail[i].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 x[i] < 0 or not x[i]:isWholeNumber then return false end end -- check that t/u values also work for z axis local t,u = x[1]:toWholeNumber, x[2]:toWholeNumber local hz = hail[i].position.z + t*hail[i].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 seen[key] then return false end -- already checked seen[key] = 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(hail[1].position, hail[1].velocity) local a2,b2,c2,d2 = abcd(hail[2].position, hail[2].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 = hail[i], hail[j] -- for all pairs of coords: for k=1,#coords do       for l=k+1,#coords do          local d1,d2 = coords[k], coords[l] -- for all prime factors in common between the velocity vectors local common1 = gcd(abs(hail1.velocity[d1]), abs(hail1.velocity[d2])) local common2 = gcd(abs(hail2.velocity[d1]), abs(hail2.velocity[d2])) 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.position[d1] % pp               local pos2d1 = hail2.position[d1] % pp                local pos1d2 = hail1.position[d2] % pp                local pos2d2 = hail2.position[d2] % 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) vResults[key] = { coord=d1, mod=pp } --print("pos", d1,"=", pos1d1, "mod", pp) key = string.format("p%s=%smod%s", d1, pos1d1, pp) pResults[key] = { 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) vResults[key] = { coord=d2, mod=pp } --print("pos", d2,"=", pos1d2, "mod", pp) key = string.format("p%s=%smod%s", d2, pos1d2, pp) pResults[key] = { 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 = vResults[k] print(string.format("Velocity %s = 0 mod %s", r.coord, r.mod)) v[r.coord] = v[r.coord] * r.mod end print("So:") for _,c in ipairs(coords) do   print(string.format("Velocity %s = * %s", c, v[c])) end for _,k in ipairs(sortedKeys(pResults)) do   local r = pResults[k] 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 = * 3 Velocity y = * 6 Velocity z = * 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 = hail[1], hail[2], hail[3]

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 = hail[1] local rows = {} local rhs = {} for i=2,5 do     local h2 = hail[i] table.insert(rows, {                    (h2.velocity[d] - h1.velocity[d]),                     (h1.velocity.x - h2.velocity.x),                     (h1.position[d] - h2.position[d]),                     (h2.position.x - h1.position.x),      }) table.insert(rhs,                  h2.position.x * h2.velocity[d] -                   h2.position[d] * h2.velocity.x -                   h1.position.x * h1.velocity[d] +                   h1.position[d] * 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 = x[1] ans.position[d] = x[2] ans.velocity.x = x[3] ans.velocity[d] = x[4] 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 modules[name] == nil then modules[name] = true modules[name] = (builders[name])(myrequire) end return modules[name] end return myrequire('day24') end)