2019/day22.m4

   1 divert(-1)dnl -*- m4 -*-
   2 # Usage: m4 [-Dfile=day22.input] day22.m4
   3 # Optionally use -Donly=[12] to skip a part
   4
   5 include(`common.m4')ifelse(common(22), `ok', `',
   6 `errprint(`Missing common initialization
   7 ')m4exit(1)')
   8
   9 include(`math64.m4')
  10
  11 define(`input', quote(translit(include(defn(`file')), `
  12 ', `,')))
  13
  14 # https://en.wikipedia.org/wiki/Montgomery_modular_multiplication
  15 define(`splitinto', `pushdef(`i', 0)_foreach(`_$0(`$1',', `)', `', chunk($2,
  16   decr(len(B))))ifelse(eval(i <= 2 * r - 1), 1, `forloop(i, eval(2 * r - 1),
  17   `_$0(`$1', 0,', `)')')popdef(`i')')
  18 define(`_splitinto', `define(`$1'i, `$2')define(`i', incr(i))')
  19 define(`combine', `rebuild(decr(len(B))forloop(eval($2 + 0),
  20   eval(2 * r - 1), `, defn(`$1'', `)'))')
  21 define(`divide', `_$0(eval((len($1) >= len(B)) * (len($1) - len(B) + 1)), $@,
  22   `$2')')
  23 define(`_divide', `define(`$4', ifelse($1, 0, 0, `substr($2, 0,
  24   $1)'))define(`$3', ifelse($1, 0, $2, `trim(substr($2, $1))'))')
  25 define(`redc', `splitinto(`T', $1)forloop_var(`i', 0, decr(r), `define(`c',
  26   0)divide(eval(defn(`T'i) * Np), `m')forloop_var(`j', 0, eval(2 * r - 1),
  27   `_$0()')')normalize(combine(`T', r), N)')
  28 define(`_redc', `define(`x', eval(defn(`T'eval(i + j)) + m * defn(`N'j) +
  29   c))divide(x, `T'eval(i + j), `c')')
  30 define(`normalize', `ifelse(lt64($1, $2), 1, $1, `sub64($1, $2)')')
  31 define(`mont', `redc(mul64($1, R2_N))')
  32
  33 # Convert each line into f(x)=ax+b linear function with non-negative a and b.
  34 # Store a and b in Montgomery form to avoid needing 64-bit divide for part 2.
  35 define(`parse', `ifelse($1, `', `', `_$0(translit($1, ` ', `,'))')')
  36 define(`_parse', `pushdef(`act', ifelse($1, cut, ``1, sub64(N, $2)'',
  37   $3, new, ``N_, N_'', ``$4, 0''))')
  38 foreach(`parse', input)
  39
  40 define(`setup', `define(`B', $1)define(`r', $2)define(`N', $3)define(`Np',
  41   $4)define(`N_', sub64(N, 1))define(`R2_N', $5)splitinto(`N',
  42   N)define(`R_N', redc(R2_N))define(`a', R_N)define(`b', 0)stack_foreach(`act',
  43   `docompose')')
  44 define(`apply', `redc(add64(redc(mul64($2, mont($1))), $3))')
  45 define(`docompose', `compose($1)')
  46 define(`compose', `_$0(a, b, mont($1), mont($2), `a', `b')')
  47 define(`_compose', `define(`$5', redc(mul64($3, $1)))define(`$6',
  48   normalize(add64(redc(mul64($3, $2)), $4), N))')
  49 define(`bits', `ifelse($1, 0, 0, $1, 1, 1, `_$0(mul64($1, 5)), eval(substr($1,
  50   decr(len($1))) & 1)')')
  51 define(`_bits', `bits(substr($1, 0, decr(len($1))))')
  52 define(`compose_n', `define(`an', a)define(`bn', b)foreach(`_$0',
  53   shift(bits($1)))')
  54 define(`_compose_n', `_compose(an, bn, an, bn, `an', `bn')ifelse($1, 1,
  55   `_compose(an, bn, a, b, `an', `bn')')')
  56
  57 # precomputed: for R=1000000==B^r and N=10007, Nprime=57 and R^2%N=9664
  58 ifelse(ifdef(`only', only, 1), 1, `
  59 setup(100, 3, 10007, 57, 9664)
  60 define(`part1', apply(2019, a, b))
  61 ')
  62
  63 # precomputed: for R=10000000000000000==B^r and N=119315717514047,
  64 # Nprime=9617 and R^2%N=76941077250887
  65 ifelse(ifdef(`only', only, 2), 2, `
  66 setup(10000, 4, 119315717514047, 9617, 76941077250887)
  67 compose_n(sub64(N_, 101741582076661))
  68 define(`part2', apply(2020, an, bn))
  69 ')
  70
  71 divert`'part1
  72 part2