// MIT License
//
// Copyright (c) 2021 Joel Gibson
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

// A program for calculating permutation and irreducible characters of the symmetric group.
//
// Author: Joel Gibson
// Last updated: 2021-09-02
//
// TODO:
//    - Decomposing modules when n >= 20 or so is still slow, and a lot of time is spent
//      in the inner product function. It is fast enough to be useful though.


// Assert that a condition is true, otherwise crash the program with a message.
function assert(cond: boolean, msg?: string) {
    if (!cond)
        throw new Error("Assertion failed" + ((msg == undefined) ? "." : ": " + msg))
}

// Memoise a function.
function memoise<K, V>(f: (k: K) => V): (k: K) => V {
    let cache = new Map<K, V>()
    return (k: K) => {
        if (!cache.has(k))
            cache.set(k, f(k))

        return cache.get(k)!
    }
}

// Sums/products of numbers/bigints.
function sum(ns: number[]): number { return ns.reduce((a, b) => a + b, 0) }
function bigsum(ns: bigint[]): bigint {
    let acc = 0n
    for (let i = 0; i < ns.length; i++)
        acc += ns[i]

    return acc
}
function bigproduct(ns: bigint[]): bigint { return ns.reduce((a, b) => a * b, 1n) }

// The factorial of n.
const factorial = memoise(function factorialHelp(n: number): bigint {
    assert(n >= 0)
    return (n == 0) ? 1n : BigInt(n) * factorial(n - 1)
})

// The greatest common divisor of two integers.
function gcd(a: number, b: number): number {
    if (b == 0)
        return a;
    return gcd(b, a % b);
}

// Return base^exp for bigints.
function bigpow(base: bigint, exp: number): bigint {
    assert(exp >= 0)
    let acc = 1n
    for (let i = 0; i < exp; i++)
        acc *= base

    return acc
}


// Mapping partitions to integer tokens and back.
//
// To use the Javascript Map type, we need some way of mapping partitions to and from primitive objects like numbers.
// We will map a partition to a pair (size of partition, index in list of partitions of that size), and further encode
// that pair as an integer using a high-bits low-bits strategy. This compound integer is called a token.
//
// The lower 6 bits of an integer records the size of a partition, in the range [0, 64).
// The higher bits form an index into an enumeration of the partitions of a given size.
// There are 1,505,499 ≈ 2^21 partitions of size 63, so we will use at most 27 bits of each number (though I suspect
// getting character values for symmetric groups larger than S_30 or so will be challenging).
//
// The partitions of n are indexed by sorting the first part in increasing order, then the second, and so on (lex order).
// For example the partitions of 6 are:
//  x    xx    xx    xx    xxx    xxx    xxx    xxxx    xxxx    xxxxx    xxxxxx
//  x    x     xx    xx    x      xx     xxx    x       xx      x
//  x    x     x     xx    x      x             x
//  x    x     x           x
//  x    x
//  x
//
//      Index:      0,     1,    2,   3,    4,   5,  6,   7,  8,  9, 10
//  Partition: 111111, 21111, 2211, 222, 3111, 321, 33, 411, 42, 51, 6
//
// For example, the partition 3111 would be encoded as 4 * N + 6, where N = 64.
//
// In order to convert to and from integers, we will keep two NxN arrays, where N=64 is the maximum size bound.
//   P(n, k) = # of partitions of n with first part equal to k.
//   R(n, k) = # of partitions of n of width <= k.
//
// These are mutually recursive: for P we have
//   P(0, 0) = 1,
//   P(n, 0) = 0           if n >= 1,
//   P(n, k) = 0           if k > n, and
//   P(n, k) = R(n - k, k) if k <= n.
// The recurrence in the last formula comes from cutting off the first row of length k, and being left with a
// partition of (n - k) which has width at most k. For R we have
//   R(n, k) = P(n, 0) + P(n, 1) + ... + P(n, k)  for 0 <= k
// so R is basically just a cumulative sum speeding up the computation.
// The number of partitions of size n is R(n, n).
// Note that R(n, n) = R(n, n+1) = ...
//
// Take the example above again, we have
//        k: 0  1  2  3  4  5  6
//  P(6, k): 0  1  3  3  2  1  1
//  R(6, k): 0  1  4  7  9  10 11
//
// An example: converting 2211 to a index. Looking at the first part 2 and the value R(6, 2 - 1) = 1, we can tell that
// the partitions with first part 2 start at index 1. Now, there is an order-preserving bijection
// {partitions of 6 with first part 2} -> {partitions of 4 of width at most 2}, by chopping off the first row, so we
// recurse the same process for 211, and so on.
namespace tok {
    // N = 64 = 2^6, so 6 is magic bitshift number, 0x3f is magic bitmask number.
    const N = 64

