malachite-nz 0.13.0

// Copyright © 2026 Mikhail Hogrefe
//
// This file is part of Malachite.
//
// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.

//! Threshold tuning, after GMP's `tune/tuneup.c`.
//!
//! For each threshold, two algorithm variants are measured head-to-head at increasing sizes and
//! the crossover is located. Like GMP, this sidesteps the "recursive thresholds must be tuned
//! simultaneously" problem by tuning bottom-up: each threshold is finalized (edited into
//! `platform_64.rs`, library rebuilt) before any threshold above it is measured. Since algorithms
//! near a crossover are nearly equal by definition, modest error in earlier levels barely
//! perturbs later ones; iterate to a fixpoint (2-3 passes) if paranoid.
//!
//! The measurement machinery (batched best-of-k interleaved timing and the crossover scan) is
//! shared with the other crates' tuners: see `malachite_base::test_util::bench::tune`.
//!
//! Usage: `cargo run --release --features bin_build -p malachite-nz -- -g tune_mul`
//! (acquire perf/bench-lock.sh first; results are garbage on a busy machine)

use malachite_base::num::basic::integers::PrimitiveInt;
use malachite_base::num::random::random_primitive_ints;
use malachite_base::random::EXAMPLE_SEED;
use malachite_base::test_util::bench::tune::{
    INPUT_SETS, find_crossover_spec, interleaved_min_pair, time_batch,
};
use malachite_nz::natural::arithmetic::add_mul::limbs_slice_add_mul_limb_same_length_in_place_left;
use malachite_nz::natural::arithmetic::div::{
    limbs_div_barrett_approx, limbs_div_barrett_approx_scratch_len,
    limbs_div_divide_and_conquer_approx, limbs_div_schoolbook_approx,
};
use malachite_nz::natural::arithmetic::div_exact::{
    limbs_modular_div_barrett, limbs_modular_div_barrett_scratch_len,
    limbs_modular_div_divide_and_conquer, limbs_modular_div_mod_barrett,
    limbs_modular_div_mod_barrett_scratch_len, limbs_modular_div_mod_divide_and_conquer,
    limbs_modular_div_mod_schoolbook, limbs_modular_div_schoolbook, limbs_modular_invert,
    limbs_modular_invert_limb, limbs_modular_invert_scratch_len,
};
use malachite_nz::natural::arithmetic::div_mod::{
    limbs_div_mod_barrett, limbs_div_mod_barrett_scratch_len, limbs_div_mod_divide_and_conquer,
    limbs_div_mod_schoolbook, limbs_div_mod_to_out, limbs_invert_approx,
    limbs_invert_approx_scratch_len, limbs_invert_basecase_approx, limbs_invert_newton_approx,
    limbs_two_limb_inverse_helper,
};
use malachite_nz::natural::arithmetic::mod_power_of_2_square::{
    limbs_square_low_basecase, limbs_square_low_divide_and_conquer, limbs_square_low_scratch_len,
};
use malachite_nz::natural::arithmetic::mul::fft::{
    mpn_mul_fft_for_tuning, mpn_square_fft_for_tuning,
};
use malachite_nz::natural::arithmetic::mul::limbs_mul_greater_to_out_basecase;
use malachite_nz::natural::arithmetic::mul::mul_low::{
    limbs_mul_low_same_length_basecase, limbs_mul_low_same_length_divide_and_conquer,
    limbs_mul_low_same_length_divide_and_conquer_scratch_len,
};
use malachite_nz::natural::arithmetic::mul::toom::{
    TUNE_PROGRAM_BUILD, limbs_mul_greater_to_out_toom_6h,
    limbs_mul_greater_to_out_toom_6h_input_sizes_valid,
    limbs_mul_greater_to_out_toom_6h_scratch_len, limbs_mul_greater_to_out_toom_8h,
    limbs_mul_greater_to_out_toom_8h_input_sizes_valid,
    limbs_mul_greater_to_out_toom_8h_scratch_len, limbs_mul_greater_to_out_toom_22,
    limbs_mul_greater_to_out_toom_22_input_sizes_valid,
    limbs_mul_greater_to_out_toom_22_scratch_len, limbs_mul_greater_to_out_toom_32,
    limbs_mul_greater_to_out_toom_32_input_sizes_valid,
    limbs_mul_greater_to_out_toom_32_scratch_len, limbs_mul_greater_to_out_toom_33,
    limbs_mul_greater_to_out_toom_33_input_sizes_valid,
    limbs_mul_greater_to_out_toom_33_scratch_len, limbs_mul_greater_to_out_toom_42,
    limbs_mul_greater_to_out_toom_42_input_sizes_valid,
    limbs_mul_greater_to_out_toom_42_scratch_len, limbs_mul_greater_to_out_toom_43,
    limbs_mul_greater_to_out_toom_43_input_sizes_valid,
    limbs_mul_greater_to_out_toom_43_scratch_len, limbs_mul_greater_to_out_toom_44,
    limbs_mul_greater_to_out_toom_44_input_sizes_valid,
    limbs_mul_greater_to_out_toom_44_scratch_len, limbs_mul_greater_to_out_toom_53,
    limbs_mul_greater_to_out_toom_53_input_sizes_valid,
    limbs_mul_greater_to_out_toom_53_scratch_len, limbs_mul_greater_to_out_toom_63,
    limbs_mul_greater_to_out_toom_63_input_sizes_valid,
    limbs_mul_greater_to_out_toom_63_scratch_len,
};
use malachite_nz::natural::arithmetic::mul::{
    limbs_mul_greater_to_out, limbs_mul_greater_to_out_scratch_len,
};
use malachite_nz::natural::arithmetic::square::{
    SQR_TOOM2_THRESHOLD, limbs_square_to_out_basecase, limbs_square_to_out_toom_2,
    limbs_square_to_out_toom_2_scratch_len, limbs_square_to_out_toom_3,
    limbs_square_to_out_toom_3_scratch_len, limbs_square_to_out_toom_4,
    limbs_square_to_out_toom_4_scratch_len, limbs_square_to_out_toom_6,
    limbs_square_to_out_toom_6_scratch_len, limbs_square_to_out_toom_8,
    limbs_square_to_out_toom_8_scratch_len,
};
use malachite_nz::natural::conversion::digits::general_digits::{
    GET_STR_PRECOMPUTE_THRESHOLD, limbs_to_digits_small_base_basecase,
};
use malachite_nz::platform::Limb;
use std::hint::black_box;

type MulFn<'a> = &'a dyn Fn(&mut [Limb], &[Limb], &[Limb], &mut [Limb]);

// A mul-shaped algorithm: validity predicate, scratch size, and the routine itself.
struct Algo<'a> {
    name: &'a str,
    valid: &'a dyn Fn(usize) -> bool,
    scratch_len: &'a dyn Fn(usize) -> usize,
    run: MulFn<'a>,
}

// Measure two mul-shaped algorithms at balanced size n x n on identical, rotating input sets.
fn measure_mul_pair(n: usize, a: &Algo, b: &Algo) -> Option<(f64, f64)> {
    if !(a.valid)(n) || !(b.valid)(n) {
        return None;
    }
    let inputs: Vec<(Vec<Limb>, Vec<Limb>)> = (0..INPUT_SETS)
        .map(|k| {
            let xs = random_primitive_ints(EXAMPLE_SEED.fork(&format!("x{k}")))
                .take(n)
                .collect();
            let ys = random_primitive_ints(EXAMPLE_SEED.fork(&format!("y{k}")))
                .take(n)
                .collect();
            (xs, ys)
        })
        .collect();
    let mut out_a = vec![0; n << 1];
    let mut out_b = vec![0; n << 1];
    let mut scratch_a = vec![0; (a.scratch_len)(n)];
    let mut scratch_b = vec![0; (b.scratch_len)(n)];
    // Warmup: fault in pages, settle the core.
    for (xs, ys) in &inputs {
        (a.run)(&mut out_a, xs, ys, &mut scratch_a);
        (b.run)(&mut out_b, xs, ys, &mut scratch_b);
    }
    let (mut i, mut j) = (0usize, 0usize);
    let times = interleaved_min_pair(
        &mut || {
            let (xs, ys) = &inputs[i & (INPUT_SETS - 1)];
            i += 1;
            (a.run)(black_box(&mut out_a), xs, ys, &mut scratch_a);
        },
        &mut || {
            let (xs, ys) = &inputs[j & (INPUT_SETS - 1)];
            j += 1;
            (b.run)(black_box(&mut out_b), xs, ys, &mut scratch_b);
        },
    );
    Some(times)
}

struct Level<'a> {
    threshold_name: &'a str,
    min_size: usize,
    max_size: usize,
    lower: Algo<'a>,
    upper: Algo<'a>,
}

fn find_crossover(c: &Level) {
    find_crossover_spec(
        c.threshold_name,
        "usize",
        c.lower.name,
        c.upper.name,
        c.min_size,
        c.max_size,
        &|n| measure_mul_pair(n, &c.lower, &c.upper),
    );
}

fn basecase_algo<'a>() -> Algo<'a> {
    Algo {
        name: "basecase",
        valid: &|_| true,
        scratch_len: &|_| 0,
        run: &|out, xs, ys, _| limbs_mul_greater_to_out_basecase(out, xs, ys),
    }
}

fn toom22_algo<'a>() -> Algo<'a> {
    Algo {
        name: "toom22",
        valid: &|n| limbs_mul_greater_to_out_toom_22_input_sizes_valid(n, n),
        scratch_len: &|n| limbs_mul_greater_to_out_toom_22_scratch_len(n, n),
        run: &|out, xs, ys, scratch| limbs_mul_greater_to_out_toom_22(out, xs, ys, scratch),
    }
}

fn toom33_algo<'a>() -> Algo<'a> {
    Algo {
        name: "toom33",
        valid: &|n| limbs_mul_greater_to_out_toom_33_input_sizes_valid(n, n),
        scratch_len: &|n| limbs_mul_greater_to_out_toom_33_scratch_len(n, n),
        run: &|out, xs, ys, scratch| limbs_mul_greater_to_out_toom_33(out, xs, ys, scratch),
    }
}

fn toom44_algo<'a>() -> Algo<'a> {
    Algo {
        name: "toom44",
        valid: &|n| limbs_mul_greater_to_out_toom_44_input_sizes_valid(n, n),
        scratch_len: &|n| limbs_mul_greater_to_out_toom_44_scratch_len(n, n),
        run: &|out, xs, ys, scratch| limbs_mul_greater_to_out_toom_44(out, xs, ys, scratch),
    }
}

fn tune_mul_toom22() {
    find_crossover(&Level {
        threshold_name: "MUL_TOOM22_THRESHOLD",
        min_size: 4,
        max_size: 1000,
        lower: basecase_algo(),
        upper: toom22_algo(),
    });
}

fn tune_mul_toom33() {
    find_crossover(&Level {
        threshold_name: "MUL_TOOM33_THRESHOLD",
        min_size: 20,
        max_size: 2000,
        lower: toom22_algo(),
        upper: toom33_algo(),
    });
}

fn tune_mul_toom44() {
    find_crossover(&Level {
        threshold_name: "MUL_TOOM44_THRESHOLD",
        min_size: 60,
        max_size: 4000,
        lower: toom33_algo(),
        upper: toom44_algo(),
    });
}

// The squaring algorithms reuse the mul-shaped `Algo` plumbing, ignoring `ys`. Validity predicates
// replicate each function's split asserts (the sqr functions have no `_input_sizes_valid` helpers).

fn sqr_basecase_algo<'a>() -> Algo<'a> {
    Algo {
        name: "sqr_basecase",
        // The basecase's stack buffer is sized by the compiled-in threshold, so it cannot be
        // measured above it; the crossover scan is capped accordingly. To scan higher, raise
        // SQR_TOOM2_THRESHOLD in platform_64.rs and rebuild.
        valid: &|n| n <= SQR_TOOM2_THRESHOLD,
        scratch_len: &|_| 0,
        run: &|out, xs, _, _| limbs_square_to_out_basecase(out, xs),
    }
}

fn sqr_toom2_algo<'a>() -> Algo<'a> {
    Algo {
        name: "sqr_toom2",
        valid: &|n| n > 1,
        scratch_len: &limbs_square_to_out_toom_2_scratch_len,
        run: &|out, xs, _, scratch| limbs_square_to_out_toom_2(out, xs, scratch),
    }
}

fn sqr_toom3_algo<'a>() -> Algo<'a> {
    Algo {
        name: "sqr_toom3",
        // n = ceil(len / 3), s = len - 2n; s must be in 1..=n.
        valid: &|len| {
            let n = len.div_ceil(3);
            len > n << 1 && len <= 3 * n
        },
        scratch_len: &limbs_square_to_out_toom_3_scratch_len,
        run: &|out, xs, _, scratch| limbs_square_to_out_toom_3(out, xs, scratch),
    }
}

fn sqr_toom4_algo<'a>() -> Algo<'a> {
    Algo {
        name: "sqr_toom4",
        // n = ceil(len / 4), s = len - 3n; s must be in 1..=n.
        valid: &|len| {
            let n = (len + 3) >> 2;
            len > 3 * n && len <= n << 2
        },
        scratch_len: &limbs_square_to_out_toom_4_scratch_len,
        run: &|out, xs, _, scratch| limbs_square_to_out_toom_4(out, xs, scratch),
    }
}

fn tune_sqr_toom2() {
    find_crossover(&Level {
        threshold_name: "SQR_TOOM2_THRESHOLD",
        min_size: 4,
        max_size: SQR_TOOM2_THRESHOLD,
        lower: sqr_basecase_algo(),
        upper: sqr_toom2_algo(),
    });
}

fn tune_sqr_toom3() {
    find_crossover(&Level {
        threshold_name: "SQR_TOOM3_THRESHOLD",
        min_size: 20,
        max_size: 2000,
        lower: sqr_toom2_algo(),
        upper: sqr_toom3_algo(),
    });
}

fn sqr_toom6_algo<'a>() -> Algo<'a> {
    Algo {
        name: "sqr_toom6",
        // n = 1 + (len - 1) / 6, s = len - 5n; needs len >= 18, s in 1..=n, and 10n + 3 <= 2 * len.
        valid: &|len| {
            if len < 18 {
                return false;
            }
            let n = 1 + (len - 1) / 6;
            len > 5 * n && len - 5 * n <= n && 10 * n + 3 <= len << 1
        },
        scratch_len: &limbs_square_to_out_toom_6_scratch_len,
        run: &|out, xs, _, scratch| limbs_square_to_out_toom_6(out, xs, scratch),
    }
}

fn sqr_toom8_algo<'a>() -> Algo<'a> {
    Algo {
        name: "sqr_toom8",
        // n = ceil(len / 8), s = len - 7n; needs len >= 40 and s in 2..=n.
        valid: &|len| {
            if len < 40 {
                return false;
            }
            let n = len.div_ceil(8);
            len > 7 * n + 1 && len - 7 * n <= n
        },
        scratch_len: &limbs_square_to_out_toom_8_scratch_len,
        run: &|out, xs, _, scratch| limbs_square_to_out_toom_8(out, xs, scratch),
    }
}

fn tune_sqr_toom4() {
    find_crossover(&Level {
        threshold_name: "SQR_TOOM4_THRESHOLD",
        min_size: 60,
        max_size: 4000,
        lower: sqr_toom3_algo(),
        upper: sqr_toom4_algo(),
    });
}

fn tune_sqr_toom6() {
    find_crossover(&Level {
        threshold_name: "SQR_TOOM6_THRESHOLD",
        min_size: 200,
        max_size: 6000,
        lower: sqr_toom4_algo(),
        upper: sqr_toom6_algo(),
    });
}

fn toom6h_algo<'a>() -> Algo<'a> {
    Algo {
        name: "toom6h",
        valid: &|n| limbs_mul_greater_to_out_toom_6h_input_sizes_valid(n, n),
        scratch_len: &|n| limbs_mul_greater_to_out_toom_6h_scratch_len(n, n),
        run: &|out, xs, ys, scratch| limbs_mul_greater_to_out_toom_6h(out, xs, ys, scratch),
    }
}

fn toom8h_algo<'a>() -> Algo<'a> {
    Algo {
        name: "toom8h",
        valid: &|n| limbs_mul_greater_to_out_toom_8h_input_sizes_valid(n, n),
        scratch_len: &|n| limbs_mul_greater_to_out_toom_8h_scratch_len(n, n),
        run: &|out, xs, ys, scratch| limbs_mul_greater_to_out_toom_8h(out, xs, ys, scratch),
    }
}

fn tune_mul_toom6h() {
    find_crossover(&Level {
        threshold_name: "MUL_TOOM6H_THRESHOLD",
        min_size: 100,
        max_size: 4000,
        lower: toom44_algo(),
        upper: toom6h_algo(),
    });
}

// toom6h's measured crossover vs toom44 (229) is below toom44's own threshold (a ~315-465 plateau
// vs toom33), suggesting toom44 may have no winning range for balanced mul, as toom4 has none for
// squaring; this measures toom6h against the real incumbent directly.
fn tune_mul_toom6h_vs_toom33() {
    find_crossover(&Level {
        threshold_name: "MUL_TOOM6H_THRESHOLD",
        min_size: 60,
        max_size: 4000,
        lower: toom33_algo(),
        upper: toom6h_algo(),
    });
}

fn tune_mul_toom8h() {
    find_crossover(&Level {
        threshold_name: "MUL_TOOM8H_THRESHOLD",
        min_size: 200,
        max_size: 8000,
        lower: toom6h_algo(),
        upper: toom8h_algo(),
    });
}

// Division algorithms all share the `mpn_sbpi1_div_qr` shape: quotient out, dividend mutated in
// place (its low limbs become the remainder), normalized divisor, precomputed two-limb inverse.
type DivAlgoFn = fn(&mut [Limb], &mut [Limb], &[Limb], Limb) -> bool;

