numerical_analysis 0.1.1

// build.rs

use std::fs::OpenOptions;
use std::io::{BufWriter, Write};
use std::path::Path;

pub fn legendre_nodes(size: usize) -> Vec<f64>
{
    use faer::Mat;

    let n = size;
    let mut matrix = Mat::<f64>::zeros(size, size);

    let a_fn = |i| (2. * i as f64 - 1.) / i as f64;
    let c_fn = |i| (i as f64 - 1.) / i as f64;
    for i in 1..=n
    {
        let a = a_fn(i);
        let an = a_fn(i + 1);
        let cn = c_fn(i + 1);

        let beta = (cn / (a * an)).sqrt();

        if i < size
        {
            matrix[(i, i - 1)] = beta;
            matrix[(i - 1, i)] = beta;
        }
    }

    let evals = matrix.eigenvalues().unwrap();
    let res: Vec<_> = evals.iter().copied().collect();
    let mut vals: Vec<_> = res.iter().map(|v| v.re).collect();
    vals.sort_by(|a, b| a.partial_cmp(b).unwrap());

    vals
}

pub fn legendre_weights(legendre_nodes: &[f64]) -> Vec<f64>
{
    let n = legendre_nodes.len();

    // Multiple simplifications were done based on the fact that xi are the roots of
    // the polynomial p_n.
    let mut res = Vec::new();
    for xi in legendre_nodes
    {
        let mut p_n2 = 1.;
        let mut p_n1 = *xi;

        for i in 2..=n
        {
            let i = i as f64;
            let next = ((2. * i - 1.) * xi * p_n1 - (i - 1.) * p_n2) / i;

            p_n2 = p_n1;
            p_n1 = next;
        }

        let n = n as f64;
        let n1 = n * p_n2;
        let d = 1. - xi * xi;
        let w = if n > 0. { (2. * d) / (n1 * n1) } else { 2. };

        res.push(w);
    }

    res
}

fn main()
{
    let out_dir = std::env::var_os("OUT_DIR").unwrap();

    let dest_path = Path::new(&out_dir).join("legendre_roots.rs");
    let _ = std::fs::write(&dest_path, "");

    let f = OpenOptions::new()
        .write(true)
        .truncate(true)
        .open(dest_path)
        .expect("unable to open file");
    let mut f = BufWriter::new(f);

    let max_node_count = 100;
    let total_length = max_node_count * (max_node_count + 1) / 2;

    write!(
        f,
        "pub (crate) const LEGENDRE_NODE_TABLES: [f64; {total_length}] = [\n"
    )
    .unwrap();

    let mut weights = Vec::new();
    for node_count in 1..=max_node_count
    {
        let values = legendre_nodes(node_count);
        weights.extend(legendre_weights(&values));
        let string = format!(
            "\t// n = {}\n{},\n",
            node_count,
            values
                .iter()
                .map(|x| {
                    let mut v = "\t\t".to_string();
                    v.push_str(&format!("{:.25}", x).to_string());
                    v
                })
                .collect::<Vec<String>>()
                .join(",\n")
        );
        write!(f, "{}", string).unwrap();
    }
    write!(f, "];\n\n").unwrap();

    write!(
        f,
        "pub (crate) const LEGENDRE_WEIGHT_TABLES: [f64; {total_length}] = [\n"
    )
    .unwrap();
    for node_count in 1..=max_node_count
    {
        let start = (node_count - 1) * node_count / 2;
        let values = &weights[start..start + node_count];
        let string = format!(
            "\t// n = {}\n{},\n",
            node_count,
            values
                .iter()
                .map(|x| {
                    let mut v = "\t\t".to_string();
                    v.push_str(&format!("{:.20}", x).to_string());
                    v
                })
                .collect::<Vec<String>>()
                .join(",\n")
        );
        write!(f, "{}", string).unwrap();
    }
    write!(f, "];").unwrap();
    println!("cargo::rerun-if-changed=build.rs");
}