type Face = number[]
type XYZ = [number, number, number]
type Polyhedron = {face: Face[], xyz: XYZ[], name?: string}

type Notation = string

// Useful constants
rt2 := Math.sqrt 2
rth := rt2/2
phi := (1 + Math.sqrt 5)/2
ihp := 1/phi
tau := 2*Math.PI

// Sadly needs to be early because we initialize the Dodecahedron as the dual of
// the icosahedron:
// Only add one direction of each, will auto reverse as well
specialDuals: Record<string, string> :=
   Tetrahedron: 'Tetrahedron',
   Cube: 'Octahdron',
   Dodecahedron: 'Icosahedron',
   Cuboctahedron: 'Rhombic dodecahedron',
   'truncated Tetrahedron': 'triakis Tetrahedron',
   'truncated Cube': 'triakis Octahedron',
   'truncated Octahedron': 'tetrakis Cube',
   Rhombicuboctahedron: 'Deltoidal icositetrahedron',
   'truncated Cuboctahedron': 'Disdyakis dodecahedron',
   'snub Cube': 'pentagonal icositetrahedron',
   Icosidodecahedron: 'Rhombic triacontahedron',
   'truncated Dodecahedron': 'triakis Icosahedron',
   'truncated Icosahedron': 'pentakis Dodecahedron',
   Rhombicosidodecahedron: 'Deltoidal hexecontahedron',
   'truncated Icosidodecahedron': 'Disdyakis triacontahedron',
   'snub Dodecahedron': 'pentagonal hexecontahedron'

// All seeds are canonical, i.e., edges tangent to unit sphere
icosahedron: Polyhedron :=
   face: [[0,1,9], [0,8,1], [0,4,8], [0,5,4], [0,9,5],
          [4,5,2], [4,2,10], [4,10,8], [8,10,6], [8,6,1],
          [1,6,7], [1,7,9], [9,7,11], [9,11,5], [5,11,2],
          [3,2,11], [3,10,2], [3,6,10], [3,7,6], [3,11,7]]
   xyz: [[0,ihp,1], [0,-ihp,1], [0,ihp,-1], [0,-ihp,-1],
         [ihp,1,0], [-ihp,1,0], [ihp,-1,0], [-ihp,-1,0],
         [1,0,ihp], [-1,0,ihp], [1,0,-ihp], [-1,0,-ihp]]
   name: 'Icosahedron'

polyCache: Record<Notation, Polyhedron> :=
   '': face: [], xyz: []
   T:
      face: [[0,2,1], [0,3,2], [0,1,3], [1,2,3]]
      xyz: [[1,1,1], [1,-1,-1], [-1,1,-1], [-1,-1,1]]
      name: 'Tetrahedron'
   O:
      face: [[0,2,1], [0,3,2], [0,4,3], [0,1,4],
             [1,5,4], [1,2,5], [2,3,5], [3,4,5]]
      xyz: [[0,0,rt2], [rt2,0,0], [0,rt2,0], [-rt2,0,0], [0,-rt2,0], [0,0,-rt2]]
      name: 'Octahedron'
   C:
      face: [[0,3,2,1], [0,5,4,3], [0,1,6,5], [1,2,7,6], [2,3,4,7], [4,5,6,7]]
      xyz: [[rth,rth,rth], [-rth,rth,rth], [-rth,-rth,rth], [rth,-rth,rth],
            [rth,-rth,-rth], [rth,rth,-rth], [-rth,rth,-rth], [-rth,-rth,-rth]]
      name: 'Cube'
   I: icosahedron
   D: geomDual icosahedron

export function generateVRML(notation: Notation): string
   outputVRML notation, generatePoly notation

function generatePoly(notation: Notation): Polyhedron
   getStandardPoly inform standardize notation

function getStandardPoly(notation: Notation): Polyhedron
   if notation in polyCache then return polyCache[notation]
   [ , op, arg, rest] := notation.match(/^(.)(\d*)(.*)$/) or []
   parent := getStandardPoly rest
   // may have created what we want by side effect
   if notation in polyCache then return polyCache[notation]
   dispatch op, arg, parent, notation  // will do the caching

// Convenience tuple maker
function ð<Tup extends unknown[]>(...arg: Tup): Tup arg