    // We only use the R array (the only trick is filling it correctly.)
    export const RData = new Int32Array(N * N)

    // Base case: P(0, 0) = R(0, 0) = 1.
    RData[0*N + 0] = 1

    // Recursive case: we already have R(n, 0) = 0, so for the other points we use
    // R(n, k) = R(n, k - 1) + P(n, k)
    for (let n = 0; n < N; n++)
        for (let k = 1; k < N; k++)
            RData[N*n + k] = RData[N*n + k - 1] + ((k <= n) ? RData[N*(n - k) + k] : 0)


    export function R(n: number, k: number): number { return (n >= 0 && k >= 0) ? RData[N*n + k] : 0 }

    /** Convert a partition (which may have trailing zeros) to a token. */
    export function parToTok(parts: readonly number[]): number {
        let size = 0
        for (let i = 0; i < parts.length; i++)
            size += parts[i]

        if (size >= N)
            throw Error("Partition was too large")

        let index = 0
        for (let i = 0, sz = size; i < parts.length && parts[i] > 0; i++) {
            index += R(sz, parts[i] - 1)
            sz -= parts[i]
        }

        return (index << 6) + size
    }

    /** Convert a token to a partition. */
    export function tokToPar(index: number): number[] {
        let parts: number[] = []
        let size = index & 0x3f
        let idx = index >>> 6

        // Here we should really be doing a binary search...
        while (size > 0) {
            // The indices for partitions with first part k are [R(n, k-1), R(n, k)).
            let k = 1
            while (k <= size && !(idx < R(size, k)))
                k++

            if (k > size)
                throw new Error("Out of bounds")

            parts.push(k)
            idx -= R(size, k - 1)
            size -= k
        }

        return parts
    }

    /** The number of partitions of a size. */
    export function countParOfSize(size: number): number { return R(size, size) }

    /** A list of the tokens of all partitions of a given size. */
    export function toksOfSize(size: number): number[] {
        let result: number[] = []
        for (let i = 0; i < RData[N*size + size]; i++)
            result[i] = (i << 6) + size

        return result
    }
}



// ------------------
// --- Partitions ---
// ------------------

// A partition is a weakly decreasing list of positive integers.
type Partition = number[];

// Check that a list is weakly decreasing, i.e. [a >= b >= ... >= d].
function isWeaklyDecreasing(ns: number[]): boolean {
    for (let i = 0; i < ns.length - 1; i++)
        if (ns[i] < ns[i + 1])
            return false;

    return true;
}

// Check a list is a partition.
function isPartition(ns: number[]): boolean {
    return isWeaklyDecreasing(ns) && ns.every(x => x > 0);
}

// Turn any weak composition (list of nonnegative integers) into a partition,
// by removing any zero parts, and sorting.
function partitionFromComposition(comp: number[]): Partition {
    return comp.filter(x => x > 0).sort((a, b) => b - a);
}

// Generate all partitions of n, in lex order.
const partitionsOf = memoise(function partitionsOfHelp(n: number): Partition[] {
    // Helper function: returns all the partitions of n which fit inside a width w strip.
    function helper(n: number, w: number): Partition[] {
        assert(n >= 0);
        assert(w >= 0);

        if (n == 0)
            return [[]];
        if (w <= 0)
            return [];

        let results: Partition[] = [];
        for (let row = Math.min(n, w); row >= 1; row--) {
            let smaller = helper(n - row, row);
            results.push(...smaller.map(part => [row].concat(part)));
        }

        return results;
    }
    return helper(n, n);
})

// Turn a partition into "group notation", for example the partition
// [4, 4, 3, 2, 2, 2, 1] => [[4, 2], [3, 1], [2, 3], [1, 1]].
function groupPartition(partition: Partition): [number, number][] {
    if (partition.length == 0)
        return [];

    let i = 0, j = 0;
    let groups: [number, number][] = [];
    while (i < partition.length) {
        for (j = i; j < partition.length; j++)
            if (partition[i] != partition[j])
                break;

        groups.push([partition[i], j - i]);
        i = j;
    }

    return groups;
}