// Measure two division algorithms dividing 2n limbs by n limbs on identical, rotating input sets.
// The dividend is refreshed from a pristine copy before each call (the copy cost is incurred
// identically by both sides).
fn measure_div_pair(n: usize, min_d: usize, a: DivAlgoFn, b: DivAlgoFn) -> Option<(f64, f64)> {
    if n < min_d {
        return None;
    }
    let inputs: Vec<(Vec<Limb>, Vec<Limb>, Limb)> = (0..INPUT_SETS)
        .map(|k| {
            let ns: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("dn{k}")))
                .take(n << 1)
                .collect();
            let mut ds: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("dd{k}")))
                .take(n)
                .collect();
            // The divisor must be normalized (highest bit set).
            ds[n - 1] |= 1 << (Limb::WIDTH - 1);
            let d_inv = limbs_two_limb_inverse_helper(ds[n - 1], ds[n - 2]);
            (ns, ds, d_inv)
        })
        .collect();
    let mut ns_a = vec![0; n << 1];
    let mut ns_b = vec![0; n << 1];
    let mut qs_a = vec![0; n];
    let mut qs_b = vec![0; n];
    // Warmup: fault in pages, settle the core.
    for (ns, ds, d_inv) in &inputs {
        ns_a.copy_from_slice(ns);
        a(&mut qs_a, &mut ns_a, ds, *d_inv);
        ns_b.copy_from_slice(ns);
        b(&mut qs_b, &mut ns_b, ds, *d_inv);
    }
    let (mut i, mut j) = (0usize, 0usize);
    let times = interleaved_min_pair(
        &mut || {
            let (ns, ds, d_inv) = &inputs[i & (INPUT_SETS - 1)];
            i += 1;
            ns_a.copy_from_slice(ns);
            a(black_box(&mut qs_a), &mut ns_a, ds, *d_inv);
        },
        &mut || {
            let (ns, ds, d_inv) = &inputs[j & (INPUT_SETS - 1)];
            j += 1;
            ns_b.copy_from_slice(ns);
            b(black_box(&mut qs_b), &mut ns_b, ds, *d_inv);
        },
    );
    Some(times)
}

// Newton inversion vs the basecase approximate inversion, at divisor length n (normalized).
fn tune_inv_newton() {
    find_crossover_spec(
        "INV_NEWTON_THRESHOLD",
        "usize",
        "invert_basecase",
        "invert_newton",
        5,
        4000,
        &|n| {
            if n < 5 {
                return None;
            }
            let inputs: Vec<Vec<Limb>> = (0..INPUT_SETS)
                .map(|k| {
                    let mut ds: Vec<Limb> =
                        random_primitive_ints(EXAMPLE_SEED.fork(&format!("inv{k}")))
                            .take(n)
                            .collect();
                    ds[n - 1] |= 1 << (Limb::WIDTH - 1);
                    ds
                })
                .collect();
            let mut is_a = vec![0; n];
            let mut is_b = vec![0; n];
            let mut scratch_a = vec![0; limbs_invert_approx_scratch_len(n)];
            let mut scratch_b = vec![0; limbs_invert_approx_scratch_len(n)];
            for ds in &inputs {
                limbs_invert_basecase_approx(&mut is_a, ds, &mut scratch_a);
                limbs_invert_newton_approx(&mut is_b, ds, &mut scratch_b);
            }
            let (mut i, mut j) = (0usize, 0usize);
            Some(interleaved_min_pair(
                &mut || {
                    let ds = &inputs[i & (INPUT_SETS - 1)];
                    i += 1;
                    black_box(limbs_invert_basecase_approx(
                        black_box(&mut is_a),
                        ds,
                        &mut scratch_a,
                    ));
                },
                &mut || {
                    let ds = &inputs[j & (INPUT_SETS - 1)];
                    j += 1;
                    black_box(limbs_invert_newton_approx(
                        black_box(&mut is_b),
                        ds,
                        &mut scratch_b,
                    ));
                },
            ))
        },
    );
}

// Divide-and-conquer division vs Barrett (MU) division, dividing 2n limbs by n limbs. The DC side
// consumes its dividend, so its per-call refresh copy is included (as any dividend-preserving
// caller would pay it); Barrett reads the dividend by reference.
fn measure_mu_pair(
    n: usize,
    dc: DivAlgoFn,
    barrett: fn(&mut [Limb], &mut [Limb], &[Limb], &[Limb], &mut [Limb]) -> bool,
    barrett_scratch: fn(usize, usize) -> usize,
) -> Option<(f64, f64)> {
    if n < 6 {
        return None;
    }
    let inputs: Vec<(Vec<Limb>, Vec<Limb>, Limb)> = (0..INPUT_SETS)
        .map(|k| {
            let ns: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("mn{k}")))
                .take(n << 1)
                .collect();
            let mut ds: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("md{k}")))
                .take(n)
                .collect();
            ds[n - 1] |= 1 << (Limb::WIDTH - 1);
            let d_inv = limbs_two_limb_inverse_helper(ds[n - 1], ds[n - 2]);
            (ns, ds, d_inv)
        })
        .collect();
    let mut ns_work = vec![0; n << 1];
    let mut qs_a = vec![0; n + 1];
    let mut qs_b = vec![0; n + 1];
    let mut rs = vec![0; n];
    let mut scratch = vec![0; barrett_scratch(n << 1, n)];
    for (ns, ds, d_inv) in &inputs {
        ns_work.copy_from_slice(ns);
        dc(&mut qs_a[..n], &mut ns_work, ds, *d_inv);
        barrett(&mut qs_b[..n], &mut rs, ns, ds, &mut scratch);
    }
    let (mut i, mut j) = (0usize, 0usize);
    Some(interleaved_min_pair(
        &mut || {
            let (ns, ds, d_inv) = &inputs[i & (INPUT_SETS - 1)];
            i += 1;
            ns_work.copy_from_slice(ns);
            dc(black_box(&mut qs_a[..n]), &mut ns_work, ds, *d_inv);
        },
        &mut || {
            let (ns, ds, _) = &inputs[j & (INPUT_SETS - 1)];
            j += 1;
            barrett(black_box(&mut qs_b[..n]), &mut rs, ns, ds, &mut scratch);
        },
    ))
}

fn tune_mu_div_qr() {
    find_crossover_spec(
        "MU_DIV_QR_THRESHOLD",
        "usize",
        "dc_div_qr",
        "barrett_div_qr",
        61,
        8000,
        &|n| {
            measure_mu_pair(
                n,
                limbs_div_mod_divide_and_conquer,
                limbs_div_mod_barrett,
                limbs_div_mod_barrett_scratch_len,
            )
        },
    );
}

fn tune_mu_divappr_q() {
    find_crossover_spec(
        "MU_DIVAPPR_Q_THRESHOLD",
        "usize",
        "dc_divappr_q",
        "barrett_divappr_q",
        61,
        10000,
        &|n| {
            measure_mu_pair(
                n,
                limbs_div_divide_and_conquer_approx,
                |qs, _rs, ns, ds, scratch| {
                    // The approx variant takes no remainder buffer; adapt to the common shape.
                    limbs_div_barrett_approx(qs, ns, ds, scratch)
                },
                limbs_div_barrett_approx_scratch_len,
            )
        },
    );
}

// The kernels for multiplying `NaturalPolynomial`s modulo a word, timed over a grid of lengths and
// modulus sizes: the full product reduced afterwards (as FLINT does for `fmpz_mod_poly`), and
// malachite-base's word kernels, by schoolbook and by Karatsuba multiplication, with the
// conversions to and from limbs included. POLY_GRID_OP chooses the operation (mul, square,
// mul_truncated, or square_truncated), and POLY_GRID_LENS and POLY_GRID_BITS (comma-separated)
// replace the default grid. The modulus has exactly the given number of bits and is odd. A batch is
// calibrated to at least 20 ms, and the best of 5 is kept; `-` means the kernel was skipped as too
// slow there.
fn tune_poly_mod_mul_grid() {
    use malachite_base::num::conversion::traits::ExactFrom;
    use malachite_base::unsigned_polynomial::arithmetic::mod_mul::{
        mod_mul_to_out_classical as word_mul_classical,
        mod_mul_to_out_karatsuba as word_mul_karatsuba,
    };
    use malachite_base::unsigned_polynomial::arithmetic::mod_mul_truncated::{
        mod_mul_truncated_to_out_classical as word_mul_truncated_classical,
        mod_mul_truncated_to_out_karatsuba as word_mul_truncated_karatsuba,
    };
    use malachite_base::unsigned_polynomial::arithmetic::mod_square::{
        mod_square_to_out_classical as word_square_classical,
        mod_square_to_out_karatsuba as word_square_karatsuba,
    };
    use malachite_base::unsigned_polynomial::arithmetic::mod_square_truncated::{
        mod_square_truncated_to_out_classical as word_square_truncated_classical,
        mod_square_truncated_to_out_karatsuba as word_square_truncated_karatsuba,
    };
    use malachite_nz::natural::Natural;
    use malachite_nz::natural_polynomial::arithmetic::mod_mul::*;
    use malachite_nz::natural_polynomial::arithmetic::mod_mul_truncated::*;
    use malachite_nz::natural_polynomial::arithmetic::mod_square::*;
    use malachite_nz::natural_polynomial::arithmetic::mod_square_truncated::*;
    use std::time::Instant;
    fn reduced(seed: &str, n: usize, m: Limb) -> Vec<Natural> {
        random_primitive_ints::<Limb>(EXAMPLE_SEED.fork(seed))
            .take(n)
            .map(|x| Natural::from(x % m))
            .collect()
    }
    fn time_one(f: &mut dyn FnMut()) -> f64 {
        let mut iters = 1u64;
        loop {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            let t = t0.elapsed().as_secs_f64();
            if t >= 0.02 {
                break;
            }
            iters *= if t < 0.002 { 10 } else { 2 };
        }
        let mut best = f64::INFINITY;
        for _ in 0..5 {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            best = best.min(t0.elapsed().as_secs_f64() / iters as f64);
        }
        best
    }
    let grid = |var: &str, default: &[u64]| -> Vec<u64> {
        std::env::var(var).map_or_else(
            |_| default.to_vec(),
            |v| v.split(',').map(|x| x.trim().parse().unwrap()).collect(),
        )
    };
    let op = std::env::var("POLY_GRID_OP").unwrap_or_else(|_| "mul".to_string());
    let lens = grid(
        "POLY_GRID_LENS",
        &[2, 4, 8, 12, 16, 24, 32, 48, 64, 100, 150, 200, 300, 500, 1000, 2000, 4000],
    );
    let all_bits = grid("POLY_GRID_BITS", &[8, 16, 32, 48, 60, 63, 64]);
    println!(
        "{:>6} {:>6} {:>12} {:>12} {:>12}  {:<10}",
        "len", "bits", "full", "classical", "karatsuba", "winner"
    );
    for &n in &lens {
        for &bits in &all_bits {
            if bits > Limb::WIDTH {
                continue;
            }
            let n = usize::try_from(n).unwrap();
            let top: Limb = 1 << (bits - 1);
            let noise = random_primitive_ints::<Limb>(EXAMPLE_SEED.fork(&format!("m{bits}")))
                .next()
                .unwrap();
            let m = top | ((top - 1) & noise) | 1;
            let m_natural = Natural::from(m);
            let xs = reduced(&format!("a{n}_{bits}"), n, m);
            let ys = reduced(&format!("b{n}_{bits}"), n, m);
            let classical_ok = (n as f64) * (n as f64) <= 2.0e9;
            let to_limbs =
                |xs: &[Natural]| -> Vec<Limb> { xs.iter().map(Limb::exact_from).collect() };
            let from_limbs =
                |xs: Vec<Limb>| -> Vec<Natural> { xs.into_iter().map(Natural::from).collect() };
            let word_mul = |classical: bool| -> Vec<Natural> {
                let (a, b) = (to_limbs(&xs), to_limbs(&ys));
                let mut out = vec![0; (n << 1) - 1];
                if classical {
                    word_mul_classical(&mut out, &a, &b, m);
                } else {
                    word_mul_karatsuba(&mut out, &a, &b, m);
                }
                from_limbs(out)
            };
            let word_square = |classical: bool| -> Vec<Natural> {
                let a = to_limbs(&xs);
                let mut out = vec![0; (n << 1) - 1];
                if classical {
                    word_square_classical(&mut out, &a, m);
                } else {
                    word_square_karatsuba(&mut out, &a, m);
                }
                from_limbs(out)
            };
            let word_mul_truncated = |classical: bool| -> Vec<Natural> {
                let (a, b) = (to_limbs(&xs), to_limbs(&ys));
                let mut out = vec![0; n];
                if classical {
                    word_mul_truncated_classical(&mut out, &a, &b, m);
                } else {
                    word_mul_truncated_karatsuba(&mut out, &a, &b, m);
                }
                from_limbs(out)
            };
            let word_square_truncated = |classical: bool| -> Vec<Natural> {
                let a = to_limbs(&xs);
                let mut out = vec![0; n];
                if classical {
                    word_square_truncated_classical(&mut out, &a, m);
                } else {
                    word_square_truncated_karatsuba(&mut out, &a, m);
                }
                from_limbs(out)
            };
            let (full, classical, karatsuba): (Box<dyn Fn()>, Box<dyn Fn()>, Box<dyn Fn()>) =
                match op.as_str() {
                    "mul" => (
                        Box::new(|| drop(black_box(mod_mul_full(&xs, &ys, &m_natural)))),
                        Box::new(|| drop(black_box(word_mul(true)))),
                        Box::new(|| drop(black_box(word_mul(false)))),
                    ),
                    "square" => (
                        Box::new(|| drop(black_box(mod_square_full(&xs, &m_natural)))),
                        Box::new(|| drop(black_box(word_square(true)))),
                        Box::new(|| drop(black_box(word_square(false)))),
                    ),
                    "mul_truncated" => (
                        Box::new(|| {
                            drop(black_box(mod_mul_truncated_full(&xs, &ys, n, &m_natural)));
                        }),
                        Box::new(|| drop(black_box(word_mul_truncated(true)))),
                        Box::new(|| drop(black_box(word_mul_truncated(false)))),
                    ),
                    "square_truncated" => (
                        Box::new(|| {
                            drop(black_box(mod_square_truncated_full(&xs, n, &m_natural)));
                        }),
                        Box::new(|| drop(black_box(word_square_truncated(true)))),
                        Box::new(|| drop(black_box(word_square_truncated(false)))),
                    ),
                    _ => panic!("unknown POLY_GRID_OP {op}"),
                };
            let t_full = time_one(&mut || full());
            let t_classical = if classical_ok {
                time_one(&mut || classical())
            } else {
                -1.0
            };
            let t_karatsuba = time_one(&mut || karatsuba());
            let names = ["full", "classical", "karatsuba"];
            let ts = [t_full, t_classical, t_karatsuba];
            let mut winner = 0;
            for k in 1..3 {
                if ts[k] >= 0.0 && ts[k] < ts[winner] {
                    winner = k;
                }
            }
            let ns = |t: f64| {
                if t < 0.0 {
                    "-".to_string()
                } else {
                    format!("{:.0}", t * 1e9)
                }
            };
            println!(
                "{:>6} {:>6} {:>12} {:>12} {:>12}  {:<10}",
                n,
                bits,
                ns(t_full),
                ns(t_classical),
                ns(t_karatsuba),
                names[winner]
            );
        }
    }
}

// Geometric evaluation of `UnsignedPolynomial<u64>`s modulo a word, timed over a grid of polynomial
// lengths and numbers of points: evaluating at each power of the ratio, several points at a time,
// and Bluestein's trick, with a middle product and with a truncated product. POLY_GRID_LENS,
// POLY_GRID_POINTS, and POLY_GRID_BITS (comma-separated) replace the default grid. The modulus has
// exactly the given number of bits and is odd, and the ratio is a unit. A batch is calibrated to at
// least 20 ms, and the best of 5 is kept.
fn tune_poly_mod_evaluate_geometric_grid() {
    use malachite_base::num::arithmetic::traits::ModInverse;
    use malachite_base::unsigned_polynomial::arithmetic::evaluate::{
        mod_evaluate_geometric_fast, mod_evaluate_geometric_fast_truncated,
        mod_evaluate_geometric_iter,
    };
    use std::time::Instant;
    fn time_one(f: &mut dyn FnMut()) -> f64 {
        let mut iters = 1u64;
        loop {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            let t = t0.elapsed().as_secs_f64();
            if t >= 0.02 {
                break;
            }
            iters *= if t < 0.002 { 10 } else { 2 };
        }
        let mut best = f64::INFINITY;
        for _ in 0..5 {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            best = best.min(t0.elapsed().as_secs_f64() / iters as f64);
        }
        best
    }
    let grid = |var: &str, default: &[u64]| -> Vec<u64> {
        std::env::var(var).map_or_else(
            |_| default.to_vec(),
            |v| v.split(',').map(|x| x.trim().parse().unwrap()).collect(),
        )
    };
    let lens = grid(
        "POLY_GRID_LENS",
        &[8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096],
    );
    let points = grid(
        "POLY_GRID_POINTS",
        &[8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096],
    );
    let all_bits = grid("POLY_GRID_BITS", &[20, 40, 63, 64]);
    println!(
        "{:>6} {:>6} {:>5} {:>12} {:>12} {:>12} {:>8}",
        "len", "points", "bits", "iter", "fast", "truncated", "iter/fast"
    );
    for &bits in &all_bits {
        let top: u64 = 1 << (bits - 1);
        let noise = random_primitive_ints::<u64>(EXAMPLE_SEED.fork(&format!("m{bits}")))
            .next()
            .unwrap();
        let m = top | ((top - 1) & noise) | 1;
        let mut q = random_primitive_ints::<u64>(EXAMPLE_SEED.fork(&format!("q{bits}")))
            .next()
            .unwrap()
            % m;
        let q_inverse = loop {
            if q != 0
                && let Some(q_inverse) = q.mod_inverse(m)
            {
                break q_inverse;
            }
            q = (q + 1) % m;
        };
        for &n in &lens {
            let n = usize::try_from(n).unwrap();
            let coefficients: Vec<u64> =
                random_primitive_ints::<u64>(EXAMPLE_SEED.fork(&format!("c{n}_{bits}")))
                    .take(n)
                    .map(|c| c % m)
                    .collect();
            for &k in &points {
                let k = usize::try_from(k).unwrap();
                let t_iter = time_one(&mut || {
                    drop(black_box(mod_evaluate_geometric_iter(
                        &coefficients,
                        q,
                        k,
                        m,
                    )));
                });
                let t_fast = time_one(&mut || {
                    drop(black_box(mod_evaluate_geometric_fast(
                        &coefficients,
                        q,
                        q_inverse,
                        k,
                        m,
                    )));
                });
                let t_truncated = time_one(&mut || {
                    drop(black_box(mod_evaluate_geometric_fast_truncated(
                        &coefficients,
                        q,
                        q_inverse,
                        k,
                        m,
                    )));
                });
                println!(
                    "{:>6} {:>6} {:>5} {:>12.0} {:>12.0} {:>12.0} {:>8.2}",
                    n,
                    k,
                    bits,
                    t_iter * 1e9,
                    t_fast * 1e9,
                    t_truncated * 1e9,
                    t_iter / t_fast
                );
            }
        }
    }
}

