pdf-odin-llm/function.odin

package pdf

import "core:math"

// PDF functions (ISO 32000-2, 7.10). Four types share one struct because
// callers only ever evaluate them, and a union would put a switch at every
// call site for no gain.

Function_Kind :: enum u8 {
    Sampled     = 0,   // type 0
    Exponential = 2,   // type 2
    Stitching   = 3,   // type 3
    PostScript  = 4,   // type 4
}

// The most inputs and outputs the library will accept. The format sets no
// limit, but a shading with 33 outputs is a malformed file, not a design.
MAX_FUNCTION_INPUTS  :: 32
MAX_FUNCTION_OUTPUTS :: 32

Function :: struct {
    kind:   Function_Kind,
    n_in:   int,
    n_out:  int,
    // Always 2*n_in long. Inputs are clamped to it before evaluation.
    domain: []f64,
    // 2*n_out long, or empty for a type 2 or 3 function that omits it.
    // Outputs are clamped to it when present.
    range:  []f64,

    // Type 0: a sampled table.
    size:          []int,     // n_in entries
    bits_per_sample: int,
    encode:        []f64,     // 2*n_in
    decode:        []f64,     // 2*n_out
    samples:       []u8,      // raw, packed at bits_per_sample
    order:         int,       // 1 linear, 3 cubic (cubic falls back to linear)

    // Type 2: C0 + x^N * (C1 - C0).
    c0: []f64,
    c1: []f64,
    exponent: f64,

    // Type 3: k sub-functions stitched across the domain.
    functions: []^Function,
    bounds:    []f64,

    // Type 4: a compiled PostScript calculator program.
    program: []Ps_Op,
}

// Ps_Op is one instruction of a type 4 function, already parsed. Branches carry
// absolute jump targets so evaluation is a flat loop rather than a tree walk,
// which keeps the interpreter's step budget easy to enforce.
Ps_Op :: struct {
    kind:  Ps_Op_Kind,
    value: f64,
    // For Jump and Jump_If_False: the index to continue at.
    target: int,
}

Ps_Op_Kind :: enum u8 {
    Push,
    // Arithmetic
    Add, Sub, Mul, Div, Idiv, Mod, Neg, Abs, Ceiling, Floor, Round, Truncate,
    Sqrt, Sin, Cos, Atan, Exp, Ln, Log, Cvi, Cvr,
    // Relational and boolean
    Eq, Ne, Gt, Ge, Lt, Le, And, Or, Xor, Not, Bitshift, True, False,
    // Stack
    Pop, Exch, Dup, Copy, Index, Roll,
    // Control
    Jump, Jump_If_False,
}

// ---------------------------------------------------------------------------
// Limits
//
// Everything a file controls -- how many sub-functions stitch together, how
// long a calculator program is, how many samples a table claims -- is bounded
// here, so that a hostile or merely broken function costs a fixed amount of
// time and memory instead of an unbounded one.
// ---------------------------------------------------------------------------

// How deeply type 3 functions may nest. This is also what stops a stitching
// function that lists itself among its own /Functions: the nested load stops
// here, and because a function is only entered into the cache once it is
// completely built, the graph that results is always acyclic.
MAX_FUNCTION_DEPTH :: 8
// Entries read from any one numeric array (/Domain, /Bounds, /Encode ...).
MAX_FUNCTION_ARRAY :: 4096
// Sub-functions of one type 3 function.
MAX_SUB_FUNCTIONS  :: 1024
// Samples one type 0 table may hold, across every input and output. 64 Mi
// samples is far beyond any real gradient and keeps size * n_out in range.
MAX_FUNCTION_SAMPLES :: 1 << 26
// Inputs a type 0 function is interpolated over. Multilinear interpolation
// visits 2^n_in corners per evaluation, so beyond this the nearest sample is
// used instead. Real sampled functions have one or two inputs; a shading's
// has exactly one.
MAX_INTERPOLATED_INPUTS :: 8

// Operands a type 4 program may stack. ISO 32000-2, 7.10.5 fixes this at 100.
MAX_PS_STACK :: 100
// Instructions in a compiled type 4 program.
MAX_PS_OPS   :: 1 << 16
// Nesting of { } in a type 4 program.
MAX_PS_DEPTH :: 32
// Instructions one evaluation may execute. A type 4 function has no loop
// construct, so a real program always finishes in fewer steps than it has
// instructions; this only ever fires on a program built to spin.
MAX_PS_STEPS :: 10000

// ---------------------------------------------------------------------------
// Loading
// ---------------------------------------------------------------------------

// function_load builds the function `obj` describes -- a dictionary for types
// 2 and 3, a stream for types 0 and 4.
//
// The result is owned by the Document and must not be freed. It is cached,
// keyed by the object it was built from, because a shading that names the same
// tint transform from every one of its patterns must build it once.
function_load :: proc(doc: ^Document, obj: Object) -> (^Function, Error) {
    return fn_load(doc, obj, 0)
}