// The symmetric group S_l acts on a partition of length l by permuting the parts.
// When the partition is written in "group notation" [[a1, b1], ..., [ar, br]], the
// order of the stabiliser of this action is b1! * ... * br!.
function stabiliserOrder(partition: Partition): bigint {
    let factorials = groupPartition(partition).map(([_, b]) => factorial(b));
    return bigproduct(factorials);
}

// For a partition of n, how many elements of the symmetric group S_n are conjugate
// to an element with the given cycle type? When the cycle type is written in
// "group notation" [[a1, b1], ..., [ar, br]], the answer is
// n! / ((a1^b1 * ... * ar^br) * (b1! * ... * br!)
const conjugacySize = memoise(function conjugacySizeHelp(cycleType: Partition): bigint {
    let n = sum(cycleType);
    let groups = groupPartition(cycleType);
    let factorials = groups.map(([_, b]) => factorial(b));
    let powers = groups.map(([a, b]) => bigpow(BigInt(a), b));

    return factorial(n) / (bigproduct(powers) * bigproduct(factorials));
})
const conjugacySizeTok = memoise(function conjugacySizeTokHelp(token: number) { return conjugacySize(tok.tokToPar(token))})

// For an element of the symmetric group with a given cycle type, what cycle type
// does the pth power of that element have? The answer can be computed independently
// for each cycle; for a cycle of size n, the power n^p breaks into gcd(n, p) cycles,
// each of length n / gcd(n, p).
function cyclePow(cycleType: Partition, p: number): Partition {
    let newCycleType: number[] = [];
    for (let cycle of cycleType) {
        let div = gcd(cycle, p);
        for (let i = 0; i < div; i++)
            newCycleType.push(cycle / div);
    }
    return partitionFromComposition(newCycleType);
}

// Turn a partition into a string
function showPartition(partition: Partition): string {
    return JSON.stringify(partition);
}

// Turn a string into a partition
function readPartition(partition: string): Partition {
    return JSON.parse(partition);
}

// Return a LaTeX expression of the partition in "group notation", so
// [4, 4, 3, 2, 2, 2, 1] => [4^2, 3^1, 2^3, 1^1].
function showPartitionCompact(partition: Partition): string {
    return "[" + groupPartition(partition).map(([a, b]) => `${a}^{${b}}`).join(",") + "]";
}

function showPartitionCompactHTML(partition: Partition): string {
    return "[" + groupPartition(partition).map(([a, b]) => `${a}<sup>${b}</sup>`).join("") + "]";
}

// A map of type (Partition -> T). Since Javascript only understands strings as
// may keys, we convert to and from strings a lot.
class PartMap<T> {
    protected store: Map<number, T> = new Map()

    has(k: Partition): boolean {
        return this.store.has(tok.parToTok(k))
    }

    get(k: Partition): T | null {
        let idx = tok.parToTok(k)
        if (!this.store.has(idx))
            return null

        return this.store.get(idx)!
    }

    set(k: Partition, v: T): void {
        this.store.set(tok.parToTok(k), v)
    }

    keys(): Partition[] {
        return [...this.store.keys()].map(tok.tokToPar)
    }
}

// -----------------------------------------
// --- Linear combinations of partitions ---
// -----------------------------------------

// An integer-linear combination of partitions.
class Lin extends PartMap<bigint> {
    get(k: Partition): bigint {
        let idx = tok.parToTok(k)
        if (!this.store.has(idx))
            return 0n

        return this.store.get(idx)!
    }

    getTok(k: number): bigint {
        let result = this.store.get(k)
        return (result === undefined) ? 0n : result
    }

    addTerm(k: readonly number[], v: bigint): void {
        let idx = tok.parToTok(k)
        let val = this.store.get(idx)
        this.store.set(idx, (val === undefined) ? v : val + v)
    }

    support(): Partition[] {
        return this.keys().filter(k => this.get(k) != 0n);
    }

    toList(): [Partition, bigint][] {
        return this.support().map(k => [k, this.get(k)]);
    }

    show(letter?: string): string {
        const output: string[] = [];
        for (const part of this.support()) {
            let coeffStr = "" + this.get(part);
            let partStr = showPartition(part);

            if (partStr == "[]") {
                output.push(coeffStr);
                continue;
            }
            if (coeffStr == "1")
                coeffStr = "";

            output.push(coeffStr + ((letter === undefined) ? '' : letter) + partStr)
        }
        if (output.length == 0)
            return "0";
        return output.join(" + ");
    }