// Middle products of `u64` coefficients modulo a word, timed over a grid of lengths and modulus
// sizes: classical and Karatsuba. The polynomial and the number of outputs both have the given
// length. POLY_GRID_LENS and POLY_GRID_BITS (comma-separated) replace the default grid. A batch is
// calibrated to at least 20 ms, and the best of 5 is kept.
fn tune_poly_mod_mul_middle_grid() {
    use malachite_base::unsigned_polynomial::arithmetic::mod_mul_middle::{
        mod_mul_middle_to_out_classical, mod_mul_middle_to_out_karatsuba,
    };
    use std::time::Instant;
    fn time_one(f: &mut dyn FnMut()) -> f64 {
        let mut iters = 1u64;
        loop {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            let t = t0.elapsed().as_secs_f64();
            if t >= 0.02 {
                break;
            }
            iters *= if t < 0.002 { 10 } else { 2 };
        }
        let mut best = f64::INFINITY;
        for _ in 0..5 {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            best = best.min(t0.elapsed().as_secs_f64() / iters as f64);
        }
        best
    }
    let grid = |var: &str, default: &[u64]| -> Vec<u64> {
        std::env::var(var).map_or_else(
            |_| default.to_vec(),
            |v| v.split(',').map(|x| x.trim().parse().unwrap()).collect(),
        )
    };
    let lens = grid(
        "POLY_GRID_LENS",
        &[8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 512, 1024],
    );
    let all_bits = grid("POLY_GRID_BITS", &[20, 40, 62, 63, 64]);
    println!(
        "{:>6} {:>5} {:>12} {:>12} {:>8}",
        "len", "bits", "classical", "karatsuba", "cl/kar"
    );
    for &bits in &all_bits {
        let top: u64 = 1 << (bits - 1);
        let noise = random_primitive_ints::<u64>(EXAMPLE_SEED.fork(&format!("m{bits}")))
            .next()
            .unwrap();
        let m = top | ((top - 1) & noise) | 1;
        for &n in &lens {
            let n = usize::try_from(n).unwrap();
            let xs: Vec<u64> = random_primitive_ints::<u64>(EXAMPLE_SEED.fork(&format!("x{n}")))
                .take(n)
                .map(|c| c % m)
                .collect();
            let ys: Vec<u64> = random_primitive_ints::<u64>(EXAMPLE_SEED.fork(&format!("y{n}")))
                .take((n << 1) - 1)
                .map(|c| c % m)
                .collect();
            let mut out = vec![0; n];
            let t_classical = time_one(&mut || {
                mod_mul_middle_to_out_classical(black_box(&mut out), &xs, &ys, m);
            });
            let t_karatsuba = time_one(&mut || {
                mod_mul_middle_to_out_karatsuba(black_box(&mut out), &xs, &ys, m);
            });
            println!(
                "{:>6} {:>5} {:>12.0} {:>12.0} {:>8.2}",
                n,
                bits,
                t_classical * 1e9,
                t_karatsuba * 1e9,
                t_classical / t_karatsuba
            );
        }
    }
}

// Times the IntegerPolynomial powering kernels over a grid of lengths, coefficient sizes, and
// exponents, through `pow_ref_with_kernel`, for choosing between them in `pow_to_out`. Coefficients
// are random with exactly `bits` bits and random signs. "default" is the current dispatch; the last
// columns compare binary exponentiation with the multinomial recurrence and with addition chains,
// and a ratio above 1 means binary exponentiation is slower. Grids are overridable through
// POW_GRID_LENS, POW_GRID_BITS, and POW_GRID_EXPS; configurations whose power would have more than
// POW_GRID_MAX_BITS bits in total (default 2^27) are skipped.
fn tune_poly_pow_grid() {
    use malachite_base::num::arithmetic::traits::{ModPowerOf2, Parity};
    use malachite_base::num::logic::traits::BitAccess;
    use malachite_nz::integer::Integer;
    use malachite_nz::integer_polynomial::arithmetic::pow::binexp::pow_to_out_binexp;
    use malachite_nz::integer_polynomial::arithmetic::pow::binomial::pow_to_out_binomial;
    use malachite_nz::integer_polynomial::arithmetic::pow::multinomial::pow_to_out_multinomial;
    use malachite_nz::integer_polynomial::arithmetic::pow::{
        pow_ref_with_kernel, pow_to_out, pow_to_out_addchains_e,
    };
    use malachite_nz::natural::Natural;
    use malachite_nz::test_util::integer_polynomial::arithmetic::pow::*;
    use std::time::Instant;
    fn dense(seed: &str, n: usize, bits: u64) -> Vec<Integer> {
        let mut limbs = random_primitive_ints::<Limb>(EXAMPLE_SEED.fork(seed));
        let mut signs = random_primitive_ints::<u8>(EXAMPLE_SEED.fork(&format!("{seed}s")));
        let limb_count = usize::try_from(bits.div_ceil(Limb::WIDTH)).unwrap();
        (0..n)
            .map(|_| {
                let xs: Vec<Limb> = (&mut limbs).take(limb_count).collect();
                let mut x = Natural::from_owned_limbs_asc(xs).mod_power_of_2(bits);
                x.set_bit(bits - 1);
                if signs.next().unwrap().odd() {
                    -Integer::from(x)
                } else {
                    Integer::from(x)
                }
            })
            .collect()
    }
    // Seconds per call: iterations are doubled (or multiplied by 10) until a batch takes at least
    // 20 ms, and the best of 3 batches is kept.
    fn time_one(f: &mut dyn FnMut()) -> f64 {
        let mut iters = 1u64;
        loop {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            let t = t0.elapsed().as_secs_f64();
            if t >= 0.02 {
                break;
            }
            iters *= if t < 0.002 { 10 } else { 2 };
        }
        let mut best = f64::INFINITY;
        for _ in 0..3 {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            best = best.min(t0.elapsed().as_secs_f64() / iters as f64);
        }
        best
    }
    let grid = |var: &str, default: &[u64]| -> Vec<u64> {
        std::env::var(var).map_or_else(
            |_| default.to_vec(),
            |v| v.split(',').map(|x| x.trim().parse().unwrap()).collect(),
        )
    };
    let lens = grid("POW_GRID_LENS", &[2, 3, 4, 6, 8, 16, 32, 64]);
    let all_bits = grid("POW_GRID_BITS", &[8, 32, 64, 256, 1024, 4096]);
    let exps = grid("POW_GRID_EXPS", &[5, 8, 16, 32, 64, 148]);
    let max_bits: u64 = std::env::var("POW_GRID_MAX_BITS").map_or(1 << 27, |v| v.parse().unwrap());
    println!(
        "{:>4} {:>5} {:>4} {:>12} {:>12} {:>12} {:>12} {:>12} {:>8} {:>8}",
        "len",
        "bits",
        "e",
        "default",
        "binexp",
        "multinom",
        "multi FLINT",
        "addchains",
        "bin/mul",
        "bin/add"
    );
    for &len in &lens {
        for &bits in &all_bits {
            for &e in &exps {
                if e.saturating_mul(e).saturating_mul(len).saturating_mul(bits) > max_bits {
                    continue;
                }
                let n = usize::try_from(len).unwrap();
                let xs = dense(&format!("p{len}_{bits}"), n, bits);
                let t_default = time_one(&mut || {
                    black_box(pow_ref_with_kernel(&xs, e, pow_to_out));
                });
                // For length 2, the binomial kernel stands in the multinomial column.
                let t_binexp = time_one(&mut || {
                    black_box(pow_ref_with_kernel(&xs, e, pow_to_out_binexp));
                });
                let t_multi = if len == 2 {
                    time_one(&mut || {
                        black_box(pow_ref_with_kernel(&xs, e, pow_to_out_binomial));
                    })
                } else {
                    time_one(&mut || {
                        black_box(pow_ref_with_kernel(&xs, e, pow_to_out_multinomial));
                    })
                };
                let t_flint = time_one(&mut || {
                    black_box(pow_ref_with_kernel(&xs, e, pow_to_out_multinomial_flint));
                });
                let t_add = time_one(&mut || {
                    black_box(pow_ref_with_kernel(&xs, e, pow_to_out_addchains_e));
                });
                println!(
                    concat!(
                        "{:>4} {:>5} {:>4} {:>12.0} {:>12.0} {:>12.0} {:>12.0} {:>12.0} ",
                        "{:>8.2} {:>8.2}"
                    ),
                    len,
                    bits,
                    e,
                    t_default * 1e9,
                    t_binexp * 1e9,
                    t_multi * 1e9,
                    t_flint * 1e9,
                    t_add * 1e9,
                    t_binexp / t_multi,
                    t_binexp / t_add
                );
            }
        }
    }
}

// Times the forms of the IntegerPolynomial multinomial powering kernel: the default choice between
// the two below, the precomputed multiples, the sum split by sign, and FLINT's single signed sum.
// Grids are overridable through POW_GRID_LENS, POW_GRID_BITS, and POW_GRID_EXPS.
fn tune_poly_pow_multinomial_grid() {
    use malachite_base::num::arithmetic::traits::{ModPowerOf2, Parity};
    use malachite_base::num::logic::traits::BitAccess;
    use malachite_nz::integer::Integer;
    use malachite_nz::integer_polynomial::arithmetic::pow::multinomial::*;
    use malachite_nz::integer_polynomial::arithmetic::pow::pow_ref_with_kernel;
    use malachite_nz::natural::Natural;
    use malachite_nz::test_util::integer_polynomial::arithmetic::pow::pow_to_out_multinomial_flint;
    use std::time::Instant;
    fn dense(seed: &str, n: usize, bits: u64) -> Vec<Integer> {
        let mut limbs = random_primitive_ints::<Limb>(EXAMPLE_SEED.fork(seed));
        let mut signs = random_primitive_ints::<u8>(EXAMPLE_SEED.fork(&format!("{seed}s")));
        let limb_count = usize::try_from(bits.div_ceil(Limb::WIDTH)).unwrap();
        (0..n)
            .map(|_| {
                let xs: Vec<Limb> = (&mut limbs).take(limb_count).collect();
                let mut x = Natural::from_owned_limbs_asc(xs).mod_power_of_2(bits);
                x.set_bit(bits - 1);
                if signs.next().unwrap().odd() {
                    -Integer::from(x)
                } else {
                    Integer::from(x)
                }
            })
            .collect()
    }
    // Seconds per call: iterations are doubled (or multiplied by 10) until a batch takes at least
    // 20 ms, and the best of 3 batches is kept.
    fn time_one(f: &mut dyn FnMut()) -> f64 {
        let mut iters = 1u64;
        loop {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            let t = t0.elapsed().as_secs_f64();
            if t >= 0.02 {
                break;
            }
            iters *= if t < 0.002 { 10 } else { 2 };
        }
        let mut best = f64::INFINITY;
        for _ in 0..3 {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            best = best.min(t0.elapsed().as_secs_f64() / iters as f64);
        }
        best
    }
    let grid = |var: &str, default: &[u64]| -> Vec<u64> {
        std::env::var(var).map_or_else(
            |_| default.to_vec(),
            |v| v.split(',').map(|x| x.trim().parse().unwrap()).collect(),
        )
    };
    let lens = grid("POW_GRID_LENS", &[3, 4, 8, 16]);
    let all_bits = grid("POW_GRID_BITS", &[8, 32, 64, 256, 1024]);
    let exps = grid("POW_GRID_EXPS", &[16, 64, 148]);
    println!(
        "{:>4} {:>5} {:>4} {:>12} {:>12} {:>12} {:>12}",
        "len", "bits", "e", "default", "multiples", "split", "FLINT"
    );
    for &len in &lens {
        for &bits in &all_bits {
            for &e in &exps {
                let n = usize::try_from(len).unwrap();
                let xs = dense(&format!("p{len}_{bits}"), n, bits);
                let t = [
                    time_one(&mut || {
                        black_box(pow_ref_with_kernel(&xs, e, pow_to_out_multinomial));
                    }),
                    time_one(&mut || {
                        black_box(pow_ref_with_kernel(
                            &xs,
                            e,
                            pow_to_out_multinomial_multiples,
                        ));
                    }),
                    time_one(&mut || {
                        black_box(pow_ref_with_kernel(&xs, e, pow_to_out_multinomial_split));
                    }),
                    time_one(&mut || {
                        black_box(pow_ref_with_kernel(&xs, e, pow_to_out_multinomial_flint));
                    }),
                ];
                println!(
                    "{:>4} {:>5} {:>4} {:>12.0} {:>12.0} {:>12.0} {:>12.0}",
                    len,
                    bits,
                    e,
                    t[0] * 1e9,
                    t[1] * 1e9,
                    t[2] * 1e9,
                    t[3] * 1e9
                );
            }
        }
    }
}

// Times binary exponentiation against addition chains for IntegerPolynomial powering, over
// exponents that are not powers of 2. A ratio above 1 means binary exponentiation is slower.
fn tune_poly_pow_chain_grid() {
    use malachite_base::num::arithmetic::traits::{ModPowerOf2, Parity};
    use malachite_base::num::logic::traits::BitAccess;
    use malachite_nz::integer::Integer;
    use malachite_nz::integer_polynomial::arithmetic::pow::binexp::pow_to_out_binexp;
    use malachite_nz::integer_polynomial::arithmetic::pow::{
        pow_ref_with_kernel, pow_to_out_addchains_e,
    };
    use malachite_nz::natural::Natural;
    use std::time::Instant;
    fn dense(seed: &str, n: usize, bits: u64) -> Vec<Integer> {
        let mut limbs = random_primitive_ints::<Limb>(EXAMPLE_SEED.fork(seed));
        let mut signs = random_primitive_ints::<u8>(EXAMPLE_SEED.fork(&format!("{seed}s")));
        let limb_count = usize::try_from(bits.div_ceil(Limb::WIDTH)).unwrap();
        (0..n)
            .map(|_| {
                let xs: Vec<Limb> = (&mut limbs).take(limb_count).collect();
                let mut x = Natural::from_owned_limbs_asc(xs).mod_power_of_2(bits);
                x.set_bit(bits - 1);
                if signs.next().unwrap().odd() {
                    -Integer::from(x)
                } else {
                    Integer::from(x)
                }
            })
            .collect()
    }
    // Seconds per call: iterations are doubled (or multiplied by 10) until a batch takes at least
    // 20 ms, and the best of 3 batches is kept.
    fn time_one(f: &mut dyn FnMut()) -> f64 {
        let mut iters = 1u64;
        loop {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            let t = t0.elapsed().as_secs_f64();
            if t >= 0.02 {
                break;
            }
            iters *= if t < 0.002 { 10 } else { 2 };
        }
        let mut best = f64::INFINITY;
        for _ in 0..3 {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            best = best.min(t0.elapsed().as_secs_f64() / iters as f64);
        }
        best
    }
    let grid = |var: &str, default: &[u64]| -> Vec<u64> {
        std::env::var(var).map_or_else(
            |_| default.to_vec(),
            |v| v.split(',').map(|x| x.trim().parse().unwrap()).collect(),
        )
    };
    let lens = grid("POW_GRID_LENS", &[3, 8, 32]);
    let all_bits = grid("POW_GRID_BITS", &[64, 1024, 4096]);
    let exps = grid(
        "POW_GRID_EXPS",
        &[5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 23, 27, 31, 47, 63, 100, 127],
    );
    println!(
        "{:>4} {:>5} {:>4} {:>12} {:>12} {:>8}",
        "len", "bits", "e", "binexp", "addchains", "bin/add"
    );
    for &len in &lens {
        for &bits in &all_bits {
            for &e in &exps {
                let n = usize::try_from(len).unwrap();
                let xs = dense(&format!("p{len}_{bits}"), n, bits);
                let t_binexp = time_one(&mut || {
                    black_box(pow_ref_with_kernel(&xs, e, pow_to_out_binexp));
                });
                let t_add = time_one(&mut || {
                    black_box(pow_ref_with_kernel(&xs, e, pow_to_out_addchains_e));
                });
                println!(
                    "{:>4} {:>5} {:>4} {:>12.0} {:>12.0} {:>8.2}",
                    len,
                    bits,
                    e,
                    t_binexp * 1e9,
                    t_add * 1e9,
                    t_binexp / t_add
                );
            }
        }
    }
}

// The Malachite side of the FFT-region mul comparison; the C sides are perf/scratch/{mul_gmp.c,
// mul_flint.c} (make mul-gmp / mul-gmp-noasm / mul-flint). Inputs use the same LCG so all four
// harnesses multiply identical operands. Times go through the full dispatch, so sizes >=
// MUL_FFT_THRESHOLD exercise the fft_small port.
fn tune_mul_fft_probe() {
    let mut n = 1024;
    while n <= 131072 {
        let mut xs = vec![0; n];
        let mut ys = vec![0; n];
        lcg_fill(&mut xs, 1);
        lcg_fill(&mut ys, 2);
        let mut out = vec![0; n << 1];
        let mut scratch = vec![0; limbs_mul_greater_to_out_scratch_len(n, n)];
        limbs_mul_greater_to_out(&mut out, &xs, &ys, &mut scratch); // warmup
        let mut f = || {
            black_box(limbs_mul_greater_to_out(
                black_box(&mut out),
                &xs,
                &ys,
                &mut scratch,
            ));
        };
        let iters = 1 + ((1u64 << 22) / n as u64);
        let mut best = f64::INFINITY;
        for _ in 0..7 {
            let t = time_batch(&mut f, iters);
            if t < best {
                best = t;
            }
        }
        println!("limbs_mul_greater_to_out n={n:<7} {best:>14.1} ns");
        n <<= 1;
    }
}

