Don't hardcode XI_DEGREES and TAU_DEGREES

Generate them with const functions instead.

Author: dalcde

ext/crates/algebra/src/algebra/combinatorics.rs

@@ -4,42 +4,47 @@ use fp::vector::{FpVector, FpVectorT};
 
 pub const MAX_XI_TAU : usize = fp::prime::MAX_MULTINOMIAL_LEN;
 
-// Generated by Mathematica:
-// "[\n    " <> # <> "\n]" &[
-//  StringJoin @@ 
-//   StringReplace[
-//    ToString /@ 
-//     Riffle[Map[If[# > 2^31, 0, #] &, 
-//       Function[p, Function[k, (p^k - 1)/(p - 1)] /@ Range[10]] /@ 
-//        Prime[Range[8]], {2}], ",\n    "], {"{" -> "[", "}" -> "]"}]]
-static XI_DEGREES : [[i32; 10]; 8] = [
-    [1, 3, 7, 15, 31, 63, 127, 255, 511, 1023],
-    [1, 4, 13, 40, 121, 364, 1093, 3280, 9841, 29524],
-    [1, 6, 31, 156, 781, 3906, 19531, 97656, 488281, 2441406],
-    [1, 8, 57, 400, 2801, 19608, 137257, 960800, 6725601, 47079208],
-    [1, 12, 133, 1464, 16105, 177156, 1948717, 21435888, 235794769, 0],
-    [1, 14, 183, 2380, 30941, 402234, 5229043, 67977560, 883708281, 0],
-    [1, 18, 307, 5220, 88741, 1508598, 25646167, 435984840, 0, 0],
-    [1, 20, 381, 7240, 137561, 2613660, 49659541, 943531280, 0, 0]
-];
-
-// Generated by Mathematica:
-// "[\n    " <> # <> "\n]" &[
-//  StringJoin @@ 
-//   StringReplace[
-//    ToString /@ 
-//     Riffle[Map[If[# > 2^31, 0, #] &, 
-//       Function[p, Function[k, 2 p^k - 1] /@ Range[10]] /@ Prime[Range[8]], {2}], ",\n    "], {"{" -> "[", "}" -> "]"}]]
-static TAU_DEGREES : [[i32; 10]; 8] = [
-    [1, 3, 7, 15, 31, 63, 127, 255, 511, 1023],
-    [1, 5, 17, 53, 161, 485, 1457, 4373, 13121, 39365],
-    [1, 9, 49, 249, 1249, 6249, 31249, 156249, 781249, 3906249],
-    [1, 13, 97, 685, 4801, 33613, 235297, 1647085, 11529601, 80707213],
-    [1, 21, 241, 2661, 29281, 322101, 3543121, 38974341, 428717761, 0],
-    [1, 25, 337, 4393, 57121, 742585, 9653617, 125497033, 1631461441, 0],
-    [1, 33, 577, 9825, 167041, 2839713, 48275137, 820677345, 0, 0],
-    [1, 37, 721, 13717, 260641, 4952197, 94091761, 1787743477, 0, 0]
-];
+macro_rules! const_for {
+    ($i:ident in $a:literal .. $b:ident $contents:block) => {
+        let mut $i = 0;
+        while $i < $b {
+            $contents;
+            $i += 1;
+        }
+    };
+}
+
+/// If p is the nth prime, then XI_DEGREES[n][i - 1] is the degree of ξ_i at the prime p divided by
+/// q, where q = 2p - 2 if p != 2 and 1 if p = 2.
+const XI_DEGREES : [[i32; MAX_XI_TAU]; NUM_PRIMES] = {
+    let mut res = [[0; 10]; 8];
+    const_for! { p_idx in 0 .. NUM_PRIMES {
+        let p = PRIMES[p_idx];
+        let mut p_to_the_i = p;
+        const_for! { x in 0 .. MAX_XI_TAU {
+            res[p_idx][x] = ((p_to_the_i - 1) / (p - 1)) as i32;
+            // At some point the powers overflow. The values are not going to be useful, so use
+            // something replacement that is suitable for const evaluation.
+            p_to_the_i = p_to_the_i.overflowing_mul(p).0;
+        }}
+    }}
+    res
+};
+
+/// If p is the nth prime, then TAU_DEGREES[n][i] is the degree of τ_i at the prime p. Its value is
+/// nonsense at the prime 2
+const TAU_DEGREES : [[i32; MAX_XI_TAU]; NUM_PRIMES] = {
+    let mut res = [[0; 10]; 8];
+    const_for! { p_idx in 0 .. NUM_PRIMES {
+        let p = PRIMES[p_idx];
+        let mut p_to_the_i: u32 = 1;
+        const_for! { x in 0 .. MAX_XI_TAU {
+            res[p_idx][x] = (2_u32.overflowing_mul(p_to_the_i).0 - 1) as i32;
+            p_to_the_i = p_to_the_i.overflowing_mul(p).0;
+        }}
+    }}
+    res
+};
 
 pub fn adem_relation_coefficient(p : ValidPrime, x : u32, y : u32, j : u32, e1 : u32, e2 : u32) -> u32{
     let pi32 = *p as i32;
@@ -473,4 +478,4 @@ mod tests {
         assert_eq!(result, expected);
     }
 
-}
+}

ext/crates/fp/src/prime.rs

@@ -472,6 +472,8 @@ mod tests {
 
 }
 
+pub const PRIMES: [u32; NUM_PRIMES] = [2, 3, 5, 7, 11, 13, 17, 19];
+
 pub const PRIME_TO_INDEX_MAP : [usize; MAX_PRIME+1] = [
     !1, !1, 0, 1, !1, 2, !1, 3, !1, !1, !1, 4, !1, 5, !1, !1, !1, 6, !1, 7
 ];