    static fromList(pairs: [Partition, bigint][]): Lin {
        let lin = new Lin();
        for (let [partition, coeff] of pairs)
            lin.set(partition, lin.get(partition) + coeff);

        return lin;
    }

    static zero(): Lin {
        return Lin.fromList([]);
    }

    static scale(scalar: bigint, lin: Lin): Lin {
        return Lin.fromList(lin.toList().map(([k, v]) => [k, scalar * v]));
    }

    /** Ensure each term divides evenly into divisor, and return result. */
    static divide(divisor: bigint, lin: Lin): Lin {
        lin.toList().forEach(([k, v]) => {
            assert(v % divisor == 0n)
        })
        return Lin.fromList(lin.toList().map(([k, v]) => [k, v / divisor]));
    }

    static add(left: Lin, right: Lin): Lin {
        return Lin.fromList(left.toList().concat(right.toList()));
    }

    static sub(left: Lin, right: Lin): Lin {
        return Lin.add(left, Lin.scale(-1n, right));
    }

    static applyDiagonal(f: (a: Partition) => bigint, lin: Lin): Lin {
        return Lin.fromList(lin.toList().map(([k, v]) => [k, f(k) * v]));
    }

    static applyBilinear(f: (a: Partition, b: Partition) => Lin, left: Lin, right: Lin): Lin {
        let pairs: [Partition, bigint][] = [];
        for (let [k1, v1] of left.toList())
            for (let [k2, v2] of right.toList())
                pairs.push(...Lin.scale(v1 * v2, f(k1, k2)).toList());

        return Lin.fromList(pairs);
    }
}


// ------------------------------------
// ---   Power sum to Schur basis   ---
// ------------------------------------

/** Iterate over the partitions outer for which outer/shape is a border strip of size length.
 * The array 'outer' which is returned is a partition, but it has trailing zeros. It will also be
 * overwritten between calls.
*/
function iterBorderStrip(shape: readonly number[], length: number, fn: (outer: readonly number[], height: number) => void): void {
    if (length <= 0)
        throw new Error(`Length must be positive for iterBorderStrip, got length=${length}`)

    // Inner will store the inner shape, but pad the array with zeros at the end so we don't need
    // to worry about out-of-bounds errors. Outer will store the inner shape plus the border strip.
    let inner: number[] = shape.slice()
    for (let i = 0; i < length; i++)
        inner.push(0)

    let outer = inner.slice()

    // A border strip has a "top" and "bottom" cell, when the partition is written in English notation.
    // So long as B and T are placed on the end of rows, and T is north-east of B, there is a valid border
    // strip from B to T, including both of them. (B may be placed in the "zeros" at the end of the partition,
    // and T may be placed in the empty space to the right of the first row, as well).
    //
    //   012345678901
    // 0 xxxxxxxxx  T             xxxxxxxxxooT    Border strip of length 14
    // 1 xxxxxx                   xxxxxxoooo
    // 2 xxxxx                    xxxxxoo
    // 3 xxxxx                    xxxxxo
    // 4 xxxxx                    xxxxxo
    // 5 xxxB                     xxxBoo
    // 6 x                        x
    // 7 x                        x
    //
    // Assuming 0-indexed rows and columns, B has coordinates (5, 3) and T has coordinates (0, 11),
    // giving a signed L1 distance of (Row(B) - Row(T)) + (Col(T) - Col(B)) + 1 = (5 - 0) + (11 - 3) + 1 = 14.
    // We use the signed L1 distance since it will be able to detect if B and T ever invert positions.
    //
    // Row(B) and Col(T) are used as the "anchor points" which determine the positions of B and T: the values Col(B)
    // is determined by Row(B), and likewise Row(T) is determined by Col(T).
    //
    // In this example, we'll find border strips of length 4 next to the partition (4, 2).
    // The strategy is to start with B and T in the top row, creating a border strip of the correct length.
    // The resulting shape is "emitted", then B is moved to the next position. This will create a signed distance
    // which is too long, so T is moved to the next position. If the signed distance is right, another shape is
    // emitted, otherwise if the signed distance is still too long, T is moved again. If too short, B is moved.
    // And so on.
    //
    // xxxxBooT   xxxxoooT    xxxxooT    xxxxoT     xxxxT    and so on.
    // xx         xxBoo       xxBoo      xxBoo      xxBoo
    // emit       too long    too long   too long   emit
    // move B     move T      move T     move T     move B

    // Initial: top row. (inner always has positive length, so inner[0] is always defined).
    let B_row = 0
    let B_col = inner[0]
    let T_col = inner[0] + length - 1
    let T_row = 0

    while (true) {
        let L1dist = (B_row - T_row) + (T_col - B_col) + 1

        // Emit if we have the correct length.
        if (L1dist == length) {
            // "Set" the strip, i.e. fill in outer with the strip between B and T. The top row goes all the way to T.
            outer[T_row] = T_col + 1
            // The other rows get extended out to 1 past their immediate upper neighbour.
            for (let r = T_row + 1; r <= B_row; r++)
                outer[r] += inner[r - 1] - inner[r] + 1

            // Emit the shape and the height. According to Stanley, the height is 1 less than the number of rows the strip touches.
            let height = B_row - T_row
            fn(outer, height)
            //console.log(`B = (${B_row}, ${B_col}), T = (${T_row}, ${T_col}), inner = (${inner}), outer = (${outer})`)

            // Set outer back to normal.
            for (let r = T_row; r <= B_row; r++)
                outer[r] = inner[r]
        }

        // We may have just emitted (==), or perhaps the strip is too short. In this case move B.
        if (L1dist <= length) {
            // If we're trying to increase B and it's at the end, we must have found all border strips.
            if (B_row + 1 == inner.length)
                return

            // Every row is a legal position for B, so we just need to increase its row, and calculate the
            // column position (place B in the row far to the right of the partition, and let it "fall" left).
            B_row += 1
            B_col = inner[B_row]
        }

        // Otherwise, the strip is too long, so we need to move T.
        if (L1dist > length) {
            // If we're moving T and it's in the 0th column, there's nowhere to go.
            if (T_col == 0)
                return

            // Otherwise, every column is a legal position for T, so we need to decrease its column, and
            // calculate the row position.
            T_col -= 1
            while (inner[T_row] > T_col) T_row += 1
        }
    }
}