// Note we now allow numeric arguments on all of the basic operations,
// kis/truncate, join/ambo, and gyro/snub. Likely some of the operations
// we are taking as composite could have reasonable numeric versions, but
// there didn't seem to be any sensible way to propagate such an argument
// to the operations in their rewrites. In other words, the numeric-limited
// operations may not be composite, or at least not in the same way. So
// we have just left them as applying throughout the polyhedron.
rawStandardizations :=
   P4$: 'C', A3$: 'O', Y3$: 'T',                       // Seed synonyms
   e: 'aa', b: 'ta', o: 'jj', m: 'kj',                 // abbreviations
   [String.raw`t(\d*)`]: 'dk$1d',
   [String.raw`a(\d*)`]: 'dj$1d',                      // dual operations
   [String.raw`s(\d*)`]: 'dg$1d',
   dd: '', rr: '', jd: 'j', gd: 'rgr',                 // absorption rules
   rd: 'dr',  // these commute; others? If so, move 'r' in to cancel w/ seed
   // Remainder are all simplifications/unique selections for seeds:
   jY: 'dA', dT: 'T', gT: 'D', jT: 'C', dC: 'O', dO: 'C',
   dI: 'D', dD: 'I', rO: 'O', rC: 'C', rI: 'I', rD: 'D',
   jO: 'jC', jI: 'jD', gO: 'gC', gI: 'gD'
standardizations :=
   (ð RegExp(pat, 'g'), rep for pat, rep in rawStandardizations)

function standardize(notation: Notation): Notation
   lastNotation .= ''
   while lastNotation != notation  // iterate in case of rdrd, e.g.
      lastNotation = notation
      for [pat, rep] of standardizations
         notation = notation.replace(pat, rep)
   notation

function orb(r: number, n: number,
             height: number, t = 0): XYZ[]
   // A regular n-gon inscribed in a horizontal circle of radius r
   // at given height, rotated t*tau/n from standard position.
   theta := tau/n
   rho := t*theta
   for i of [0...n]
      [r*Math.cos(rho + i*theta), r*Math.sin(rho + i*theta), height] as XYZ

operator add(a: XYZ, b: XYZ)
   accumulate copy(a), b

ngonalNames := 3: 'triangular', 4: 'square', 5: 'pentagonal', 6: 'hexagonal',
   7: 'heptagonal', 8: 'octagonal', 9: 'enneagonal', 10: 'decagonal',
   12: 'dodecagonal'
type SpecialNgon = keyof typeof ngonalNames
function ngonal(n: number)
   ngonalNames[n as SpecialNgon] ?? n+`-gonal`

seeds :=
   P: (n: number) => // Prism
         unless n then n = 3
         theta := tau/n
         halfEdge := Math.sin theta/2
         xyz := orb(1, n, halfEdge) ++ orb(1, n, -halfEdge)
         face := [[n-1..0], [n...2*n]] // top and bottom
         for i of [0...n]
            ip1 := (i+1)%n  // next vertex around
            face.push [i, ip1, ip1+n, i+n]
         {face, xyz, name: ngonal(n)+` prism`}
   A: (n: number) => // Antiprism
         unless n then n = 4
         theta := tau/n
         halfHeight .= Math.sqrt
            1 - 4/(4 + 2*Math.cos(theta/2) - 2*Math.cos(theta))
         faceRadius .= Math.sqrt 1-halfHeight*halfHeight
         // Scale to put edge midpoints on unit sphere
         f := mag [halfHeight, faceRadius*Math.cos(theta/2), 0]
         halfHeight /= f
         faceRadius /= f
         xyz := orb(faceRadius, n, halfHeight)
	        ++ orb(faceRadius, n, -halfHeight, 0.5)
         face := [[n-1..0], [n...2*n]] // top and bottom
         for i of [0...n]
            face.push [i, (i+1)%n, i+n]
            face.push [i, i+n, n + (n+i-1)%n]
         {face, xyz, name: ngonal(n)+` antiprism`}
   Y: (n: number) => // pYramid
         unless n then n = 4
         // Canonical solution by Intelligenti Pauca and Ed Pegg, see
         // https://math.stackexchange.com/questions/2286628/canonical-pyramid-polynomials
         theta := tau/n
         c := Math.cos theta/2
         baseRadius := Math.sqrt 2/(c*(1+c))
         depth := Math.sqrt (1-c)/(1+c)
         height := 2*Math.sqrt 1/(1 - c*c)
         xyz := orb baseRadius, n, depth
         edgeMid2 := xyz[0] add xyz[1]
         xyz.push [0, 0, depth-height]
         face := ([i, (i+1)%n, n] for i of [0...n])
         face.unshift [n-1..0]
         {face, xyz, name: ngonal(n)+` pyramid`}
