oxinum 0.1.3

OxiNum — Pure-Rust arbitrary-precision math (GMP/MPFR-free via dashu)
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

//! Solve a 3x3 linear system exactly over `RBig` using Cramer's rule.
//!
//! Demonstrates exact (no-rounding) arithmetic through the `oxinum` facade.
//!
//! Run:
//!     cargo run --example exact_rational_linear_solve
//!
//! The example solves
//!     2x +  y -  z =  8
//!    -3x -  y + 2z = -11
//!    -2x +  y + 2z = -3
//! whose unique solution is (x, y, z) = (2, 3, -1).

use oxinum::{IBig, RBig};

fn det_3x3(m: &[[RBig; 3]; 3]) -> RBig {
    &m[0][0] * &(&(&m[1][1] * &m[2][2]) - &(&m[1][2] * &m[2][1]))
        - &m[0][1] * &(&(&m[1][0] * &m[2][2]) - &(&m[1][2] * &m[2][0]))
        + &m[0][2] * &(&(&m[1][0] * &m[2][1]) - &(&m[1][1] * &m[2][0]))
}

fn main() {
    let r = |n: i64, d: i64| RBig::from_parts_signed(IBig::from(n), IBig::from(d));
    let a = [
        [r(2, 1), r(1, 1), r(-1, 1)],
        [r(-3, 1), r(-1, 1), r(2, 1)],
        [r(-2, 1), r(1, 1), r(2, 1)],
    ];
    let b = [r(8, 1), r(-11, 1), r(-3, 1)];

    let det_a = det_3x3(&a);
    assert_ne!(
        det_a,
        RBig::from_parts_signed(IBig::from(0), IBig::from(1)),
        "system is singular"
    );

    // Cramer's rule: x_i = det(A_i) / det(A), where A_i replaces column i with b.
    let mut solution: [RBig; 3] = [
        RBig::from_parts_signed(IBig::from(0), IBig::from(1)),
        RBig::from_parts_signed(IBig::from(0), IBig::from(1)),
        RBig::from_parts_signed(IBig::from(0), IBig::from(1)),
    ];

    for i in 0..3 {
        let mut a_i = a.clone();
        for (row, b_val) in b.iter().enumerate() {
            a_i[row][i] = b_val.clone();
        }
        let det_i = det_3x3(&a_i);
        solution[i] = &det_i / &det_a;
    }

    println!("System:");
    println!("   2x +  y -  z =  8");
    println!("  -3x -  y + 2z = -11");
    println!("  -2x +  y + 2z = -3");
    println!();
    println!(
        "Exact solution: x = {}, y = {}, z = {}",
        solution[0], solution[1], solution[2]
    );

    // Verify A * solution == b exactly.
    for (row_idx, row) in a.iter().enumerate() {
        let lhs =
            &(&(&row[0] * &solution[0]) + &(&row[1] * &solution[1])) + &(&row[2] * &solution[2]);
        let rhs = &b[row_idx];
        assert_eq!(&lhs, rhs, "row {row_idx}: lhs={lhs}, rhs={rhs}");
    }
    println!("Verification: A * x == b (exact)");
}