/** Return a list of partitions with the border strip added. */
function borderStrips(shape: readonly number[], length: number): number[][] {
    let result: number[][] = []
    iterBorderStrip(shape, length, parts => result.push(parts.slice()))
    return result
}

function powersumToSchur(partition: Partition): Lin {
    let acc = Lin.fromList([[[], 1n]])

    // Multiplying the partition in reverse order (smaller parts first) seems a bit faster.
    for (let i = partition.length - 1; i >= 0; i--) {
        let prod = new Lin()
        acc.toList().map(([parts, coeff]) => {
            iterBorderStrip(parts, partition[i], (shape, height) => {
                prod.addTerm(shape, ((height % 2 == 0) ? 1n : -1n) * coeff)
            })
        })

        acc = prod
    }

    return acc
}

// To generate the character table for the specht modules, we just need to transpose the
// turn-power-sums-into-schur data.
const characterTablePowersum = memoise(function characterTablePowersumHelp(n: number): CharacterTable {
    let start = performance.now()
    let table = new CharacterTable(n);

    for (let col of table.partitions) {
        let colentries = powersumToSchur(col)
        for (let row of table.partitions)
            table.get(row).set(col, colentries.get(row))
    }

    console.log(`Time to generate irreducible character table for S${n}: ${performance.now() - start}ms`)
    return table;
})

// ------------------------------------
// --- Power sum to monomial basis  ---
// ------------------------------------

// Multiply the monomial symmetric function m_(r) with an augmented monomial symmetric function.
// The product M_r * M(l1, ..., lk) is equal to the sum
//    M(l1 + r, ..., lk) + ... + M(l1, ..., lk + r) + M(r, l1, ..., lk),
// where we need to remember to sort the compositions on the right side back into being partitions.
function multRowWithAugmented(rowPartition: Partition, partition: Partition): Lin {
    if (partition.length == 1)
        [rowPartition, partition] = [partition, rowPartition];
    assert(rowPartition.length == 1);
    let row = rowPartition[0];
    let compositions = partition.map(
        (part, i) => partition.slice(0, i).concat([part + row]).concat(partition.slice(i + 1)));
    compositions.push(partition.concat([row]));
    return Lin.fromList(compositions.map(comp => [partitionFromComposition(comp), 1n]));
}

// Expand a power sum symmetric function into an augmented monomial symmetric function.
function powerToAugmentedMonomial(partition: Partition): Lin {
    return partition
        .map(part => Lin.fromList([[[part], 1n]]))
        .reduce((a, b) => Lin.applyBilinear(multRowWithAugmented, a, b), Lin.fromList([[[], 1n]]));
}