type SeedOp = keyof typeof seeds

// Syntactic sugar to deal with weird TypeScript typing:
operator æ<T>(A: T[], i: number) A.at(i) as T
function wordsof(s: string) 1 + (s.match(/\s+/g)?.length ?? 0)
specialJoins :=
   Tetrahedron: 'Cube',
   Cube: 'Rhombic dodecahedron', Octahedron: 'Rhombic dodecahedron',
   Dodecahedron: 'Rhombic triacontahedron', Icosahedron: 'Rhombic triacontahedron',
   'Rhombic dodecahedron': 'Deltoidal icositetrahedron',
   'Rhombic triacontahedron': 'Deltoidal hexecontahedron',
   'pentakis Dodecahedron': 'Rhombic enneacontahedron'
type SpecialJoin = keyof typeof specialJoins
specialKis :=
   'Rhombic dodecahedron': 'Disdyakis dodecahedron'
   'Rhombic triacontahedron': 'Disdyakis triacontahedron'

type SpecialKis = keyof typeof specialKis
kisWords:= 3: 'triakis', 4: 'tetrakis', 5: 'pentakis', 6: 'hexakis', 7: 'heptakis'
type KisNumber = keyof typeof kisWords

function kisjoin(P: Polyhedron, notation: string,
                 digits: string, join: boolean): Polyhedron
   // kis and join are closely related operations. Both of them add a
   // pyramid on a selection of faces; join then further deletes any
   // _original_ edge bordered by two _new_ triangles, producing a quad.
   // Faces are selected by their numbers of sides, using the given digits.
   // If there are none, all faces are used. Otherwise, the digits are turned
   // into a list of numbers by breaking after every digit except as needed
   // to prevent leading 0s or isolated 1s or 2s (since no face has one or
   // two sides); this way you can list any subset of the numbers 3 - 32,
   // which is plenty.
   // The operation is then applied just to faces with the numbers of edges on
   // the list. e.g. k3412 will add pyramids to the triangles, quads, and
   // dodecagon faces.
   allowed := parseSides digits
   // first collect a directory from face indices to new vertex numbers
   nextVertex .= P.xyz.length
   newVixes :=
      for each f of P.face
         !digits or f.length is in allowed ? nextVertex++ : 0
   if nextVertex is P.xyz.length then return P  // nothing to do
   xyz := P.xyz ++ faceCenters(P).filter (f,ix) => newVixes[ix]
   face: Face[] := []
   for each f, ix of P.face
      v := newVixes[ix]
      if v is 0
         face.push f.slice()
         continue
      // Add the pyramid, possibly eliding edges:
      for each w, jx of f
         pw := f æ jx-1
         neighbor .= 0
         if join
            neighbor = P.face.findIndex (g, gx) =>
                          gx !== ix and w is in g and pw is in g
         if join and newVixes[neighbor]  // elide this edge
            if pw < w // avoid adding same face twice
               face.push [v, pw, newVixes[neighbor], w]
         else face.push [v, pw, w]
   // Nomenclature
   let name: string|undefined
   if P.name and wordsof(P.name) < 3 // don't go hog-wild with naming
      if join
         unless digits
            name = specialJoins[P.name as SpecialJoin] ?? 'joined ' + P.name
      else
         size .= P.face[0].length
         unless P.face.every (f) => f.length is size
            size = 0
         if !size and digits then size = Number(digits) or 0
         // very special case
         if size is 5 and P.name is 'pentagonal antiprism'
            name = 'Icosahedron'
         else if (!digits or Number(digits) is size) and size in kisWords
            name = (specialKis[P.name as SpecialKis]
                    ?? kisWords[size as KisNumber] + ' ' + P.name)
   // Cheaty super special case
   if notation is 'jk5djP5'
      name = 'dual elongated pentagonal orthobirotunda'
   // Done, fix up the vertices a bit
   adjustXYZ {face, xyz, name}, notation, 3