// The kernels for multiplying `NaturalPolynomial`s modulo $2^k$, timed over a grid of lengths and
// powers: the full product reduced afterwards (as FLINT does), and the products with low halves of
// coefficient products, by schoolbook and by Karatsuba multiplication. POLY_GRID_OP chooses the
// operation (mul, square, mul_truncated, or square_truncated), and POLY_GRID_LENS and
// POLY_GRID_POWS (comma-separated) replace the default grid. A batch is calibrated to at least 20
// ms, and the best of 5 is kept; `-` means the kernel was skipped as too slow there.
fn tune_poly_mod_power_of_2_mul_grid() {
    use malachite_base::num::arithmetic::traits::ModPowerOf2;
    use malachite_base::num::basic::traits::Zero;
    use malachite_base::num::conversion::traits::ExactFrom;
    use malachite_base::unsigned_polynomial::arithmetic::mod_power_of_2_mul::{
        mod_power_of_2_mul_to_out_classical as word_mul_classical,
        mod_power_of_2_mul_to_out_karatsuba as word_mul_karatsuba,
    };
    use malachite_base::unsigned_polynomial::arithmetic::mod_power_of_2_mul_truncated::{
        mod_power_of_2_mul_truncated_to_out_classical as word_mul_truncated_classical,
        mod_power_of_2_mul_truncated_to_out_karatsuba as word_mul_truncated_karatsuba,
    };
    use malachite_base::unsigned_polynomial::arithmetic::mod_power_of_2_square::{
        mod_power_of_2_square_to_out_classical as word_square_classical,
        mod_power_of_2_square_to_out_karatsuba as word_square_karatsuba,
    };
    use malachite_base::unsigned_polynomial::arithmetic::mod_power_of_2_square_truncated::{
        mod_power_of_2_square_truncated_to_out_classical as word_square_truncated_classical,
        mod_power_of_2_square_truncated_to_out_karatsuba as word_square_truncated_karatsuba,
    };
    use malachite_nz::natural::Natural;
    use malachite_nz::natural_polynomial::arithmetic::mod_power_of_2_mul::*;
    use malachite_nz::natural_polynomial::arithmetic::mod_power_of_2_mul_truncated::*;
    use malachite_nz::natural_polynomial::arithmetic::mod_power_of_2_square::*;
    use malachite_nz::natural_polynomial::arithmetic::mod_power_of_2_square_truncated::*;
    use std::time::Instant;
    fn reduced(seed: &str, n: usize, pow: u64) -> Vec<Natural> {
        let mut limbs = random_primitive_ints::<Limb>(EXAMPLE_SEED.fork(seed));
        let limb_count = usize::try_from(pow.div_ceil(Limb::WIDTH)).unwrap();
        (0..n)
            .map(|_| {
                let xs: Vec<Limb> = (&mut limbs).take(limb_count).collect();
                Natural::from_owned_limbs_asc(xs).mod_power_of_2(pow)
            })
            .collect()
    }
    fn time_one(f: &mut dyn FnMut()) -> f64 {
        let mut iters = 1u64;
        loop {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            let t = t0.elapsed().as_secs_f64();
            if t >= 0.02 {
                break;
            }
            iters *= if t < 0.002 { 10 } else { 2 };
        }
        let mut best = f64::INFINITY;
        for _ in 0..5 {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            best = best.min(t0.elapsed().as_secs_f64() / iters as f64);
        }
        best
    }
    let grid = |var: &str, default: &[u64]| -> Vec<u64> {
        std::env::var(var).map_or_else(
            |_| default.to_vec(),
            |v| v.split(',').map(|x| x.trim().parse().unwrap()).collect(),
        )
    };
    let op = std::env::var("POLY_GRID_OP").unwrap_or_else(|_| "mul".to_string());
    let lens = grid(
        "POLY_GRID_LENS",
        &[2, 4, 8, 12, 16, 24, 32, 48, 64, 100, 150, 200, 300, 500, 1000, 2000, 4000],
    );
    let pows = grid(
        "POLY_GRID_POWS",
        &[8, 16, 32, 48, 64, 96, 128, 200, 300, 500, 1000, 2000, 5000, 10000, 20000],
    );
    println!(
        "{:>6} {:>6} {:>12} {:>12} {:>12}  {:<10}",
        "len", "pow", "full", "classical", "karatsuba", "winner"
    );
    for &n in &lens {
        for &pow in &pows {
            if n as f64 * pow as f64 > 4.1e7 {
                continue;
            }
            let n = usize::try_from(n).unwrap();
            let xs = reduced(&format!("a{n}_{pow}"), n, pow);
            let ys = reduced(&format!("b{n}_{pow}"), n, pow);
            // Rough costs, in limb operations, of the quadratic and the Karatsuba kernels, to skip
            // the cells where they would take seconds.
            let k = pow.div_ceil(Limb::WIDTH) as f64;
            let nf = n as f64;
            let classical_ok = nf * nf * k * k <= 2.0e9;
            let karatsuba_ok = nf.powf(1.585) * k.powf(1.585) <= 2.0e9;
            // For words, the low-half kernels are malachite-base's word kernels; the conversions to
            // and from limbs are timed with them, since a dispatcher would pay for them too.
            let word = pow <= Limb::WIDTH;
            let to_limbs =
                |xs: &[Natural]| -> Vec<Limb> { xs.iter().map(Limb::exact_from).collect() };
            let from_limbs =
                |xs: Vec<Limb>| -> Vec<Natural> { xs.into_iter().map(Natural::from).collect() };
            let word_mul = |classical: bool| -> Vec<Natural> {
                let (a, b) = (to_limbs(&xs), to_limbs(&ys));
                let mut out = vec![0; (n << 1) - 1];
                if classical {
                    word_mul_classical(&mut out, &a, &b, pow);
                } else {
                    word_mul_karatsuba(&mut out, &a, &b, pow);
                }
                from_limbs(out)
            };
            let word_square = |classical: bool| -> Vec<Natural> {
                let a = to_limbs(&xs);
                let mut out = vec![0; (n << 1) - 1];
                if classical {
                    word_square_classical(&mut out, &a, pow);
                } else {
                    word_square_karatsuba(&mut out, &a, pow);
                }
                from_limbs(out)
            };
            let word_mul_truncated = |classical: bool| -> Vec<Natural> {
                let (a, b) = (to_limbs(&xs), to_limbs(&ys));
                let mut out = vec![0; n];
                if classical {
                    word_mul_truncated_classical(&mut out, &a, &b, pow);
                } else {
                    word_mul_truncated_karatsuba(&mut out, &a, &b, pow);
                }
                from_limbs(out)
            };
            let word_square_truncated = |classical: bool| -> Vec<Natural> {
                let a = to_limbs(&xs);
                let mut out = vec![0; n];
                if classical {
                    word_square_truncated_classical(&mut out, &a, pow);
                } else {
                    word_square_truncated_karatsuba(&mut out, &a, pow);
                }
                from_limbs(out)
            };
            let (full, classical, karatsuba): (Box<dyn Fn()>, Box<dyn Fn()>, Box<dyn Fn()>) =
                match op.as_str() {
                    "mul" if word => (
                        Box::new(|| drop(black_box(mod_power_of_2_mul_full(&xs, &ys, pow)))),
                        Box::new(|| drop(black_box(word_mul(true)))),
                        Box::new(|| drop(black_box(word_mul(false)))),
                    ),
                    "square" if word => (
                        Box::new(|| drop(black_box(mod_power_of_2_square_full(&xs, pow)))),
                        Box::new(|| drop(black_box(word_square(true)))),
                        Box::new(|| drop(black_box(word_square(false)))),
                    ),
                    "mul_truncated" if word => (
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_mul_truncated_full(
                                &xs, &ys, n, pow,
                            )));
                        }),
                        Box::new(|| drop(black_box(word_mul_truncated(true)))),
                        Box::new(|| drop(black_box(word_mul_truncated(false)))),
                    ),
                    "square_truncated" if word => (
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_square_truncated_full(&xs, n, pow)));
                        }),
                        Box::new(|| drop(black_box(word_square_truncated(true)))),
                        Box::new(|| drop(black_box(word_square_truncated(false)))),
                    ),
                    "mul" => (
                        Box::new(|| drop(black_box(mod_power_of_2_mul_full(&xs, &ys, pow)))),
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_mul_low_classical(&xs, &ys, pow)));
                        }),
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_mul_low_karatsuba(&xs, &ys, pow)));
                        }),
                    ),
                    "square" => (
                        Box::new(|| drop(black_box(mod_power_of_2_square_full(&xs, pow)))),
                        Box::new(|| drop(black_box(mod_power_of_2_square_low_classical(&xs, pow)))),
                        Box::new(|| drop(black_box(mod_power_of_2_square_low_karatsuba(&xs, pow)))),
                    ),
                    "mul_truncated" => (
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_mul_truncated_full(
                                &xs, &ys, n, pow,
                            )));
                        }),
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_mul_truncated_low_classical(
                                &xs, &ys, n, pow,
                            )));
                        }),
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_mul_truncated_low_karatsuba(
                                &xs, &ys, n, pow,
                            )));
                        }),
                    ),
                    "square_truncated" => (
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_square_truncated_full(&xs, n, pow)));
                        }),
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_square_truncated_low_classical(
                                &xs, n, pow,
                            )));
                        }),
                        Box::new(|| {
                            drop(black_box(mod_power_of_2_square_truncated_low_karatsuba(
                                &xs, n, pow,
                            )));
                        }),
                    ),
                    _ => panic!("unknown POLY_GRID_OP {op}"),
                };
            let t_full = time_one(&mut || full());
            let t_classical = if classical_ok {
                time_one(&mut || classical())
            } else {
                -1.0
            };
            let t_karatsuba = if karatsuba_ok {
                time_one(&mut || karatsuba())
            } else {
                -1.0
            };
            let names = ["full", "classical", "karatsuba"];
            let ts = [t_full, t_classical, t_karatsuba];
            let mut winner = 0;
            for k in 1..3 {
                if ts[k] >= 0.0 && ts[k] < ts[winner] {
                    winner = k;
                }
            }
            let ns = |t: f64| {
                if t < 0.0 {
                    "-".to_string()
                } else {
                    format!("{:.0}", t * 1e9)
                }
            };
            println!(
                "{:>6} {:>6} {:>12} {:>12} {:>12}  {:<10}",
                n,
                pow,
                ns(t_full),
                ns(t_classical),
                ns(t_karatsuba),
                names[winner]
            );
            let _ = Natural::ZERO;
        }
    }
}

// The Malachite side of the polynomial-multiplication shootout; the C side is
// perf/scratch/poly_mul_flint.c (make poly-mul-flint). The grid of lengths and coefficient sizes,
// the dense inputs (every coefficient has exactly `bits` bits and a random sign), and the timing (a
// batch calibrated to at least 20 ms, best of 5) mirror the C harness, so the tables compare line
// for line. There is one more column, the small-prime FFT kernel, which FLINT's harness reaches
// only through its dispatcher; `-` means the kernel declined.
fn tune_poly_mul_grid() {
    use malachite_base::num::arithmetic::traits::{ModPowerOf2, Parity};
    use malachite_base::num::basic::traits::Zero;
    use malachite_base::num::logic::traits::BitAccess;
    use malachite_nz::integer::Integer;
    use malachite_nz::integer_polynomial::arithmetic::mul::classical::mul_to_out_classical;
    use malachite_nz::integer_polynomial::arithmetic::mul::karatsuba::mul_to_out_karatsuba;
    use malachite_nz::integer_polynomial::arithmetic::mul::kronecker::mul_to_out_kronecker;
    use malachite_nz::integer_polynomial::arithmetic::mul::mul_greater_to_out;
    use malachite_nz::integer_polynomial::arithmetic::mul::schonhage_strassen::*;
    use malachite_nz::integer_polynomial::arithmetic::mul_middle::classical::*;
    use malachite_nz::integer_polynomial::arithmetic::mul_middle::fft::mul_middle_to_out_fft;
    use malachite_nz::integer_polynomial::arithmetic::mul_middle::kronecker::*;
    use malachite_nz::integer_polynomial::arithmetic::mul_middle::mul_middle_to_out;
    use malachite_nz::integer_polynomial::arithmetic::mul_middle::schonhage_strassen::*;
    use malachite_nz::integer_polynomial::arithmetic::mul_truncated::classical::*;
    use malachite_nz::integer_polynomial::arithmetic::mul_truncated::karatsuba::*;
    use malachite_nz::integer_polynomial::arithmetic::mul_truncated::kronecker::*;
    use malachite_nz::integer_polynomial::arithmetic::mul_truncated::mul_truncated_to_out;
    use malachite_nz::integer_polynomial::arithmetic::mul_truncated::schonhage_strassen::*;
    use malachite_nz::integer_polynomial::arithmetic::square::classical::square_to_out_classical;
    use malachite_nz::integer_polynomial::arithmetic::square::karatsuba::square_to_out_karatsuba;
    use malachite_nz::integer_polynomial::arithmetic::square::kronecker::square_to_out_kronecker;
    use malachite_nz::integer_polynomial::arithmetic::square::schonhage_strassen::*;
    use malachite_nz::integer_polynomial::arithmetic::square::square_to_out;
    use malachite_nz::integer_polynomial::arithmetic::square_truncated::classical::*;
    use malachite_nz::integer_polynomial::arithmetic::square_truncated::karatsuba::*;
    use malachite_nz::integer_polynomial::arithmetic::square_truncated::kronecker::*;
    use malachite_nz::integer_polynomial::arithmetic::square_truncated::schonhage_strassen::*;
    use malachite_nz::integer_polynomial::arithmetic::square_truncated::square_truncated_to_out;
    use malachite_nz::natural::Natural;
    use std::time::Instant;
    fn dense(seed: &str, n: usize, bits: u64) -> Vec<Integer> {
        let mut limbs = random_primitive_ints::<Limb>(EXAMPLE_SEED.fork(seed));
        let mut signs = random_primitive_ints::<u8>(EXAMPLE_SEED.fork(&format!("{seed}s")));
        let limb_count = usize::try_from(bits.div_ceil(Limb::WIDTH)).unwrap();
        (0..n)
            .map(|_| {
                let xs: Vec<Limb> = (&mut limbs).take(limb_count).collect();
                let mut x = Natural::from_owned_limbs_asc(xs).mod_power_of_2(bits);
                x.set_bit(bits - 1);
                if signs.next().unwrap().odd() {
                    -Integer::from(x)
                } else {
                    Integer::from(x)
                }
            })
            .collect()
    }
    // Seconds per call: iterations are doubled (or multiplied by 10) until a batch takes at least
    // 20 ms, and the best of 5 batches is kept.
    fn time_one(f: &mut dyn FnMut()) -> f64 {
        let mut iters = 1u64;
        loop {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            let t = t0.elapsed().as_secs_f64();
            if t >= 0.02 {
                break;
            }
            iters *= if t < 0.002 { 10 } else { 2 };
        }
        let mut best = f64::INFINITY;
        for _ in 0..5 {
            let t0 = Instant::now();
            for _ in 0..iters {
                f();
            }
            best = best.min(t0.elapsed().as_secs_f64() / iters as f64);
        }
        best
    }
    // POLY_GRID_LENS and POLY_GRID_BITS (comma-separated) replace the default grid, for zooming in
    // on a crossover; the C harness has no such override, so the joined comparison needs the
    // defaults.
    let grid = |var: &str, default: &[u64]| -> Vec<u64> {
        std::env::var(var).map_or_else(
            |_| default.to_vec(),
            |v| v.split(',').map(|x| x.trim().parse().unwrap()).collect(),
        )
    };
    // POLY_GRID_OP chooses the operation: mul (the default, and the only one the C harness
    // measures), square, mul_truncated, square_truncated, or mul_middle.
    let op = std::env::var("POLY_GRID_OP").unwrap_or_else(|_| "mul".to_string());
    let lens: Vec<usize> = grid(
        "POLY_GRID_LENS",
        &[8, 16, 32, 48, 64, 75, 100, 200, 500, 1000, 2000, 5000, 10000],
    )
    .into_iter()
    .map(|n| usize::try_from(n).unwrap())
    .collect();
    let bitss = grid(
        "POLY_GRID_BITS",
        &[10, 64, 200, 500, 800, 1500, 2000, 4000, 8000, 20000],
    );
    println!(
        "{:>6} {:>6} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12}  {:<10} {:>8} {:>8}",
        "len",
        "bits",
        "classical",
        "karatsuba",
        "KS",
        "SS",
        "fft",
        "dispatcher",
        "winner",
        "SS/best",
        "disp/best"
    );
    for &n in &lens {
        for &bits in &bitss {
            if n as f64 * bits as f64 > 4.1e7 {
                continue;
            }
            // The shapes: full products of two length-n factors, products truncated to n
            // coefficients, and the middle product of a length-(2n - 1) factor and a length-n one
            // (the coefficients from n - 1 up to 2n - 1), as Newton iterations use it.
            let (xs_len, out_len, nlo) = match op.as_str() {
                "mul" | "square" => (n, (n << 1) - 1, 0),
                "mul_truncated" | "square_truncated" => (n, n, 0),
                "mul_middle" => ((n << 1) - 1, n, n - 1),
                _ => panic!("unknown POLY_GRID_OP {op}"),
            };
            let nhi = nlo + out_len;
            let xs = dense(&format!("a{n}_{bits}"), xs_len, bits);
            let ys = dense(&format!("b{n}_{bits}"), n, bits);
            let mut out = vec![Integer::ZERO; out_len];
            type Kernel<'a> = &'a dyn Fn(&mut [Integer], &[Integer], &[Integer]);
            let (classical, karatsuba, kronecker, ss, dispatcher): (
                Kernel,
                Option<Kernel>,
                Kernel,
                Kernel,
                Kernel,
            ) = match op.as_str() {
                "mul" => (
                    &mul_to_out_classical,
                    Some(&mul_to_out_karatsuba),
                    &mul_to_out_kronecker,
                    &mul_to_out_schonhage_strassen,
                    &mul_greater_to_out,
                ),
                "square" => (
                    &|o, x, _| square_to_out_classical(o, x),
                    Some(&|o, x, _| square_to_out_karatsuba(o, x)),
                    &|o, x, _| square_to_out_kronecker(o, x),
                    &|o, x, _| square_to_out_schonhage_strassen(o, x),
                    &|o, x, _| square_to_out(o, x),
                ),
                "mul_truncated" => (
                    &mul_truncated_to_out_classical,
                    Some(&mul_truncated_to_out_karatsuba),
                    &mul_truncated_to_out_kronecker,
                    &mul_truncated_to_out_schonhage_strassen,
                    &mul_truncated_to_out,
                ),
                "square_truncated" => (
                    &|o, x, _| square_truncated_to_out_classical(o, x),
                    Some(&|o, x, _| square_truncated_to_out_karatsuba(o, x)),
                    &|o, x, _| square_truncated_to_out_kronecker(o, x),
                    &|o, x, _| square_truncated_to_out_schonhage_strassen(o, x),
                    &|o, x, _| square_truncated_to_out(o, x),
                ),
                _ => (
                    &|o, x, y| mul_middle_to_out_classical(o, x, y, nlo, nhi),
                    None,
                    &|o, x, y| mul_middle_to_out_kronecker(o, x, y, nlo, nhi),
                    &|o, x, y| mul_middle_to_out_schonhage_strassen(o, x, y, nlo, nhi),
                    &|o, x, y| mul_middle_to_out(o, x, y, nlo, nhi),
                ),
            };
            let squaring = op.starts_with("square");
            let fft: Kernel = &|o, x, y| {
                mul_middle_to_out_fft(o, x, if squaring { x } else { y }, nlo, nhi);
            };
            let mut run = |f: Kernel| -> f64 {
                time_one(&mut || f(black_box(&mut out), black_box(&xs), black_box(&ys)))
            };
            let tc = if n <= 200 { run(classical) } else { -1.0 };
            let tk = match karatsuba {
                Some(f) if n <= 2000 => run(f),
                _ => -1.0,
            };
            let tks = run(kronecker);
            let tss = run(ss);
            let fft_accepts = {
                let mut scratch = vec![Integer::ZERO; out_len];
                mul_middle_to_out_fft(
                    &mut scratch,
                    &xs,
                    if squaring { &xs } else { &ys },
                    nlo,
                    nhi,
                )
            };
            let tfft = if fft_accepts { run(fft) } else { -1.0 };
            let td = run(dispatcher);
            let names = ["classical", "karatsuba", "KS", "SS", "fft"];
            let ts = [tc, tk, tks, tss, tfft];
            let mut winner = 0;
            let mut best = f64::INFINITY;
            let mut best_non_ss = f64::INFINITY;
            for (k, &t) in ts.iter().enumerate() {
                if t < 0.0 {
                    continue;
                }
                if t < best {
                    best = t;
                    winner = k;
                }
                if k != 3 && t < best_non_ss {
                    best_non_ss = t;
                }
            }
            let ns = |t: f64| {
                if t < 0.0 {
                    "-".to_string()
                } else {
                    format!("{:.0}", t * 1e9)
                }
            };
            println!(
                "{:>6} {:>6} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12}  {:<10} {:>8.2} {:>8.2}",
                n,
                bits,
                ns(tc),
                ns(tk),
                ns(tks),
                ns(tss),
                ns(tfft),
                ns(td),
                names[winner],
                tss / best_non_ss,
                td / best
            );
        }
    }
}