@(private)
fn_load :: proc(doc: ^Document, obj: Object, depth: int) -> (^Function, Error) {
    if doc == nil {
        return nil, .Invalid_Argument
    }
    if depth > MAX_FUNCTION_DEPTH {
        add_repair(doc, .Recursion_Limit_Hit, 0)
        return nil, .Recursion_Limit
    }

    o := resolve(doc, obj)
    key := object_identity(o)
    if key != 0 {
        if cached, ok := doc.fn_cache[key]; ok {
            return cached, .None
        }
    }

    d, is_dict := to_dict(o)
    if !is_dict {
        return nil, .Bad_Function
    }
    stm, _ := to_stream(o)

    kind: Function_Kind
    switch dict_int(doc, d, "FunctionType", -1) {
    case 0: kind = .Sampled
    case 2: kind = .Exponential
    case 3: kind = .Stitching
    case 4: kind = .PostScript
    case:   return nil, .Bad_Function
    }

    f := new(Function, doc.allocator)
    f.kind   = kind
    f.domain = fn_num_array(doc, d, "Domain")
    if len(f.domain) < 2 {
        // /Domain is required, but a function that omits it is still usable:
        // every type but 0 has an obvious default and a shading feeding it
        // will clamp to [0 1] anyway.
        f.domain = fn_make_pair(doc, 0, 1)
        add_repair(doc, .Bad_Resource, 0)
    }
    f.n_in  = min(len(f.domain) / 2, MAX_FUNCTION_INPUTS)
    f.range = fn_num_array(doc, d, "Range")
    f.n_out = min(len(f.range) / 2, MAX_FUNCTION_OUTPUTS)

    err: Error
    switch kind {
    case .Sampled:     err = fn_load_sampled(doc, f, d, stm)
    case .Exponential: err = fn_load_exponential(doc, f, d)
    case .Stitching:   err = fn_load_stitching(doc, f, d, depth)
    case .PostScript:  err = fn_load_postscript(doc, f, stm)
    }
    if err != .None {
        return nil, err
    }
    if f.n_in < 1 || f.n_out < 1 {
        return nil, .Bad_Function
    }

    // Entered only now that it is complete: a half-built function in the cache
    // is how a self-referential /Functions array would become a cycle.
    if key != 0 {
        if doc.fn_cache == nil {
            doc.fn_cache = make(map[uintptr]^Function, 16, doc.allocator)
        }
        doc.fn_cache[key] = f
    }
    return f, .None
}

@(private)
fn_num_array :: proc(doc: ^Document, d: ^Dict, key: string) -> []f64 {
    a := dict_array(doc, d, key)
    n := array_len(a)
    if n <= 0 {
        return nil
    }
    n = min(n, MAX_FUNCTION_ARRAY)
    out := make([]f64, n, doc.allocator)
    for i in 0 ..< n {
        // A non-finite bound would propagate into every value the function
        // produces and from there into the geometry.
        out[i] = clamp_finite(array_num(doc, a, i, 0), 0)
    }
    return out
}

@(private)
fn_make_pair :: proc(doc: ^Document, lo, hi: f64) -> []f64 {
    out := make([]f64, 2, doc.allocator)
    out[0], out[1] = lo, hi
    return out
}

@(private)
fn_make_single :: proc(doc: ^Document, v: f64) -> []f64 {
    out := make([]f64, 1, doc.allocator)
    out[0] = v
    return out
}

// ---------------------------------------------------------------------------
// Type 0: sampled
// ---------------------------------------------------------------------------