specialGyro :=
   Cube: 'pentagonal icositetrahedron', Octahedron: 'pentagonal icositetrahedron',
   Icosahedron: 'pentagonal hexecontahedron',
   Dodecahedron: 'pentagonal hexecontahedron'
type SpecialGyro = keyof typeof specialGyro

// how enums ought to work?
FromCenter := Symbol()
AlongEdge := Symbol()
type Gyway = typeof FromCenter | typeof AlongEdge

function gyropel(P: Polyhedron, notation: string,
                 digits: string, ...ways: Gyway[]): Polyhedron
   // gyro and propellor are closely related operations. Both of them add new
   // vertices one third of the way along each edge of each face selected
   // by the digits argument (see kisjoin). They then differ in what edges
   // are drawn to the new vertices. In gyro, another new vertex is added
   // at the center of each face and connected to each of them; in propellor,
   // they are just connected in sequence. For completeness, we also allow
   // both at the same time, which is equivalent to propellor followed by kis
   // just on the new rotated faces.
   fromCenter := FromCenter is in ways
   alongEdge := AlongEdge is in ways
   unless fromCenter or alongEdge then return P  // nothing to do
   allowed := parseSides digits
   // first collect a directory from directed edges to new vertex numbers
   xyz := P.xyz.slice()
   startV := xyz.length
   edgeV: number[][] := []
   for each f of P.face
      if digits and f.length is not in allowed then continue
      for each v, ix of f
         pv := f æ ix-1
         (edgeV[pv] ??= [])[v] = xyz.length
         xyz.push lerp xyz[pv], xyz[v], 1/3
   if xyz.length is startV then return P  // nothing to do

   // Now revisit each face, accumulating the new faces.
   face: Face[] := []
   centers: XYZ[] := fromCenter ? faceCenters P : []
   for each f, fx of P.face
      if digits and f.length is not in allowed
         // Just collect all of the vertices around
         newFace: Face := []
         for each v, ix of f
            pv := f æ ix-1
            reverseV := edgeV[v]?[pv] ?? -1
            if reverseV >= 0 then newFace.push reverseV
            newFace.push v
         face.push newFace
         continue
      centerV := xyz.length
      if fromCenter then xyz.push centers[fx]
      aroundOutside: Face .= []
      for each v, ix of f
         pv := f æ ix-1
         ppv := f æ ix-2
         firstNew := edgeV[ppv][pv]
         newSection := [firstNew]
         reverseV := edgeV[pv][ppv] ?? -1
         if reverseV >= 0 then newSection.push reverseV
         newSection.push pv
         secondNew := edgeV[pv][v]
         if alongEdge
            newSection.push secondNew
            face.push newSection
            if fromCenter then face.push ð centerV, firstNew, secondNew
            else aroundOutside.push firstNew
         else if fromCenter
            newSection.push secondNew, centerV
            face.push newSection
         else aroundOutside ++= newSection
      if aroundOutside.length then face.push aroundOutside
   let name: string|undefined
   nw .= 0
   if ((fromCenter xor alongEdge)
       and P.name and !digits and (nw = wordsof(P.name)) < 3)
      if alongEdge
         if nw is 1 then name = 'propello' + P.name
         else name = 'propellorized '+P.name
      else
         if P.name in specialGyro then name = specialGyro[P.name as SpecialGyro]
         else if nw is 1 then name = 'gyro' + P.name
         else name = 'gyrated '+P.name
   // Done, fix up the vertices a bit
   adjustXYZ {face, xyz, name}, notation, 3

function parseSides(digits: string): number[]
   unless digits return []
   tooSmall := ['1', '2']
   last := digits.length - 1
   return := []
   current .= ''
   pos .= 0
   while pos <= last
      current += digits[pos++]
      nextDigit := digits[pos]
      if (current is in tooSmall
          or nextDigit is '0'
          or pos == last and nextDigit is in tooSmall)
         continue
      return.value.push parseInt current
      current = ''

