1  // Copyright (c) 2019-2025 Alpha-Delta Network Inc.
  2  // This file is part of the deltavm library.
  3  
  4  // Licensed under the Apache License, Version 2.0 (the "License");
  5  // you may not use this file except in compliance with the License.
  6  // You may obtain a copy of the License at:
  7  
  8  // http://www.apache.org/licenses/LICENSE-2.0
  9  
 10  // Unless required by applicable law or agreed to in writing, software
 11  // distributed under the License is distributed on an "AS IS" BASIS,
 12  // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 13  // See the License for the specific language governing permissions and
 14  // limitations under the License.
 15  
 16  use crate::fft::{DensePolynomial, domain::*};
 17  use deltavm_curves::bls12_377::{Fr, G1Projective};
 18  use deltavm_fields::{FftField, Field, One, Zero};
 19  use deltavm_utilities::rand::{TestRng, Uniform};
 20  use rand::Rng;
 21  
 22  #[test]
 23  fn vanishing_polynomial_evaluation() {
 24      let rng = &mut TestRng::default();
 25  
 26      for coeffs in 0..10 {
 27          let domain = EvaluationDomain::<Fr>::new(coeffs).unwrap();
 28          let z = domain.vanishing_polynomial();
 29          for _ in 0..100 {
 30              let point: Fr = rng.r#gen();
 31              assert_eq!(z.evaluate(point), domain.evaluate_vanishing_polynomial(point))
 32          }
 33      }
 34  }
 35  
 36  #[test]
 37  fn vanishing_polynomial_vanishes_on_domain() {
 38      for coeffs in 0..1000 {
 39          let domain = EvaluationDomain::<Fr>::new(coeffs).unwrap();
 40          let z = domain.vanishing_polynomial();
 41          for point in domain.elements() {
 42              assert!(z.evaluate(point).is_zero())
 43          }
 44      }
 45  }
 46  
 47  #[test]
 48  fn size_of_elements() {
 49      for coeffs in 1..10 {
 50          let size = 1 << coeffs;
 51          let domain = EvaluationDomain::<Fr>::new(size).unwrap();
 52          let domain_size = domain.size();
 53          assert_eq!(domain_size, domain.elements().count());
 54      }
 55  }
 56  
 57  #[test]
 58  fn elements_contents() {
 59      for coeffs in 1..10 {
 60          let size = 1 << coeffs;
 61          let domain = EvaluationDomain::<Fr>::new(size).unwrap();
 62          for (i, element) in domain.elements().enumerate() {
 63              assert_eq!(element, domain.group_gen.pow([i as u64]));
 64          }
 65      }
 66  }
 67  
 68  /// Test that lagrange interpolation for a random polynomial at a random
 69  /// point works.
 70  #[test]
 71  fn non_systematic_lagrange_coefficients_test() {
 72      let mut rng = TestRng::default();
 73  
 74      for domain_dim in 1..10 {
 75          let domain_size = 1 << domain_dim;
 76          let domain = EvaluationDomain::<Fr>::new(domain_size).unwrap();
 77          // Get random pt + lagrange coefficients
 78          let rand_pt = Fr::rand(&mut rng);
 79          let lagrange_coeffs = domain.evaluate_all_lagrange_coefficients(rand_pt);
 80  
 81          // Sample the random polynomial, evaluate it over the domain and the random
 82          // point.
 83          let rand_poly = DensePolynomial::<Fr>::rand(domain_size - 1, &mut rng);
 84          let poly_evals = domain.fft(rand_poly.coeffs());
 85          let actual_eval = rand_poly.evaluate(rand_pt);
 86  
 87          // Do lagrange interpolation, and compare against the actual evaluation
 88          let mut interpolated_eval = Fr::zero();
 89          for i in 0..domain_size {
 90              interpolated_eval += lagrange_coeffs[i] * poly_evals[i];
 91          }
 92          assert_eq!(actual_eval, interpolated_eval);
 93      }
 94  }
 95  
 96  /// Test that lagrange coefficients for a point in the domain is correct
 97  #[test]
 98  fn systematic_lagrange_coefficients_test() {
 99      // This runs in time O(N^2) in the domain size, so keep the domain dimension
100      // low. We generate lagrange coefficients for each element in the domain.
101      for domain_dim in 1..5 {
102          let domain_size = 1 << domain_dim;
103          let domain = EvaluationDomain::<Fr>::new(domain_size).unwrap();
104          let all_domain_elements: Vec<Fr> = domain.elements().collect();
105          for (i, element) in all_domain_elements.iter().enumerate().take(domain_size) {
106              let lagrange_coeffs = domain.evaluate_all_lagrange_coefficients(*element);
107              for (j, &coeff) in lagrange_coeffs.iter().enumerate().take(domain_size) {
108                  // Lagrange coefficient for the evaluation point, which should be 1
109                  if i == j {
110                      assert_eq!(coeff, Fr::one());
111                  } else {
112                      assert_eq!(coeff, Fr::zero());
113                  }
114              }
115          }
116      }
117  }
118  
119  #[test]
120  fn test_fft_correctness() {
121      // Tests that the ffts output the correct result.
122      // This assumes a correct polynomial evaluation at point procedure.
123      // It tests consistency of FFT/IFFT, and coset_fft/coset_ifft,
124      // along with testing that each individual evaluation is correct.
125  
126      // Runs in time O(degree^2)
127      let log_degree = 5;
128      let degree = 1 << log_degree;
129      let rand_poly = DensePolynomial::<Fr>::rand(degree - 1, &mut TestRng::default());
130  
131      for log_domain_size in log_degree..(log_degree + 2) {
132          let domain_size = 1 << log_domain_size;
133          let domain = EvaluationDomain::<Fr>::new(domain_size).unwrap();
134          let poly_evals = domain.fft(&rand_poly.coeffs);
135          let poly_coset_evals = domain.coset_fft(&rand_poly.coeffs);
136          for (i, x) in domain.elements().enumerate() {
137              let coset_x = Fr::multiplicative_generator() * x;
138  
139              assert_eq!(poly_evals[i], rand_poly.evaluate(x));
140              assert_eq!(poly_coset_evals[i], rand_poly.evaluate(coset_x));
141          }
142  
143          let rand_poly_from_subgroup = DensePolynomial::from_coefficients_vec(domain.ifft(&poly_evals));
144          let rand_poly_from_coset = DensePolynomial::from_coefficients_vec(domain.coset_ifft(&poly_coset_evals));
145  
146          assert_eq!(rand_poly, rand_poly_from_subgroup, "degree = {degree}, domain size = {domain_size}");
147          assert_eq!(rand_poly, rand_poly_from_coset, "degree = {degree}, domain size = {domain_size}");
148      }
149  }
150  
151  #[test]
152  fn test_roots_of_unity() {
153      // Tests that the roots of unity result is the same as domain.elements()
154      let max_degree = 10;
155      for log_domain_size in 0..max_degree {
156          let domain_size = 1 << log_domain_size;
157          let domain = EvaluationDomain::<Fr>::new(domain_size).unwrap();
158          let actual_roots = domain.roots_of_unity(domain.group_gen);
159          for &value in &actual_roots {
160              assert!(domain.evaluate_vanishing_polynomial(value).is_zero());
161          }
162          let expected_roots_elements = domain.elements();
163          for (expected, &actual) in expected_roots_elements.zip(&actual_roots) {
164              assert_eq!(expected, actual);
165          }
166          assert_eq!(actual_roots.len(), domain_size / 2);
167      }
168  }
169  
170  #[test]
171  #[cfg(not(feature = "serial"))]
172  fn parallel_fft_consistency() {
173      // This implements the Cooley-Turkey FFT, derived from libfqfft
174      // The libfqfft implementation uses pseudocode from [CLRS 2n Ed, pp. 864].
175      fn serial_radix2_fft(a: &mut [Fr], omega: Fr, log_n: u32) {
176          #[inline]
177          pub(crate) fn bitreverse(mut n: u32, l: u32) -> u32 {
178              let mut r = 0;
179              for _ in 0..l {
180                  r = (r << 1) | (n & 1);
181                  n >>= 1;
182              }
183              r
184          }
185          use core::convert::TryFrom;
186          let n = u32::try_from(a.len()).expect("cannot perform FFTs larger on vectors of len > (1 << 32)");
187          assert_eq!(n, 1 << log_n);
188  
189          // swap coefficients in place
190          for k in 0..n {
191              let rk = bitreverse(k, log_n);
192              if k < rk {
193                  a.swap(rk as usize, k as usize);
194              }
195          }
196  
197          let mut m = 1;
198          for _i in 1..=log_n {
199              // w_m is 2^i-th root of unity
200              let w_m = omega.pow([(n / (2 * m)) as u64]);
201  
202              let mut k = 0;
203              while k < n {
204                  // w = w_m^j at the start of every loop iteration
205                  let mut w = Fr::one();
206                  for j in 0..m {
207                      let mut t = a[(k + j + m) as usize];
208                      t *= w;
209                      let mut tmp = a[(k + j) as usize];
210                      tmp -= t;
211                      a[(k + j + m) as usize] = tmp;
212                      a[(k + j) as usize] += t;
213                      w *= &w_m;
214                  }
215  
216                  k += 2 * m;
217              }
218  
219              m *= 2;
220          }
221      }
222  
223      fn serial_radix2_ifft(a: &mut [Fr], omega: Fr, log_n: u32) {
224          serial_radix2_fft(a, omega.inverse().unwrap(), log_n);
225          let domain_size_inv = Fr::from(a.len() as u64).inverse().unwrap();
226          for coeff in a.iter_mut() {
227              *coeff *= domain_size_inv;
228          }
229      }
230  
231      fn serial_radix2_coset_fft(a: &mut [Fr], omega: Fr, log_n: u32) {
232          let coset_shift = Fr::multiplicative_generator();
233          let mut cur_pow = Fr::one();
234          for coeff in a.iter_mut() {
235              *coeff *= cur_pow;
236              cur_pow *= coset_shift;
237          }
238          serial_radix2_fft(a, omega, log_n);
239      }
240  
241      fn serial_radix2_coset_ifft(a: &mut [Fr], omega: Fr, log_n: u32) {
242          serial_radix2_ifft(a, omega, log_n);
243          let coset_shift = Fr::multiplicative_generator().inverse().unwrap();
244          let mut cur_pow = Fr::one();
245          for coeff in a.iter_mut() {
246              *coeff *= cur_pow;
247              cur_pow *= coset_shift;
248          }
249      }
250  
251      fn test_consistency<R: Rng>(rng: &mut R, max_coeffs: u32) {
252          for _ in 0..5 {
253              for log_d in 0..max_coeffs {
254                  let d = 1 << log_d;
255  
256                  let expected_poly = (0..d).map(|_| Fr::rand(rng)).collect::<Vec<_>>();
257                  let mut expected_vec = expected_poly.clone();
258                  let mut actual_vec = expected_vec.clone();
259  
260                  let domain = EvaluationDomain::new(d).unwrap();
261  
262                  serial_radix2_fft(&mut expected_vec, domain.group_gen, log_d);
263                  domain.fft_in_place(&mut actual_vec);
264                  assert_eq!(expected_vec, actual_vec);
265  
266                  serial_radix2_ifft(&mut expected_vec, domain.group_gen, log_d);
267                  domain.ifft_in_place(&mut actual_vec);
268                  assert_eq!(expected_vec, actual_vec);
269                  assert_eq!(expected_vec, expected_poly);
270  
271                  serial_radix2_coset_fft(&mut expected_vec, domain.group_gen, log_d);
272                  domain.coset_fft_in_place(&mut actual_vec);
273                  assert_eq!(expected_vec, actual_vec);
274  
275                  serial_radix2_coset_ifft(&mut expected_vec, domain.group_gen, log_d);
276                  domain.coset_ifft_in_place(&mut actual_vec);
277                  assert_eq!(expected_vec, actual_vec);
278              }
279          }
280      }
281  
282      let rng = &mut TestRng::default();
283  
284      test_consistency(rng, 10);
285  }
286  
287  #[test]
288  fn fft_composition() {
289      fn test_fft_composition<F: FftField, T: crate::fft::DomainCoeff<F> + Uniform + core::fmt::Debug + Eq, R: Rng>(
290          rng: &mut R,
291          max_coeffs: usize,
292      ) {
293          for coeffs in 0..max_coeffs {
294              let coeffs = 1 << coeffs;
295  
296              let domain = EvaluationDomain::new(coeffs).unwrap();
297  
298              let mut v = vec![];
299              for _ in 0..coeffs {
300                  v.push(T::rand(rng));
301              }
302              // Fill up with zeros.
303              v.resize(domain.size(), T::zero());
304              let mut v2 = v.clone();
305  
306              domain.ifft_in_place(&mut v2);
307              domain.fft_in_place(&mut v2);
308              assert_eq!(v, v2, "ifft(fft(.)) != iden");
309  
310              domain.fft_in_place(&mut v2);
311              domain.ifft_in_place(&mut v2);
312              assert_eq!(v, v2, "fft(ifft(.)) != iden");
313  
314              domain.coset_ifft_in_place(&mut v2);
315              domain.coset_fft_in_place(&mut v2);
316              assert_eq!(v, v2, "coset_fft(coset_ifft(.)) != iden");
317  
318              domain.coset_fft_in_place(&mut v2);
319              domain.coset_ifft_in_place(&mut v2);
320              assert_eq!(v, v2, "coset_ifft(coset_fft(.)) != iden");
321          }
322      }
323  
324      let rng = &mut TestRng::default();
325  
326      test_fft_composition::<Fr, Fr, _>(rng, 10);
327      test_fft_composition::<Fr, G1Projective, _>(rng, 10);
328  }
329  
330  #[test]
331  fn evaluate_over_domain() {
332      let rng = &mut TestRng::default();
333  
334      for domain_size in (1..10).map(|i| 2usize.pow(i)) {
335          let domain = EvaluationDomain::<Fr>::new(domain_size).unwrap();
336          for degree in [domain_size - 2, domain_size - 1, domain_size + 10] {
337              let p = DensePolynomial::rand(degree, rng);
338              assert_eq!(
339                  p.evaluate_over_domain_by_ref(domain).evaluations,
340                  domain.elements().map(|e| p.evaluate(e)).collect::<Vec<_>>()
341              );
342          }
343      }
344  }