@(private)
fn_load_sampled :: proc(doc: ^Document, f: ^Function, d: ^Dict, stm: ^Stream) -> Error {
    if stm == nil {
        return .Bad_Function
    }
    // /Range is what says how many outputs a sample has; without it the
    // sample table cannot be read at all.
    if f.n_out < 1 {
        return .Bad_Function
    }

    switch dict_int(doc, d, "BitsPerSample", 0) {
    case 1, 2, 4, 8, 12, 16, 24, 32:
        f.bits_per_sample = int(dict_int(doc, d, "BitsPerSample", 8))
    case:
        return .Bad_Function
    }

    sizes := dict_array(doc, d, "Size")
    n := min(array_len(sizes), MAX_FUNCTION_INPUTS)
    if n < 1 {
        return .Bad_Function
    }
    // /Size and /Domain must agree on how many inputs there are. Trusting the
    // shorter of the two keeps every later index in range.
    f.n_in = min(f.n_in, n)
    f.size = make([]int, f.n_in, doc.allocator)
    total := i64(1)
    for i in 0 ..< f.n_in {
        s := array_int(doc, sizes, i, 0)
        if s < 1 {
            // A zero or negative dimension would make the table empty and
            // every index into it invalid.
            s = 1
            add_repair(doc, .Bad_Resource, 0)
        }
        // A single dimension larger than the whole table may be is not a
        // table this library will read; nor is one whose dimensions multiply
        // out past the limit.
        if s > MAX_FUNCTION_SAMPLES {
            return .Bad_Function
        }
        f.size[i] = int(s)
        total *= s
        if total > MAX_FUNCTION_SAMPLES {
            return .Bad_Function
        }
    }
    if total * i64(f.n_out) > MAX_FUNCTION_SAMPLES {
        return .Bad_Function
    }

    f.order = int(dict_int(doc, d, "Order", 1))

    f.encode = fn_num_array(doc, d, "Encode")
    if len(f.encode) < 2 * f.n_in {
        // Default /Encode is [0 Size0-1 0 Size1-1 ...] (Table 38).
        enc := make([]f64, 2 * f.n_in, doc.allocator)
        for i in 0 ..< f.n_in {
            enc[2 * i]     = 0
            enc[2 * i + 1] = f64(f.size[i] - 1)
        }
        f.encode = enc
    }

    f.decode = fn_num_array(doc, d, "Decode")
    if len(f.decode) < 2 * f.n_out {
        f.decode = f.range
    }

    data, derr := stream_data(doc, stm)
    if derr != .None {
        return .Bad_Function
    }
    f.samples = data

    have := i64(len(data)) * 8
    need := total * i64(f.n_out) * i64(f.bits_per_sample)
    if have < need {
        // Short sample tables do occur, usually from a truncated stream. The
        // missing samples read as zero rather than the whole function failing.
        add_repair(doc, .Bad_Resource, 0, stm.ref)
    }
    return .None
}

// fn_sample reads sample number `ord` (already including the output index) out
// of the packed table. Out-of-range ordinals read as 0, which is what makes a
// truncated table safe to evaluate.
@(private)
fn_sample :: proc(f: ^Function, ord: i64) -> f64 {
    if ord < 0 {
        return 0
    }
    bps := f.bits_per_sample
    bit := ord * i64(bps)
    if bit + i64(bps) > i64(len(f.samples)) * 8 {
        return 0
    }
    if bps == 8 {
        return f64(f.samples[bit >> 3])
    }
    v := u64(0)
    for k in 0 ..< bps {
        b := bit + i64(k)
        byte_val := f.samples[b >> 3]
        bit_val := (byte_val >> (7 - u8(b & 7))) & 1
        v = v << 1 | u64(bit_val)
    }
    return f64(v)
}

@(private)
fn_eval_sampled :: proc(f: ^Function, in_: []f64, out: []f64) {
    if len(f.size) < f.n_in || f.n_in < 1 {
        return
    }

    // Position within the table, split into the sample below and the fraction
    // between it and the next one.
    base:  [MAX_FUNCTION_INPUTS]int
    frac:  [MAX_FUNCTION_INPUTS]f64
    for i in 0 ..< f.n_in {
        x := in_[i] if i < len(in_) else 0
        e := fn_interp(x, f.domain[2 * i], f.domain[2 * i + 1],
                          f.encode[2 * i], f.encode[2 * i + 1])
        e = clamp(clamp_finite(e, 0), 0, f64(f.size[i] - 1))
        lo := int(math.floor(e))
        if lo > f.size[i] - 1 { lo = f.size[i] - 1 }
        if lo < 0 { lo = 0 }
        base[i] = lo
        frac[i] = e - f64(lo)
    }

    // Strides through the flattened table: the first input varies fastest
    // (ISO 32000-2, 7.10.2).
    stride: [MAX_FUNCTION_INPUTS]i64
    s := i64(1)
    for i in 0 ..< f.n_in {
        stride[i] = s
        s *= i64(f.size[i])
    }

    max_val := f64(u64(1) << u64(f.bits_per_sample)) - 1
    if f.bits_per_sample == 32 {
        max_val = 4294967295
    }

    // /Order 3 asks for a cubic spline. Linear interpolation is used instead:
    // the difference is invisible in the gradients these tables encode, and
    // the format explicitly allows a reader to substitute a lower order.
    corners := 1 << uint(f.n_in) if f.n_in <= MAX_INTERPOLATED_INPUTS else 1

    for j in 0 ..< f.n_out {
        if j >= len(out) {
            break
        }
        acc := 0.0
        for corner in 0 ..< corners {
            weight := 1.0
            idx := i64(0)
            for i in 0 ..< f.n_in {
                bit := (corner >> uint(i)) & 1
                coord := base[i]
                if bit == 1 {
                    weight *= frac[i]
                    if coord + 1 <= f.size[i] - 1 {
                        coord += 1
                    }
                } else {
                    weight *= 1 - frac[i]
                }
                idx += i64(coord) * stride[i]
            }
            if weight == 0 {
                continue
            }
            acc += weight * fn_sample(f, idx * i64(f.n_out) + i64(j))
        }

        lo, hi := 0.0, 1.0
        if 2 * j + 1 < len(f.decode) {
            lo, hi = f.decode[2 * j], f.decode[2 * j + 1]
        }
        out[j] = lo + acc * (hi - lo) / max_val
    }
}