// The Malachite side of the small-kernel gap analysis; the C side is perf/scratch/small_gmp.c (make
// small-gmp / small-gmp-noasm). Sizes, seeds, and the best-of-batches loop mirror the C harness
// exactly, so the tables compare directly.
fn tune_small_kernel_probe() {
    fn best_of(f: &mut dyn FnMut(), iters: u64) -> f64 {
        let mut best = f64::INFINITY;
        for _ in 0..9 {
            let t = time_batch(f, iters);
            if t < best {
                best = t;
            }
        }
        best
    }
    for n in [1usize, 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64] {
        let mut xs = vec![0; n];
        let mut ys = vec![0; n];
        lcg_fill(&mut xs, 1);
        lcg_fill(&mut ys, 2);
        let mut out = vec![0; n << 1];
        let mut scratch = vec![0; limbs_mul_greater_to_out_scratch_len(n, n)];
        let iters = 1 + ((1u64 << 22) / (n * n + 16) as u64);
        let t_dispatch = best_of(
            &mut || {
                black_box(limbs_mul_greater_to_out(
                    black_box(&mut out),
                    &xs,
                    &ys,
                    &mut scratch,
                ));
            },
            iters,
        );
        let t_basecase = best_of(
            &mut || {
                limbs_mul_greater_to_out_basecase(black_box(&mut out), &xs, &ys);
            },
            iters,
        );
        println!("mul n={n:<3} dispatch {t_dispatch:>10.1} ns   basecase {t_basecase:>10.1} ns");
    }
    for n in [2usize, 4, 8, 16, 32, 64] {
        let mut ns = vec![0; n << 1];
        let mut ds = vec![0; n];
        lcg_fill(&mut ns, 3);
        lcg_fill(&mut ds, 4);
        ds[n - 1] |= 1 << (Limb::WIDTH - 1);
        let mut qs = vec![0; n + 1];
        let mut rs = vec![0; n];
        let iters = 1 + ((1u64 << 22) / (n * n + 16) as u64);
        let t = best_of(
            &mut || {
                limbs_div_mod_to_out(black_box(&mut qs), &mut rs, &ns, &ds);
            },
            iters,
        );
        println!("div 2n/n n={n:<3} {t:>10.1} ns");
    }
    for n in [16usize, 64, 256] {
        let mut xs = vec![0; n];
        let mut ys = vec![0; n];
        lcg_fill(&mut xs, 5);
        lcg_fill(&mut ys, 6);
        // The golden-ratio constant for the limb width, 2^`Limb::WIDTH` / phi: the top half of the
        // 64-bit one is the 32-bit one, so taking the leading `Limb::WIDTH` bits gives the right
        // value either way. Truncating instead would clear the top bit on 32-bit limbs, making the
        // multiplier narrower than a limb and the benchmark unrepresentative.
        let z = (0x9E3779B97F4A7C15u64 >> (u64::WIDTH - Limb::WIDTH)) as Limb;
        let iters = 1 + ((1u64 << 22) / n as u64);
        let t = best_of(
            &mut || {
                black_box(limbs_slice_add_mul_limb_same_length_in_place_left(
                    black_box(&mut xs),
                    &ys,
                    z,
                ));
            },
            iters,
        );
        println!("addmul_1 n={n:<3} {t:>10.1} ns");
    }
}

// The Malachite side of the gcd comparison; the C side is perf/scratch/gcd_gmp.c (make gcd-gmp /
// gcd-gmp-noasm). Identical operands via the shared LCG.
fn tune_gcd_probe() {
    use malachite_base::num::arithmetic::traits::Gcd;
    use malachite_nz::natural::Natural;
    for n in [10usize, 25, 50, 100, 200, 400, 800, 1600, 3200, 6400] {
        let mut xs = vec![0; n];
        let mut ys = vec![0; n];
        lcg_fill(&mut xs, 1);
        lcg_fill(&mut ys, 2);
        let x = Natural::from_owned_limbs_asc(xs);
        let y = Natural::from_owned_limbs_asc(ys);
        let mut best = f64::INFINITY;
        let iters = 1 + ((1u64 << 24) / (n * n) as u64);
        for _ in 0..7 {
            let mut f = || {
                black_box((&x).gcd(&y));
            };
            let t = time_batch(&mut f, iters);
            if t < best {
                best = t;
            }
        }
        println!("gcd n={n:<5} {best:>12.1} ns");
    }
}

// Extended-gcd timing at n-limb operand sizes, for the GCDEXT_DC / MATRIX22_STRASSEN sweeps.
fn tune_xgcd_probe() {
    use malachite_base::num::arithmetic::traits::ExtendedGcd;
    use malachite_nz::natural::Natural;
    for n in [25usize, 50, 100, 200, 400, 800, 1600, 3200] {
        let mut xs = vec![0; n];
        let mut ys = vec![0; n];
        lcg_fill(&mut xs, 1);
        lcg_fill(&mut ys, 2);
        let x = Natural::from_owned_limbs_asc(xs);
        let y = Natural::from_owned_limbs_asc(ys);
        let mut best = f64::INFINITY;
        let iters = 1 + ((1u64 << 23) / (n * n) as u64);
        for _ in 0..7 {
            let mut f = || {
                black_box((&x).extended_gcd(&y));
            };
            let t = time_batch(&mut f, iters);
            if t < best {
                best = t;
            }
        }
        println!("xgcd n={n:<5} {best:>12.1} ns");
    }
}

// Matches the `lcg_fill` in the C harnesses, so all sides see identical operands.
fn lcg_fill(p: &mut [Limb], seed: u64) {
    let mut s = 0x9E3779B97F4A7C15u64 ^ seed;
    for x in p.iter_mut() {
        s = s
            .wrapping_mul(6364136223846793005)
            .wrapping_add(1442695040888963407);
        *x = s as Limb;
    }
}

fn fft_mul_algo<'a>() -> Algo<'a> {
    Algo {
        name: "fft",
        valid: &|_| true,
        scratch_len: &|_| 0,
        run: &|out, xs, ys, _| mpn_mul_fft_for_tuning(out, xs, ys),
    }
}

fn tune_mul_fft() {
    find_crossover(&Level {
        threshold_name: "MUL_FFT_THRESHOLD",
        min_size: 400,
        max_size: 8000,
        lower: toom8h_algo(),
        upper: fft_mul_algo(),
    });
}

fn fft_sqr_algo<'a>() -> Algo<'a> {
    Algo {
        name: "fft_sqr",
        valid: &|_| true,
        scratch_len: &|_| 0,
        run: &|out, xs, _, _| mpn_square_fft_for_tuning(out, xs),
    }
}

// SQR_FFT_THRESHOLD is currently derived (SQR_FFT_MODF_THRESHOLD * 10, frozen at 11700 in
// square.rs); this measures where the FFT square actually overtakes toom8, to inform replacing the
// derivation with a measured constant.
fn tune_sqr_fft() {
    find_crossover(&Level {
        threshold_name: "SQR_FFT_THRESHOLD",
        min_size: 400,
        max_size: 16000,
        lower: sqr_toom8_algo(),
        upper: fft_sqr_algo(),
    });
}

// Both FFT crossovers landed below the toom8 thresholds, so toom8h/toom8 may be squeezed out; these
// measure the FFT against the real incumbents at those sizes.
fn tune_mul_fft_vs_toom6h() {
    find_crossover(&Level {
        threshold_name: "MUL_FFT_THRESHOLD",
        min_size: 64,
        max_size: 8000,
        lower: toom6h_algo(),
        upper: fft_mul_algo(),
    });
}

fn tune_sqr_fft_vs_toom6() {
    find_crossover(&Level {
        threshold_name: "SQR_FFT_THRESHOLD",
        min_size: 64,
        max_size: 16000,
        lower: sqr_toom6_algo(),
        upper: fft_sqr_algo(),
    });
}

// Correctness sweep for the FFT at small sizes (the `l <= LG_BLK_SZ` small-transform branch is only
// reachable below ~400 limbs, which production never did while MUL_FFT_THRESHOLD was 1500):
// compares the FFT against the standard dispatch on many random inputs.
fn tune_fft_small_check() {
    let mut mismatches = 0;
    for n in 64..=1024 {
        for k in 0..4u32 {
            let xs: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("fx{n}_{k}")))
                .take(n)
                .collect();
            let ys: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("fy{n}_{k}")))
                .take(n)
                .collect();
            let mut out_ref = vec![0; n << 1];
            let mut scratch = vec![0; limbs_mul_greater_to_out_scratch_len(n, n)];
            let xs2 = xs.clone();
            let ys2 = ys.clone();
            let mul_ok = std::panic::catch_unwind(move || {
                let mut out = vec![0; xs2.len() << 1];
                mpn_mul_fft_for_tuning(&mut out, &xs2, &ys2);
                out
            });
            limbs_mul_greater_to_out(&mut out_ref, &xs, &ys, &mut scratch);
            match mul_ok {
                Err(_) => {
                    println!("MUL PANIC at n={n} set {k}");
                    mismatches += 1;
                }
                Ok(out) if out != out_ref => {
                    println!("MUL MISMATCH at n={n} set {k}");
                    mismatches += 1;
                }
                _ => {}
            }
            let xs2 = xs.clone();
            let sqr_ok = std::panic::catch_unwind(move || {
                let mut out = vec![0; xs2.len() << 1];
                mpn_square_fft_for_tuning(&mut out, &xs2);
                out
            });
            let mut sq_ref = vec![0; n << 1];
            limbs_mul_greater_to_out(&mut sq_ref, &xs, &xs, &mut scratch);
            match sqr_ok {
                Err(_) => {
                    println!("SQR PANIC at n={n} set {k}");
                    mismatches += 1;
                }
                Ok(out) if out != sq_ref => {
                    println!("SQR MISMATCH at n={n} set {k}");
                    mismatches += 1;
                }
                _ => {}
            }
        }
    }
    println!("fft small-size check: {mismatches} failures over sizes 64..=1024 x4 input sets");
}

// An algorithm entry for the unbalanced-multiplication crossovers, parameterized over (xs_len,
// ys_len) rather than a single balanced size.
struct UnbalancedAlgo<'a> {
    name: &'a str,
    valid: &'a dyn Fn(usize, usize) -> bool,
    scratch_len: &'a dyn Fn(usize, usize) -> usize,
    run: &'a dyn Fn(&mut [Limb], &[Limb], &[Limb], &mut [Limb]),
}

// Like measure_mul_pair, but for distinct operand lengths.
fn measure_unbalanced_pair(
    xs_len: usize,
    ys_len: usize,
    a: &UnbalancedAlgo,
    b: &UnbalancedAlgo,
) -> Option<(f64, f64)> {
    if xs_len < ys_len || !(a.valid)(xs_len, ys_len) || !(b.valid)(xs_len, ys_len) {
        return None;
    }
    let inputs: Vec<(Vec<Limb>, Vec<Limb>)> = (0..INPUT_SETS)
        .map(|k| {
            let xs = random_primitive_ints(EXAMPLE_SEED.fork(&format!("ux{k}")))
                .take(xs_len)
                .collect();
            let ys = random_primitive_ints(EXAMPLE_SEED.fork(&format!("uy{k}")))
                .take(ys_len)
                .collect();
            (xs, ys)
        })
        .collect();
    let out_len = xs_len + ys_len;
    let mut out_a = vec![0; out_len];
    let mut out_b = vec![0; out_len];
    let mut scratch_a = vec![0; (a.scratch_len)(xs_len, ys_len)];
    let mut scratch_b = vec![0; (b.scratch_len)(xs_len, ys_len)];
    for (xs, ys) in &inputs {
        (a.run)(&mut out_a, xs, ys, &mut scratch_a);
        (b.run)(&mut out_b, xs, ys, &mut scratch_b);
    }
    assert_eq!(out_a, out_b);
    let (mut i, mut j) = (0usize, 0usize);
    Some(interleaved_min_pair(
        &mut || {
            let (xs, ys) = &inputs[i & (INPUT_SETS - 1)];
            i += 1;
            (a.run)(black_box(&mut out_a), xs, ys, &mut scratch_a);
        },
        &mut || {
            let (xs, ys) = &inputs[j & (INPUT_SETS - 1)];
            j += 1;
            (b.run)(black_box(&mut out_b), xs, ys, &mut scratch_b);
        },
    ))
}

macro_rules! unbalanced_algo {
    ($name: literal, $valid: ident, $scratch: ident, $run: ident) => {
        UnbalancedAlgo {
            name: $name,
            valid: &|x, y| $valid(x, y),
            scratch_len: &|x, y| $scratch(x, y),
            run: &|out, xs, ys, scratch| $run(out, xs, ys, scratch),
        }
    };
}

// Each interior threshold selects between two Toom variants within an aspect-ratio band of the
// unbalanced dispatch (see limbs_mul_greater_to_out); the scan sweeps ys_len with xs_len at a
// representative aspect ratio inside that band.
fn tune_unbalanced_interior(
    threshold_name: &str,
    lower: &UnbalancedAlgo,
    upper: &UnbalancedAlgo,
    aspect: &dyn Fn(usize) -> usize,
    max_size: usize,
) {
    find_crossover_spec(
        threshold_name,
        "usize",
        lower.name,
        upper.name,
        15,
        max_size,
        &|y| measure_unbalanced_pair(aspect(y), y, lower, upper),
    );
}

fn tune_mul_toom32_to_toom43() {
    // dispatch band: 7 * y / 6 <= x < 3 * y / 2; representative aspect 4 / 3
    tune_unbalanced_interior(
        "MUL_TOOM32_TO_TOOM43_THRESHOLD",
        &unbalanced_algo!(
            "toom32",
            limbs_mul_greater_to_out_toom_32_input_sizes_valid,
            limbs_mul_greater_to_out_toom_32_scratch_len,
            limbs_mul_greater_to_out_toom_32
        ),
        &unbalanced_algo!(
            "toom43",
            limbs_mul_greater_to_out_toom_43_input_sizes_valid,
            limbs_mul_greater_to_out_toom_43_scratch_len,
            limbs_mul_greater_to_out_toom_43
        ),
        &|y| y * 4 / 3,
        1200,
    );
}

fn tune_mul_toom32_to_toom53() {
    // dispatch band: 3 * y / 2 <= x < 7 * y / 4; toom32's validity (2x < 3(y + 1)) only overlaps
    // the band's bottom edge, so the comparison runs at x = 3y/2 + 1
    tune_unbalanced_interior(
        "MUL_TOOM32_TO_TOOM53_THRESHOLD",
        &unbalanced_algo!(
            "toom32",
            limbs_mul_greater_to_out_toom_32_input_sizes_valid,
            limbs_mul_greater_to_out_toom_32_scratch_len,
            limbs_mul_greater_to_out_toom_32
        ),
        &unbalanced_algo!(
            "toom53",
            limbs_mul_greater_to_out_toom_53_input_sizes_valid,
            limbs_mul_greater_to_out_toom_53_scratch_len,
            limbs_mul_greater_to_out_toom_53
        ),
        &|y| y * 3 / 2 + 1,
        2000,
    );
}

fn tune_mul_toom42_to_toom53() {
    // dispatch band: 7 * y / 4 <= x < 11 * y / 6; representative aspect 9 / 5
    tune_unbalanced_interior(
        "MUL_TOOM42_TO_TOOM53_THRESHOLD",
        &unbalanced_algo!(
            "toom42",
            limbs_mul_greater_to_out_toom_42_input_sizes_valid,
            limbs_mul_greater_to_out_toom_42_scratch_len,
            limbs_mul_greater_to_out_toom_42
        ),
        &unbalanced_algo!(
            "toom53",
            limbs_mul_greater_to_out_toom_53_input_sizes_valid,
            limbs_mul_greater_to_out_toom_53_scratch_len,
            limbs_mul_greater_to_out_toom_53
        ),
        &|y| y * 9 / 5,
        3000,
    );
}

fn tune_mul_toom42_to_toom63() {
    // dispatch band: 11 * y / 6 <= x < 5 * y / 2; representative aspect 21 / 10
    tune_unbalanced_interior(
        "MUL_TOOM42_TO_TOOM63_THRESHOLD",
        &unbalanced_algo!(
            "toom42",
            limbs_mul_greater_to_out_toom_42_input_sizes_valid,
            limbs_mul_greater_to_out_toom_42_scratch_len,
            limbs_mul_greater_to_out_toom_42
        ),
        &unbalanced_algo!(
            "toom63",
            limbs_mul_greater_to_out_toom_63_input_sizes_valid,
            limbs_mul_greater_to_out_toom_63_scratch_len,
            limbs_mul_greater_to_out_toom_63
        ),
        &|y| y * 21 / 10,
        3000,
    );
}

fn mullo_basecase_algo<'a>() -> Algo<'a> {
    Algo {
        name: "mullo_basecase",
        valid: &|_| true,
        scratch_len: &|_| 0,
        run: &|out, xs, ys, _| limbs_mul_low_same_length_basecase(out, xs, ys),
    }
}

fn mullo_dc_algo<'a>() -> Algo<'a> {
    Algo {
        name: "mullo_dc",
        valid: &|n| n >= 8,
        scratch_len: &|n| limbs_mul_low_same_length_divide_and_conquer_scratch_len(n),
        run: &|out, xs, ys, scratch| {
            limbs_mul_low_same_length_divide_and_conquer(out, xs, ys, scratch);
        },
    }
}

fn tune_mullo_dc() {
    find_crossover(&Level {
        threshold_name: "MULLO_DC_THRESHOLD",
        min_size: 8,
        max_size: 500,
        lower: mullo_basecase_algo(),
        upper: mullo_dc_algo(),
    });
}

