factorization-coords.myr

   1 use std
   2
   3 use t
   4 use yakmo
   5
   6 use "types"
   7 use "kc"
   8 use "roots"
   9 use "zzz"
  10
  11 pkg =
  12         /*
  13            The input is mc = [ (1, g1), (2, g2), ..., (m, gm) ][:], where
  14
  15              gk = Δ^{γ_k}(y) for y in N_{-} .
  16
  17            Note that the first r coordinates are not coordinates
  18            assigned to the quiver vertices. This is because they turn
  19            out to be rational functions of solely edge coordinates.
  20
  21            The output is [ `X(i1, t1), `X(i2, t2), ..., `X(im, tm) ][:],
  22            which give factorization coordinates for x in N_{+}, such
  23            that
  24
  25              y = π_{-}(x) = w0^{-1} [ x w0^{-1} ]_{+} wo
  26
  27            Note that, due to the π_{-}(x) stuff, this map is actually
  28            (ΨΦΨ)(u). Give minor coordinates for u, receive factorization
  29            coordinates for (ΨΦΨ)(u).
  30          */
  31         const interior_fc_from_mc : (kc : KC_type, mc : minor_coordinate[:] -> std.result(factorization_coordinate[:], byte[:]))
  32         const interior_fc_from_mc_and_undo_psi_phi_psi : (kc : KC_type, mc : minor_coordinate[:] -> std.result(factorization_coordinate[:], byte[:]))
  33
  34         /*
  35            The inverse: given [ `X(i1, v1), ... `X(im, vm) ][:], produce
  36            [ g1, g2, ..., gm ].
  37
  38            Note that, due to the π_{-}(x) stuff, this map is actually
  39            (ΨΦΨ)^{-1}(u). Give factorization coordinates for u, receive
  40            minor coordinates for (ΨΦΨ)^{-1}(u).
  41          */
  42         const interior_mc_from_fc : (kc : KC_type, fc : factorization_coordinate[:]  -> std.result(minor_coordinate[:], byte[:]))
  43
  44         /*
  45            A degenerate case - minor coordinate are exactly
  46            factorization coordinates on edges.
  47
  48            Note that this means that if mc[k] = (j, v), that j is used as
  49
  50              v = χ_{ω_{i_j}}(h)
  51          */
  52         const edge_mc_from_fc : (kc : KC_type, fc : factorization_coordinate[:]  -> std.result(minor_coordinate[:], byte[:]))
  53         const edge_fc_from_mc : (kc : KC_type, mc : minor_coordinate[:]  -> std.result(factorization_coordinate[:], byte[:]))
  54
  55         /*
  56            A specialized, degenerate case. Given H = [ `H(i1, v1), ...,
  57            `H(in, vn) ], return
  58
  59                v = χ_{ω_{j}}(H)
  60
  61            Which may be 1, or may be a product of vjs.
  62          */
  63         const chi : (fc : factorization_coordinate[:], j : weyl_reflection -> std.option(yakmo.ratfunc#))
  64
  65         /*
  66            Another specialized, degenerate case. Given H = [ [.k=k1,
  67            .c=c1], ..., [.k=kn, .c=cn] ], and k, return the ci
  68            corresponding to ki = k.
  69          */
  70         const delta_gamma_k : (h : minor_coordinate[:], k : int -> yakmo.ratfunc#)
  71
  72         /* The Φ and Ψ maps. */
  73         const psi : (kc : KC_type, fc : factorization_coordinate[:] -> std.result(factorization_coordinate[:], byte[:]))
  74         const phi : (kc : KC_type, fc : factorization_coordinate[:] -> std.result(factorization_coordinate[:], byte[:]))
  75         const phi_inv : (kc : KC_type, fc : factorization_coordinate[:] -> std.result(factorization_coordinate[:], byte[:]))
  76
  77         /* ι : given (y, h, x), return (X, Y, H) such that yhx = XHY */
  78         const iota : (kc : KC_type, y : factorization_coordinate[:], h : factorization_coordinate[:], x : factorization_coordinate[:] -> std.result((factorization_coordinate[:], factorization_coordinate[:], factorization_coordinate[:]), byte[:]))
  79
  80         /*
  81            compute w(h) = w h w^-1, where w \in W and h is given by
  82            factorization coordinates. TODO: remove apply_word, it's
  83            superceded by apply_wk.
  84
  85            apply_wk computes w_k h w_k⁻¹. required_sign is a little
  86            special: it controls whether the "numerator", "denominator",
  87            or "entire fraction" is returned (scare quotes because the
  88            exponents are treated abstractly).
  89          */
  90         const apply_word : (kc : KC_type, w : weyl_reflection[:], h : factorization_coordinate[:] -> std.result(factorization_coordinate[:], byte[:]))
  91         const apply_wk : (kc : KC_type, k : int, h : factorization_coordinate[:], required_sign : int -> std.result(factorization_coordinate[:], byte[:]))
  92
  93         /* Convenience methods for fiddling with H */
  94         const multiply_h_mc_by_sG : (kc : KC_type, h : minor_coordinate[:] -> std.result(minor_coordinate[:], byte[:]))
  95         const apply_wo_of_inv : (kc : KC_type, h : minor_coordinate[:] -> std.result(minor_coordinate[:], byte[:]))
  96
  97         /* Take a jumbled mess of factorization coords, put it in Y H X (or other) format */
  98         type canonical_ordering = union
  99                 `Ordering_YHX
 100                 `Ordering_XYH
 101                 `Ordering_XHY
 102         ;;
 103         const canonicalize_factorization_coords : (kc : KC_type, fc : factorization_coordinate[:]#, ord : canonical_ordering -> std.result(void, byte[:]))
 104
 105         /*
 106            Given an alpha constructed from a quiver (missing the first
 107            $rank coordinates, fill them in
 108          */
 109         const fill_in_alpha : (kc : KC_type, alpha : conf_3_star# -> std.result(void, byte[:]))
 110         const fill_in_u_from_h01_h12_h20 : (kc : KC_type, u : minor_coordinate[:]#, h01 : minor_coordinate[:], h12 : minor_coordinate[:], h02 : minor_coordinate[:] -> std.result(void, byte[:]))
 111
 112         /*
 113            The "other" Ψ maps. These are the ones that split Conf₄* into
 114            two copies of Conf₃*
 115          */
 116         const psi_02 : (kc : KC_type, alpha : conf_4_star# -> std.result((conf_3_star#, conf_3_star#), byte[:]))
 117         const psi_13 : (kc : KC_type, alpha : conf_4_star# -> std.result((conf_3_star#, conf_3_star#), byte[:]))
 118
 119         /* Perform non-quiver rotation */
 120         const rotate : (kc : KC_type, alpha : conf_3_star# -> std.result(conf_3_star#, byte[:]))
 121 ;;
 122
 123 const interior_fc_from_mc = {kc, mc
 124         var c = get_coxeter_elt(kc)
 125         var iforward = get_w0(kc)
 126         if mc.len != iforward.len
 127                 -> `std.Err std.fmt("expected {} minor coordinates, got {}", iforward.len, mc.len)
 128         ;;
 129         for var k = 0; k < mc.len; ++k
 130                 if mc[k].k != k + 1
 131                         -> `std.Err std.fmt("expected minor coordinate {}  to be  Δ^{{γ_{}}}, was declared as Δ^{{γ_{}}}", k + 1, k + 1, mc[k].k)
 132                 ;;
 133         ;;
 134
 135         var irev : weyl_reflection[:] = std.slalloc(iforward.len)
 136         for var k = 0; k < iforward.len; ++k
 137                 irev[k] = iforward[iforward.len - 1 - k]
 138         ;;
 139         auto irev
 140
 141         var rank : int = get_n(kc)
 142         var A : int[:][:] = get_cartan_matrix(kc)
 143
 144         var fc : factorization_coordinate[:] = [][:]
 145
 146         for var k = 0; k < iforward.len; ++k
 147                 /*
 148                    Construct t_{k + 1}. This is given by (in the
 149                    following notation, everything is 1-based):
 150
 151                      den = delta^{w_{k + 0} ω_{i_k}}(y) · delta^{w_{k + 1} ω_{i_k}}(y)
 152                      num = prod_{j != i_k} ( delta^{w_k ω_j}(y) )^{-A_{j, i_k}}
 153                      t_k = num / den
 154
 155                    Since our input is minor coordinates of y, computing
 156                    these deltas is simply a matter of figuring out which
 157                    γ_k each thing refers to.
 158                  */
 159
 160                 var num : yakmo.ratfunc# = yakmo.rid()
 161                 auto num
 162                 for var j = 0; j < rank; ++j
 163                         /* Simplify w_k ω_j to some chamber weight */
 164                         var jw : weyl_reflection = (j + 1 : weyl_reflection)
 165                         var which_k = simplify_chamber_weight(irev[:irev.len - k], jw, irev.len)
 166                         if which_k <= 0
 167                                 continue
 168                         ;;
 169
 170                         /*
 171                            Now that we know the chamber weight, pull out
 172                            the correct coordinate, raise it to
 173                            -A_{j,i_k}, and multiply the result into the
 174                            numerator.
 175                          */
 176                         var this_delta_of_y = mc[which_k - 1].c
 177                         for var l = 0; l > A[j][iforward[k] - 1]; --l
 178                                 yakmo.rmul_ip(num, this_delta_of_y)
 179                         ;;
 180                 ;;
 181
 182                 var den : yakmo.ratfunc# = t.dup(mc[k].c)
 183                 auto den
 184                 if irev.len - k - 1 > 0
 185                         var other_k = simplify_chamber_weight(irev[:irev.len - k - 1], iforward[k], irev.len)
 186                         if other_k > 0
 187                                 yakmo.rmul_ip(den, mc[other_k - 1].c)
 188                         ;;
 189                 ;;
 190
 191                 match yakmo.div_maybe(num, den)
 192                 | `std.None: -> `std.Err std.fmt("t_{} = [ {} ] / [ {} ] is not a Laurent polynomial: damn", k + 1, num, den)
 193                 | `std.Some x: std.slpush(&fc, `X(iforward[k], x))
 194                 ;;
 195         ;;
 196
 197         -> `std.Ok fc
 198 }
 199
 200 /*
 201    Given w_j ω_l = (s_i1 s_i2 ⋯ s_im) ω_l, with w_j a prefix of i^{-1},
 202    elements from the end of w_j, repeatedly, following the rule
 203
 204      for i != j,     s_i ω_j = ω_j .
 205
 206    The result will either be 1) a shorter prefix of i^{-1}, or 2) empty.
 207    If it's a shorter prefix of length k, then we've found γ_k, so we
 208    return k. If it's empty, we should just return 0 (we don't
 209    particularly care about which ω_j we're using in that case).
 210
 211    To ensure that w_j is really a prefix of i^{-1}, the input is simply
 212    i^{-1} and the length of that prefix.
 213  */
 214 const simplify_chamber_weight = {w : weyl_reflection[:], which_omega : weyl_reflection, full_len : std.size
 215         while w.len > 0 && w[w.len - 1] != which_omega
 216                 w = w[:w.len - 1]
 217         ;;
 218
 219         if w.len == 0
 220                 -> -1
 221         ;;
 222
 223         -> (full_len - w.len : int) + 1
 224 }
 225
 226 const interior_mc_from_fc = {kc, fc
 227         var bkl : yakmo.Z#[:][:] = [][:]
 228         match get_bkl(kc)
 229         | `std.Err e: -> `std.Err std.fmt("get_bkl({}): {}", kc, e)
 230         | `std.Ok x: bkl = x
 231         ;;
 232
 233         var mc : minor_coordinate[:] = [][:]
 234         var w0 : weyl_reflection[:] = get_w0(kc)
 235
 236         /*
 237            First, make sure the factorization coordinates follow our
 238            idea of the longest word.
 239          */
 240         if fc.len != w0.len
 241                 -> `std.Err std.fmt("expected {} factorization coordinates, got {}", w0.len, fc.len)
 242         ;;
 243
 244         var fcv = std.slalloc(fc.len)
 245         for var j = 0; j < fc.len; ++j
 246                 match fc[j]
 247                 | `X (i, l):
 248                         if i != w0[j]
 249                                 -> `std.Err std.fmt("factorization coordinate {} should be X({}, _), was X({}, _)", j + 1, w0[j], i)
 250                         ;;
 251                         fcv[j] = l
 252                 | `Y (i, l): -> `std.Err std.fmt("factorization coordinate {} should be X({}, _), was Y({}, _)", j + 1, w0[j], i)
 253                 | `H (i, l): -> `std.Err std.fmt("factorization coordinate {} should be X({}, _), was H({}, _)", j + 1, w0[j], i)
 254                 ;;
 255         ;;
 256
 257
 258         for var k = 0; k < fc.len; ++k
 259                 var gk : yakmo.ratfunc# = yakmo.rid()
 260                 for var l = k; l < fc.len; l++
 261                         /* gk *= (fcv[l])^bkl(k,l) */
 262                         var power = (t.dup(bkl[k][l]) : std.bigint#)
 263                         var factor = fcv[l]
 264
 265                         if power.sign < 0
 266                                 while !std.bigiszero(power)
 267                                         match yakmo.div_maybe(gk, factor)
 268                                         | `std.None: -> `std.Err std.fmt("cannot compute g_{}; could not divide {} / {}", k + 1, gk, factor)
 269                                         | `std.Some r:
 270                                                 __dispose__(gk)
 271                                                 gk = r
 272                                         ;;
 273                                         std.bigaddi(power, 1)
 274                                 ;;
 275                         else
 276                                 while !std.bigiszero(power)
 277                                         yakmo.rmul_ip(gk, fcv[l])
 278                                         std.bigaddi(power, -1)
 279                                 ;;
 280                         ;;
 281                         std.bigfree(power)
 282                 ;;
 283                 yakmo.reduceratfunc(gk)
 284                 std.slpush(&mc, [ .k = k + 1, .c = gk ])
 285         ;;
 286
 287         /*
 288            Don't __dispose__; that will recurse and kill the variables,
 289            which actually belong to the input fc.
 290          */
 291         std.slfree(fcv)
 292
 293         -> `std.Ok mc
 294 }
 295
 296 const edge_mc_from_fc = {kc : KC_type, fc : factorization_coordinate[:]
 297         var rank = get_n(kc)
 298         if fc.len > rank
 299                 /* We're allowed to omit `H(_, 1), it's the identity */
 300                 -> `std.Err std.fmt("expected {} factorization coordinates, got {}", rank, fc.len)
 301         ;;
 302         var mc : minor_coordinate[:] = [][:]
 303         for var k = 1; k <= rank; ++k
 304                 std.slpush(&mc, [ .k = k, .c = yakmo.rid() ])
 305         ;;
 306         for var k = 0; k < fc.len; ++k
 307                 match fc[k]
 308                 | `H(i, val):
 309                         /*
 310                            `H(coxeter[j], rf) =~= [ .k = j, rf ]
 311
 312                            Also, the minor coordinate with .k = j is at index j.
 313                          */
 314                         var j = inv_coxeter(kc, i)
 315                         yakmo.rmul_ip(mc[j].c, val)
 316                 | _: -> `std.Err std.fmt("expected factorization coordinate {} to be H(_ ,_), got {}", k + 1, fc[k])
 317                 ;;
 318         ;;
 319
 320         -> `std.Ok mc
 321 }
 322
 323 const edge_fc_from_mc = {kc : KC_type, mc : minor_coordinate[:]
 324         var rank = get_n(kc)
 325         var c = std.try(get_coxeter_elt(kc))
 326         if mc.len != rank
 327                 -> `std.Err std.fmt("expected {} minor coordinates, got {}", rank, mc.len)
 328         ;;
 329         var fc = [][:]
 330         for var k = 0; k < rank; ++k
 331                 /* `H(coxeter[j], rf) =~= [ .k = j, rf ] */
 332                 std.slpush(&fc, `H(c[k], delta_gamma_k(mc, k + 1)))
 333         ;;
 334
 335         -> `std.Ok fc
 336 }
 337
 338 /*
 339    Given
 340
 341      [ `X(i1, v1), `X(i2, v2), ... `X(im, vm) ][:] = x in N_{+} ∩ G_0 w_0,
 342
 343    return
 344
 345      [ `Y(im, vm), `Y(i{m-1}, v{m-1}), ... `Y(i1, v1) ][:] = y in N_{-} ∩ w_0 G_0,
 346
 347    such that
 348
 349      y = [ x w0 ]_{-}
 350  */
 351 const phi = {kc : KC_type, fc : factorization_coordinate[:]
 352         var i : weyl_reflection[:] = get_w0(kc)
 353         var ireverse : weyl_reflection[:] = std.slalloc(i.len)
 354         for var k = 0; k < i.len; ++k
 355                 ireverse[k] = i[i.len - 1 - k]
 356         ;;
 357         auto ireverse
 358         var ret : factorization_coordinate[:] = std.slalloc(fc.len)
 359
 360         /* Ensure that fc may actually be something in the set we care about */
 361         if fc.len != i.len
 362                 -> `std.Err std.fmt("fc.len = {}, expected {}", fc.len, i.len)
 363         ;;
 364         for var j = 0; j < fc.len; ++j
 365                 match fc[j]
 366                 | `X (k, l):
 367                         if k != i[j]
 368                                 -> `std.Err std.fmt("coordinate {} was X({}, _), should have been X({}, _)", j + 1, k, i[j])
 369                         ;;
 370                         ret[j] = `X(k, t.dup(l))
 371                 | `Y (k, _): -> `std.Err std.fmt("coordinate {} was Y({}, _), should have been X({}, _)", j + 1, k, i[j])
 372                 | `H (k, _): -> `std.Err std.fmt("coordinate {} was H({}, _), should have been X({}, _)", j + 1, k, i[j])
 373                 ;;
 374         ;;
 375
 376         /*
 377            Consider (ret, i), representing (initially) the product x w0.
 378            Look at the last factorization element of x and the first
 379            Weyl reflection of w0.
 380
 381            They had better be X(i, _) and s_i for the same i. So use identity
 382
 383              X(i, t) s_i = Y(i, 1/t)
 384
 385            to transform the right end of x and the left end of the word.
 386
 387            Rem: we're actually using w0^-1 here, so that it lines up
 388            with the factorization of x correctly. That's okay, because
 389            w0^-1 * w0 * w0 = w0, and w0^2 is an element of H, which is
 390            ignored by [·]_{-}. Therefore
 391
 392              [ x w0 ]_{-} = [x w0^{-1} w0^2 ]_{-} = [ x w0^{-1} ]_{-}.
 393          */
 394         for var j = 0; j < i.len; ++j
 395                 match ret[ret.len - 1]
 396                 | `X(xj, val):
 397                         if xj != ireverse[j]
 398                                 -> `std.Err std.fmt("coordinates, word = {}, {}, this doesn't work[1]", ret, i[j:])
 399                         ;;
 400
 401                         /* Quick inversion */
 402                         if val.num.terms.len == 0
 403                                 -> `std.Err std.fmt("coordinates, word = {}, {}, this doesn't work[2]", ret, i[j:])
 404                         ;;
 405                         std.swap(&val.num, &val.den)
 406                         ret[ret.len - 1] = `Y(xj, val)
 407                 | _: -> `std.Err std.fmt("coordinates, word = {}, {}, this doesn't work[3]", ret, i[j:])
 408                 ;;
 409
 410                 match canonicalize_factorization_coords(kc, &ret, `Ordering_YHX)
 411                 | `std.Ok void:
 412                 | `std.Err e: -> `std.Err std.fmt("canonicalize: {}", e)
 413                 ;;
 414         ;;
 415
 416         /*
 417            We now have something in form Y H X. We had better have
 418            exactly i.len Y-shaped elements, because we created one for
 419            each element of i. So let's kill everything after that.
 420          */
 421         match ret[i.len - 1]
 422         | `Y(_, _):
 423         | _: -> `std.Err std.fmt("not enough elements of Y in decomposition: expected {}, got {}", i.len, ret)
 424         ;;
 425
 426         while ret.len > i.len
 427                 match ret[i.len]
 428                 | `X(_, xval): __dispose__(xval)
 429                 | `H(_, hval): __dispose__(hval)
 430                 | `Y(_, _): -> `std.Err std.fmt("too many elements of Y in decomposition: expected {}, got {}", i.len, ret)
 431                 ;;
 432                 std.sldel(&ret, i.len)
 433         ;;
 434
 435         -> `std.Ok ret
 436 }
 437
 438 /*
 439    Given
 440
 441      [ `Y(im, vm), ..., `Y(i2, v2), `Y(i1, v1) ][:] = y in N_{-} ∩ w_0 G_0,
 442
 443    return
 444
 445      [ `X(i1, v1), `X(i2, v2), ..., `X(im, vm) ][:] = x in N_{+} ∩ G_0 w_0,
 446
 447    such that Φ(x) = y. We do this by computing w₀( [w₀⁻¹ y]_{-} ).
 448  */
 449 const phi_inv = {kc : KC_type, fc : factorization_coordinate[:]
 450         var i : weyl_reflection[:] = get_w0(kc)
 451         var ireverse : weyl_reflection[:] = std.slalloc(i.len)
 452         for var k = 0; k < i.len; ++k
 453                 ireverse[k] = i[i.len - 1 - k]
 454         ;;
 455         auto ireverse
 456         var ret : factorization_coordinate[:] = std.slalloc(fc.len)
 457
 458         /* Ensure that fc may actually be something in the set we care about */
 459         if fc.len != i.len
 460                 -> `std.Err std.fmt("fc.len = {}, expected {}", fc.len, i.len)
 461         ;;
 462         for var j = 0; j < fc.len; ++j
 463                 match fc[j]
 464                 | `Y (k, l):
 465                         if k != ireverse[j]
 466                                 -> `std.Err std.fmt("coordinate {} was Y({}, _), should have been Y({}, _)", j + 1, k, ireverse[j])
 467                         ;;
 468                         ret[j] = `Y(k, t.dup(l))
 469                 | `X (k, _): -> `std.Err std.fmt("coordinate {} was X({}, _), should have been X({}, _)", j + 1, k, i[j])
 470                 | `H (k, _): -> `std.Err std.fmt("coordinate {} was H({}, _), should have been X({}, _)", j + 1, k, i[j])
 471                 ;;
 472         ;;
 473
 474         for var j = 0; j < i.len; ++j
 475                 match ret[ret.len - 1]
 476                 | `Y(yj, val):
 477                         if yj != i[j]
 478                                 -> `std.Err std.fmt("coordinates, word = {}, {}, this doesn't work[1]", ret, ireverse[j:])
 479                         ;;
 480
 481                         /* Quick inversion */
 482                         if val.num.terms.len == 0
 483                                 -> `std.Err std.fmt("coordinates, word = {}, {}, this doesn't work[2]", ret, ireverse[j:])
 484                         ;;
 485                         std.swap(&val.num, &val.den)
 486                         ret[ret.len - 1] = `X(yj, val)
 487                 | _: -> `std.Err std.fmt("coordinates, word = {}, {}, this doesn't work[5]", ret, ireverse[j:])
 488                 ;;
 489
 490                 canonicalize_factorization_coords(kc, &ret, `Ordering_XHY)
 491         ;;
 492
 493         match ret[i.len - 1]
 494         | `X(_, _):
 495         | _: -> `std.Err std.fmt("not enough elements of X in decomposition: expected {}, got {}", i.len, ret)
 496         ;;
 497
 498         while ret.len > i.len
 499                 match ret[i.len]
 500                 | `Y(_, xval): __dispose__(xval)
 501                 | `H(_, hval): __dispose__(hval)
 502                 | `X(_, _): -> `std.Err std.fmt("too many elements of X in decomposition: expected {}, got {}", i.len, ret)
 503                 ;;
 504                 std.sldel(&ret, i.len)
 505         ;;
 506
 507         -> `std.Ok ret
 508
 509 }
 510
 511 const iota = { kc : KC_type, y, h, x
 512         /* concatenate */
 513         var total : factorization_coordinate[:] = [][:]
 514         for q : [ y, h, x ][:]
 515                 for var j = 0; j < q.len; ++j
 516                         match q[j]
 517                         | `Y(k, v): std.slpush(&total, `Y(k, t.dup(v)))
 518                         | `H(k, v): std.slpush(&total, `H(k, t.dup(v)))
 519                         | `X(k, v): std.slpush(&total, `X(k, t.dup(v)))
 520                         ;;
 521                 ;;
 522         ;;
 523
 524         match canonicalize_factorization_coords(kc, &total, `Ordering_XYH)
 525         | `std.Ok void:
 526         | `std.Err e: -> `std.Err std.fmt("canonicalize_factorization_coords: {}", e)
 527         ;;
 528
 529         /*
 530            It's possible that canonicalize didn't put things into pure X
 531            Y H form, for example if dividing by zero would happen. So we
 532            have to be picking when splitting this up into components.
 533          */
 534         var xout : factorization_coordinate[:] = [][:]
 535         var yout : factorization_coordinate[:] = [][:]
 536         var hout : factorization_coordinate[:] = [][:]
 537         var j = 0
 538
 539         while j < total.len
 540                 match total[j]
 541                 | `X(i, v): std.slpush(&xout, `X(i, v))
 542                 | _: break
 543                 ;;
 544                 j++
 545         ;;
 546
 547         while j < total.len
 548                 match total[j]
 549                 | `X(_, _): -> `std.Err std.fmt("canonicalization failed (trailing X)")
 550                 | `Y(i, v): std.slpush(&yout, `Y(i, v))
 551                 | _: break
 552                 ;;
 553                 j++
 554         ;;
 555
 556         while j < total.len
 557                 match total[j]
 558                 | `X(_, _): -> `std.Err std.fmt("canonicalization failed (trailing X)")
 559                 | `Y(_, _): -> `std.Err std.fmt("canonicalization failed (trailing Y)")
 560                 | `H(i, v): std.slpush(&hout, `H(i, v))
 561                 | _: break
 562                 ;;
 563                 j++
 564         ;;
 565
 566         -> `std.Ok (xout, yout, hout)
 567 }
 568
 569 const chi = {fc : factorization_coordinate[:], j : weyl_reflection
 570         var ret = yakmo.rid()
 571         for var k = 0; k < fc.len; ++k
 572                 match fc[k]
 573                 | `H(j2, x):
 574                         if j2 == j
 575                                 yakmo.rmul_ip(ret, x)
 576                         ;;
 577                 | _:
 578                         __dispose__(ret)
 579                         -> `std.None
 580                 ;;
 581         ;;
 582
 583         -> `std.Some ret
 584 }
 585
 586 const delta_gamma_k = {h : minor_coordinate[:], k : int
 587         for var j = 0; j < h.len; ++j
 588                 if h[j].k == k
 589                         -> t.dup(h[j].c)
 590                 ;;
 591         ;;
 592
 593         -> yakmo.gid()
 594 }
 595
 596 /* Given a list of X, Y, H, 1) swap the order, and 2) swap between X and Y. */
 597 const psi = {kc, fc
 598         if fc.len == 0
 599                 -> `std.Ok [][:]
 600         ;;
 601
 602         var ret : factorization_coordinate[:] = std.slalloc(fc.len)
 603
 604         for var j = 0; j < fc.len; ++j
 605                 var k = fc.len - 1 - j
 606                 match fc[j]
 607                 | `H (i, l): ret[k] = `H(i, t.dup(l))
 608                 | `X (i, l): ret[k] = `Y(i, t.dup(l))
 609                 | `Y (i, l): ret[k] = `X(i, t.dup(l))
 610                 ;;
 611         ;;
 612
 613         /* TODO: maybe unnecessary */
 614         for var j = 0; j < ret.len; ++j
 615                 match ret[j]
 616                 | `H (i, rf): yakmo.reduceratfunc(rf)
 617                 | `X (i, rf): yakmo.reduceratfunc(rf)
 618                 | `Y (i, rf): yakmo.reduceratfunc(rf)
 619                 ;;
 620         ;;
 621
 622         -> `std.Ok ret
 623 }
 624
 625 const apply_word = {kc, w, h_orig
 626         var A : int[:][:] = get_cartan_matrix(kc)
 627         var h = std.sldup(h_orig)
 628         for var k = 0; k < h_orig.len; ++k
 629                 h[k] = t.dup(h_orig[k])
 630         ;;
 631
 632         /* Make sure we collapse all duplicates */
 633         match canonicalize_factorization_coords(kc, &h, `Ordering_YHX)
 634         | `std.Ok void:
 635         | `std.Err e:
 636                 auto (e : t.doomed_str)
 637                 -> `std.Err std.fmt("initial canonicalize: {}", e)
 638         ;;
 639
 640         for var k = 0; k < h.len; ++k
 641                 match h[k]
 642                 | `H(_, _):
 643                 | _: -> `std.Err std.fmt("did not expect to apply Weyl action to anything but H")
 644                         /*
 645                            We could handle this: just construct
 646
 647                                 xi(-1) yi(1) xi(-1) g xi(1) yi(-1) xi(1),
 648
 649                            canonicalize it, and ship it back. But that's
 650                            actually very slow.
 651                          */
 652                 ;;
 653         ;;
 654
 655         /*
 656            Using Lusztig's calculus, it is straight-forward to show that for
 657
 658                w = si,      H = h1^a1 h2^a2 ⋯ hr^ar,      [ A_{ij} ] the Cartan matrix,
 659
 660            the Weyl action is
 661
 662                w(H) = H hi^( \prod_k ak^A_{ki} )
 663
 664            therefore we can apply each si to h by multiplying around the
 665            various components of h
 666          */
 667         for var j = w.len - 1; j >= 0; --j
 668                 var i : weyl_reflection = w[j]
 669
 670                 /* First, find `H(i, _) so that we can mess with the variable */
 671                 var found = false
 672                 var ai = yakmo.rid()
 673                 for var k = 0; k < h.len; ++k
 674                         match h[k]
 675                         | `H(imaybe, amaybe):
 676                                 if imaybe == i
 677                                         found = true
 678                                         __dispose__(ai)
 679                                         ai = amaybe
 680                                         break
 681                                 ;;
 682                         | _:
 683                         ;;
 684                 ;;
 685
 686                 if !found
 687                         std.slpush(&h, `H(i, ai))
 688                 ;;
 689
 690                 /* Now build the \prod_j aj^A_{ji} factor */
 691                 var factor : yakmo.ratfunc# = yakmo.rid()
 692                 auto factor
 693                 for var l = 0; l < h.len; ++l
 694                         match h[l]
 695                         | `H(k, ak):
 696                                 var Aki = A[k - 1][i - 1]
 697                                 if Aki < 0
 698                                         for var z = 0; z > Aki; --z
 699                                                 yakmo.rmul_ip(factor, ak)
 700                                         ;;
 701                                 elif Aki > 0
 702                                         match yakmo.finv(ak)
 703                                         | `std.None:
 704                                         | `std.Some aki:
 705                                                 for var z = 0; z < Aki; ++z
 706                                                         yakmo.rmul_ip(factor, aki)
 707                                                 ;;
 708                                         ;;
 709                                 ;;
 710                         | _:
 711                         ;;
 712                 ;;
 713                 yakmo.rmul_ip(ai, factor)
 714                 yakmo.reduceratfunc(ai)
 715         ;;
 716
 717         -> `std.Ok h
 718 }
 719
 720 /* Compute w_k h w_k⁻¹. Note that to get w_k = w₀, use k = 1. */
 721 const apply_wk = {kc, k, h_orig, required_sign
 722         var wk_H_action_exps = get_wk_H_action_exps(kc)
 723         var h = std.sldup(h_orig)
 724         var rank = get_n(kc)
 725         if k <= 0
 726                 -> `std.Err std.fmt("k = {} is wrong. (For w₀, use k=1)", k)
 727         ;;
 728         for var l = 0; l < h_orig.len; ++l
 729                 h[l] = `H((l + 1 : weyl_reflection), yakmo.rid())
 730         ;;
 731
 732         for var i = 0; i < rank; ++i
 733                 for var j = 0; j < rank; ++j
 734                         var e = wk_H_action_exps[k - 1][i][j]
 735                         var chiv : yakmo.ratfunc#
 736                         match chi(h_orig, (j + 1 : weyl_reflection))
 737                         | `std.Some chivv: chiv = chivv
 738                         | `std.None: chiv = yakmo.rid()
 739                         ;;
 740                         auto chiv
 741
 742                         var a : yakmo.ratfunc# = yakmo.rid()
 743                         var found = false
 744                         for var l = 0; l < h.len; ++l
 745                                 match h[l]
 746                                 | `H(imaybe, amaybe):
 747                                         if imaybe == (i + 1 : weyl_reflection)
 748                                                 __dispose__(a)
 749                                                 a = amaybe
 750                                                 found = true
 751                                                 break
 752                                         ;;
 753                                 | _:
 754                                 ;;
 755                         ;;
 756                         if !found
 757                                 std.slpush(&h, `H((i + 1 : weyl_reflection), a))
 758                         ;;
 759
 760                         if e > 0 && required_sign >= 0
 761                                 for var z = 0; z < e; ++z
 762                                         yakmo.rmul_ip(a, chiv)
 763                                 ;;
 764                         elif e < 0 && required_sign <= 0
 765                                 match yakmo.finv(chiv)
 766                                 | `std.Some hoji:
 767                                         for var z = 0; z > e; --z
 768                                                 yakmo.rmul_ip(a, hoji)
 769                                         ;;
 770                                 | `std.None:
 771                                 ;;
 772                         ;;
 773                 ;;
 774         ;;
 775
 776         -> `std.Ok h
 777 }
 778
 779 const canonicalize_factorization_coords = {kc : KC_type, fc : factorization_coordinate[:]#, ord : canonical_ordering
 780         /*
 781            Use the identities
 782
 783              X(i, s) Y(j, t) = Y(j, t) X(i, s)                            for i != j
 784
 785              X(i, s) Y(i, t) = Y(i, t/(1 + st)) H(i, 1 + st) X(i, 1 + st)
 786
 787              H(i, s) Y(j, t) = Y(j, t s^{-A_{i,j}}) H(i, s)
 788
 789              Y(j, t) H(i, s) = H(i, s) Y(j, s^{A_{i,j}})
 790
 791              H(i, s) X(j, t) = X(j, t s^{A_{i,j}}) H(i, s)
 792
 793              X(j, t) H(i, s) = H(i, s) X(j, t s^{-A_{i,j}})
 794
 795            To move the components around until nothing more can be done.
 796
 797            Sometimes, we may find ourselves attempting to compute
 798            infinity accidentally.
 799
 800            Ex 1: s_1 h_1(foo) s_1^-1. This expands out to
 801
 802                     X(1, -1) Y(1, 1) X(1, -1) H(1, foo), X(1, 1) Y(1, -1) X(1, 1)
 803
 804                The H is "supposed" to move through the Xs and Ys
 805                (direction doesn't really matter), which then cancel out
 806                after creating a few more Hs. However, the naive, greedy
 807                algorithm immediately runs into problems trying to switch
 808                the first X and Y, since that involves division by zero.
 809
 810                As a workaround, just ignore division by zero. This means
 811                the result might not always be in perfect Y H X form.
 812         */
 813         var A : int[:][:] = get_cartan_matrix(kc)
 814         var one : yakmo.ratfunc# = yakmo.rid()
 815         auto one
 816
 817         /*
 818            This typemess is so that I can keep all the ordering
 819            information on the guard line. Splitting things into ordering
 820            functions and such just doesn't work with my head.
 821          */
 822         type YH_ordering = union
 823                 `Y_LT_H
 824                 `H_LT_Y
 825         ;;
 826         type YX_ordering = union
 827                 `Y_LT_X
 828                 `X_LT_Y
 829         ;;
 830         type HX_ordering = union
 831                 `H_LT_X
 832                 `X_LT_H
 833         ;;
 834
 835         var yh, yx, hx
 836         match ord
 837         | `Ordering_YHX:
 838                 yh = `Y_LT_H
 839                 yx = `Y_LT_X
 840                 hx = `H_LT_X
 841         | `Ordering_XHY:
 842                 yh = `H_LT_Y
 843                 yx = `X_LT_Y
 844                 hx = `X_LT_H
 845         | `Ordering_XYH:
 846                 yh = `Y_LT_H
 847                 yx = `X_LT_Y
 848                 hx = `X_LT_H
 849         ;;
 850
 851         var need_another_pass = true
 852         while need_another_pass
 853                 need_another_pass = false
 854                 for var j = fc#.len - 2; j >= 0 && j + 1 < fc#.len; --j
 855                         match (fc#[j + 0], fc#[j + 1], yh, yx, hx)
 856                         | (`X(xi, xval), `Y(yj, yval), _, `Y_LT_X, _):
 857                                 if xi != yj
 858                                         fc#[j + 0] = `Y(yj, yval)
 859                                         fc#[j + 1] = `X(xi, xval)
 860                                         need_another_pass = true
 861                                 else
 862                                         var one_plus_st : yakmo.ratfunc# = yakmo.rmul(xval, yval)
 863                                         yakmo.gadd_ip(one_plus_st, one)
 864                                         auto one_plus_st
 865                                         match yakmo.finv(one_plus_st)
 866                                         | `std.None: /* Ex 1? */
 867                                         | `std.Some opsti:
 868                                                 need_another_pass = true
 869                                                 var t_over = yakmo.rmul(yval, opsti)
 870                                                 yakmo.reduceratfunc(t_over)
 871                                                 fc#[j + 0] = `Y(yj, t_over)
 872                                                 var s_over = yakmo.rmul(xval, opsti)
 873                                                 fc#[j + 1] = `X(xi, s_over)
 874                                                 yakmo.reduceratfunc(s_over)
 875                                                 std.slput(fc, j + 1, `H(xi, t.dup(one_plus_st)))
 876                                                 __dispose__(xval)
 877                                                 __dispose__(yval)
 878                                         ;;
 879                                 ;;
 880                         | (`Y(yi, yval), `X(xj, xval), _, `X_LT_Y, _):
 881                                 if yi != xj
 882                                         fc#[j + 0] = `X(xj, xval)
 883                                         fc#[j + 1] = `Y(yi, yval)
 884                                         need_another_pass = true
 885                                 else
 886                                         var one_plus_st : yakmo.ratfunc# = yakmo.rmul(xval, yval)
 887                                         yakmo.gadd_ip(one_plus_st, one)
 888                                         auto one_plus_st
 889                                         match yakmo.finv(one_plus_st)
 890                                         | `std.None: /* Ex 1? */
 891                                         | `std.Some opsti:
 892                                                 need_another_pass = true
 893                                                 var t_over = yakmo.rmul(xval, opsti)
 894                                                 yakmo.reduceratfunc(t_over)
 895                                                 fc#[j + 0] = `X(xj, t_over)
 896                                                 var s_over = yakmo.rmul(yval, opsti)
 897                                                 fc#[j + 1] = `Y(yi, s_over)
 898                                                 yakmo.reduceratfunc(s_over)
 899                                                 std.slput(fc, j + 1, `H(xj, t.dup(opsti)))
 900                                                 __dispose__(xval)
 901                                                 __dispose__(yval)
 902                                         ;;
 903                                 ;;
 904                         | (`X(xj, xval), `H(hi, hval), _, _, `H_LT_X):
 905                                 need_another_pass = true
 906                                 /* A is 0-based, xj and hi are 1-based */
 907                                 var Aij = A[hi - 1][xj - 1]
 908                                 if Aij > 0
 909                                         match yakmo.finv(hval)
 910                                         | `std.None: -> `std.Err std.fmt("cannot invert {} while simplifying {}", hval, fc#)
 911                                         | `std.Some hvali:
 912                                                 for var k = 0; k < Aij; ++k
 913                                                         yakmo.rmul_ip(xval, hvali)
 914                                                 ;;
 915                                                 yakmo.reduceratfunc(xval)
 916                                                 __dispose__(hvali)
 917                                         ;;
 918                                 elif Aij < 0
 919                                         for var k = 0; k > Aij; --k
 920                                                 yakmo.rmul_ip(xval, hval)
 921                                         ;;
 922                                         yakmo.reduceratfunc(xval)
 923                                 ;;
 924                                 fc#[j + 0] = `H(hi, hval)
 925                                 fc#[j + 1] = `X(xj, xval)
 926                         | (`H(hi, hval), `X(xj, xval), _, _, `X_LT_H):
 927                                 need_another_pass = true
 928                                 /* A is 0-based, xj and hi are 1-based */
 929                                 var Aij = A[hi - 1][xj - 1]
 930                                 if Aij < 0
 931                                         match yakmo.finv(hval)
 932                                         | `std.None: -> `std.Err std.fmt("cannot invert {} while simplifying {}", hval, fc#)
 933                                         | `std.Some hvali:
 934                                                 for var k = 0; k > Aij; --k
 935                                                         yakmo.rmul_ip(xval, hvali)
 936                                                 ;;
 937                                                 yakmo.reduceratfunc(xval)
 938                                                 __dispose__(hvali)
 939                                         ;;
 940                                 elif Aij > 0
 941                                         for var k = 0; k < Aij; ++k
 942                                                 yakmo.rmul_ip(xval, hval)
 943                                         ;;
 944                                         yakmo.reduceratfunc(xval)
 945                                 ;;
 946                                 fc#[j + 0] = `X(xj, xval)
 947                                 fc#[j + 1] = `H(hi, hval)
 948                         | (`H(hi, hval), `Y(yj, yval), `Y_LT_H, _, _):
 949                                 need_another_pass = true
 950                                 var Aij = A[hi - 1][yj - 1]
 951                                 if Aij > 0
 952                                         match yakmo.finv(hval)
 953                                         | `std.None: -> `std.Err std.fmt("cannot invert {} while simplifying {}", hval, fc#)
 954                                         | `std.Some hvali:
 955                                                 for var k = 0; k < Aij; ++k
 956                                                         yakmo.rmul_ip(yval, hvali)
 957                                                 ;;
 958                                                 yakmo.reduceratfunc(yval)
 959                                                 __dispose__(hvali)
 960                                         ;;
 961                                 elif Aij < 0
 962                                         for var k = 0; k > Aij; --k
 963                                                 yakmo.rmul_ip(yval, hval)
 964                                         ;;
 965                                         yakmo.reduceratfunc(yval)
 966                                 ;;
 967                                 fc#[j + 0] = `Y(yj, yval)
 968                                 fc#[j + 1] = `H(hi, hval)
 969                         | (`Y(yj, yval), `H(hi, hval), `H_LT_Y, _, _):
 970                                 need_another_pass = true
 971                                 var Aij = A[hi - 1][yj - 1]
 972                                 if Aij < 0
 973                                         match yakmo.finv(hval)
 974                                         | `std.None: -> `std.Err std.fmt("cannot invert {} while simplifying {}", hval, fc#)
 975                                         | `std.Some hvali:
 976                                                 for var k = 0; k > Aij; --k
 977                                                         yakmo.rmul_ip(yval, hvali)
 978                                                 ;;
 979                                                 yakmo.reduceratfunc(yval)
 980                                                 __dispose__(hvali)
 981                                         ;;
 982                                 elif Aij > 0
 983                                         for var k = 0; k < Aij; ++k
 984                                                 yakmo.rmul_ip(yval, hval)
 985                                         ;;
 986                                         yakmo.reduceratfunc(yval)
 987                                 ;;
 988                                 fc#[j + 0] = `H(hi, hval)
 989                                 fc#[j + 1] = `Y(yj, yval)
 990                         | (`H(hi, hval1), `H(hj, hval2), _, _, _):
 991                                 if hi > hj
 992                                         fc#[j + 0] = `H(hj, hval2)
 993                                         fc#[j + 1] = `H(hi, hval1)
 994                                         need_another_pass = true
 995                                 elif hi == hj
 996                                         yakmo.rmul_ip(hval1, hval2)
 997                                         __dispose__(hval2)
 998                                         std.sldel(fc, j + 1)
 999                                         if std.eq(hval1.num, hval1.den)
1000                                                 __dispose__(hval1)
1001                                                 std.sldel(fc, j)
1002                                         ;;
1003                                         need_another_pass = true
1004                                 ;;
1005                         | (`X(xi, xval1), `X(xj, xval2), _, _, _):
1006                                 if xi == xj
1007                                         yakmo.gadd_ip(xval1, xval2)
1008                                         __dispose__(xval2)
1009                                         std.sldel(fc, j + 1)
1010                                         if yakmo.eq_gid(xval1)
1011                                                 __dispose__(xval1)
1012                                                 std.sldel(fc, j)
1013                                         ;;
1014                                         need_another_pass = true
1015                                 ;;
1016                         | (`Y(yi, yval1), `Y(yj, yval2), _, _, _):
1017                                 if yi == yj
1018                                         yakmo.gadd_ip(yval1, yval2)
1019                                         __dispose__(yval2)
1020                                         std.sldel(fc, j + 1)
1021                                         if yakmo.eq_gid(yval1)
1022                                                 __dispose__(yval1)
1023                                                 std.sldel(fc, j)
1024                                         ;;
1025                                         need_another_pass = true
1026                                 ;;
1027                         | (_, `H(hi, hval), _, _, _):
1028                                 if std.eq(hval.num, hval.den)
1029                                         __dispose__(hval)
1030                                         std.sldel(fc, j + 1)
1031                                         need_another_pass = true
1032                                 ;;
1033                         | (`H(hi, hval), _, _, _, _):
1034                                 if std.eq(hval.num, hval.den)
1035                                         __dispose__(hval)
1036                                         std.sldel(fc, j)
1037                                         need_another_pass = true
1038                                 ;;
1039                         | _:
1040                         ;;
1041
1042                         //for var j = 0; j < fc#.len; ++j
1043                         //        match fc#[j]
1044                         //        | `H (i, rf): yakmo.reduceratfunc(rf)
1045                         //        | `X (i, rf): yakmo.reduceratfunc(rf)
1046                         //        | `Y (i, rf): yakmo.reduceratfunc(rf)
1047                         //        ;;
1048                         //;;
1049                 ;;
1050         ;;
1051
1052         -> `std.Ok void
1053 }
1054
1055 const fill_in_alpha = {kc : KC_type, alpha : conf_3_star#
1056         var rank = get_n(kc)
1057         var w0 = get_w0(kc)
1058
1059         if alpha.u.len == w0.len
1060                 -> `std.Ok void
1061         elif alpha.u.len != w0.len - rank
1062                 -> `std.Err std.fmt("alpha.u is not complete, and it's not missing exactly r elements")
1063         ;;
1064
1065         var h1_fc
1066         match edge_fc_from_mc(kc, alpha.h1)
1067         | `std.Ok q: h1_fc = q
1068         | `std.Err e: -> `std.Err std.fmt("edge_fc_from_mc(h1): {}", e)
1069         ;;
1070         auto h1_fc
1071
1072         var h2_fc
1073         match edge_fc_from_mc(kc, alpha.h2)
1074         | `std.Ok q: h2_fc = q
1075         | `std.Err e: -> `std.Err std.fmt("edge_fc_from_mc(h2): {}", e)
1076         ;;
1077         auto h2_fc
1078
1079         var h3_fc
1080         match edge_fc_from_mc(kc, alpha.h3)
1081         | `std.Ok q: h3_fc = q
1082         | `std.Err e: -> `std.Err std.fmt("edge_fc_from_mc(h3): {}", e)
1083         ;;
1084         auto h3_fc
1085
1086         var h1_h3_fc : factorization_coordinate[:] = [][:]
1087         for var j = 0; j < h1_fc.len; ++j
1088                 match h1_fc[j]
1089                 | `H(i, v): std.slpush(&h1_h3_fc, `H(i, t.dup(v)))
1090                 | _: -> `std.Err std.fmt("h1 has non-H elements")
1091                 ;;
1092         ;;
1093         for var j = 0; j < h3_fc.len; ++j
1094                 match h3_fc[j]
1095                 | `H(i, v): std.slpush(&h1_h3_fc, `H(i, t.dup(v)))
1096                 | _: -> `std.Err std.fmt("h3 has non-H elements")
1097                 ;;
1098         ;;
1099         match canonicalize_factorization_coords(kc, &h1_h3_fc, `Ordering_YHX)
1100         | `std.Ok void:
1101         | `std.Err e: -> `std.Err std.fmt("canonicalize_factorization_coords( h1 h3 ): {}", e)
1102         ;;
1103         auto h1_h3_fc
1104
1105         var wo_h1_h3_fc
1106         match apply_wk(kc, 1, h1_h3_fc, 0)
1107         | `std.Ok q: wo_h1_h3_fc = q
1108         | `std.Err e: -> `std.Err std.fmt("w₀(h1 h3): {}", e)
1109         ;;
1110         auto wo_h1_h3_fc
1111
1112         var mash = [][:]
1113         for var j = 0; j < wo_h1_h3_fc.len; ++j
1114                 match wo_h1_h3_fc[j]
1115                 | `H(i, v): std.slpush(&mash, `H(i, t.dup(v)))
1116                 | _: -> `std.Err std.fmt("w₀(h1 h3) has non-H elements")
1117                 ;;
1118         ;;
1119         for var j = 0; j < h2_fc.len; ++j
1120                 match h2_fc[j]
1121                 | `H(i, v): std.slpush(&mash, `H(i, t.dup(v)))
1122                 | _: -> `std.Err std.fmt("h2 has non-H elements")
1123                 ;;
1124         ;;
1125         match canonicalize_factorization_coords(kc, &mash, `Ordering_YHX)
1126         | `std.Ok void:
1127         | `std.Err e: -> `std.Err std.fmt("canonicalize_factorization_coords( w₀(h1 h3) h2 ): {}", e)
1128         ;;
1129         auto mash
1130
1131         for var j = 0; j < mash.len; ++j
1132                 match mash[j]
1133                 | `H(i, v):
1134                         match yakmo.finv(v)
1135                         | `std.Some vi:
1136                                 yakmo.reduceratfunc(vi)
1137                                 __dispose__(v)
1138                                 mash[j] = `H(i, vi)
1139                         | `std.None: -> `std.Err std.fmt("cannot invert w₀(h1 h3) h2")
1140                         ;;
1141                 | _: -> `std.Err std.fmt("w₀(h1 h3) h2 contains non-H elements")
1142                 ;;
1143         ;;
1144
1145         match edge_mc_from_fc(kc, mash)
1146         | `std.Ok mc:
1147                 for var j = 0; j < alpha.u.len; ++j
1148                         std.slpush(&mc, alpha.u[j])
1149                 ;;
1150                 std.slfree(alpha.u)
1151                 alpha.u = mc
1152                 -> `std.Ok void
1153         | `std.Err e: -> `std.Err std.fmt("edge_mc_from_fc(mash): {}", e)
1154         ;;
1155 }
1156
1157 const fill_in_u_from_h01_h12_h20 = {kc : KC_type, u : minor_coordinate[:]#, h01 : minor_coordinate[:], h12 : minor_coordinate[:], h20 : minor_coordinate[:]
1158         var alpha = [
1159                 .h1 = h01,
1160                 .h2 = h12,
1161                 .h3 = h20,
1162                 .u = u#,
1163         ]
1164
1165         match fill_in_alpha(kc, &alpha)
1166         | `std.Ok void:
1167         | `std.Err e: -> `std.Err std.fmt("fill_in_alpha: {}", e)
1168         ;;
1169
1170         u# = alpha.u
1171
1172         -> `std.Ok void
1173 }
1174
1175 /* Given h in minor coordinates, compute w₀(h⁻¹) */
1176 const apply_wo_of_inv = {kc : KC_type, h : minor_coordinate[:]
1177         var h_fc
1178         match edge_fc_from_mc(kc, h)
1179         | `std.Err e: -> `std.Err std.fmt("edge_fc_from_mc({}, h): {}", kc, e)
1180         | `std.Ok y: h_fc = y
1181         ;;
1182         auto h_fc
1183
1184         var h_inv : factorization_coordinate[:] = std.slalloc(h_fc.len)
1185         for var k = 0; k < h_fc.len; ++k
1186                 match h_fc[k]
1187                 | `H(i, c):
1188                         match yakmo.finv(c)
1189                         | `std.Some ci: h_inv[k] = `H(i, ci)
1190                         | `std.None: -> `std.Err std.fmt("cannot invert h")
1191                         ;;
1192                 | _: -> `std.Err std.fmt("h contained non-H elements")
1193                 ;;
1194         ;;
1195         auto h_inv
1196
1197         var wo_h_inv
1198         match apply_wk(kc, 1, h_inv, 0)
1199         | `std.Ok y: wo_h_inv = y
1200         | `std.Err e: -> `std.Err std.fmt("apply_wk({}, 1, h_inv, 0): {}", e)
1201         ;;
1202         auto wo_h_inv
1203
1204         var wo_h_inv_mc
1205         match edge_mc_from_fc(kc, wo_h_inv)
1206         | `std.Ok y: wo_h_inv_mc = y
1207         | `std.Err e: -> `std.Err std.fmt("edge_mc_from_fc({}, wo(h^-1)): {}", kc, e)
1208         ;;
1209
1210         -> `std.Ok wo_h_inv_mc
1211 }
1212
1213 /* Given h in minor coordinates, compute h sG = h w₀ w₀ */
1214 const multiply_h_mc_by_sG = {kc : KC_type, h : minor_coordinate[:]
1215         var sG
1216         match get_sG(kc)
1217         | `std.Err e: -> `std.Err std.fmt("get_sG: {}", e)
1218         | `std.Ok sGG: sG = sGG
1219         ;;
1220
1221         var h_fc
1222         match edge_fc_from_mc(kc, h)
1223         | `std.Err e: -> `std.Err std.fmt("edge_fc_from_mc({}, h): {}", kc, e)
1224         | `std.Ok y: h_fc = y
1225         ;;
1226         for var k = 0; k < sG.len; ++k
1227                 std.slpush(&h_fc, t.dup(sG[k]))
1228         ;;
1229
1230         match canonicalize_factorization_coords(kc, &h_fc, `Ordering_YHX)
1231         | `std.Ok void:
1232         | `std.Err e: -> `std.Err std.fmt("canonicalize_factorization_coords: {}", e)
1233         ;;
1234         auto h_fc
1235
1236         match edge_mc_from_fc(kc, h_fc)
1237         | `std.Err e: -> `std.Err std.fmt("edge_mc_from_fc({}, h sG) {}", kc, e)
1238         | `std.Ok y: -> `std.Ok y
1239         ;;
1240 }
1241
1242 /* Return ( rot(α₁₂₀) = rot(g₁ N, g₂ N, g₀ sG N),  α₀₂₃ = (g₀ N, g₂ N, g₃) ) */
1243 const psi_02 = {kc : KC_type, alpha : conf_4_star#
1244         /*
1245            α₁₂₀ (in terms of α = (g₀ N, g₁ N, g₂ N, g₃ N) ) is (g₁ N, g₂ N, g₀ sG N)
1246             - h₁ = [w₀⁻¹ g₁⁻¹ g₂]₀ = h₁₂
1247             - h₂ = [w₀⁻¹ g₂⁻¹ g₀ sG]₀ = h₂₀ sG = w₀(h₀₂⁻¹) sG sG = w₀(h₀₂⁻¹)
1248             - h₃ = [w₀⁻¹ sG⁻¹ g₀⁻¹ g₁]₀ = sG h₀₁ = h₀₁ sG
1249             - u = u
1250          */
1251         var wo_h02_inv
1252         match apply_wo_of_inv(kc, alpha.h02)
1253         | `std.Err e: -> `std.Err std.fmt("w₀(h₀₂⁻¹): {}", e)
1254         | `std.Ok h: wo_h02_inv = h
1255         ;;
1256
1257         var h01_sG
1258         match multiply_h_mc_by_sG(kc, alpha.h01)
1259         | `std.Err e: -> `std.Err std.fmt("h₀₁ sG: {}", e)
1260         | `std.Ok h: h01_sG = h
1261         ;;
1262
1263         var alpha_120_base : conf_3_star = [
1264                 .h1 = t.dup(alpha.h12),
1265                 .h2 = wo_h02_inv,
1266                 .h3 = h01_sG,
1267                 .u = t.dup(alpha.u),
1268         ]
1269         var alpha_120 = auto std.mk(alpha_120_base)
1270
1271         var alpha_012 : conf_3_star#
1272         match rotate(kc, alpha_120)
1273         | `std.Err e: -> `std.Err std.fmt("rotate(α₁₂₀): {}", e)
1274         | `std.Ok a: alpha_012 = a
1275         ;;
1276
1277         /*
1278             α₀₂₃ (in terms of α as above):
1279              - h1 = [w₀⁻¹ g₀⁻¹ g₂]₀ = h₀₂
1280              - h2 = [w₀⁻¹ g₂⁻¹ g₃]₀ = h₂₃
1281              - h3 = [w₀⁻¹ g₃⁻¹ g₀]₀ = h₃₀ = w₀(h₀₃⁻¹) sG
1282              - u = v
1283          */
1284         var h30_sans_sG
1285         match apply_wo_of_inv(kc, alpha.h03)
1286         | `std.Err e: -> `std.Err std.fmt("w₀(h₀₃⁻¹): {}", e)
1287         | `std.Ok h: h30_sans_sG = h
1288         ;;
1289         auto h30_sans_sG
1290         var h30
1291         match multiply_h_mc_by_sG(kc, h30_sans_sG)
1292         | `std.Err e: -> `std.Err std.fmt("h₃₀ sG: {}", e)
1293         | `std.Ok h: h30 = h
1294         ;;
1295         var alpha_023 : conf_3_star = [
1296                 .h1 = t.dup(alpha.h02),
1297                 .h2 = t.dup(alpha.h23),
1298                 .h3 = h30,
1299                 .u = t.dup(alpha.v),
1300         ]
1301
1302         -> `std.Ok (alpha_012, std.mk(alpha_023))
1303 }
1304
1305 /* Return ( α₁₂₃ = (g₁ N, g₂ N, g₃), rot(α₁₃₀) = rot(g₁ sG N, g₃ N, g₀ N) ) */
1306 const psi_13 = {kc : KC_type, alpha : conf_4_star#
1307         var w0 = get_w0(kc)
1308         if alpha.u.len != w0.len || alpha.v.len != w0.len
1309                 -> `std.Err std.fmt("u ({}) and/or v ({}) are too short ({})", alpha.u.len, alpha.v.len, w0.len)
1310         ;;
1311
1312         /* Compute the components for ι: k = w₀(h₀₂ h₁₂⁻¹), Φ⁻¹(v) and u. */
1313
1314         /* v */
1315         var v_fc
1316         match interior_fc_from_mc_and_undo_psi_phi_psi(kc, alpha.v)
1317         | `std.Ok q: v_fc = q
1318         | `std.Err e: -> `std.Err std.fmt("interior_fc_from_mc_and_undo_ΨΦΨ(v): {}", e)
1319         ;;
1320
1321         var phi_inverse_v
1322         match phi_inv(kc, v_fc)
1323         | `std.Err e: -> `std.Err std.fmt("Φ⁻¹(v): {}", e)
1324         | `std.Ok q: phi_inverse_v = q
1325         ;;
1326         auto phi_inverse_v
1327
1328         /* u */
1329         var u_fc
1330         match interior_fc_from_mc_and_undo_psi_phi_psi(kc, alpha.u)
1331         | `std.Ok q: u_fc = q
1332         | `std.Err e: -> `std.Err std.fmt("interior_fc_from_mc_and_undo_ΨΦΨ(u): {}", e)
1333         ;;
1334         auto u_fc
1335
1336         /* k */
1337         var h02_fc
1338         match edge_fc_from_mc(kc, alpha.h02)
1339         | `std.Ok fc: h02_fc = fc
1340         | `std.Err e: -> `std.Err std.fmt("h₀₂ broken: {}", e)
1341         ;;
1342         auto h02_fc
1343
1344         var h12_fc
1345         match edge_fc_from_mc(kc, alpha.h12)
1346         | `std.Ok fc: h12_fc = fc
1347         | `std.Err e: -> `std.Err std.fmt("h₁₂ broken: {}", e)
1348         ;;
1349         auto h12_fc
1350
1351         var h12_h02i = [][:]
1352         for var k = 0; k < h12_fc.len; ++k
1353                 std.slpush(&h12_h02i, t.dup(h12_fc[k]))
1354         ;;
1355
1356         for var k = 0; k < h02_fc.len; ++k
1357                 match h02_fc[k]
1358                 | `H(j, v):
1359                         match yakmo.finv(v)
1360                         | `std.Some vi: std.slpush(&h12_h02i, `H(j, vi))
1361                         | `std.None: -> `std.Err std.fmt("cannot invert h₀₂")
1362                         ;;
1363                 | _: -> `std.Err std.fmt("h₀₂ contains non-H elements")
1364                 ;;
1365         ;;
1366
1367         match canonicalize_factorization_coords(kc, &h12_h02i, `Ordering_YHX)
1368         | `std.Ok void:
1369         | `std.Err e: -> `std.Err std.fmt("canonicalize(h₁₂ h₀₂⁻¹): {}", e)
1370         ;;
1371         auto h12_h02i
1372
1373         var k
1374         match apply_wk(kc, 1, h12_h02i, 0)
1375         | `std.Ok kk: k = kk
1376         | `std.Err e: -> `std.Err std.fmt("w₀(h₁₂ h₀₂⁻¹): {}", e)
1377         ;;
1378
1379         /* The results of ι: c, d, l */
1380         var phi_inverse_c, d, l
1381         match iota(kc, u_fc, k, phi_inverse_v)
1382         | `std.Ok ii: (phi_inverse_c, d, l) = ii
1383         | `std.Err e: -> `std.Err std.fmt("iota(u, k, Φ⁻¹(v)): {}", e)
1384         ;;
1385         auto phi_inverse_c
1386         auto d
1387         auto l
1388
1389         var c
1390         match phi(kc, phi_inverse_c)
1391         | `std.Ok q: c = q
1392         | `std.Err e: -> `std.Err std.fmt("Φ(phi_inverse_c): {}", e)
1393         ;;
1394         auto c
1395
1396         var c_mc
1397         match interior_mc_from_fc_and_do_psi_phi_psi(kc, c)
1398         | `std.Ok q: c_mc = q
1399         | `std.Err e: -> `std.Err std.fmt("interior_mc_from_fc_and_do_psi_phi_psi(c): {}", e)
1400         ;;
1401
1402         var d_mc
1403         match interior_mc_from_fc_and_do_psi_phi_psi(kc, d)
1404         | `std.Ok q: d_mc = q
1405         | `std.Err e: -> `std.Err std.fmt("interior_mc_from_fc_and_do_psi_phi_psi(d): {}", e)
1406         ;;
1407
1408         /*
1409            α₁₂₃ (in terms of α = (g₀ N, g₁ N, g₂ N, g₃ N) ) is (g₁ N, g₂ N, g₃ N)
1410             - h₁ = [w₀⁻¹ g₁⁻¹ g₂]₀ = h₁₂
1411             - h₂ = [w₀⁻¹ g₂⁻¹ g₃]₀ = h₂₃
1412             - h₃ = [w₀⁻¹ g₃⁻¹ g₁]₀ = h₃₁ = w₀(h₁₃⁻¹) sG
1413             - u = c
1414          */
1415         var wo_l
1416         match apply_wk(kc, 1, l, 0)
1417         | `std.Ok q: wo_l = q
1418         | `std.Err e: -> `std.Err std.fmt("w₀(l): {}", e)
1419         ;;
1420         auto wo_l
1421
1422         var h03_fc
1423         match edge_fc_from_mc(kc, alpha.h03)
1424         | `std.Ok q: h03_fc = q
1425         | `std.Err e: -> `std.Err std.fmt("edge_fc_from_mc(h₀₃): {}", e)
1426         ;;
1427
1428         var h13_fc = [][:]
1429         for var k = 0; k < wo_l.len; ++k
1430                 match wo_l[k]
1431                 | `H(j, v): std.slpush(&h13_fc, `H(j, t.dup(v)))
1432                 | _: -> `std.Err std.fmt("w₀(l) contained non-H")
1433                 ;;
1434         ;;
1435         for var k = 0; k < h03_fc.len; ++k
1436                 match h03_fc[k]
1437                 | `H(j, v): std.slpush(&h13_fc, `H(j, t.dup(v)))
1438                 | _: -> `std.Err std.fmt("h₀₃ contained non-H")
1439                 ;;
1440         ;;
1441         match canonicalize_factorization_coords(kc, &h13_fc, `Ordering_YHX)
1442         | `std.Ok void:
1443         | `std.Err e: -> `std.Err std.fmt("h₁₃ is broken")
1444         ;;
1445         auto h13_fc
1446
1447         var h13
1448         match edge_mc_from_fc(kc, h13_fc)
1449         | `std.Ok q: h13 = q
1450         | `std.Err e: -> `std.Err std.fmt("edge_mc_from_fc(h₁₃): {}", e)
1451         ;;
1452         auto h13
1453
1454         for var j = 0; j < h13.len; ++j
1455                 yakmo.reduceratfunc(h13[j].c)
1456         ;;
1457
1458         var wo_h13i
1459         match apply_wo_of_inv(kc, h13)
1460         | `std.Ok q: wo_h13i = q
1461         | `std.Err e: -> `std.Err std.fmt("w₀(h₁₃⁻¹): {}", e)
1462         ;;
1463         auto wo_h13i
1464
1465         var h31
1466         match multiply_h_mc_by_sG(kc, wo_h13i)
1467         | `std.Ok q: h31 = q
1468         | `std.Err e: -> `std.Err std.fmt("w₀(h₁₃⁻¹) sG: {}", e)
1469         ;;
1470
1471         for var j = 0; j < h31.len; ++j
1472                 yakmo.reduceratfunc(h31[j].c)
1473         ;;
1474
1475         var alpha_123 : conf_3_star = [
1476                 .h1 = t.dup(alpha.h12),
1477                 .h2 = t.dup(alpha.h23),
1478                 .h3 = h31,
1479                 .u = c_mc,
1480         ]
1481
1482         /*
1483            α₁₃₀ (in terms of α) is (g₁ sG N, g₃ N, g₀ N)
1484             - h₁ = [w₀⁻¹ sG g₁⁻¹ g₃]₀ = h₁₃ sG
1485             - h₂ = [w₀⁻¹ g₃⁻¹ g₀]₀ = h₃₀
1486             - h₃ = [w₀⁻¹ g₀⁻¹ g₁]₀ = h₀₁
1487             - u = d
1488          */
1489         var h13_sG
1490         match multiply_h_mc_by_sG(kc, h13)
1491         | `std.Ok q: h13_sG = q
1492         | `std.Err e: -> `std.Err std.fmt("h₁₃ sG: {}", e)
1493         ;;
1494
1495         var wo_h03i
1496         match apply_wo_of_inv(kc, alpha.h03)
1497         | `std.Err e: -> `std.Err std.fmt("w₀(h₀₃⁻¹): {}", e)
1498         | `std.Ok h: wo_h03i = h
1499         ;;
1500         auto wo_h03i
1501
1502         var h30
1503         match multiply_h_mc_by_sG(kc, wo_h03i)
1504         | `std.Err e: -> `std.Err std.fmt("w₀(h₀₃⁻¹) sG: {}", e)
1505         | `std.Ok h: h30 = h
1506         ;;
1507
1508         var h01_sG
1509         match multiply_h_mc_by_sG(kc, alpha.h01)
1510         | `std.Ok q: h01_sG = q
1511         | `std.Err e: -> `std.Err std.fmt("h₀₁ sG: {}", e)
1512         ;;
1513
1514         var alpha_130_base : conf_3_star = [
1515                 .h1 = h13_sG,
1516                 .h2 = h30,
1517                 .h3 = h01_sG,
1518                 .u = d_mc,
1519         ]
1520         var alpha_130 = auto std.mk(alpha_130_base)
1521
1522         var alpha_013 : conf_3_star#
1523         match rotate(kc, alpha_130)
1524         | `std.Err e: -> `std.Err std.fmt("rotate(α₁₃₀): {}", e)
1525         | `std.Ok a: alpha_013 = a
1526         ;;
1527
1528         -> `std.Ok (std.mk(alpha_123), alpha_013)
1529 }
1530
1531 /*
1532    Like interior_fc_from_mc, but undo ΨΦΨ as well.
1533
1534    This isn't the default because Zickert's paper treats ΨΦΨ as a
1535    primitive.
1536  */
1537 const interior_fc_from_mc_and_undo_psi_phi_psi = {kc, y
1538         var psi_phi_psi_y_fc
1539         match interior_fc_from_mc(kc, y)
1540         | `std.Ok q: psi_phi_psi_y_fc = q
1541         | `std.Err e: -> `std.Err std.fmt("interior_fc_from_mc(y): {}", e)
1542         ;;
1543         auto psi_phi_psi_y_fc
1544
1545         var phi_psi_y_fc
1546         match psi(kc, psi_phi_psi_y_fc)
1547         | `std.Ok q: phi_psi_y_fc = q
1548         | `std.Err e: -> `std.Err std.fmt("Ψ( ΨΦΨ(y) ): {}", e)
1549         ;;
1550         auto phi_psi_y_fc
1551
1552         var psi_y_fc
1553         match phi_inv(kc, phi_psi_y_fc)
1554         | `std.Ok q: psi_y_fc = q
1555         | `std.Err e: -> `std.Err std.fmt("Φ⁻¹( ΦΨ(y) ): {}", e)
1556         ;;
1557         auto psi_y_fc
1558
1559         var y_fc
1560         match psi(kc, psi_y_fc)
1561         | `std.Ok q: y_fc = q
1562         | `std.Err e: -> `std.Err std.fmt("Ψ( Ψ(y) ): {}", e)
1563         ;;
1564
1565         -> `std.Ok y_fc
1566 }
1567
1568 /*
1569    Like interior_mc_from_fc, but undo (ΨΦΨ)⁻¹ as well.
1570
1571    This isn't the default because Zickert's paper treats ΨΦΨ as a
1572    primitive.
1573  */
1574 const interior_mc_from_fc_and_do_psi_phi_psi = {kc, y
1575         var psi_y
1576         match psi(kc, y)
1577         | `std.Ok q: psi_y = q
1578         | `std.Err e: -> `std.Err std.fmt("Ψ(y): {}", e)
1579         ;;
1580         auto psi_y
1581
1582         var phi_psi_y
1583         match phi(kc, psi_y)
1584         | `std.Ok q: phi_psi_y = q
1585         | `std.Err e: -> `std.Err std.fmt("Φ( Ψ(y) ): {}", e)
1586         ;;
1587         auto phi_psi_y
1588
1589         var psi_phi_psi_y
1590         match psi(kc, phi_psi_y)
1591         | `std.Ok q: psi_phi_psi_y = q
1592         | `std.Err e: -> `std.Err std.fmt("Ψ( ΦΨ(y) ): {}", e)
1593         ;;
1594         auto psi_phi_psi_y
1595
1596         var y_mc
1597         match interior_mc_from_fc(kc, psi_phi_psi_y)
1598         | `std.Ok q: y_mc = q
1599         | `std.Err e: -> `std.Err std.fmt("interior_mc_from_fc(ΨΦΨ(y)): {}", e)
1600         ;;
1601
1602         -> `std.Ok y_mc
1603 }
1604
1605 const rotate = {kc : KC_type, alpha : conf_3_star#
1606         var c = std.try(get_coxeter_elt(kc))
1607         var sigmaG : weyl_reflection[:] = [][:]
1608         match get_sigmaG_permutation(kc)
1609         | `std.Err e: -> `std.Err std.fmt("cannot compute σG for {}: {}", kc, e)
1610         | `std.Ok perm: sigmaG = perm
1611         ;;
1612
1613         /*
1614            For the hi, factorization coordinates are exactly minor
1615            coordinates. For u, there's some work to do.
1616          */
1617         var fc_h1 : factorization_coordinate[:] = [][:]
1618         var fc_h2 : factorization_coordinate[:] = [][:]
1619         var fc_h3 : factorization_coordinate[:] = [][:]
1620         var fc_u  : factorization_coordinate[:] = [][:]
1621
1622         match edge_fc_from_mc(kc, alpha.h1)
1623         | `std.Err e:
1624                 auto (e : t.doomed_str)
1625                 -> `std.Err std.fmt("edge_fc_from_mc(alpha.h1 = {}): {}", alpha.h1, e)
1626         | `std.Ok y: fc_h1 = y
1627         ;;
1628         canonicalize_factorization_coords(kc, &fc_h1, `Ordering_YHX)
1629         auto fc_h1
1630
1631         match edge_fc_from_mc(kc, alpha.h2)
1632         | `std.Err e:
1633                 auto (e : t.doomed_str)
1634                 -> `std.Err std.fmt("edge_fc_from_mc(alpha.h2 = {}): {}", alpha.h2, e)
1635         | `std.Ok y: fc_h2 = y
1636         ;;
1637         canonicalize_factorization_coords(kc, &fc_h2, `Ordering_YHX)
1638         auto fc_h2
1639
1640         match edge_fc_from_mc(kc, alpha.h3)
1641         | `std.Err e:
1642                 auto (e : t.doomed_str)
1643                 -> `std.Err std.fmt("edge_fc_from_mc(alpha.h3 = {}): {}", alpha.h3, e)
1644         | `std.Ok y: fc_h3 = y
1645         ;;
1646         canonicalize_factorization_coords(kc, &fc_h3, `Ordering_YHX)
1647         auto fc_h3
1648
1649         match interior_fc_from_mc(kc, alpha.u)
1650         | `std.Err e:
1651                 auto (e : t.doomed_str)
1652                 -> `std.Err std.fmt("interior_fc_from_mc(alpha.u = {}): {}", alpha.u, e)
1653         | `std.Ok y: fc_u = y
1654         ;;
1655         auto fc_u
1656
1657 for var j = 0; j < fc_u.len; ++j
1658         match fc_u[j]
1659         | `H (i, rf): yakmo.reduceratfunc(rf)
1660         | `X (i, rf): yakmo.reduceratfunc(rf)
1661         | `Y (i, rf): yakmo.reduceratfunc(rf)
1662         ;;
1663 ;;
1664
1665         /* Demand that h1, h2, h3 are in correct, increasing order */
1666         var rank = get_n(kc)
1667         if fc_h1.len != fc_h2.len || fc_h2.len != fc_h3.len || fc_h3.len != rank
1668                 -> `std.Err std.fmt("alpha.hi do not all have the same length")
1669         ;;
1670
1671         for var j = 0; j < fc_h1.len; ++j
1672                 match (fc_h1[j], fc_h2[j], fc_h3[j])
1673                 | (`H(i1, _), `H(i2, _), `H(i3, _)):
1674                         if i1 != j + 1 || i2 != j + 1 || i3 != j + 1
1675                                 -> `std.Err std.fmt("alpha.hi are not in order")
1676                         ;;
1677                 | _: -> `std.Err std.fmt("alpha.hi are not elements of maximal torus (j = {})", j)
1678                 ;;
1679         ;;
1680
1681         /*
1682            First, construct (Φ Ψ)(u). Since our tranlation between
1683            minor and factorization coordinates necessarily conjugates by
1684            ΨΦΨ, we actually compute (Ψ Φ) and trust the conjugation to
1685            take care of it.
1686          */
1687         var half_rotated_fc_u : factorization_coordinate[:] = [][:]
1688         match phi(kc, fc_u)
1689         | `std.Err e:
1690                 auto (e : t.doomed_str)
1691                 -> `std.Err std.fmt("Ψ(u): {}", e)
1692         | `std.Ok phi_u:
1693                 match psi(kc, phi_u)
1694                 | `std.Err e:
1695                         auto (e : t.doomed_str)
1696                         -> `std.Err std.fmt("(ΦΨ)[(ΨΦΨ)(u)]: {}", e)
1697                 | `std.Ok psi_phi_u: half_rotated_fc_u = psi_phi_u
1698                 ;;
1699                 auto phi_u
1700         ;;
1701 for var j = 0; j < half_rotated_fc_u.len; ++j
1702         match half_rotated_fc_u[j]
1703         | `X(_, rf): yakmo.reduceratfunc(rf)
1704         | `Y(_, rf): yakmo.reduceratfunc(rf)
1705         | `H(_, rf): yakmo.reduceratfunc(rf)
1706         ;;
1707 ;;
1708
1709         /*
1710            At this point, conjugate by [w0(h1) h2]^{-1}. Specifically,
1711            jam all the `H(i, ___)s that make up this element onto the
1712            left of half_rotated_fc_u, their inverses on the right, call
1713            canonicalize_factorization_coords. All those Hs had better
1714            bubble through the Ys and annihilate themselves.
1715          */
1716         match psi(kc, half_rotated_fc_u)
1717         | `std.Err e:
1718                 auto (e : t.doomed_str)
1719                 -> `std.Err std.fmt("Ψ(half_rotated_u): {}", e)
1720         | `std.Ok x:
1721                 __dispose__(half_rotated_fc_u)
1722                 half_rotated_fc_u = x
1723         ;;
1724
1725         /*
1726            Pull out all the ratfunc#s from fc_hi. They had better be in
1727            the right order, because we checked above.
1728          */
1729         var h1_rf : yakmo.ratfunc#[:] = std.slalloc((rank : std.size))
1730         var h2_rf : yakmo.ratfunc#[:] = std.slalloc((rank : std.size))
1731         for var j = 0; j < rank; ++j
1732                 match fc_h1[j]
1733                 | `H(_, rf): h1_rf[j] = rf
1734                 | _:
1735                 ;;
1736
1737                 match fc_h2[j]
1738                 | `H(_, rf): h2_rf[j] = rf
1739                 | _:
1740                 ;;
1741         ;;
1742
1743         var pre_conjugation_length = half_rotated_fc_u.len
1744         for var j = 0; j < fc_h1.len; ++j
1745                 /* h1_{σG[j]}^{+1} goes on the left, ^{-1} on the right */
1746                 var sigmaj = sigmaG[j + 1] - 1
1747                 std.slput(&half_rotated_fc_u, 0, `H(j + 1, t.dup(h1_rf[sigmaj])))
1748                 match yakmo.finv(h1_rf[sigmaj])
1749                 | `std.Some v: std.slpush(&half_rotated_fc_u, `H(j +  1, v))
1750                 | `std.None: -> `std.Err std.fmt("cannot conjugate by [w0(h1) h2]^{{-1}}: cannot invert {}", fc_h1[sigmaj])
1751                 ;;
1752
1753                 /* h2_j^{-1} goes on the left, ^{+1} on the right */
1754                 std.slpush(&half_rotated_fc_u, `H(j + 1, t.dup(h2_rf[j])))
1755                 match yakmo.finv(h2_rf[j])
1756                 | `std.Some v: std.slput(&half_rotated_fc_u, 0, `H(j + 1, v))
1757                 | `std.None: -> `std.Err std.fmt("cannot conjugate by [w0(h1) h2]^{{-1}}: cannot invert {}", fc_h2[j])
1758                 ;;
1759         ;;
1760         canonicalize_factorization_coords(kc, &half_rotated_fc_u, `Ordering_YHX)
1761         var post_conjugation_length = half_rotated_fc_u.len
1762         if pre_conjugation_length != post_conjugation_length
1763                 -> `std.Err std.fmt("conjugating by [w0(h1) h2]^{{-1}} did not go as expected: result is {}", half_rotated_fc_u)
1764         ;;
1765         std.slfree(h1_rf)
1766         std.slfree(h2_rf)
1767
1768
1769         match psi(kc, half_rotated_fc_u)
1770         | `std.Err e:
1771                 auto (e : t.doomed_str)
1772                 -> `std.Err std.fmt("Ψ(half_rotated_u): {}", e)
1773         | `std.Ok x:
1774                 __dispose__(half_rotated_fc_u)
1775                 half_rotated_fc_u = x
1776         ;;
1777
1778         /* Finish off the rotation with another hit of ΨΦ */
1779         for var j = 0; j < half_rotated_fc_u.len; ++j
1780                 match half_rotated_fc_u[j]
1781                 | `X(_, rf): yakmo.reduceratfunc(rf)
1782                 | `Y(_, rf): yakmo.reduceratfunc(rf)
1783                 | `H(_, rf): yakmo.reduceratfunc(rf)
1784                 ;;
1785         ;;
1786         var full_rotated_fc_u : factorization_coordinate[:] = [][:]
1787         match phi(kc, half_rotated_fc_u)
1788         | `std.Err e:
1789                 auto (e : t.doomed_str)
1790                 -> `std.Err std.fmt("Ψ(half_rotated_u): {}", e)
1791         | `std.Ok round_2_phi_u:
1792                 for var j = 0; j < round_2_phi_u.len; ++j
1793                         match round_2_phi_u[j]
1794                         | `X(_, rf): yakmo.reduceratfunc(rf)
1795                         | `Y(_, rf): yakmo.reduceratfunc(rf)
1796                         | `H(_, rf): yakmo.reduceratfunc(rf)
1797                         ;;
1798                 ;;
1799                 match psi(kc, round_2_phi_u)
1800                 | `std.Err e:
1801                         auto (e : t.doomed_str)
1802                         -> `std.Err std.fmt("(ΦΨ)[half_rotated_u]: {}", e)
1803                 | `std.Ok psi_phi_u: full_rotated_fc_u = psi_phi_u
1804                 ;;
1805                 auto round_2_phi_u
1806         ;;
1807         auto half_rotated_fc_u
1808
1809         /* Now, take new_u back to minor coordinates */
1810         var rotated_mc_u = [][:]
1811         match interior_mc_from_fc(kc, full_rotated_fc_u)
1812         | `std.Err e:
1813                 auto (e : t.doomed_str)
1814                 -> `std.Err std.fmt("interior_mc_from_fc({}, full_rotated_u = {}): {}", kc, full_rotated_fc_u, e)
1815         | `std.Ok m:
1816                 for var j = 0; j < m.len; ++j
1817                         yakmo.reduceratfunc(m[j].c)
1818                 ;;
1819                 rotated_mc_u = m
1820                 __dispose__(full_rotated_fc_u)
1821         ;;
1822
1823         var ret_h1 = std.sldup(alpha.h3)
1824         var ret_h2 = std.sldup(alpha.h1)
1825         var ret_h3 = std.sldup(alpha.h3)
1826         for var k = 0; k < alpha.h1.len; ++k
1827                 ret_h1[k] = [ .k = alpha.h3[k].k, .c = t.dup(alpha.h3[k].c) ]
1828                 ret_h2[k] = [ .k = alpha.h3[k].k, .c = t.dup(alpha.h1[k].c) ]
1829                 ret_h3[k] = [ .k = alpha.h3[k].k, .c = t.dup(alpha.h2[k].c) ]
1830         ;;
1831         var retbase : conf_3_star = [
1832                 .h1 = ret_h1,
1833                 .h2 = ret_h2,
1834                 .h3 = ret_h3,
1835                 .u = rotated_mc_u,
1836         ]
1837
1838         -> `std.Ok std.mk(retbase)
1839 }