numerica 1.2.1

Open-source math library for exact and floating point computations
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170

# Numerica

<p align="center">
<a href="https://symbolica.io"><img alt="Symbolica website" src="https://img.shields.io/static/v1?label=symbolica&message=website&color=orange&style=flat-square"></a>
  <a href="https://zulip.symbolica.io"><img alt="Zulip Chat" src="https://img.shields.io/static/v1?label=zulip&message=discussions&color=blue&style=flat-square"></a>
    <a href="https://github.com/symbolica-dev/numerica"><img alt="Numerica repository" src="https://img.shields.io/static/v1?label=github&message=development&color=green&style=flat-square&logo=github"></a>
    <a href="https://app.codecov.io/gh/symbolica-dev/numerica"><img alt="Codecov" src="https://img.shields.io/codecov/c/github/symbolica-dev/numerica?token=W5GTATIVZI&style=flat-square"></a>
</p>

Numerica is an open-source mathematics library for Rust, that provides high-performance number types, such as error-tracking floats and finite field elements.

It provides
- Abstractions over rings, Euclidean domains, fields and floats
- High-performance Integer with automatic up-and-downgrading to arbitrary precision types
- Rational numbers with reconstruction algorithms
- Fast finite field arithmetic
- Error-tracking floating point types
- Generic dual numbers for automatic (higher-order) differentiation
- Matrix operations and linear system solving
- Numerical integration using Vegas algorithm with discrete layer support

For operations on symbols, check out the parent project [Symbolica](https://symbolica.io).


# Examples


### Solving a linear system

Solve a linear system over the rationals:

```rust
let a = Matrix::from_linear(
    vec![
        1.into(), 2.into(), 3.into(),
        4.into(), 5.into(), 16.into(),
        7.into(), 8.into(), 9.into(),
    ],
    3, 3, Q,
)
.unwrap();

let b = Matrix::new_vec(vec![1.into(), 2.into(), 3.into()], Q);

let r = a.solve(&b).unwrap();
assert_eq!(r.into_vec(), [(-1, 3), (2, 3), (0, 1)]);
```
Solution: $(-1/3, 2/3, 0)$.

Solve over the finite field $\mathbb{Z}_7$:

```rust
let z_7 = Zp::new(7);

let a = Matrix::from_linear(
    vec![
        z_7.to_element(5), z_7.to_element(8),
        z_7.to_element(2), z_7.to_element(1),
    ],
    2, 2, z_7,
)
.unwrap();

let r = a.inv();
assert_eq!(r.into_vec(), [z_7.to_element(4), z_7.to_element(3), z_7.to_element(2), z_7.to_element(0)]);
```

### Error-tracking floating points

Wrap `f64` and `Float` in an `ErrorPropagatingFloat` to propagate errors through
your computation. For example, a number with 60 accurate digits only has 10 remaining after the following operations:

```rust
let a = ErrorPropagatingFloat::new(Float::with_val(200, 1e-50), 60.);
let r = (a.exp() - a.one()) / a; // large cancellation
assert_eq!(format!("{r}"), "1.000000000");
assert_eq!(r.get_precision(), Some(10.205999132796238));
```

### Automatic differentiation with dual numbers

Create a dual number that fits your needs (supports multiple variables and higher-order differentiation).
Here, we create a simple dual number in three variables:

```rust
create_hyperdual_single_derivative!(Dual, 3);

fn main() {
    let x = Dual::<Rational>::new_variable(0, (1, 1).into());
    let y = Dual::new_variable(1, (2, 1).into());
    let z = Dual::new_variable(2, (3, 1).into());

    let t3 = x * y * z;

    println!("{}", t3.inv());
}
```
It yields `(1/6)+(-1/6)*ε0+(-1/12)*ε1+(-1/18)*ε2`.

The multiplication table is computed and unrolled at compile time for maximal performance.

### Solve integer relations

Solve 

$$
-32.0177 c_1 + 3.1416 c_2 + 2.7183 c_3 = 0
$$

over the integers using PSLQ:

```rust
let result = Integer::solve_integer_relation(
    &[F64::from(-32.0177), F64::from(3.1416), F64::from(2.7183)],
    F64::from(1e-4),
    1,
    Some(Integer::from(100000u64)),
    None,
)
.unwrap();

assert_eq!(result, &[1, 5, 6]);
```

Or via LLL basis reduction:

```rust
let v1 = Vector::new(vec![1.into(), 0.into(), 0.into(), 31416.into()], Z);
let v2 = Vector::new(vec![0.into(), 1.into(), 0.into(), 27183.into()], Z);
let v3 = Vector::new(vec![0.into(), 0.into(), 1.into(), (-320177).into()], Z);

let basis = Vector::basis_reduction(&[v1, v2, v3], (3, 4).into());

assert_eq!(basis[0].into_vec(), [5, 6, 1, 1]);
```

### Numerical integration

Integrate $sin(x \pi) + y$ using the Vegas algorithm and using random numbers suitable for Monte Carlo integration:

```rust
let f = |x: &[f64]| (x[0] * std::f64::consts::PI).sin() + x[1];
let mut grid = Grid::Continuous(ContinuousGrid::new(2, 128, 100, None, false));
let mut rng = MonteCarloRng::new(0, 0);
let mut sample = Sample::new();
for iteration in 1..20 {
     for _ in 0..10_000 {
         grid.sample(&mut rng, &mut sample);
    
         if let Sample::Continuous(_cont_weight, xs) = &sample {
             grid.add_training_sample(&sample, f(xs)).unwrap();
         }
     }
    
     grid.update(1.5, 1.5);
    
     println!(
         "Integral at iteration {}: {}",
         iteration,
         grid.get_statistics().format_uncertainty()
     );
}
```




## Development

Follow the development and discussions on [Zulip](https://reform.zulipchat.com)!