transforms :=
   k: (P: Polyhedron, notation: string, digits: string) =>  // kis[n]
      kisjoin P, notation, digits, false
   j: (P: Polyhedron, notation: string, digits: string) =>  // join
      kisjoin P, notation, digits, true
   g: (P: Polyhedron, notation: string, digits: string) =>  // gyro
      gyropel P, notation, digits, FromCenter
   p: (P: Polyhedron, notation: string, digits: string) =>  // propellor
      gyropel P, notation, digits, AlongEdge
   f: (P: Polyhedron, notation: string, digits: string) =>  // fan [new? name?]
      gyropel P, notation, digits, AlongEdge, FromCenter
   r: (P: Polyhedron) =>                                    // reverse (mirror)
      face: (f.toReversed() for each f of P.face)
      xyz: (scale copy(v), -1 for each v of P.xyz)
      name: P.name  // should we record that it's mirrored?
   d: (P: Polyhedron, notation: string, digits: string) =>  // dual
      if digits
         console.error `Ignoring ${digits} arg of d in ${notation}`
      // Create a "fake" of P and adjust it (so we don't disturb the
      // cached P), which will create the dual
      parentNotation := notation[1+digits.length..]
      adjustXYZ {face: P.face.slice(), P.xyz, P.name}, parentNotation, 1
      polyCache.`d${parentNotation}`
   c: (P: Polyhedron, notation: string, digits: string) =>  // canonicalize
      face: P.face.slice()
      xyz: canonicalXYZ P, notation, Number(digits) or 10
      name: P.name
   x: approxCanonicalize  // iterative direct adjustment algorithm
type TransformOp = keyof typeof transforms

function dispatch(op: string, digits: string,
                  P: Polyhedron, mynotation: string): Polyhedron
   // Note mynotation starts with op!
   return .= P
   if op in seeds
      return = seeds[op as SeedOp] Number(digits) or 0
   else if op in transforms
      return = transforms[op as TransformOp] P, mynotation, digits
   else
      console.error `Unknown operation ${op}${digits} in ${mynotation}.`
      return = polyCache.T
   polyCache[mynotation] = return.value

function topoDual(P: Polyhedron): Polyhedron
   // Note this maintains correspondence between V and F indices, but
   // places every vertex at the origin! Don't use without geometrizing
   // in some way.
   face:
      for v of [0...P.xyz.length]
         infaces :=  // gather labeled list of faces contining v
            for f, index of P.face
               unless f.includes v continue
               ð f, index
         start := infaces[0][1];
         current .= start
         newface := []
         do
           verts := P.face[current]
           preV := verts æ verts.indexOf(v)-1
           nextIx := infaces.findIndex ([face, label]) =>
              label !== current and face.includes preV
           current = infaces[nextIx][1]
           newface.push current
           if newface.length > infaces.length
              console.error 'In topoDual: Malformed polyhedron', P
              break
         until current is start
         newface
   xyz:
      Array(P.face.length).fill([0,0,0]) // warning, every vertex is ===
   name: dualName P.name


function geomDual(P: Polyhedron): Polyhedron
   return := topoDual P
   return.value.xyz = approxDualVertices P

function approxDualVertices(P: Polyhedron): XYZ[]
   P.face.map (f) => approxDualVertex f, P.xyz

operator dot same (/) (v: number[], w: number[])
   v.reduce (l,r,i) => l + r*w[i], 0

function approxDualVertex(f: Face, v: XYZ[]): XYZ
   // For each edge of f, there is a plane containing it perpendicular
   // to the line joining the origin to its nearest approach to the origin.
   // This function returns the point closest to being on all of those planes
   // (in the least-squares sense).
   // This method seems to work OK when the neighborhood of f is convex,
   // and very poorly otherwise. So it probably would not provide any better
   // canonicalization than other methods of approximating the dual.
   // It might be better to replace this with the approximate intersection of the
   // reciprocal planes to each of the vertices of the face...
   normals := (tangentPoint(v[f æ i-1], v[f[i]]) for i of [0...f.length])
   sqlens := normals.map mag2
   columns := (normals.map(&[i]) for i of [0..2])
   target := (columns[i] dot sqlens for i of [0..2]) as XYZ
   CMsource := (for c of [0..2]
        (columns[r] dot columns[c] for r of [0..2])) as [XYZ, XYZ, XYZ]
   cramerD := det ...CMsource
   if Math.abs(cramerD) < 1e-6
      console.error `Face ${f} of ${v.map (p) => '['+p+']'} ill conditioned`
      return [0, 0, 0]
   [ det(target,CMsource[1],CMsource[2])/cramerD,
     det(CMsource[0],target,CMsource[2])/cramerD,
     det(CMsource[0],CMsource[1],target)/cramerD ]