// ---------------------------------------------------------------------------
// Type 2: exponential interpolation
// ---------------------------------------------------------------------------

@(private)
fn_load_exponential :: proc(doc: ^Document, f: ^Function, d: ^Dict) -> Error {
    f.c0 = fn_num_array(doc, d, "C0")
    f.c1 = fn_num_array(doc, d, "C1")
    // C0 defaults to [0.0] and C1 to [1.0] (ISO 32000-2, Table 40).
    if len(f.c0) == 0 {
        f.c0 = fn_make_single(doc, 0)
    }
    if len(f.c1) == 0 {
        f.c1 = fn_make_single(doc, 1)
    }
    f.exponent = clamp_finite(dict_num(doc, d, "N", 1), 1)

    n := min(len(f.c0), len(f.c1))
    n = min(n, MAX_FUNCTION_OUTPUTS)
    if n < 1 {
        return .Bad_Function
    }
    // /Range is optional here; when it is absent the number of outputs is the
    // length of C0 and C1.
    if f.n_out == 0 || f.n_out > n {
        f.n_out = n
    }
    // A type 2 function takes exactly one input (7.10.3).
    f.n_in = 1
    return .None
}

@(private)
fn_eval_exponential :: proc(f: ^Function, in_: []f64, out: []f64) {
    x := len(in_) > 0 ? in_[0] : 0
    p := x
    if f.exponent != 1 {
        p = math.pow(x, f.exponent)
        // A negative base under a fractional exponent is not a real number;
        // the format says the domain must exclude it, so a file that does not
        // gets zero rather than a NaN travelling into the raster.
        if !is_finite(p) {
            p = 0
        }
    }
    for j in 0 ..< f.n_out {
        if j >= len(out) {
            break
        }
        a := f.c0[j] if j < len(f.c0) else 0
        b := f.c1[j] if j < len(f.c1) else 1
        out[j] = a + p * (b - a)
    }
}

// ---------------------------------------------------------------------------
// Type 3: stitching
// ---------------------------------------------------------------------------

@(private)
fn_load_stitching :: proc(doc: ^Document, f: ^Function, d: ^Dict, depth: int) -> Error {
    subs := dict_array(doc, d, "Functions")
    k := min(array_len(subs), MAX_SUB_FUNCTIONS)
    if k < 1 {
        return .Bad_Function
    }
    f.functions = make([]^Function, k, doc.allocator)
    n_out := 0
    for i in 0 ..< k {
        sub, serr := fn_load(doc, array_get_raw(subs, i), depth + 1)
        if serr != .None {
            // One broken sub-function leaves a hole that evaluates to zero
            // rather than discarding the whole stitch.
            add_repair(doc, .Bad_Resource, 0)
            continue
        }
        f.functions[i] = sub
        if n_out == 0 {
            n_out = sub.n_out
        }
    }
    f.bounds = fn_num_array(doc, d, "Bounds")
    f.encode = fn_num_array(doc, d, "Encode")
    if f.n_out == 0 {
        // /Range is optional for a type 3 function; the sub-functions say how
        // many outputs there are.
        f.n_out = min(n_out, MAX_FUNCTION_OUTPUTS)
    }
    // A type 3 function takes exactly one input (7.10.4).
    f.n_in = 1
    return .None
}

@(private)
fn_eval_stitching :: proc(f: ^Function, in_: []f64, out: []f64) {
    k := len(f.functions)
    if k == 0 {
        return
    }
    x := len(in_) > 0 ? in_[0] : 0
    d0, d1 := f.domain[0], f.domain[1]

    // Subdomain i is [Bounds[i-1], Bounds[i]), so an input exactly on a bound
    // belongs to the later sub-function (ISO 32000-2, 7.10.4).
    i := 0
    for i < k - 1 {
        b := f.bounds[i] if i < len(f.bounds) else d1
        if x < b {
            break
        }
        i += 1
    }

    lo := d0
    if i > 0 && i - 1 < len(f.bounds) {
        lo = f.bounds[i - 1]
    }
    hi := d1
    if i < k - 1 && i < len(f.bounds) {
        hi = f.bounds[i]
    }

    e0, e1 := 0.0, 1.0
    if 2 * i + 1 < len(f.encode) {
        e0, e1 = f.encode[2 * i], f.encode[2 * i + 1]
    }

    sub_in: [1]f64
    sub_in[0] = fn_interp(x, lo, hi, e0, e1)
    function_eval(f.functions[i], sub_in[:], out)
}

// ---------------------------------------------------------------------------
// Type 4: the PostScript calculator
// ---------------------------------------------------------------------------