// Expand a power sum symmetric function into monomial symmetric function.
function powerToMonomial(partition: Partition): Lin {
    return Lin.applyDiagonal(stabiliserOrder, powerToAugmentedMonomial(partition));
}

// ----------------------------------------
// --- Characters and character tables  ---
// ----------------------------------------

// A character of the symmetric group is a linear combination of partitions. A character table
// is a map from partitions to such linear combinations.
class CharacterTable extends PartMap<Lin> {
    readonly partitions: Partition[]

    constructor(readonly n: number) {
        super();
        this.partitions = partitionsOf(n);
        for (let partition of this.partitions)
            this.set(partition, Lin.zero());
    }

    get(partition: Partition): Lin {
        let result = super.get(partition);
        assert(result != null);
        return <Lin>result;
    }
}

// To generate the character table for the permutation modules, we just need to transpose the
// turn-power-sums-into-monomial data.
function characterTablePermutation(n: number): CharacterTable {
    let start = performance.now()
    let table = new CharacterTable(n);

    for (let col of table.partitions)
        for (let [row, coeff] of powerToMonomial(col).toList())
            table.set(row, Lin.add(table.get(row), Lin.fromList([[col, coeff]])))

    console.log(`Time to generate permutation character table for S${n}: ${performance.now() - start}ms`)
    return table;
}

// Inner product on class functions in the symmetric group on n letters.
function innerProduct(n: number, left: Lin, right: Lin): bigint {
    // Iterating over the maps inside of left and right using the iterator protocol is bad: so much garbage that GC freaks out.
    return bigsum(tok.toksOfSize(n).map(token => left.getTok(token) * right.getTok(token) * conjugacySizeTok(token))) / factorial(n)
}

// Character table for the Specht modules (we will memoise this).
const characterTableSpecht = memoise(function characterTableSpechtHelp(n: number): CharacterTable {
    let permTable = characterTablePermutation(n);
    let partitions = partitionsOf(n);
    for (let i = 0; i < partitions.length; i++) {
        // This row of permTable is currently a Specht module. We look at everything below it, compute
        // the inner product, and subtract off where needed.
        let spechtMod = permTable.get(partitions[i]);
        for (let j = i + 1; j < partitions.length; j++) {
            let mod = permTable.get(partitions[j]);
            let multiplicity = innerProduct(n, spechtMod, mod);
            permTable.set(partitions[j], Lin.sub(mod, Lin.scale(multiplicity, spechtMod)));
        }
    }

    return permTable;
})


// Tensor product on class functions in the symmetric group on n letters.
function tensorProduct(n: number, left: Lin, right: Lin): Lin {
    let character = Lin.zero();
    for (let part of partitionsOf(n))
        character.set(part, left.get(part) * right.get(part));
    return character;
}

// Decompose a character into a linear combination of irreducible characters of S_n.
function decomposeCharacter(n: number, character: Lin): Lin {
    let irrCharacters = characterTablePowersum(n);
    return Lin.fromList(partitionsOf(n).map(partition => [partition, innerProduct(n, irrCharacters.get(partition), character)]))
}

// Return the trivial character for S_n.
function trivialCharacter(n: number): Lin {
    return Lin.fromList(partitionsOf(n).map(part => [part, 1n]));
}

// Given a character 𝜒, return the character (g -> 𝜒(g^p)).
function characterPow(n: number, character: Lin, p: number): Lin {
    return Lin.fromList(partitionsOf(n).map(part => [part, character.get(cyclePow(part, p))]));
}

// Compute the (r + 1) exterior powers 𝛬^0(V), 𝛬^1(V), ..., 𝛬^r(V)
function exteriorPowers(n: number, r: number, character: Lin): Lin[] {
    let exteriors = [trivialCharacter(n)];

    // Now apply the recurrence e_i = (1/i)(e_{i-1} p_1 - e_{i-2} p_2 + ... + (-1)^(i-1) e_0 p_i)
    for (let i = 1; i <= r; i++) {
        let altSum = Lin.zero();
        for (let j = 1; j <= i; j++) {
            let term = tensorProduct(n, exteriors[i - j], characterPow(n, character, j));
            altSum = Lin.add(altSum, Lin.scale((j % 2 == 0) ? -1n : 1n, term));
        }
        exteriors[i] = Lin.divide(BigInt(i), altSum);
    }

    return exteriors;
}