// Above MULLO_MUL_N_THRESHOLD, the low product is computed as a full multiplication (whose high
// half is discarded); the harness's output buffer has 2n limbs, so both run in place. (The
// balanced-measure harness does not compare outputs, so the differing high halves are fine.)
fn tune_mullo_mul_n() {
    find_crossover(&Level {
        threshold_name: "MULLO_MUL_N_THRESHOLD",
        min_size: 300,
        max_size: 30000,
        lower: mullo_dc_algo(),
        upper: Algo {
            name: "full_mul",
            valid: &|_| true,
            scratch_len: &|n| limbs_mul_greater_to_out_scratch_len(n, n),
            run: &|out, xs, ys, scratch| {
                black_box(limbs_mul_greater_to_out(out, xs, ys, scratch));
            },
        },
    });
}

// SQRLO_DC_THRESHOLD is capped by SQRLO_DC_THRESHOLD_LIMIT (500); the scan respects the limit.
fn tune_sqrlo_dc() {
    find_crossover(&Level {
        threshold_name: "SQRLO_DC_THRESHOLD",
        min_size: 15,
        max_size: 499,
        lower: Algo {
            name: "sqrlo_basecase",
            valid: &|_| true,
            scratch_len: &|_| 0,
            run: &|out, xs, _, _| limbs_square_low_basecase(out, xs),
        },
        upper: Algo {
            name: "sqrlo_dc",
            valid: &|n| n >= 8,
            scratch_len: &|n| limbs_square_low_scratch_len(n),
            run: &|out, xs, _, scratch| limbs_square_low_divide_and_conquer(out, xs, scratch),
        },
    });
}

// Measure two Hensel (modular/bdiv) division algorithms dividing 2n limbs by an odd n-limb divisor.
// `full_q` selects the Q-only shape, whose quotient has the dividend's full length.
fn measure_bdiv_pair(
    n: usize,
    full_q: bool,
    a: &dyn Fn(&mut [Limb], &mut [Limb], &[Limb], Limb),
    b: &dyn Fn(&mut [Limb], &mut [Limb], &[Limb], Limb),
) -> Option<(f64, f64)> {
    if n < 2 {
        return None;
    }
    let inputs: Vec<(Vec<Limb>, Vec<Limb>, Limb)> = (0..INPUT_SETS)
        .map(|k| {
            let ns: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("bn{k}")))
                .take(n << 1)
                .collect();
            let mut ds: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("bd{k}")))
                .take(n)
                .collect();
            ds[0] |= 1;
            let d_inv = limbs_modular_invert_limb(ds[0]).wrapping_neg();
            (ns, ds, d_inv)
        })
        .collect();
    let qs_len = if full_q { n << 1 } else { n };
    let mut ns_a = vec![0; n << 1];
    let mut ns_b = vec![0; n << 1];
    let mut qs_a = vec![0; qs_len];
    let mut qs_b = vec![0; qs_len];
    for (ns, ds, d_inv) in &inputs {
        ns_a.copy_from_slice(ns);
        a(&mut qs_a, &mut ns_a, ds, *d_inv);
        ns_b.copy_from_slice(ns);
        b(&mut qs_b, &mut ns_b, ds, *d_inv);
    }
    let (mut i, mut j) = (0usize, 0usize);
    Some(interleaved_min_pair(
        &mut || {
            let (ns, ds, d_inv) = &inputs[i & (INPUT_SETS - 1)];
            i += 1;
            ns_a.copy_from_slice(ns);
            a(black_box(&mut qs_a), &mut ns_a, ds, *d_inv);
        },
        &mut || {
            let (ns, ds, d_inv) = &inputs[j & (INPUT_SETS - 1)];
            j += 1;
            ns_b.copy_from_slice(ns);
            b(black_box(&mut qs_b), &mut ns_b, ds, *d_inv);
        },
    ))
}

fn tune_dc_bdiv_qr() {
    find_crossover_spec(
        "DC_BDIV_QR_THRESHOLD",
        "usize",
        "bdiv_qr_schoolbook",
        "bdiv_qr_dc",
        4,
        1500,
        &|n| {
            measure_bdiv_pair(
                n,
                false,
                &|qs, ns, ds, d_inv| {
                    limbs_modular_div_mod_schoolbook(qs, ns, ds, d_inv);
                },
                &|qs, ns, ds, d_inv| {
                    limbs_modular_div_mod_divide_and_conquer(qs, ns, ds, d_inv);
                },
            )
        },
    );
}

fn tune_dc_bdiv_q() {
    find_crossover_spec(
        "DC_BDIV_Q_THRESHOLD",
        "usize",
        "bdiv_q_schoolbook",
        "bdiv_q_dc",
        4,
        1500,
        &|n| {
            measure_bdiv_pair(
                n,
                true,
                &|qs, ns, ds, d_inv| limbs_modular_div_schoolbook(qs, ns, ds, d_inv),
                &|qs, ns, ds, d_inv| limbs_modular_div_divide_and_conquer(qs, ns, ds, d_inv),
            )
        },
    );
}

// Measure DC-vs-Barrett Hensel division at 2n / n (odd divisor). The DC side consumes its dividend
// (refresh copy included); the Barrett side reads it directly and manages its own inverse
// internally.
fn measure_mu_bdiv_pair(
    n: usize,
    full_q: bool,
    dc: &dyn Fn(&mut [Limb], &mut [Limb], &[Limb], Limb),
    barrett: &dyn Fn(&mut [Limb], &mut [Limb], &[Limb], &[Limb], &mut [Limb]),
    barrett_scratch: fn(usize, usize) -> usize,
) -> Option<(f64, f64)> {
    if n < 4 {
        return None;
    }
    let inputs: Vec<(Vec<Limb>, Vec<Limb>, Limb)> = (0..INPUT_SETS)
        .map(|k| {
            let ns: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("mbn{k}")))
                .take(n << 1)
                .collect();
            let mut ds: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("mbd{k}")))
                .take(n)
                .collect();
            ds[0] |= 1;
            let d_inv = limbs_modular_invert_limb(ds[0]).wrapping_neg();
            (ns, ds, d_inv)
        })
        .collect();
    let qs_len = if full_q { n << 1 } else { n };
    let mut ns_work = vec![0; n << 1];
    let mut qs_a = vec![0; qs_len];
    let mut qs_b = vec![0; qs_len];
    let mut rs = vec![0; n];
    let mut scratch = vec![0; barrett_scratch(n << 1, n)];
    for (ns, ds, d_inv) in &inputs {
        ns_work.copy_from_slice(ns);
        dc(&mut qs_a, &mut ns_work, ds, *d_inv);
        barrett(&mut qs_b, &mut rs, ns, ds, &mut scratch);
    }
    let (mut i, mut j) = (0usize, 0usize);
    Some(interleaved_min_pair(
        &mut || {
            let (ns, ds, d_inv) = &inputs[i & (INPUT_SETS - 1)];
            i += 1;
            ns_work.copy_from_slice(ns);
            dc(black_box(&mut qs_a), &mut ns_work, ds, *d_inv);
        },
        &mut || {
            let (ns, ds, _) = &inputs[j & (INPUT_SETS - 1)];
            j += 1;
            barrett(black_box(&mut qs_b), &mut rs, ns, ds, &mut scratch);
        },
    ))
}

fn tune_mu_bdiv_qr() {
    find_crossover_spec(
        "MU_BDIV_QR_THRESHOLD",
        "usize",
        "dc_bdiv_qr",
        "barrett_bdiv_qr",
        50,
        20000,
        &|n| {
            measure_mu_bdiv_pair(
                n,
                false,
                &|qs, ns, ds, d_inv| {
                    limbs_modular_div_mod_divide_and_conquer(qs, ns, ds, d_inv);
                },
                &|qs, rs, ns, ds, scratch| {
                    limbs_modular_div_mod_barrett(qs, rs, ns, ds, scratch);
                },
                limbs_modular_div_mod_barrett_scratch_len,
            )
        },
    );
}

fn tune_mu_bdiv_q() {
    find_crossover_spec(
        "MU_BDIV_Q_THRESHOLD",
        "usize",
        "dc_bdiv_q",
        "barrett_bdiv_q",
        50,
        8000,
        &|n| {
            measure_mu_bdiv_pair(
                n,
                true,
                &|qs, ns, ds, d_inv| limbs_modular_div_divide_and_conquer(qs, ns, ds, d_inv),
                &|qs, _rs, ns, ds, scratch| limbs_modular_div_barrett(qs, ns, ds, scratch),
                limbs_modular_div_barrett_scratch_len,
            )
        },
    );
}

// Times limbs_modular_invert across sizes, for the BINV_NEWTON_THRESHOLD rebuild-per-candidate
// sweep (the threshold is compiled into the Newton recursion).
fn tune_binv_probe() {
    for n in [400usize, 800, 1600, 3200, 6400] {
        let mut ds = vec![0; n];
        lcg_fill(&mut ds, 9);
        ds[0] |= 1;
        let mut is = vec![0; n];
        let mut scratch = vec![0; limbs_modular_invert_scratch_len(n)];
        let mut best = f64::INFINITY;
        let iters = 1 + ((1u64 << 21) / n as u64);
        for _ in 0..7 {
            let mut f = || {
                limbs_modular_invert(black_box(&mut is), &ds, &mut scratch);
            };
            let t = time_batch(&mut f, iters);
            if t < best {
                best = t;
            }
        }
        println!("binv n={n:<5} {best:>12.1} ns");
    }
}

// Basecase vs divide-and-conquer digit parsing, sweeping the digit count (the threshold is in
// digits, not limbs).
fn tune_from_digits_dc() {
    use malachite_nz::natural::conversion::digits::general_digits::{
        limbs_compute_power_table, limbs_digits_power_table_scratch_len_for_tuning,
        limbs_from_digits_small_base_basecase, limbs_from_digits_small_base_divide_and_conquer,
    };
    find_crossover_spec(
        "FROM_DIGITS_DIVIDE_AND_CONQUER_THRESHOLD",
        "usize",
        "parse_basecase",
        "parse_dc",
        500,
        40000,
        &|digit_count| {
            let ds: Vec<u8> = random_primitive_ints(EXAMPLE_SEED.fork("fd"))
                .take(digit_count)
                .map(|b: u8| b % 10)
                .collect();
            let len = digit_count / 19 + 2;
            let mut power_table_memory =
                vec![0; limbs_digits_power_table_scratch_len_for_tuning(len)];
            let (power_len, powers) =
                limbs_compute_power_table(&mut power_table_memory, len, 10, None);
            let mut out_a = vec![0; len + 1];
            let mut out_b = vec![0; len + 1];
            let mut scratch = vec![0; len + 1 + (Limb::WIDTH as usize)];
            let mut best_a = f64::INFINITY;
            let mut best_b = f64::INFINITY;
            let iters = 1 + ((1u64 << 22) / (digit_count * 12) as u64);
            for _ in 0..7 {
                let mut f = || {
                    black_box(limbs_from_digits_small_base_basecase(
                        black_box(&mut out_a),
                        &ds,
                        10,
                    ));
                };
                let t = time_batch(&mut f, iters);
                if t < best_a {
                    best_a = t;
                }
            }
            for _ in 0..7 {
                let mut f = || {
                    black_box(limbs_from_digits_small_base_divide_and_conquer(
                        black_box(&mut out_b),
                        &ds,
                        10,
                        &powers,
                        power_len,
                        &mut scratch,
                    ));
                };
                let t = time_batch(&mut f, iters);
                if t < best_b {
                    best_b = t;
                }
            }
            Some((best_a, best_b))
        },
    );
}

// The measured Barrett-bdiv crossovers vs DC (120/139) fall below DC_BDIV_*_THRESHOLD (218), so DC
// may have no winning range; these measure Barrett against the real incumbent there.
fn tune_mu_bdiv_qr_vs_schoolbook() {
    find_crossover_spec(
        "MU_BDIV_QR_THRESHOLD",
        "usize",
        "bdiv_qr_schoolbook",
        "barrett_bdiv_qr",
        10,
        2000,
        &|n| {
            measure_mu_bdiv_pair(
                n,
                false,
                &|qs, ns, ds, d_inv| {
                    limbs_modular_div_mod_schoolbook(qs, ns, ds, d_inv);
                },
                &|qs, rs, ns, ds, scratch| {
                    limbs_modular_div_mod_barrett(qs, rs, ns, ds, scratch);
                },
                limbs_modular_div_mod_barrett_scratch_len,
            )
        },
    );
}

fn tune_mu_bdiv_q_vs_schoolbook() {
    find_crossover_spec(
        "MU_BDIV_Q_THRESHOLD",
        "usize",
        "bdiv_q_schoolbook",
        "barrett_bdiv_q",
        10,
        2000,
        &|n| {
            measure_mu_bdiv_pair(
                n,
                true,
                &|qs, ns, ds, d_inv| limbs_modular_div_schoolbook(qs, ns, ds, d_inv),
                &|qs, _rs, ns, ds, scratch| limbs_modular_div_barrett(qs, ns, ds, scratch),
                limbs_modular_div_barrett_scratch_len,
            )
        },
    );
}

// A shootout table over the mod-by-limb kernels: times each applicable kernel at each size for
// normalized and unnormalized (< 2^(W - 2)) divisors, from which the MOD_1* thresholds and
// MOD_1_1P_METHOD are read off manually (the dispatch conditions are multi-way).
fn tune_mod_1_shootout() {
    use malachite_nz::natural::arithmetic::mod_op::{
        limbs_mod_limb_any_leading_zeros_1, limbs_mod_limb_any_leading_zeros_2,
        limbs_mod_limb_at_least_1_leading_zero, limbs_mod_limb_at_least_2_leading_zeros,
        limbs_mod_limb_small_normalized, limbs_mod_limb_small_unnormalized,
    };
    type ModFn = fn(&[Limb], Limb) -> Limb;
    let norm_kernels: [(&str, ModFn); 3] = [
        ("small_norm", |ns, d| {
            limbs_mod_limb_small_normalized::<DoubleLimb, Limb>(ns, d)
        }),
        ("any_lz_1", |ns, d| {
            limbs_mod_limb_any_leading_zeros_1::<DoubleLimb, Limb>(ns, d)
        }),
        ("any_lz_2", |ns, d| {
            limbs_mod_limb_any_leading_zeros_2::<DoubleLimb, Limb>(ns, d)
        }),
    ];
    let unnorm_kernels: [(&str, ModFn); 5] = [
        ("small_unnorm", |ns, d| {
            limbs_mod_limb_small_unnormalized::<DoubleLimb, Limb>(ns, d)
        }),
        ("any_lz_1", |ns, d| {
            limbs_mod_limb_any_leading_zeros_1::<DoubleLimb, Limb>(ns, d)
        }),
        ("any_lz_2", |ns, d| {
            limbs_mod_limb_any_leading_zeros_2::<DoubleLimb, Limb>(ns, d)
        }),
        ("one_lz", |ns, d| {
            limbs_mod_limb_at_least_1_leading_zero::<DoubleLimb, Limb>(ns, d)
        }),
        ("two_lz", |ns, d| {
            limbs_mod_limb_at_least_2_leading_zeros::<DoubleLimb, Limb>(ns, d)
        }),
    ];
    for (label, d_mask, kernels) in [
        ("normalized", !(Limb::MAX >> 1), &norm_kernels[..]),
        ("unnormalized (2 leading zeros)", 0, &unnorm_kernels[..]),
    ] {
        println!("### {label} divisors");
        for n in [2usize, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 128] {
            let mut line = format!("n={n:<4}");
            for (name, f) in kernels {
                let inputs: Vec<(Vec<Limb>, Limb)> = (0..INPUT_SETS)
                    .map(|k| {
                        let ns: Vec<Limb> =
                            random_primitive_ints(EXAMPLE_SEED.fork(&format!("m1n{k}")))
                                .take(n)
                                .collect();
                        let mut d: Limb =
                            random_primitive_ints(EXAMPLE_SEED.fork(&format!("m1d{k}")))
                                .next()
                                .unwrap();
                        d >>= 2;
                        d |= d_mask | 1;
                        (ns, d)
                    })
                    .collect();
                let mut best = f64::INFINITY;
                let iters = 1 + ((1u64 << 20) / n as u64);
                for _ in 0..7 {
                    let mut i = 0;
                    let mut g = || {
                        let (ns, d) = &inputs[i & (INPUT_SETS - 1)];
                        i += 1;
                        black_box(f(ns, *d));
                    };
                    let t = time_batch(&mut g, iters);
                    if t < best {
                        best = t;
                    }
                }
                line.push_str(&format!("  {name} {best:>7.1}"));
            }
            println!("{line}");
        }
    }
}

// Times limbs_invert_approx across sizes, for the INV_MULMOD_BNM1_THRESHOLD rebuild-per-candidate
// sweep.
fn tune_invert_probe() {
    for n in [50usize, 100, 200, 400, 800, 1600] {
        let mut ds = vec![0; n];
        lcg_fill(&mut ds, 11);
        ds[n - 1] |= 1 << (Limb::WIDTH - 1);
        let mut is = vec![0; n];
        let mut scratch = vec![0; limbs_invert_approx_scratch_len(n)];
        let mut best = f64::INFINITY;
        let iters = 1 + ((1u64 << 20) / n as u64);
        for _ in 0..7 {
            let mut f = || {
                black_box(limbs_invert_approx(black_box(&mut is), &ds, &mut scratch));
            };
            let t = time_batch(&mut f, iters);
            if t < best {
                best = t;
            }
        }
        println!("invert n={n:<5} {best:>12.1} ns");
    }
}

fn tune_dc_div_qr() {
    find_crossover_spec(
        "DC_DIV_QR_THRESHOLD",
        "usize",
        "schoolbook",
        "divide_and_conquer",
        6,
        500,
        &|n| {
            measure_div_pair(
                n,
                6,
                limbs_div_mod_schoolbook,
                limbs_div_mod_divide_and_conquer,
            )
        },
    );
}

fn tune_dc_divappr_q() {
    find_crossover_spec(
        "DC_DIVAPPR_Q_THRESHOLD",
        "usize",
        "schoolbook_approx",
        "divide_and_conquer_approx",
        6,
        1000,
        &|n| {
            measure_div_pair(
                n,
                6,
                limbs_div_schoolbook_approx,
                limbs_div_divide_and_conquer_approx,
            )
        },
    );
}

// toom4 appears to have no winning range on this machine (toom6 overtakes it below the toom3/toom4
// crossover), so the effective ladder is toom3 -> toom6; this measures that crossover directly.
fn tune_sqr_toom6_vs_toom3() {
    find_crossover(&Level {
        threshold_name: "SQR_TOOM6_THRESHOLD",
        min_size: 60,
        max_size: 6000,
        lower: sqr_toom3_algo(),
        upper: sqr_toom6_algo(),
    });
}