@(private)
fn_load_postscript :: proc(doc: ^Document, f: ^Function, stm: ^Stream) -> Error {
    if stm == nil || f.n_out < 1 {
        // Without /Range there is no way to know how many of the values left
        // on the stack are outputs.
        return .Bad_Function
    }
    src, derr := stream_data(doc, stm)
    if derr != .None || len(src) == 0 {
        return .Bad_Function
    }

    // Compiled into scratch memory first: the Document's arena cannot release
    // the intermediate blocks a growing array leaves behind.
    ops := make([dynamic]Ps_Op, 0, 64, context.allocator)
    defer delete(ops)

    p := Ps_Lexer{src = src}
    if err := ps_compile(&p, &ops); err != .None {
        return err
    }
    if len(ops) == 0 {
        return .Bad_Function
    }

    prog := make([]Ps_Op, len(ops), doc.allocator)
    copy(prog, ops[:])
    f.program = prog
    return .None
}

@(private)
Ps_Token_Kind :: enum u8 {
    End,
    Number,
    Lbrace,
    Rbrace,
    Operator,
    If,
    Ifelse,
    Unknown,
}

@(private)
Ps_Token :: struct {
    kind:  Ps_Token_Kind,
    op:    Ps_Op_Kind,
    value: f64,
}

@(private)
Ps_Lexer :: struct {
    src: []byte,
    pos: int,
}

@(private)
ps_is_space :: proc "contextless" (c: byte) -> bool {
    return c == ' ' || c == '\t' || c == '\r' || c == '\n' || c == '\f' || c == 0
}

@(private)
ps_next :: proc(p: ^Ps_Lexer) -> Ps_Token {
    for p.pos < len(p.src) {
        c := p.src[p.pos]
        if ps_is_space(c) {
            p.pos += 1
            continue
        }
        // A comment runs to the end of the line, as everywhere else in PDF.
        if c == '%' {
            for p.pos < len(p.src) && p.src[p.pos] != '\n' && p.src[p.pos] != '\r' {
                p.pos += 1
            }
            continue
        }
        break
    }
    if p.pos >= len(p.src) {
        return Ps_Token{kind = .End}
    }

    c := p.src[p.pos]
    if c == '{' {
        p.pos += 1
        return Ps_Token{kind = .Lbrace}
    }
    if c == '}' {
        p.pos += 1
        return Ps_Token{kind = .Rbrace}
    }

    start := p.pos
    for p.pos < len(p.src) {
        b := p.src[p.pos]
        if ps_is_space(b) || b == '{' || b == '}' || b == '%' {
            break
        }
        p.pos += 1
    }
    word := string(p.src[start:p.pos])
    if len(word) == 0 {
        // Cannot happen -- the delimiters above are all consumed as tokens --
        // but a zero-length word would spin the parser, so it ends it.
        p.pos = len(p.src)
        return Ps_Token{kind = .End}
    }

    if c == '+' || c == '-' || c == '.' || (c >= '0' && c <= '9') {
        obj, ok := parse_number(word)
        if ok {
            v, _ := to_number(obj)
            return Ps_Token{kind = .Number, value = clamp_finite(v, 0)}
        }
        return Ps_Token{kind = .Unknown}
    }
    return ps_operator_token(word)
}

@(private)
ps_operator_token :: proc(word: string) -> Ps_Token {
    k: Ps_Op_Kind
    switch word {
    case "add":      k = .Add
    case "sub":      k = .Sub
    case "mul":      k = .Mul
    case "div":      k = .Div
    case "idiv":     k = .Idiv
    case "mod":      k = .Mod
    case "neg":      k = .Neg
    case "abs":      k = .Abs
    case "ceiling":  k = .Ceiling
    case "floor":    k = .Floor
    case "round":    k = .Round
    case "truncate": k = .Truncate
    case "sqrt":     k = .Sqrt
    case "sin":      k = .Sin
    case "cos":      k = .Cos
    case "atan":     k = .Atan
    case "exp":      k = .Exp
    case "ln":       k = .Ln
    case "log":      k = .Log
    case "cvi":      k = .Cvi
    case "cvr":      k = .Cvr
    case "eq":       k = .Eq
    case "ne":       k = .Ne
    case "gt":       k = .Gt
    case "ge":       k = .Ge
    case "lt":       k = .Lt
    case "le":       k = .Le
    case "and":      k = .And
    case "or":       k = .Or
    case "xor":      k = .Xor
    case "not":      k = .Not
    case "bitshift": k = .Bitshift
    case "true":     k = .True
    case "false":    k = .False
    case "pop":      k = .Pop
    case "exch":     k = .Exch
    case "dup":      k = .Dup
    case "copy":     k = .Copy
    case "index":    k = .Index
    case "roll":     k = .Roll
    case "if":       return Ps_Token{kind = .If}
    case "ifelse":   return Ps_Token{kind = .Ifelse}
    case:            return Ps_Token{kind = .Unknown}
    }
    return Ps_Token{kind = .Operator, op = k}
}

// ps_compile reads the whole program. The outer braces are expected but not
// required: a stream that forgot them still compiles, because the alternative
// is refusing a function every other viewer runs.
@(private)
ps_compile :: proc(p: ^Ps_Lexer, ops: ^[dynamic]Ps_Op) -> Error {
    save := p.pos
    t := ps_next(p)
    if t.kind == .Lbrace {
        return ps_compile_block(p, ops, 0)
    }
    p.pos = save
    return ps_compile_block(p, ops, 0)
}