// Compute the (r + 1) symmetric powers S^0(V), S^1(V), ..., S^r(V)
function symmetricPowers(n: number, r: number, character: Lin): Lin[] {
    let symmetrics = [trivialCharacter(n)];

    // Now apply the recurrence e_i = (1/i)(e_{i-1} p_1 - e_{i-2} p_2 + ... + (-1)^(i-1) e_0 p_i)
    for (let i = 1; i <= r; i++) {
        let sum = Lin.zero();
        for (let j = 1; j <= i; j++) {
            let term = tensorProduct(n, symmetrics[i - j], characterPow(n, character, j));
            sum = Lin.add(sum, term);
        }
        symmetrics[i] = Lin.divide(BigInt(i), sum);
    }

    return symmetrics;
}

function tensorPowers(n: number, r: number, character: Lin): Lin[] {
    let tensors = [trivialCharacter(n)];
    for (let i = 1; i <= r; i++)
        tensors[i] = tensorProduct(n, tensors[i-1], character);
    return tensors;
}



// Tests
function checkOrthonormal(n: number) {
    let table = characterTablePowersum(n)
    let partitions = partitionsOf(n)
    for (let i = 0; i < partitions.length; i++) {
        let part1 = partitions[i]
        let prod = innerProduct(n, table.get(part1), table.get(part1))
        if (prod != 1n)
            throw new Error(`Inner product (s[${part1}], s[${part1}]) = ${prod}, should have been 1.`)

        for (let j = 0; j < i; j++) {
            let part2 = partitions[j]
            let prod = innerProduct(n, table.get(part1), table.get(part2))
            if (prod != 0n)
                throw new Error(`Inner product (s[${part1}], s[${part2}]) = ${prod}, should have been 0.`)
        }
    }
}

// Run self-tests before continuing (takes about 20ms).
for (let n = 0; n < 10; n++)
    checkOrthonormal(n)





// The rest of this file is just interfacing to HTML.
// The programmer should ensure that katex is loaded externally (eg through a <script> tag).
declare var katex: any



/** Render the a character table densely. Each tbody->tr has a data-partition attribute recording
 * the JSON encoding of the partition. We do this using string-concatenation rather than making nodes,
 * because it is a ton of data for large groups. */
function denseCharacterTableHTML(letter: string, table: CharacterTable) {
    let parts = table.partitions
    let partsReverse = table.partitions.slice().reverse()

    let tableHTML = '<table>'
    tableHTML += '<thead><tr><td></td>'
    for (let part of partsReverse)
        tableHTML += "<th>" + showPartitionCompactHTML(part) + "</th>"
    tableHTML += '</tr><tr><td>#</td>'
    for (let part of partsReverse)
        tableHTML += "<td>" + conjugacySize(part) + "</td>"
    tableHTML += '</tr>'
    tableHTML += '</thead>'
    tableHTML += '<tbody>'
    for (let part of parts) {
        tableHTML += `<tr class="character" data-partition=${showPartition(part)}><td>` + letter + showPartitionCompactHTML(part) + '</td>'
        for (let colpart of partsReverse)
            tableHTML += '<td>' + table.get(part).get(colpart) + '</td>'
        tableHTML += '</tr>'
    }
    '</table>'

    return tableHTML
}

function sparseCharacterTableHTML(letter: string, table: CharacterTable) {
    let parts = table.partitions

    let tableHTML = '<table>'
    tableHTML += '<tbody>'
    for (let firstPart = table.n; firstPart >= 1; firstPart--) {
        tableHTML += '<tr><td></td>'
        for (let part of parts) {
            if (part[0] != firstPart)
                continue

            tableHTML += `<td class="character" data-partition=${showPartition(part)}>` + letter + showPartitionCompactHTML(part) + '</td>'
        }
        tableHTML += '</tr>'
    }
    '</table>'

    return tableHTML
}

/** Convert a piece of HTML, which should parse as a single DOM node, to that DOM node object. */
function toDOMNode(html: string): HTMLElement {
    return <HTMLElement>document.createRange().createContextualFragment(html).firstChild!
}