fn tune_sqr_toom8() {
    find_crossover(&Level {
        threshold_name: "SQR_TOOM8_THRESHOLD",
        min_size: 400,
        max_size: 10000,
        lower: sqr_toom6_algo(),
        upper: sqr_toom8_algo(),
    });
}

// ---------------------------------------------------------------------------------------------
// Experimental carry-propagation kernels for limbs_add_same_length_to_out (mpn_add_n analog). These
// exist to compare codegen idioms; the winner gets promoted into natural/arithmetic/add.rs. All are
// #[inline(never)] so they're separately measurable and visible to `cargo asm`.

use malachite_nz::platform::DoubleLimb;

// Variant A: the current library idiom (wrapping_add + comparisons).
#[inline(never)]
fn add_n_current(out: &mut [Limb], xs: &[Limb], ys: &[Limb]) -> bool {
    let mut carry = 0;
    for (out, (&x, &y)) in out.iter_mut().zip(xs.iter().zip(ys.iter())) {
        let result_no_carry = x.wrapping_add(y);
        let result = result_no_carry.wrapping_add(carry);
        carry = Limb::from((result_no_carry < x) || (result < result_no_carry));
        *out = result;
    }
    carry != 0
}

// Variant B: DoubleLimb (u128) accumulator.
#[inline(never)]
fn add_n_double_limb(out: &mut [Limb], xs: &[Limb], ys: &[Limb]) -> bool {
    let mut carry = 0;
    for (out, (&x, &y)) in out.iter_mut().zip(xs.iter().zip(ys.iter())) {
        let sum = DoubleLimb::from(x) + DoubleLimb::from(y) + DoubleLimb::from(carry);
        *out = sum as Limb;
        carry = (sum >> Limb::WIDTH) as Limb;
    }
    carry != 0
}

// Variant C: overflowing_add pair (LLVM's uaddo idiom).
#[inline(never)]
fn add_n_overflowing(out: &mut [Limb], xs: &[Limb], ys: &[Limb]) -> bool {
    let mut carry = false;
    for (out, (&x, &y)) in out.iter_mut().zip(xs.iter().zip(ys.iter())) {
        let (sum, c1) = x.overflowing_add(y);
        let (sum, c2) = sum.overflowing_add(Limb::from(carry));
        carry = c1 | c2;
        *out = sum;
    }
    carry
}

// Variant D: overflowing_add pair, 4x unrolled via as_chunks.
#[inline(never)]
fn add_n_overflowing_x4(out: &mut [Limb], xs: &[Limb], ys: &[Limb]) -> bool {
    let mut carry = false;
    let (out_blocks, out_rem) = out.as_chunks_mut::<4>();
    let (xs_blocks, xs_rem) = xs.as_chunks::<4>();
    let (ys_blocks, ys_rem) = ys.as_chunks::<4>();
    for ((o, x), y) in out_blocks.iter_mut().zip(xs_blocks).zip(ys_blocks) {
        for i in 0..4 {
            let (sum, c1) = x[i].overflowing_add(y[i]);
            let (sum, c2) = sum.overflowing_add(Limb::from(carry));
            carry = c1 | c2;
            o[i] = sum;
        }
    }
    for ((o, &x), &y) in out_rem.iter_mut().zip(xs_rem.iter()).zip(ys_rem.iter()) {
        let (sum, c1) = x.overflowing_add(y);
        let (sum, c2) = sum.overflowing_add(Limb::from(carry));
        carry = c1 | c2;
        *o = sum;
    }
    carry
}

// Variant E: DoubleLimb accumulator, 4x unrolled.
#[inline(never)]
fn add_n_double_limb_x4(out: &mut [Limb], xs: &[Limb], ys: &[Limb]) -> bool {
    let mut carry = 0;
    let (out_blocks, out_rem) = out.as_chunks_mut::<4>();
    let (xs_blocks, xs_rem) = xs.as_chunks::<4>();
    let (ys_blocks, ys_rem) = ys.as_chunks::<4>();
    for ((o, x), y) in out_blocks.iter_mut().zip(xs_blocks).zip(ys_blocks) {
        for i in 0..4 {
            let sum = DoubleLimb::from(x[i]) + DoubleLimb::from(y[i]) + DoubleLimb::from(carry);
            o[i] = sum as Limb;
            carry = (sum >> Limb::WIDTH) as Limb;
        }
    }
    for ((o, &x), &y) in out_rem.iter_mut().zip(xs_rem.iter()).zip(ys_rem.iter()) {
        let sum = DoubleLimb::from(x) + DoubleLimb::from(y) + DoubleLimb::from(carry);
        *o = sum as Limb;
        carry = (sum >> Limb::WIDTH) as Limb;
    }
    carry != 0
}

// ---------------------------------------------------------------------------------------------
// Experimental shift kernels for limbs_shl_to_out (mpn_lshift analog).

// Variant A: the current library idiom — remaining_bits is loop-carried, serializing iterations.
#[inline(never)]
fn shl_to_out_current(out: &mut [Limb], xs: &[Limb], bits: u64) -> Limb {
    let cobits = Limb::WIDTH - bits;
    let mut remaining_bits = 0;
    for (out, x) in out[..xs.len()].iter_mut().zip(xs.iter()) {
        *out = (x << bits) | remaining_bits;
        remaining_bits = x >> cobits;
    }
    remaining_bits
}

// Variant B: windows form — each output limb depends only on two input limbs, so iterations are
// independent and LLVM is free to unroll/vectorize.
#[inline(never)]
fn shl_to_out_windows(out: &mut [Limb], xs: &[Limb], bits: u64) -> Limb {
    let len = xs.len();
    let cobits = Limb::WIDTH - bits;
    out[0] = xs[0] << bits;
    for (o, w) in out[1..len].iter_mut().zip(xs.windows(2)) {
        *o = (w[1] << bits) | (w[0] >> cobits);
    }
    xs[len - 1] >> cobits
}

// Variant C: windows form, manually 4x unrolled.
#[inline(never)]
fn shl_to_out_windows_x4(out: &mut [Limb], xs: &[Limb], bits: u64) -> Limb {
    let len = xs.len();
    let cobits = Limb::WIDTH - bits;
    out[0] = xs[0] << bits;
    let (o_blocks, o_rem) = out[1..len].as_chunks_mut::<4>();
    let mut i = 0;
    for o in o_blocks {
        for j in 0..4 {
            o[j] = (xs[i + j + 1] << bits) | (xs[i + j] >> cobits);
        }
        i += 4;
    }
    for o in o_rem {
        *o = (xs[i + 1] << bits) | (xs[i] >> cobits);
        i += 1;
    }
    xs[len - 1] >> cobits
}

fn tune_shl() {
    type ShlFn = fn(&mut [Limb], &[Limb], u64) -> Limb;
    let variants: [(&str, ShlFn); 3] = [
        ("current", shl_to_out_current),
        ("windows", shl_to_out_windows),
        ("windows_x4", shl_to_out_windows_x4),
    ];
    // Correctness cross-check before timing.
    for n in [1, 2, 5, 17, 100] {
        for bits in [1, 7, 31, 63] {
            let xs: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork("cs"))
                .take(n)
                .collect();
            let mut reference = vec![0; n];
            let ref_carry = variants[0].1(&mut reference, &xs, bits);
            for (name, f) in &variants[1..] {
                let mut out = vec![0; n];
                let carry = f(&mut out, &xs, bits);
                assert_eq!(
                    (out, carry),
                    (reference.clone(), ref_carry),
                    "variant {name} disagrees at n={n}, bits={bits}"
                );
            }
        }
    }
    println!("all variants agree; timing (ns/call, ns/limb):");
    for n in [16usize, 64, 256, 1024, 4096] {
        println!("n = {n}:");
        for (name, f) in &variants {
            let inputs: Vec<Vec<Limb>> = (0..INPUT_SETS)
                .map(|k| {
                    random_primitive_ints(EXAMPLE_SEED.fork(&format!("s{k}")))
                        .take(n)
                        .collect()
                })
                .collect();
            let mut out_a = vec![0; n];
            let mut out_b = vec![0; n];
            for xs in &inputs {
                shl_to_out_current(&mut out_a, xs, 13);
                f(&mut out_b, xs, 13);
            }
            let (mut i, mut j) = (0usize, 0usize);
            let (t_base, t) = interleaved_min_pair(
                &mut || {
                    let xs = &inputs[i & (INPUT_SETS - 1)];
                    i += 1;
                    black_box(shl_to_out_current(black_box(&mut out_a), xs, 13));
                },
                &mut || {
                    let xs = &inputs[j & (INPUT_SETS - 1)];
                    j += 1;
                    black_box(f(black_box(&mut out_b), xs, 13));
                },
            );
            println!(
                "  {name:>12}: {t:>9.1} ns  {:>6.3} ns/limb  (vs current {:>5.2}x)",
                t / n as f64,
                t_base / t,
            );
        }
    }
}

// ---------------------------------------------------------------------------------------------
// div_mod_by_preinversion correction-step shootout. GMP's asm divrem_1 keeps the rare second
// quotient correction off the critical path as a cold branch (~13 insns/limb); LLVM if-converts
// Malachite's version into a long branchless csel chain (~22 insns/limb) on a loop that is
// inherently serial. These variants test whether restructuring recovers the difference.

use malachite_base::num::conversion::traits::{JoinHalves, SplitInHalf};
use malachite_nz::natural::arithmetic::div_mod::{div_mod_by_preinversion, limbs_invert_limb};

// MP_BASES_BIG_BASE_10, the largest power of 10 fitting in a `Limb` (private to the library;
// redeclared here for the shootout): 10^19 for 64-bit limbs, 10^9 for 32-bit limbs.
#[cfg(not(feature = "32_bit_limbs"))]
const BIG_BASE_10: Limb = 0x8ac7230489e80000;
#[cfg(feature = "32_bit_limbs")]
const BIG_BASE_10: Limb = 0x3b9aca00;

// Variant B: GMP-shaped straight-line corrections (first adjustment unconditional on r > q_low,
// second as plain if), letting LLVM choose the lowering.
#[inline]
fn div_mod_preinv_gmp_shape(n_high: Limb, n_low: Limb, d: Limb, d_inv: Limb) -> (Limb, Limb) {
    let (mut q_high, q_low) = (DoubleLimb::from(n_high) * DoubleLimb::from(d_inv))
        .wrapping_add(DoubleLimb::join_halves(n_high.wrapping_add(1), n_low))
        .split_in_half();
    let mut r = n_low.wrapping_sub(q_high.wrapping_mul(d));
    if r > q_low {
        q_high = q_high.wrapping_sub(1);
        r = r.wrapping_add(d);
    }
    if r >= d {
        q_high = q_high.wrapping_add(1);
        r -= d;
    }
    (q_high, r)
}

#[cold]
#[inline(never)]
const fn divrem_second_fixup(q_high: Limb, r: Limb, d: Limb) -> (Limb, Limb) {
    (q_high.wrapping_add(1), r - d)
}

// Variant C: like B, but the rare second correction is outlined into a cold function, forcing a
// real branch and keeping the hot dependency chain short.
#[inline]
fn div_mod_preinv_cold_fixup(n_high: Limb, n_low: Limb, d: Limb, d_inv: Limb) -> (Limb, Limb) {
    let (mut q_high, q_low) = (DoubleLimb::from(n_high) * DoubleLimb::from(d_inv))
        .wrapping_add(DoubleLimb::join_halves(n_high.wrapping_add(1), n_low))
        .split_in_half();
    let mut r = n_low.wrapping_sub(q_high.wrapping_mul(d));
    if r > q_low {
        q_high = q_high.wrapping_sub(1);
        r = r.wrapping_add(d);
    }
    if r >= d {
        return divrem_second_fixup(q_high, r, d);
    }
    (q_high, r)
}

// Variant D: current's nested first adjustment (which beat the GMP shape) plus the cold-outlined
// second correction.
#[inline]
fn div_mod_preinv_hybrid(n_high: Limb, n_low: Limb, d: Limb, d_inv: Limb) -> (Limb, Limb) {
    let (mut q_high, q_low) = (DoubleLimb::from(n_high) * DoubleLimb::from(d_inv))
        .wrapping_add(DoubleLimb::join_halves(n_high.wrapping_add(1), n_low))
        .split_in_half();
    let mut r = n_low.wrapping_sub(q_high.wrapping_mul(d));
    if r > q_low {
        let (r_plus_d, overflow) = r.overflowing_add(d);
        if overflow {
            q_high = q_high.wrapping_sub(1);
            r = r_plus_d;
        }
    } else if r >= d {
        return divrem_second_fixup(q_high, r, d);
    }
    (q_high, r)
}

#[inline(never)]
fn divrem_loop_hybrid(qs: &mut [Limb], ns: &[Limb], d: Limb, d_inv: Limb) -> Limb {
    let mut r = 0;
    for (q, &n) in qs.iter_mut().zip(ns.iter()).rev() {
        (*q, r) = div_mod_preinv_hybrid(r, n, d, d_inv);
    }
    r
}

#[inline(never)]
fn divrem_loop_current(qs: &mut [Limb], ns: &[Limb], d: Limb, d_inv: Limb) -> Limb {
    let mut r = 0;
    for (q, &n) in qs.iter_mut().zip(ns.iter()).rev() {
        (*q, r) = div_mod_by_preinversion(r, n, d, d_inv);
    }
    r
}

#[inline(never)]
fn divrem_loop_gmp_shape(qs: &mut [Limb], ns: &[Limb], d: Limb, d_inv: Limb) -> Limb {
    let mut r = 0;
    for (q, &n) in qs.iter_mut().zip(ns.iter()).rev() {
        (*q, r) = div_mod_preinv_gmp_shape(r, n, d, d_inv);
    }
    r
}

#[inline(never)]
fn divrem_loop_cold_fixup(qs: &mut [Limb], ns: &[Limb], d: Limb, d_inv: Limb) -> Limb {
    let mut r = 0;
    for (q, &n) in qs.iter_mut().zip(ns.iter()).rev() {
        (*q, r) = div_mod_preinv_cold_fixup(r, n, d, d_inv);
    }
    r
}

fn tune_divrem() {
    type DivremFn = fn(&mut [Limb], &[Limb], Limb, Limb) -> Limb;
    let variants: [(&str, DivremFn); 4] = [
        ("current", divrem_loop_current),
        ("gmp_shape", divrem_loop_gmp_shape),
        ("cold_fixup", divrem_loop_cold_fixup),
        ("hybrid", divrem_loop_hybrid),
    ];
    let d = BIG_BASE_10;
    let d_inv = limbs_invert_limb::<DoubleLimb, Limb>(d);
    // Correctness cross-check on many random inputs before timing.
    for k in 0..200 {
        let ns: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(&format!("dr{k}")))
            .take(17)
            .collect();
        let mut q_ref = vec![0; 17];
        let r_ref = variants[0].1(&mut q_ref, &ns, d, d_inv);
        for (name, f) in &variants[1..] {
            let mut q = vec![0; 17];
            let r = f(&mut q, &ns, d, d_inv);
            assert_eq!(
                (q, r),
                (q_ref.clone(), r_ref),
                "variant {name} disagrees, seed {k}"
            );
        }
    }
    println!("all variants agree; timing (ns/call, ns/limb):");
    for n in [16usize, 64, 256, 1024] {
        println!("n = {n}:");
        for (name, f) in &variants {
            let inputs: Vec<Vec<Limb>> = (0..INPUT_SETS)
                .map(|k| {
                    random_primitive_ints(EXAMPLE_SEED.fork(&format!("dv{k}")))
                        .take(n)
                        .collect()
                })
                .collect();
            let mut qs_a = vec![0; n];
            let mut qs_b = vec![0; n];
            for xs in &inputs {
                divrem_loop_current(&mut qs_a, xs, d, d_inv);
                f(&mut qs_b, xs, d, d_inv);
            }
            let (mut i, mut j) = (0usize, 0usize);
            let (t_base, t) = interleaved_min_pair(
                &mut || {
                    let xs = &inputs[i & (INPUT_SETS - 1)];
                    i += 1;
                    black_box(divrem_loop_current(black_box(&mut qs_a), xs, d, d_inv));
                },
                &mut || {
                    let xs = &inputs[j & (INPUT_SETS - 1)];
                    j += 1;
                    black_box(f(black_box(&mut qs_b), xs, d, d_inv));
                },
            );
            println!(
                "  {name:>12}: {t:>9.1} ns  {:>6.3} ns/limb  (vs current {:>5.2}x)",
                t / n as f64,
                t_base / t,
            );
        }
    }
}

// ---------------------------------------------------------------------------------------------
// GET_STR_PRECOMPUTE_THRESHOLD: at what size does building a power table + divide-and-conquer beat
// the O(n^2) basecase? The power-table construction is timed inside the candidate, since the real
// dispatch pays it on every conversion.

fn tune_get_str_precompute() {
    use malachite_nz::natural::conversion::digits::general_digits::{
        digits_in_base_per_limb_for_tuning, get_chars_per_limb, limbs_compute_power_table,
        limbs_digits_power_table_scratch_len_for_tuning, limbs_to_digits_small_base_basecase,
        limbs_to_digits_small_base_divide_and_conquer_for_tuning,
        limbs_to_digits_small_base_divide_and_conquer_scratch_len_for_tuning,
    };
    const BASE: u64 = 10;
    // The basecase asserts xs_len < GET_STR_PRECOMPUTE_THRESHOLD (its stack buffers are sized by
    // it), so the scan is capped just below the compiled-in value; lower the constant and rebuild
    // to scan higher.
    let max_size = if TUNE_PROGRAM_BUILD {
        // GET_STR_THRESHOLD_LIMIT: the basecase's lifted buffer bound under TUNE_PROGRAM_BUILD.
        150
    } else {
        GET_STR_PRECOMPUTE_THRESHOLD - 1
    };
    find_crossover(&Level {
        threshold_name: "GET_STR_PRECOMPUTE_THRESHOLD",
        min_size: 4,
        max_size,
        lower: Algo {
            name: "basecase",
            valid: &|_| true,
            scratch_len: &|_| 0,
            run: &|out_limbs, xs, _, _| {
                // out is digit bytes; reuse the limb out buffer as raw space via a local. The
                // basecase writes u8 digits; we keep a thread-local-free local buffer per call
                // shape by transmuting sizes — simplest is a fixed buffer.
                let mut digits = [0u8; 64 * 20];
                let mut xs_copy = [0; 64];
                let n = xs.len();
                xs_copy[..n].copy_from_slice(xs);
                black_box(limbs_to_digits_small_base_basecase(
                    &mut digits[..n * 20],
                    0,
                    &xs_copy[..n],
                    BASE,
                ));
                // Touch out_limbs so the Algo signature stays uniform.
                black_box(&out_limbs[0]);
            },
        },
        upper: Algo {
            name: "powtab+dc",
            valid: &|_| true,
            scratch_len: &|_| 0,
            run: &|out_limbs, xs, _, _| {
                let n = xs.len();
                let mut digits = [0u8; 64 * 20];
                let mut xs_copy = [0; 64];
                xs_copy[..n].copy_from_slice(xs);
                let mut power_table_memory =
                    vec![0; limbs_digits_power_table_scratch_len_for_tuning(n)];
                let digits_len = digits_in_base_per_limb_for_tuning(n, BASE);
                let len = 1 + usize::try_from(digits_len).unwrap() / get_chars_per_limb(BASE);
                let (power_len, powers) =
                    limbs_compute_power_table(&mut power_table_memory, len, BASE, None);
                let mut scratch = vec![
                        0;
                        limbs_to_digits_small_base_divide_and_conquer_scratch_len_for_tuning(n)
                    ];
                black_box(limbs_to_digits_small_base_divide_and_conquer_for_tuning(
                    &mut digits[..n * 20],
                    &mut xs_copy[..n],
                    BASE,
                    &powers,
                    power_len,
                    &mut scratch,
                ));
                black_box(&out_limbs[0]);
            },
        },
    });
}