// ps_compile_block emits instructions until the '}' that closes the current
// block, which it consumes. Conditionals become jumps with absolute targets:
//
//     {A} if           ->  jump_if_false end;  A;  end:
//     {A} {B} ifelse   ->  jump_if_false b;    A;  jump end;  b: B;  end:
//
// so a slot for the branch is reserved before the block body is emitted and
// patched once the operator that follows the body says what it is.
@(private)
ps_compile_block :: proc(p: ^Ps_Lexer, ops: ^[dynamic]Ps_Op, depth: int) -> Error {
    if depth > MAX_PS_DEPTH {
        return .Bad_Function
    }
    for {
        if len(ops) > MAX_PS_OPS {
            return .Bad_Function
        }
        t := ps_next(p)
        switch t.kind {
        case .End:
            // End of stream: a missing '}' is not worth rejecting a program
            // that is otherwise complete.
            return .None
        case .Rbrace:
            return .None
        case .Number:
            append(ops, Ps_Op{kind = .Push, value = t.value})
        case .Operator:
            append(ops, Ps_Op{kind = t.op})
        case .Lbrace:
            branch := len(ops)
            append(ops, Ps_Op{kind = .Jump})
            ps_compile_block(p, ops, depth + 1) or_return

            save := p.pos
            t2 := ps_next(p)
            #partial switch t2.kind {
            case .If:
                ops[branch] = Ps_Op{kind = .Jump_If_False, target = len(ops)}
            case .Lbrace:
                skip := len(ops)
                append(ops, Ps_Op{kind = .Jump})
                ps_compile_block(p, ops, depth + 1) or_return
                t3 := ps_next(p)
                if t3.kind != .Ifelse {
                    return .Bad_Function
                }
                ops[branch] = Ps_Op{kind = .Jump_If_False, target = skip + 1}
                ops[skip]   = Ps_Op{kind = .Jump, target = len(ops)}
            case:
                // A block that no `if` or `ifelse` consumes never runs in
                // PostScript either, so it is jumped over rather than
                // executed.
                ops[branch] = Ps_Op{kind = .Jump, target = len(ops)}
                p.pos = save
            }
        case .If, .Ifelse, .Unknown:
            // An `if` with no block before it, or a word this calculator has
            // no operator for. Skipping it keeps the rest of the program
            // running, which is what a partly-understood tint transform needs.
        }
    }
}

// Ps_Value is one operand. The type 4 language distinguishes booleans from
// numbers -- `and`, `or`, `xor` and `not` are logical on the first and bitwise
// on the second -- so the tag travels with the value.
@(private)
Ps_Value :: struct {
    v:       f64,
    is_bool: bool,
}

@(private)
Ps_Stack :: struct {
    v:  [MAX_PS_STACK]Ps_Value,
    sp: int,
    // Set when a push found the stack full. Evaluation stops, leaving whatever
    // has been computed so far rather than discarding it.
    overflow: bool,
}

@(private)
ps_push :: proc(s: ^Ps_Stack, val: f64, is_bool := false) {
    if s.sp >= MAX_PS_STACK {
        s.overflow = true
        return
    }
    s.v[s.sp] = Ps_Value{v = clamp_finite(val, 0), is_bool = is_bool}
    s.sp += 1
}

// ps_pop returns zero on an empty stack. An underflowing program is malformed,
// but a defined result costs nothing and a crash costs the page.
@(private)
ps_pop :: proc(s: ^Ps_Stack) -> Ps_Value {
    if s.sp <= 0 {
        return Ps_Value{}
    }
    s.sp -= 1
    return s.v[s.sp]
}

@(private)
ps_int :: proc(x: f64) -> i32 {
    if !is_finite(x) {
        return 0
    }
    y := math.trunc(x)
    if y > 2147483647 {
        return 2147483647
    }
    if y < -2147483648 {
        return -2147483648
    }
    return i32(y)
}