function showDomTable(n: number, table: CharacterTable, letter: string): HTMLElement {
    let E = function(el: string, ...children: (string | HTMLElement)[]): HTMLElement {
        let $el = document.createElement(el);
        for (let $child of children) {
            if (typeof $child == 'string') {
                let $text = document.createTextNode($child);
                $el.appendChild($text);
            } else{
                $el.appendChild($child);
            }
        }
        return $el;
    }
    let M = function(text: string): HTMLSpanElement {
        let $span = E('span');
        $span.innerHTML = katex.renderToString(text);
        return $span;
    }


    let selectedPartitions: string[] = [];

    let $product = E('div');

    function productContents(): HTMLElement[] {
        if (selectedPartitions.length == 0)
            return [E('p', "Select some characters to tensor.")];

        let selectedCharacters = selectedPartitions.map(partitionString => table.get(readPartition(partitionString)));
        let chi = selectedCharacters.reduce((a, b) => tensorProduct(n, a, b));

        let interestingCharacters: [string, Lin][] = [
            ["\\chi", chi],
            ...exteriorPowers(n, 4, chi).slice(2).map((lin, r) => <[string, Lin]>[`\\wedge^{${r+2}} \\chi`, lin]),
            ...symmetricPowers(n, 4, chi).slice(2).map((lin, r) => <[string, Lin]>[`S^{${r+2}} \\chi`, lin]),
            ...tensorPowers(n, 4, chi).slice(2).map((lin, r) => <[string, Lin]>[`\\otimes^{${r+2}} \\chi`, lin])
        ];

        return [
            E('p',
                "Selected ",
                M(`\\chi = ${selectedPartitions.map(s => letter + s).join(` \\otimes `)}`)),
            E('h3',
                "Decomposition of ",
                M(`\\chi`),
                " and its exterior, symmetric and tensor powers"),
            createDecompositionTable(interestingCharacters)
        ]
    }

    function updateProduct() {
        $product.innerHTML = ''
        $product.append(...productContents())
    }

    function createDecompositionTable(characters: [string, Lin][]): HTMLElement {
        let irreducibles = characterTablePowersum(n);
        return E('table',
            // Header: the names and dimensions of each character.
            E('thead',
                E('tr',
                    E('td'),
                    ...characters.map(([name, _]) => E('th', M(name)))),
                E('tr',
                    E('td', M('\\dim')),
                    ...characters.map(([_, lin]) => E('td', "" + lin.get(irreducibles.partitions[irreducibles.partitions.length - 1]))))),

            // Body: multiplicities.
            E('tbody',
                ...irreducibles.keys().map(partition => E('tr',
                    E('td',
                        M(`s${showPartitionCompact(partition)}`)),
                    ...characters.map(
                        ([_, lin]) => E('td', "" + innerProduct(n, irreducibles.get(partition), lin)))
                )))
        )
    }

    let $table = toDOMNode((n <= 20) ? denseCharacterTableHTML(letter, table) : sparseCharacterTableHTML(letter, table))
    $table.addEventListener('click', function(event) {
        // Find an element with a data-partition attribute.
        for (let el = <HTMLElement>event.target; el != this; el = el.parentElement!) {
            if (el.dataset.partition !== undefined) {
                let idx = selectedPartitions.indexOf(el.dataset.partition);
                if (idx >= 0) {
                    selectedPartitions.splice(idx, 1);
                    el.classList.remove('selected');
                } else {
                    selectedPartitions.push(el.dataset.partition);
                    el.classList.add('selected');
                }
                updateProduct()
                return
            }
        }
    })

    updateProduct();
    return E('div',
        E('h3',
            `Character table for the ${(letter == 's') ? 'Specht' : 'permutation'} modules of `,
            M(`S_{${n}}`)),
        $table,
        $product);
}

function doComputation() {
    let checked = <HTMLInputElement>document.querySelector('input[name="module-kind"]:checked');
    let n = (<HTMLInputElement>document.getElementById("order")!).valueAsNumber;
    let $tablePlace = <HTMLDivElement>document.getElementById("tablePlace");

    while ($tablePlace.firstChild)
        $tablePlace.removeChild($tablePlace.firstChild);

    let characterTable = (checked.value == "specht") ? characterTablePowersum(n) : characterTablePermutation(n);
    let letter = (checked.value == "specht") ? "s" : "h";

    $tablePlace.appendChild(showDomTable(n, characterTable, letter));
}

let $computationForm = <HTMLFormElement>document.getElementById('computationForm');
$computationForm.addEventListener('submit', ev => ev.preventDefault());
$computationForm.addEventListener('change', ev => {ev.preventDefault(); doComputation();});
doComputation();