function det(a: XYZ, b: XYZ, c:XYZ)
   a[0]*b[1]*c[2] + a[1]*b[2]*c[0] + a[2]*b[0]*c[1]
      - a[2]*b[1]*c[0] - a[1]*b[0]*c[2] - a[0]*b[2]*c[1]

operator sub(a: XYZ, b: XYZ)
   diminish copy(a), b

function tangentPoint(v: XYZ, w: XYZ) // closest point on vw to origin
   d := w sub v
   v sub scale d, d dot v / mag2 d

function faceCenters(P: Polyhedron): XYZ[]
   for each face of P.face
      scale face.reduce((ctr,v) => accumulate(ctr, P.xyz[v]), [0,0,0]),
         1/face.length

function dualName(p: string|undefined)
   unless 'Octahedron' in specialDuals
      // one-time reversal of all special duals
      specialDuals[dual] = poly for poly, dual in specialDuals

   unless p return undefined
   if p in specialDuals
      return specialDuals[p]
   words := p.split(' ')
   if words[0] is 'dual' return words[1..].join ' '
   if words.length is 2
      switch words[1]
         'prism'
            return words[0] + ' bipyramid'
         'bipyramid'
            return words[0] + ' prism'
         'antiprism'
            return words[0] + ' trapezohedron'
         'trapezohedron'
            return words[0] + ' antiprism'
         'pyramid'
            return p  // self-dual
   return 'dual ' + p

function adjustXYZ(P: Polyhedron, notation: string, iterations = 1): Polyhedron
   dualNotation := 'd' + notation
   D .= topoDual P
   if dualNotation in polyCache
      console.error 'Creating', notation, '_after_ its dual'
      D = polyCache[dualNotation]
   for iter of [1..iterations]
      D.xyz = reciprocalC P
      P.xyz = reciprocalC D
   polyCache[dualNotation] = D
   P

function reciprocalC(P: Polyhedron): XYZ[]
   return := faceCenters P
   for each v of return.value
      scale v, 1/mag2 v

function canonicalXYZ(P: Polyhedron, notation: string, iterations: number): XYZ[]
   dualNotation := 'd' + notation
   D .= topoDual P
   if dualNotation in polyCache
      console.error 'Creating', notation, '_after_ its dual'
      D = polyCache[dualNotation]
   tempP := Object.assign({}, P)  // algorithm is read-only on original data
   if iterations < 1 then iterations = 1
   for iter of [1..iterations]
      D.xyz = reciprocalN tempP
      tempP.xyz = reciprocalN D
   polyCache[dualNotation] = D
   tempP.xyz

operator cross tighter (*) (v: XYZ, w: XYZ): XYZ
   [ v[1]*w[2] - v[2]*w[1], v[2]*w[0] - v[0]*w[2], v[0]*w[1] - v[1]*w[0] ]

function reciprocalN(P: Polyhedron): XYZ[]
   for each f of P.face
      centroid := ð 0, 0, 0
      normal := ð 0, 0, 0
      meanEdgeRadius .= 0
      for each vlabel, ix of f
         v := P.xyz[vlabel]
         accumulate centroid, v
         pv := P.xyz[f æ ix-1]
         ppv := P.xyz[f æ ix-2]
         // original doc says unit normal but below isn't. Didn't help to try, tho
         nextNormal := (pv sub ppv) cross (v sub pv)
         // or instead, just chop down big ones: (didn't work either)
         // magNext := mag nextNormal
         // if magNext > 1 then scale nextNormal, 1/magNext
         accumulate normal, nextNormal
         meanEdgeRadius += mag tangentPoint pv, v
      scale centroid, 1/f.length
      scale normal, 1/mag normal
      meanEdgeRadius /= f.length
      scale normal, centroid dot normal
      scale normal, 1/mag2 normal  // invert in unit sphere
      scale normal, (1 + meanEdgeRadius)/2

operator dist(a: XYZ, b: XYZ) mag a sub b