@(private)
ps_exec :: proc(prog: []Ps_Op, s: ^Ps_Stack) {
    pc := 0
    for steps := 0; steps < MAX_PS_STEPS; steps += 1 {
        if pc < 0 || pc >= len(prog) || s.overflow {
            return
        }
        op := prog[pc]
        pc += 1

        switch op.kind {
        case .Push:
            ps_push(s, op.value)

        case .Add:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v + b.v)
        case .Sub:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v - b.v)
        case .Mul:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v * b.v)
        case .Div:
            b, a := ps_pop(s), ps_pop(s)
            // Division by zero is undefined in PostScript; zero is the value
            // every other reader substitutes.
            ps_push(s, b.v != 0 ? a.v / b.v : 0)
        case .Idiv:
            b, a := ps_pop(s), ps_pop(s)
            ib, ia := ps_int(b.v), ps_int(a.v)
            // -2147483648 idiv -1 overflows a two's-complement divide.
            ps_push(s, ib != 0 && !(ia == -2147483648 && ib == -1) ? f64(ia / ib) : 0)
        case .Mod:
            b, a := ps_pop(s), ps_pop(s)
            ib, ia := ps_int(b.v), ps_int(a.v)
            ps_push(s, ib != 0 && !(ia == -2147483648 && ib == -1) ? f64(ia % ib) : 0)
        case .Neg:
            a := ps_pop(s)
            ps_push(s, -a.v)
        case .Abs:
            a := ps_pop(s)
            ps_push(s, abs(a.v))
        case .Ceiling:
            a := ps_pop(s)
            ps_push(s, math.ceil(a.v))
        case .Floor:
            a := ps_pop(s)
            ps_push(s, math.floor(a.v))
        case .Round:
            a := ps_pop(s)
            // PostScript rounds a halfway case towards positive infinity.
            ps_push(s, math.floor(a.v + 0.5))
        case .Truncate:
            a := ps_pop(s)
            ps_push(s, math.trunc(a.v))
        case .Sqrt:
            a := ps_pop(s)
            ps_push(s, a.v > 0 ? math.sqrt(a.v) : 0)
        case .Sin:
            a := ps_pop(s)
            ps_push(s, math.sin(a.v * math.PI / 180))
        case .Cos:
            a := ps_pop(s)
            ps_push(s, math.cos(a.v * math.PI / 180))
        case .Atan:
            den, num := ps_pop(s), ps_pop(s)
            deg := math.atan2(num.v, den.v) * 180 / math.PI
            // `atan` answers in degrees between 0 and 360, never negative.
            if deg < 0 {
                deg += 360
            }
            ps_push(s, deg)
        case .Exp:
            e, base := ps_pop(s), ps_pop(s)
            r := math.pow(base.v, e.v)
            ps_push(s, is_finite(r) ? r : 0)
        case .Ln:
            a := ps_pop(s)
            ps_push(s, a.v > 0 ? math.ln(a.v) : 0)
        case .Log:
            a := ps_pop(s)
            ps_push(s, a.v > 0 ? math.log10(a.v) : 0)
        case .Cvi:
            a := ps_pop(s)
            ps_push(s, math.trunc(a.v))
        case .Cvr:
            a := ps_pop(s)
            ps_push(s, a.v)

        case .Eq:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v == b.v ? 1 : 0, true)
        case .Ne:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v != b.v ? 1 : 0, true)
        case .Gt:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v > b.v ? 1 : 0, true)
        case .Ge:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v >= b.v ? 1 : 0, true)
        case .Lt:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v < b.v ? 1 : 0, true)
        case .Le:
            b, a := ps_pop(s), ps_pop(s)
            ps_push(s, a.v <= b.v ? 1 : 0, true)
        case .And:
            b, a := ps_pop(s), ps_pop(s)
            if a.is_bool && b.is_bool {
                ps_push(s, (a.v != 0 && b.v != 0) ? 1 : 0, true)
            } else {
                ps_push(s, f64(ps_int(a.v) & ps_int(b.v)))
            }
        case .Or:
            b, a := ps_pop(s), ps_pop(s)
            if a.is_bool && b.is_bool {
                ps_push(s, (a.v != 0 || b.v != 0) ? 1 : 0, true)
            } else {
                ps_push(s, f64(ps_int(a.v) | ps_int(b.v)))
            }
        case .Xor:
            b, a := ps_pop(s), ps_pop(s)
            if a.is_bool && b.is_bool {
                ps_push(s, ((a.v != 0) != (b.v != 0)) ? 1 : 0, true)
            } else {
                ps_push(s, f64(ps_int(a.v) ~ ps_int(b.v)))
            }
        case .Not:
            a := ps_pop(s)
            if a.is_bool {
                ps_push(s, a.v == 0 ? 1 : 0, true)
            } else {
                ps_push(s, f64(~ps_int(a.v)))
            }
        case .Bitshift:
            sh, a := ps_pop(s), ps_pop(s)
            n := ps_int(sh.v)
            bits := u32(ps_int(a.v))
            // A shift of 32 or more is undefined in C and in Odin; the value
            // it must produce here is zero.
            switch {
            case n >= 32 || n <= -32: ps_push(s, 0)
            case n >= 0:              ps_push(s, f64(i32(bits << u32(n))))
            case:                     ps_push(s, f64(i32(bits >> u32(-n))))
            }
        case .True:
            ps_push(s, 1, true)
        case .False:
            ps_push(s, 0, true)

        case .Pop:
            ps_pop(s)
        case .Exch:
            b, a := ps_pop(s), ps_pop(s)
            if s.sp + 2 <= MAX_PS_STACK {
                s.v[s.sp] = b
                s.v[s.sp + 1] = a
                s.sp += 2
            } else {
                s.overflow = true
            }
        case .Dup:
            if s.sp > 0 {
                top := s.v[s.sp - 1]
                ps_push(s, top.v, top.is_bool)
            } else {
                ps_push(s, 0)
            }
        case .Copy:
            n := int(ps_int(ps_pop(s).v))
            if n <= 0 {
                break
            }
            if n > s.sp || s.sp + n > MAX_PS_STACK {
                s.overflow = true
                break
            }
            src := s.sp - n
            for i in 0 ..< n {
                s.v[s.sp + i] = s.v[src + i]
            }
            s.sp += n
        case .Index:
            n := int(ps_int(ps_pop(s).v))
            if n < 0 || n >= s.sp {
                // Indexing past the bottom of the stack is malformed; zero
                // keeps the program running.
                ps_push(s, 0)
                break
            }
            val := s.v[s.sp - 1 - n]
            ps_push(s, val.v, val.is_bool)
        case .Roll:
            j := int(ps_int(ps_pop(s).v))
            n := int(ps_int(ps_pop(s).v))
            if n <= 0 || n > s.sp {
                break
            }
            // Rotation by three reversals, so no scratch stack is needed.
            k := j %% n
            if k != 0 {
                lo := s.sp - n
                ps_reverse(s.v[lo:lo + n - k])
                ps_reverse(s.v[lo + n - k:lo + n])
                ps_reverse(s.v[lo:lo + n])
            }

        case .Jump:
            pc = op.target
        case .Jump_If_False:
            cond := ps_pop(s)
            if cond.v == 0 {
                pc = op.target
            }
        }
    }
}

