Day 12 turned into a learning journey into Haskell's ad-hoc polymophism and typeclasses.
Part 1
I started with the position and velocity vectors represented as V3 Integer
values from Data.Linear
. A Planet
was a record that bundled the position and velocity vectors together, and the four satellites were bundled together into Set
of such records.
import Linear (V3(..), (^+^), (^-^))
import qualified Data.Set as S
type Vec = V3 Integer
data Planet = Planet { _pos :: Vec, _vel :: Vec} deriving (Show, Eq, Ord)
type Planets = S.Set Planet
The simulation was pretty straightforward. Each simulationStep
applies the changes due to gravity, then applies the changes due to velocity. applyGravity
turns out to be O(n2) as the gravity change for each planet depends on the gravitational effect of each other planet. That's why there's both the map
and fold
in there.
simulate = iterate simulationStep
simulationStep planets = planets''
where planets' = applyGravity planets
planets'' = applyVelocity planets'
applyGravity planets = S.map (applyGravityHere planets) planets
applyGravityHere planets here = S.foldl' updateGravity here planets
updateGravity here there = here { _vel = vel'}
where vel = _vel here
vel' = vel ^+^ gravity ((_pos there) ^-^ (_pos here))
gravity (V3 x y z) = V3 (signum x) (signum y) (signum z)
applyVelocity = S.map applyVelocityHere
applyVelocityHere here = here {_pos = (_pos here) ^+^ (_vel here)}
The full simulation is just an infinite iterate
over the simulationStep
s. I find the final answer by drop
ing the first 1000 states (the first of which is the zeroth, starting state) than taking the systemEnergy
of the next one.
part1 :: Planets -> Integer
part1 = systemEnergy . head . drop 1000 . simulate
The energy calculation for each planet uses absSum
to add up a vector's components, then the systemEnergy
is another fold
over the set of planets.
absSum (V3 x y z) = (abs x) + (abs y) + (abs z)
potentalEnergy planet = absSum $ _pos planet
kineticEnergy planet = absSum $ _vel planet
totalEnergy planet = (potentalEnergy planet) * (kineticEnergy planet)
systemEnergy = (S.foldl' (+) 0) . (S.map totalEnergy)
Part 2
Simplifying part 2 relies on two observations.
- The motion of each moon in each dimension is independent of its motion in the other dimensions. That allows us to simulate each dimension independently.
- The simulation equations are symmetical in time, which means the first repeated state is the initial state.
These observations mean that I an find the period of each dimension separately, and the period of the whole system is the least common multiple of these three periods.
But simulating in one dimension is just about the same as simulating in three dimensions. All the vector calculations are the same, apart from the few low-level functions that explicitly refer to the vector's components. But Haskell's type system means that if a function is defined to operate over V3 Integers
, it won't operate over V1 Integer
s.
One way to do this is to continue to use V3 Integer
s, but just two of the components to zero. Another (more interesting) way to do it is to use Haskell's polymorphism features to ensure that all the simulation functions operate over both types of vector.
I did this by defining my own typeclass, NVec
, and defining the operations required to support it. I had to use ^+^^
rather than ^+^
to avoid confusion between my operator and the one from the Data.Linear
module. nvZero
, nvSignum
and nvAbsSum
are the low-level operations that are used in the simulation.
Given the typeclass, I can define how V1 Integer
and V3 Integer
implement that class, and give the definitions of the required functions.
class (Eq a) => NVec a where
(^+^^) :: a -> a -> a
(^-^^) :: a -> a -> a
nvZero :: a
nvSignum :: a -> a
nvAbsSum :: a -> Integer
instance NVec (V1 Integer) where
x ^+^^ y = x ^+^ y
x ^-^^ y = x ^-^ y
nvZero = V1 0
nvSignum (V1 x) = V1 (signum x)
nvAbsSum (V1 x) = abs x
instance NVec (V3 Integer) where
x ^+^^ y = x ^+^ y
x ^-^^ y = x ^-^ y
nvZero = V3 0 0 0
nvSignum (V3 x y z) = V3 (signum x) (signum y) (signum z)
nvAbsSum (V3 x y z) = (abs x) + (abs y) + (abs z)
Given this class, I can now define the data types and functions in terms of this class. For example,
data Planet a = Planet { _pos :: a , _vel :: a} deriving (Show, Eq, Ord)
type Planets a = V.Vector (Planet a)
simulate :: (NVec a) => Planets a -> [Planets a]
simulate = iterate simulationStep
applyGravity :: (NVec a) => Planets a -> Planets a
applyGravity planets = V.map (applyGravityHere planets) planets
(Another change I made was in the type of the container for the planets. The state only repeats when each planet is in its original state. But a Set
is unordered, so it could be that planet a ends up looking like planet b used to, while b now looks like a. The sets representing these state would be identical, even if the planets themselves aren't. I got around that by using a Vector
(fixed-size one-dimenional array) of planets; that allows me to distinguish the planets by their position in the Vector
.)
Once the class was defined, I needed functions to split a single vector of three-dimensional planets into a list of three vectors of one-dimensional planets. These functions need not be polymorphic.
unzipPlanets :: V.Vector (Planet (V3 Integer)) -> [V.Vector (Planet (V1 Integer))]
unzipPlanets = dimensionSlice . V.map unzipPlanet
unzipPlanet :: Planet (V3 Integer) -> [Planet (V1 Integer)]
unzipPlanet planet = map mkPlanet posVecs
where posVecs = unzipVec $ _pos planet
unzipVec :: V3 Integer -> [V1 Integer]
unzipVec (V3 x y z) = [V1 x, V1 y, V1 z]
dimensionSlice :: (NVec a) => V.Vector [Planet a] -> [Planets a]
dimensionSlice slicedPlanets = [sliceDim d | d <- [0..2]]
where sliceDim d = V.map (!!d) slicedPlanets
Code
The code is available (and on Github).