// adapted from Polyhedronisme https://levskaya.github.io/polyhedronisme/
function approxCanonicalize(P: Polyhedron, notation: string,
                            digits: String): Polyhedron
   THRESHOLD := 1e-6
   // A difficulty is that the planarizing sometimes has the effect of
   // "folding over" neighboring faces, in which case all is lost.
   // Keeping the weight of the edge smoothing high compared to planarizing
   // seems to help with that.
   EDGE_SMOOTH_FACTOR := 0.5
   PLANARIZE_FACTOR := 0.1
   edge := edges P
   xyz := P.xyz.map copy
   V := xyz.length
   normalizeEdges xyz, edge
   for iter of [1..Number(digits) or 10]
      start := xyz.map copy
      smoothEdgeDists xyz, edge, EDGE_SMOOTH_FACTOR
      normalizeEdges xyz, edge
      planarize xyz, P.face, PLANARIZE_FACTOR
      normalizeEdges xyz, edge
      if Math.max(...(xyz[i] dist start[i] for i of [0...V])) < THRESHOLD
         break
   {face: P.face.slice(), xyz, P.name}

type Edge = [number, number]
function edges(P:Polyhedron): Edge[]
   return: Edge[] := []
   for each f of P.face
      for each v, ix of f
         pv := f æ ix-1
         if pv < v then return.value.push ð pv, v

function normalizeEdges(xyz: XYZ[], edge: Edge[]): void
   // Adjusts xyz so that edge tangentpoints have centroid at origin and
   // mean radius 1
   edgeP .= edge.map ([a,b]) => tangentPoint xyz[a], xyz[b]
   edgeCentroid := centroid edgeP
   xyz.forEach (pt) => diminish pt, edgeCentroid
   edgeScale := 1/(mean edge.map ([a,b]) => mag tangentPoint xyz[a], xyz[b])
   xyz.forEach (pt) => scale pt, edgeScale

function centroid(xyz: XYZ[]): XYZ
   scale xyz.reduce(accumulate, ð 0,0,0), 1/xyz.length

function smoothEdgeDists(xyz: XYZ[], edge: Edge[], fudge: number): void
   // Attempts in the most straightforward way possible to reduce the
   // variance of the radii of the edgepoints
   V := xyz.length
   adj := (ð 0,0,0 for i of [1..V])
   edgeDistsStart := edge.map ([a,b]) => mag tangentPoint xyz[a], xyz[b]
   for each [a,b] of edge
      t := tangentPoint xyz[a], xyz[b]
      scale t, (1 - mag t)/2
      accumulate adj[a], t
      accumulate adj[b], t
   for i of [0...V]
      accumulate xyz[i], scale adj[i], fudge
   edgeDistsEnd := edge.map ([a,b]) => mag tangentPoint xyz[a], xyz[b]

function summary(data: number[])
   [Math.min(...data), Math.max(...data), mean(data),
    mean(data.map((x) => Math.abs(1-x)))]

function planarize(xyz: XYZ[], face: Face[], fudge: number): void
   V := xyz.length
   adj := (ð 0,0,0 for i of [1..V])
   for each f of face
      if f.length is 3 then continue  // triangles always planar
      fxyz := (xyz[v] for each v of f)
      c := centroid fxyz
      n := meanNormal fxyz
      if c dot n < 0 then scale n, -1
      for each v of f
         accumulate adj[v], scale copy(n), n dot (c sub xyz[v])
   for i of [0...V]
      accumulate xyz[i], scale adj[i], fudge

function meanNormal(xyz: XYZ[]): XYZ
   mNormal :=  ð 0,0,0
   [v1, v2] .= xyz.slice(-2);
   for each v3 of xyz
      nextNormal := (v2 sub v1) cross (v3 sub v2)
      magNext := mag nextNormal
      // reduce influence of long edges? (didn't seem to help in brief testing)
      // if magNext > 1 then scale nextNormal, 1/Math.sqrt magNext
      accumulate mNormal, nextNormal
      [v1, v2] = [v2, v3] // shift over one
   scale mNormal, 1/mag mNormal

function mean(a: number[])
   m .= 0
   m += e for each e of a
   m / a.length

// arithmetic on 3-vectors
function accumulate(basket: XYZ, egg: XYZ)
   basket[0] += egg[0]
   basket[1] += egg[1]
   basket[2] += egg[2]
   basket

function diminish(basket: XYZ, egg: XYZ)
   basket[0] -= egg[0]
   basket[1] -= egg[1]
   basket[2] -= egg[2]
   basket

function copy(a: XYZ)
   ð a[0], a[1], a[2]

