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

// Copyright 2023 RISC Zero, Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

//! Polynomial utilities (currently only those used in polynomial evaluation).

use alloc::vec;

use crate::field::ExtElem;

/// Evaluate a polynomial whose coefficients are in the extension field at a
/// point.
pub fn poly_eval<E: ExtElem>(coeffs: &[E], x: E) -> E {
    let mut mul = E::ONE;
    let mut tot = E::ZERO;
    for i in 0..coeffs.len() {
        tot += coeffs[i] * mul;
        mul *= x;
    }
    tot
}

/// General purpose polynomial interpolation.
///
/// Given the goal value f(x) at a set of evalation points x, compute
/// coefficients.
pub fn poly_interpolate<E: ExtElem>(out: &mut [E], x: &[E], fx: &[E], size: usize) {
    // Special case the very easy ones
    if size == 1 {
        out[0] = fx[0];
        return;
    }
    if size == 2 {
        out[1] = (fx[1] - fx[0]) * (x[1] - x[0]).inv();
        out[0] = fx[0] - out[1] * x[0];
        return;
    }
    // Compute ft = product of (x - x_i) for all i
    let mut ft = vec![E::ZERO; size + 1];
    ft[0] = E::ONE;
    for i in 0..size {
        for j in (0..i + 1).rev() {
            let value = ft[j];
            ft[j + 1] += value;
            ft[j] *= -x[i];
        }
    }
    // Clear output
    for i in 0..size {
        out[i] = E::ZERO;
    }
    for i in 0..size {
        // Compute fr = ft / (x - x_i)
        let mut fr = ft.clone();
        poly_divide(&mut fr, x[i]);
        // Evaluate at x[i]
        let fr_xi = poly_eval(&fr, x[i]);
        // Compute multiplier (fx[i] / fr_xi)
        let mul = fx[i] * fr_xi.inv();
        // Multiply into output
        for j in 0..size {
            out[j] += mul * fr[j];
        }
    }
}

/// In-place polynomial division.
///
/// Take the coefficients in P, and divide by (X - z) for some z, return the
/// remainder.
pub fn poly_divide<E: ExtElem>(p: &mut [E], z: E) -> E {
    let mut cur = E::ZERO;
    for i in (0..p.len()).rev() {
        let next = z * cur + p[i];
        p[i] = cur;
        cur = next;
    }
    cur
}