// Probe for GET_STR_DC_THRESHOLD grid search: times the full powtab+dc conversion at fixed sizes.
// The DC threshold is a compiled-in constant controlling recursion leaf size, so the driver
// (perf/tune.sh or a loop) rebuilds with each candidate value and compares these numbers.
fn tune_get_str_dc_probe() {
    use malachite_nz::natural::conversion::digits::general_digits::{
        digits_in_base_per_limb_for_tuning, get_chars_per_limb, limbs_compute_power_table,
        limbs_digits_power_table_scratch_len_for_tuning,
        limbs_to_digits_small_base_divide_and_conquer_for_tuning,
        limbs_to_digits_small_base_divide_and_conquer_scratch_len_for_tuning,
    };
    const BASE: u64 = 10;
    for n in [32usize, 48, 63] {
        let inputs: Vec<Vec<Limb>> = (0..INPUT_SETS)
            .map(|k| {
                random_primitive_ints(EXAMPLE_SEED.fork(&format!("g{k}")))
                    .take(n)
                    .collect()
            })
            .collect();
        let run = |xs: &[Limb]| {
            let mut digits = [0u8; 64 * 20];
            let mut xs_copy = [0; 64];
            xs_copy[..n].copy_from_slice(xs);
            let mut power_table_memory =
                vec![0; limbs_digits_power_table_scratch_len_for_tuning(n)];
            let digits_len = digits_in_base_per_limb_for_tuning(n, BASE);
            let len = 1 + usize::try_from(digits_len).unwrap() / get_chars_per_limb(BASE);
            let (power_len, powers) =
                limbs_compute_power_table(&mut power_table_memory, len, BASE, None);
            let mut scratch =
                vec![0; limbs_to_digits_small_base_divide_and_conquer_scratch_len_for_tuning(n)];
            black_box(limbs_to_digits_small_base_divide_and_conquer_for_tuning(
                &mut digits[..n * 20],
                &mut xs_copy[..n],
                BASE,
                &powers,
                power_len,
                &mut scratch,
            ));
        };
        // Control series: the basecase at a fixed size, which does not depend on
        // GET_STR_DC_THRESHOLD. Reporting the ratio dc/control makes numbers comparable across
        // rebuilds and runs, canceling out core-scheduling and frequency effects.
        let control_inputs: Vec<Vec<Limb>> = (0..INPUT_SETS)
            .map(|k| {
                random_primitive_ints(EXAMPLE_SEED.fork(&format!("c{k}")))
                    .take(20)
                    .collect()
            })
            .collect();
        let control = |xs: &[Limb]| {
            let mut digits = [0u8; 20 * 20];
            black_box(limbs_to_digits_small_base_basecase(
                &mut digits,
                0,
                xs,
                BASE,
            ));
        };
        for xs in &inputs {
            run(xs);
        }
        for xs in &control_inputs {
            control(xs);
        }
        let mut i = 0usize;
        let mut j = 0usize;
        let (t, t_control) = interleaved_min_pair(
            &mut || {
                let xs = &inputs[i & (INPUT_SETS - 1)];
                i += 1;
                run(xs);
            },
            &mut || {
                let xs = &control_inputs[j & (INPUT_SETS - 1)];
                j += 1;
                control(xs);
            },
        );
        println!(
            "dc_probe n={n}: {t:.1} ns, control {t_control:.1} ns, ratio {:.3}",
            t / t_control
        );
    }
}

// ---------------------------------------------------------------------------------------------
// Allocating-shl shootout: does avoiding the vec![0; n] zero-init pass via Vec::extend pay off, or
// does the iterator plumbing defeat vectorization / the TrustedLen specialization?

// Variant A: zero-init then overwrite (current limbs_shl shape). One "wasted" memset pass, but the
// fill loop is a plain slice loop with reliable codegen.
#[inline(never)]
fn shl_vec_zeroinit(xs: &[Limb], bits: u64) -> Vec<Limb> {
    let len = xs.len();
    let cobits = Limb::WIDTH - bits;
    let mut out = vec![0; len + 1];
    out[0] = xs[0] << bits;
    for (o, w) in out[1..len].iter_mut().zip(xs.windows(2)) {
        *o = (w[1] << bits) | (w[0] >> cobits);
    }
    let carry = xs[len - 1] >> cobits;
    if carry == 0 {
        out.pop();
    } else {
        *out.last_mut().unwrap() = carry;
    }
    out
}

// Variant B: no zero-init; elements are written exactly once via extend. Relies on the TrustedLen
// specialization and on LLVM dissolving the Windows/Map iterator chain.
#[inline(never)]
fn shl_vec_extend(xs: &[Limb], bits: u64) -> Vec<Limb> {
    let len = xs.len();
    let cobits = Limb::WIDTH - bits;
    let mut out = Vec::with_capacity(len + 1);
    out.push(xs[0] << bits);
    out.extend(xs.windows(2).map(|w| (w[1] << bits) | (w[0] >> cobits)));
    let carry = xs[len - 1] >> cobits;
    if carry != 0 {
        out.push(carry);
    }
    out
}

fn tune_shl_alloc() {
    type ShlVecFn = fn(&[Limb], u64) -> Vec<Limb>;
    let variants: [(&str, ShlVecFn); 2] =
        [("zeroinit", shl_vec_zeroinit), ("extend", shl_vec_extend)];
    // Correctness cross-check before timing.
    for n in [1, 2, 5, 17, 100] {
        for bits in [1, 7, 31, 63] {
            let xs: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork("ca"))
                .take(n)
                .collect();
            let reference = variants[0].1(&xs, bits);
            for (name, f) in &variants[1..] {
                assert_eq!(
                    f(&xs, bits),
                    reference,
                    "variant {name} disagrees at n={n}, bits={bits}"
                );
            }
        }
    }
    println!("all variants agree; timing includes alloc+drop (ns/call, ns/limb):");
    for n in [16usize, 64, 256, 1024, 4096] {
        println!("n = {n}:");
        for (name, f) in &variants {
            let inputs: Vec<Vec<Limb>> = (0..INPUT_SETS)
                .map(|k| {
                    random_primitive_ints(EXAMPLE_SEED.fork(&format!("a{k}")))
                        .take(n)
                        .collect()
                })
                .collect();
            for xs in &inputs {
                black_box(shl_vec_zeroinit(xs, 13));
                black_box(f(xs, 13));
            }
            let (mut i, mut j) = (0usize, 0usize);
            let (t_base, t) = interleaved_min_pair(
                &mut || {
                    let xs = &inputs[i & (INPUT_SETS - 1)];
                    i += 1;
                    black_box(shl_vec_zeroinit(black_box(xs), 13));
                },
                &mut || {
                    let xs = &inputs[j & (INPUT_SETS - 1)];
                    j += 1;
                    black_box(f(black_box(xs), 13));
                },
            );
            println!(
                "  {name:>10}: {t:>9.1} ns  {:>6.3} ns/limb  (vs zeroinit {:>5.2}x)",
                t / n as f64,
                t_base / t,
            );
        }
    }
}

fn tune_add() {
    type AddFn = fn(&mut [Limb], &[Limb], &[Limb]) -> bool;
    let variants: [(&str, AddFn); 5] = [
        ("current", add_n_current),
        ("double_limb", add_n_double_limb),
        ("overflowing", add_n_overflowing),
        ("overflowing_x4", add_n_overflowing_x4),
        ("double_limb_x4", add_n_double_limb_x4),
    ];
    // Correctness cross-check before timing.
    for n in [1, 3, 4, 7, 17, 100] {
        let xs: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork("cx"))
            .take(n)
            .collect();
        let ys: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork("cy"))
            .take(n)
            .collect();
        let mut reference = vec![0; n];
        let ref_carry = variants[0].1(&mut reference, &xs, &ys);
        for (name, f) in &variants[1..] {
            let mut out = vec![0; n];
            let carry = f(&mut out, &xs, &ys);
            assert_eq!(
                (out, carry),
                (reference.clone(), ref_carry),
                "variant {name} disagrees at n={n}"
            );
        }
    }
    println!("all variants agree; timing (ns/call, ns/limb):");
    for n in [16usize, 64, 256, 1024, 4096] {
        println!("n = {n}:");
        // Measure each variant interleaved with the current idiom as the common baseline.
        for (name, f) in &variants {
            let inputs: Vec<(Vec<Limb>, Vec<Limb>)> = (0..INPUT_SETS)
                .map(|k| {
                    let xs = random_primitive_ints(EXAMPLE_SEED.fork(&format!("x{k}")))
                        .take(n)
                        .collect();
                    let ys = random_primitive_ints(EXAMPLE_SEED.fork(&format!("y{k}")))
                        .take(n)
                        .collect();
                    (xs, ys)
                })
                .collect();
            let mut out_a = vec![0; n];
            let mut out_b = vec![0; n];
            for (xs, ys) in &inputs {
                add_n_current(&mut out_a, xs, ys);
                f(&mut out_b, xs, ys);
            }
            let (mut i, mut j) = (0usize, 0usize);
            let (t_base, t) = interleaved_min_pair(
                &mut || {
                    let (xs, ys) = &inputs[i & (INPUT_SETS - 1)];
                    i += 1;
                    black_box(add_n_current(black_box(&mut out_a), xs, ys));
                },
                &mut || {
                    let (xs, ys) = &inputs[j & (INPUT_SETS - 1)];
                    j += 1;
                    black_box(f(black_box(&mut out_b), xs, ys));
                },
            );
            println!(
                "  {name:>16}: {t:>9.1} ns  {:>6.3} ns/limb  (vs current {:>5.2}x)",
                t / n as f64,
                t_base / t,
            );
        }
    }
}

// Polynomial evaluation: Horner's rule against divide and conquer. The crossover depends on the
// size of the point as much as on the length, so no single length threshold fits; this prints a
// grid instead. For each (coefficient bits, point bits) pair it gives the ratio of the divide-and-
// conquer time to the Horner time at a range of lengths, and the first length from which divide and
// conquer wins at every length measured. The `EVALUATE_DIVIDE_AND_CONQUER_*` constants in
// `integer_polynomial::arithmetic::evaluate` are read off it.
#[allow(clippy::print_stdout)]
fn tune_evaluate() {
    use malachite_base::num::arithmetic::traits::{ModPowerOf2, Parity};
    use malachite_base::num::logic::traits::BitAccess;
    use malachite_nz::integer::Integer;
    use malachite_nz::integer_polynomial::arithmetic::evaluate::{
        evaluate_divide_and_conquer, evaluate_horner,
    };
    use malachite_nz::natural::Natural;
    let random_integer = |seed: &str, bits: u64| -> Integer {
        let limbs: Vec<Limb> = random_primitive_ints(EXAMPLE_SEED.fork(seed))
            .take(usize::try_from(bits.div_ceil(Limb::WIDTH)).unwrap())
            .collect();
        let mut n = Natural::from_owned_limbs_asc(limbs).mod_power_of_2(bits);
        n.set_bit(bits - 1);
        let negative = random_primitive_ints::<u8>(EXAMPLE_SEED.fork(&format!("{seed}s")))
            .next()
            .unwrap()
            .odd();
        if negative {
            -Integer::from(n)
        } else {
            Integer::from(n)
        }
    };
    let lens =
        [2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048];
    for coefficient_bits in [8u64, 64, 256, 2048] {
        for point_bits in [8u64, 64, 128, 256, 512, 1024, 2048, 8192] {
            let mut line = format!("c {coefficient_bits:>5} x {point_bits:>5}:");
            let mut first_win = None;
            let mut wins = 0;
            for &len in &lens {
                let sets = 4;
                let inputs: Vec<(Vec<Integer>, Integer)> = (0..sets)
                    .map(|k| {
                        let cs = (0..len)
                            .map(|i| random_integer(&format!("c{k}_{i}"), coefficient_bits))
                            .collect();
                        (cs, random_integer(&format!("x{k}"), point_bits))
                    })
                    .collect();
                let (mut i, mut j) = (0usize, 0usize);
                let (th, td) = interleaved_min_pair(
                    &mut || {
                        let (cs, x) = &inputs[i % sets];
                        i += 1;
                        black_box(evaluate_horner(black_box(cs), x));
                    },
                    &mut || {
                        let (cs, x) = &inputs[j % sets];
                        j += 1;
                        black_box(evaluate_divide_and_conquer(black_box(cs), x));
                    },
                );
                let ratio = td / th;
                line.push_str(&format!(" {len}:{ratio:.2}"));
                if ratio < 1.0 {
                    wins += 1;
                    first_win.get_or_insert(len);
                } else {
                    wins = 0;
                    first_win = None;
                }
                if wins >= 3 && ratio < 0.8 {
                    break;
                }
                // Horner is quadratic in the length; stop before a single call takes too long.
                if th > 2.0e8 {
                    break;
                }
            }
            println!(
                "{line}  => {}",
                first_win.map_or("none".to_string(), |l| l.to_string())
            );
        }
    }
}

/// Dispatch a tuning run by key. Keys mirror the bottom-up tuning order; after each level, write
/// the suggested value into platform_64.rs and rebuild before tuning the next level (perf/tune.sh
/// automates this).
pub fn tune(key: &str) {
    match key {
        "add" => tune_add(),
        "evaluate" => tune_evaluate(),
        "shl" => tune_shl(),
        "shl_alloc" => tune_shl_alloc(),
        "get_str_precompute" => tune_get_str_precompute(),
        "get_str_dc_probe" => tune_get_str_dc_probe(),
        "divrem" => tune_divrem(),
        "mul_toom22" => tune_mul_toom22(),
        "mul_toom33" => tune_mul_toom33(),
        "mul_toom44" => tune_mul_toom44(),
        "mul_toom6h" => tune_mul_toom6h(),
        "mul_toom6h_vs_toom33" => tune_mul_toom6h_vs_toom33(),
        "mul_toom8h" => tune_mul_toom8h(),
        "mul_toom32_to_toom43" => tune_mul_toom32_to_toom43(),
        "mul_toom32_to_toom53" => tune_mul_toom32_to_toom53(),
        "mul_toom42_to_toom53" => tune_mul_toom42_to_toom53(),
        "mul_toom42_to_toom63" => tune_mul_toom42_to_toom63(),
        "dc_div_qr" => tune_dc_div_qr(),
        "dc_divappr_q" => tune_dc_divappr_q(),
        "mullo_dc" => tune_mullo_dc(),
        "mullo_mul_n" => tune_mullo_mul_n(),
        "sqrlo_dc" => tune_sqrlo_dc(),
        "dc_bdiv_qr" => tune_dc_bdiv_qr(),
        "dc_bdiv_q" => tune_dc_bdiv_q(),
        "mu_bdiv_qr" => tune_mu_bdiv_qr(),
        "mu_bdiv_q" => tune_mu_bdiv_q(),
        "binv_probe" => tune_binv_probe(),
        "from_digits_dc" => tune_from_digits_dc(),
        "mu_bdiv_qr_vs_schoolbook" => tune_mu_bdiv_qr_vs_schoolbook(),
        "mu_bdiv_q_vs_schoolbook" => tune_mu_bdiv_q_vs_schoolbook(),
        "mod_1_shootout" => tune_mod_1_shootout(),
        "invert_probe" => tune_invert_probe(),
        "mul_fft_probe" => tune_mul_fft_probe(),
        "small_kernel_probe" => tune_small_kernel_probe(),
        "poly_mul_grid" => tune_poly_mul_grid(),
        "poly_mod_power_of_2_mul_grid" => tune_poly_mod_power_of_2_mul_grid(),
        "poly_mod_mul_grid" => tune_poly_mod_mul_grid(),
        "poly_mod_evaluate_geometric_grid" => tune_poly_mod_evaluate_geometric_grid(),
        "poly_mod_mul_middle_grid" => tune_poly_mod_mul_middle_grid(),
        "poly_pow_grid" => tune_poly_pow_grid(),
        "poly_pow_multinomial_grid" => tune_poly_pow_multinomial_grid(),
        "poly_pow_chain_grid" => tune_poly_pow_chain_grid(),
        "gcd_probe" => tune_gcd_probe(),
        "xgcd_probe" => tune_xgcd_probe(),
        "mul_fft" => tune_mul_fft(),
        "sqr_fft" => tune_sqr_fft(),
        "mul_fft_vs_toom6h" => tune_mul_fft_vs_toom6h(),
        "sqr_fft_vs_toom6" => tune_sqr_fft_vs_toom6(),
        "fft_small_check" => tune_fft_small_check(),
        "inv_newton" => tune_inv_newton(),
        "mu_div_qr" => tune_mu_div_qr(),
        "mu_divappr_q" => tune_mu_divappr_q(),
        "sqr_toom2" => tune_sqr_toom2(),
        "sqr_toom3" => tune_sqr_toom3(),
        "sqr_toom4" => tune_sqr_toom4(),
        "sqr_toom6" => tune_sqr_toom6(),
        "sqr_toom6_vs_toom3" => tune_sqr_toom6_vs_toom3(),
        "sqr_toom8" => tune_sqr_toom8(),
        "mul" => {
            tune_mul_toom22();
            tune_mul_toom33();
            tune_mul_toom44();
            println!();
            println!("NOTE: levels above toom22 were measured with the COMPILED-IN lower");
            println!("thresholds. Apply the suggestions to platform_64.rs bottom-up, rebuild,");
            println!("and re-run until stable (perf/tune.sh does this).");
        }
        _ => panic!("Invalid tune key: {key}"),
    }
}