@(private)
ps_reverse :: proc(vals: []Ps_Value) {
    for i, j := 0, len(vals) - 1; i < j; i, j = i + 1, j - 1 {
        vals[i], vals[j] = vals[j], vals[i]
    }
}

@(private)
fn_eval_postscript :: proc(f: ^Function, in_: []f64, out: []f64) {
    s: Ps_Stack
    for i in 0 ..< f.n_in {
        ps_push(&s, i < len(in_) ? in_[i] : 0)
    }
    ps_exec(f.program, &s)

    // The outputs are the topmost n_out values, bottom-most first. A program
    // that left fewer behind gets zeros for the rest.
    for j in 0 ..< f.n_out {
        if j >= len(out) {
            break
        }
        idx := s.sp - f.n_out + j
        out[j] = idx >= 0 && idx < s.sp ? s.v[idx].v : 0
    }
}

// ---------------------------------------------------------------------------
// Evaluation
// ---------------------------------------------------------------------------

// function_eval evaluates `fn` at `in_` and writes n_out values into `out`.
//
// Inputs are clamped to /Domain and outputs to /Range, both of which the
// format requires (ISO 32000-2, 7.10.1) and neither of which is optional in
// practice: a shading whose tint escapes its range paints garbage. `out` may
// be shorter than n_out, in which case the extra outputs are dropped; it is
// always fully written, so a nil or broken function reads as zeros rather than
// as whatever the caller's buffer held.
function_eval :: proc(fn: ^Function, in_: []f64, out: []f64) {
    for i in 0 ..< len(out) {
        out[i] = 0
    }
    if fn == nil || fn.n_out < 1 {
        return
    }

    clamped: [MAX_FUNCTION_INPUTS]f64
    n := min(fn.n_in, MAX_FUNCTION_INPUTS)
    for i in 0 ..< n {
        x := i < len(in_) ? clamp_finite(in_[i], 0) : 0
        if 2 * i + 1 < len(fn.domain) {
            lo, hi := fn.domain[2 * i], fn.domain[2 * i + 1]
            if lo > hi {
                lo, hi = hi, lo
            }
            x = clamp(x, lo, hi)
        }
        clamped[i] = x
    }

    switch fn.kind {
    case .Sampled:     fn_eval_sampled(fn, clamped[:n], out)
    case .Exponential: fn_eval_exponential(fn, clamped[:n], out)
    case .Stitching:   fn_eval_stitching(fn, clamped[:n], out)
    case .PostScript:  fn_eval_postscript(fn, clamped[:n], out)
    }

    for j in 0 ..< min(fn.n_out, len(out)) {
        out[j] = clamp_finite(out[j], 0)
        if 2 * j + 1 < len(fn.range) {
            lo, hi := fn.range[2 * j], fn.range[2 * j + 1]
            if lo > hi {
                lo, hi = hi, lo
            }
            out[j] = clamp(out[j], lo, hi)
        }
    }
}

// fn_interp maps x from [xmin xmax] onto [ymin ymax]. A degenerate source
// interval -- which /Bounds and /Encode arrays out of real files do produce --
// maps to the start of the destination rather than dividing by zero.
@(private)
fn_interp :: proc(x, xmin, xmax, ymin, ymax: f64) -> f64 {
    if xmax == xmin {
        return ymin
    }
    return ymin + (x - xmin) * (ymax - ymin) / (xmax - xmin)
}