function mag2(a: XYZ)
   a[0]*a[0] + a[1]*a[1] + a[2]*a[2]

function mag(a: XYZ)
   Math.sqrt mag2 a

function normalize(v: XYZ)
   scale v, 1/mag v

function unit

function scale(subject: XYZ, by: number)
   subject[0] *= by
   subject[1] *= by
   subject[2] *= by
   subject

function lerp(start: XYZ, end: XYZ, howFar: number)
   basket := scale copy(start), 1-howFar
   accumulate basket, scale copy(end), howFar

// Feedback

function inform(x: string)
   $('input[name="inform"]').val(x)
   x


// VRML97 generation

function outputVRML(notation: Notation, P: Polyhedron): string
   shortDescrip := P.name or `A ${P.face.length}-hedron`
   ```
#VRML V2.0 utf8
Group { children [
   WorldInfo { # Generated by GTW's reimplementation of GWH's Conway script.
      title "${notation}   ${stats P}"
      info "Generated by GTW's Conway-notation script inspired by GWH's.
By using this script, you agree that this image is released
into the public domain, although you are requested to cite
George Hart's Encyclopedia of Polyhedra as the source." }
   Background {
      groundColor [ .2 .5 1 ] # light blue
      skyColor    [ .2 .5 1 ] }
   NavigationInfo { type [ "EXAMINE" ] }
   DirectionalLight {direction -.5 -1  1 intensity 0.75}
   DirectionalLight {direction  .5  1 -1 intensity 0.75}
   Viewpoint { position 0 0 4.5 description "${shortDescrip}" }
   ${polyVRML P, colorScheme()}
   Shape {
      appearance Appearance {
         material Material {
            diffuseColor 0 0 0 } }
      geometry IndexedLineSet {
         coord ${useVerts}
         coordIndex [
            ${edgeIndices P} ] } }
   ${showDual() ? polyVRML geomDual(P), '0.5 0.5 0.5' : ''} ] }
   ```

function stats(P: Polyhedron): string
   edges := P.face.reduce((e,f) => e + f.length, 0) / 2
   `[${P.face.length} faces, ${edges} edges, ${P.xyz.length} vertices]`

useVerts := 'USE verts'
function defVerts(v: XYZ[])
   `DEF verts Coordinate {
         point [
            ${v.map(.join ' ').join(",\n            ")} ] }`

function polyVRML(P: Polyhedron, color: string): string
   facePartition := color ? [P.face]
      : P.face.reduce ((parts, f) =>
                          (parts[f.length] ??= []).push f
                          parts),
         [] as Face[][]
   emittedCoords .= ''
   shapes :=
      for part of facePartition
         unless part continue
         `
   Shape {
      appearance Appearance {
         material Material {
            diffuseColor ${color or colorBySides part[0].length} } }
      geometry IndexedFaceSet {
         ccw FALSE
         coord ${emittedCoords ? useVerts : (emittedCoords = defVerts P.xyz)}
         coordIndex [
            ${part.map(.join ', ').join(", -1,\n            ")}, -1 ] }}`
   shapes.join "\n"

function colorScheme
   button := document.getElementsByName('color')[0] as HTMLInputElement
   button.checked ? '1 1 1' : ''

function showDual
   false

faceColors: Record<number, string> :=
   3: '0.9 0.3 0.3' // red
   4: '0.4 0.4 1.0' // blue
   5: '0.2 0.9 0.3' // green
   6: '0.9 0.9 0.2' // yellow
   7: '0.5 0.25 0.25' // brown
   8: '0.8 0.2 0.8' // magenta
   9: '0.5 0.2 0.8' // purple
   10: '0.1 0.9 0.9' // cyan
   12: '1.0 0.6 0.1' // orange

function colorBySides(n: number)
   if n in faceColors
      return faceColors[n]
   '0.5 0.5 0.5' // gray

function filtmap<T,U>(A: T[], m: (e:T, i: number, arr: T[]) => U)
   A.map(m).filter (e) => !!e

function edgeIndices(P: Polyhedron)
   sep := ",\n            "
   filtmap(P.face, (thisf) =>
         filtmap(thisf, (v, ix, f) =>
            preV := f æ ix-1
            preV < v ? `${preV}, ${v}, -1` : '').join sep)
      .join sep