///| Q = 8380417, NTT 的模数，也是多项式环的模数

///| F = (R^2 / N) mod Q = 41978
/// F is used in INTT to normalize coefficients and maintain Montgomery domain
/// F 用于 INTT 中归一化系数并保持蒙哥马利域
/// 其中：R = 2^32 mod Q (蒙哥马利参数), N = 256 (多项式次数)
const F : Int64 = 41978 // (R ^ {2} / N) % Q

///| Twiddle factors for NTT in Montgomery domain
/// NTT 的旋转因子，已转换到蒙哥马利域
/// 
/// zetas[i] = r^{bit_reverse(i, 8)} * R mod Q 
/// where / 其中:
/// - r is root of unity(原根), r^{2N} = 1 mod Q, r = 1753
/// - R = 2^32 mod Q = 4193792 (Montgomery parameter / 蒙哥马利参数)
/// - bit_reverse(i, 8): 8-bit bit reversal of i / i 的 8 位比特翻转
/// 
/// Uses centered representation: zetas[i] ∈ [-Q/2, Q/2)
/// 使用中心化表示：zetas[i] 的取值范围为 [-Q/2, Q/2)
/// 
/// Note: zetas[0] is not used in the algorithm, set to R for consistency
/// 注意：zetas[0] 在算法中未使用，设为 R 以保持一致性
let zetas = [
  4193792, 25847, -2608894, -518909, 237124, -777960, -876248, 466468, 1826347, 2353451,
  -359251, -2091905, 3119733, -2884855, 3111497, 2680103, 2725464, 1024112, -1079900,
  3585928, -549488, -1119584, 2619752, -2108549, -2118186, -3859737, -1399561, -3277672,
  1757237, -19422, 4010497, 280005, 2706023, 95776, 3077325, 3530437, -1661693, -3592148,
  -2537516, 3915439, -3861115, -3043716, 3574422, -2867647, 3539968, -300467, 2348700,
  -539299, -1699267, -1643818, 3505694, -3821735, 3507263, -2140649, -1600420, 3699596,
  811944, 531354, 954230, 3881043, 3900724, -2556880, 2071892, -2797779, -3930395,
  -1528703, -3677745, -3041255, -1452451, 3475950, 2176455, -1585221, -1257611, 1939314,
  -4083598, -1000202, -3190144, -3157330, -3632928, 126922, 3412210, -983419, 2147896,
  2715295, -2967645, -3693493, -411027, -2477047, -671102, -1228525, -22981, -1308169,
  -381987, 1349076, 1852771, -1430430, -3343383, 264944, 508951, 3097992, 44288,
  -1100098, 904516, 3958618, -3724342, -8578, 1653064, -3249728, 2389356, -210977,
  759969, -1316856, 189548, -3553272, 3159746, -1851402, -2409325, -177440, 1315589,
  1341330, 1285669, -1584928, -812732, -1439742, -3019102, -3881060, -3628969, 3839961,
  2091667, 3407706, 2316500, 3817976, -3342478, 2244091, -2446433, -3562462, 266997,
  2434439, -1235728, 3513181, -3520352, -3759364, -1197226, -3193378, 900702, 1859098,
  909542, 819034, 495491, -1613174, -43260, -522500, -655327, -3122442, 2031748,
  3207046, -3556995, -525098, -768622, -3595838, 342297, 286988, -2437823, 4108315,
  3437287, -3342277, 1735879, 203044, 2842341, 2691481, -2590150, 1265009, 4055324,
  1247620, 2486353, 1595974, -3767016, 1250494, 2635921, -3548272, -2994039, 1869119,
  1903435, -1050970, -1333058, 1237275, -3318210, -1430225, -451100, 1312455, 3306115,
  -1962642, -1279661, 1917081, -2546312, -1374803, 1500165, 777191, 2235880, 3406031,
  -542412, -2831860, -1671176, -1846953, -2584293, -3724270, 594136, -3776993, -2013608,
  2432395, 2454455, -164721, 1957272, 3369112, 185531, -1207385, -3183426, 162844,
  1616392, 3014001, 810149, 1652634, -3694233, -1799107, -3038916, 3523897, 3866901,
  269760, 2213111, -975884, 1717735, 472078, -426683, 1723600, -1803090, 1910376,
  -1667432, -1104333, -260646, -3833893, -2939036, -2235985, -420899, -2286327, 183443,
  -976891, 1612842, -3545687, -554416, 3919660, -48306, -1362209, 3937738, 1400424,
  -846154, 1976782,
]

///| Forward Number Theoretic Transform (NTT)
/// 正向数论变换 (NTT)
///
/// Performs in-place Cooley-Tukey NTT using Montgomery arithmetic.
/// 使用蒙哥马利算术执行原地 Cooley-Tukey NTT。
/// 
/// Algorithm details / 算法细节:
/// - Decimation-in-frequency (DIF) butterfly operations / 频域抽取蝶形运算
/// - No modular reduction after additions/subtractions / 加减法后不进行模约简
/// - Output is in bit-reversed order / 输出为比特翻转顺序
/// - All coefficients remain in Montgomery domain / 所有系数保持在蒙哥马利域中
///
///
/// Arguments / 参数:
/// - a: input/output coefficient array (length 256) / 输入/输出系数数组（长度256）
pub fn ntt(a : Array[Int]) -> Unit {
  let mut k = 0 // Index for twiddle factors / 旋转因子索引

  // Iterate through all NTT levels (log2(N) = 8 levels)
  // 遍历所有 NTT 层级（log2(N) = 8 层）
  for len = 128; len > 0; len = len >> 1 {
    // Process all butterfly groups at current level
    // 处理当前层级的所有蝶形运算组
    for start = 0; start < N; start = start + (len << 1) {
      k = k + 1
      let zeta = zetas[k].to_int64() // Current twiddle factor / 当前旋转因子

      // Perform butterfly operations within current group
      // 在当前组内执行蝶形运算
      for j = start; j < start + len; j = j + 1 {
        // Butterfly operation: (a[j], a[j+len]) -> (a[j]+t, a[j]-t)
        // 蝶形运算：(a[j], a[j+len]) -> (a[j]+t, a[j]-t)
        // where t = zeta * a[j+len] (Montgomery multiplication)
        // 其中 t = zeta * a[j+len] (蒙哥马利乘法)
        let t = montgomery_reduce(zeta * a[j + len].to_int64())
        a[j + len] = a[j] - t
        a[j] = a[j] + t
      }
    }
  }
}

///| Inverse Number Theoretic Transform (INTT)
/// 逆向数论变换 (INTT)
///
/// Performs in-place inverse NTT with Montgomery factor multiplication.
/// 执行原地逆向 NTT 并乘以蒙哥马利因子。
///
/// Algorithm details / 算法细节:
/// - Decimation-in-time (DIT) butterfly operations / 时域抽取蝶形运算
/// - Uses negative twiddle factors / 使用负旋转因子
/// - Final multiplication by F = (R^2/N) mod Q / 最终乘以 F = (R^2/N) mod Q
/// - No modular reduction during butterfly operations / 蝶形运算期间不进行模约简
///
/// The factor F serves three purposes / 因子 F 有三个作用:
/// 1. Divide by N to normalize the transform / 除以 N 来归一化变换
/// 2. Multiply by R to keep result in Montgomery domain / 乘以 R 保持结果在蒙哥马利域
/// 3. Multiply by R again to compensate for Montgomery reduction / 再乘以 R 补偿蒙哥马利约简
///
/// Complexity / 复杂度: O(N log N) where N = 256
///
/// Arguments / 参数:
/// - a: input/output coefficient array (length 256) / 输入/输出系数数组（长度256）
///
/// Preconditions / 前置条件:
/// - Input is in bit-reversed order (output of NTT) / 输入为比特翻转顺序（NTT的输出）
/// - Input coefficients < Q in absolute value / 输入系数绝对值 < Q
///
/// Postconditions / 后置条件:
/// - Output is in standard order / 输出为标准顺序
/// - Output coefficients < Q in absolute value / 输出系数绝对值 < Q
/// - Result remains in Montgomery domain / 结果保持在蒙哥马利域中
pub fn inv_ntt_to_mont(a : Array[Int]) -> Unit {
  let mut k = 256 // Start from end of twiddle factors array / 从旋转因子数组末尾开始

  // Iterate through all INTT levels in reverse order
  // 以相反顺序遍历所有 INTT 层级
  for len = 1; len < N; len = len << 1 {
    // Process all butterfly groups at current level
    // 处理当前层级的所有蝶形运算组
    for start = 0; start < N; start = start + (len << 1) {
      k -= 1
      let zeta = -zetas[k].to_int64() // Negative twiddle factor for INTT / INTT 使用负旋转因子

      // Perform inverse butterfly operations within current group
      // 在当前组内执行逆向蝶形运算
      for j = start; j < start + len; j = j + 1 {
        let t = a[j]
        a[j] = t + a[j + len] // Addition / 加法
        a[j + len] = t - a[j + len] // Subtraction / 减法
        // Montgomery multiplication with negative twiddle factor
        // 与负旋转因子的蒙哥马利乘法
        a[j + len] = montgomery_reduce(zeta * a[j + len].to_int64())
      }
    }
  }

  // Final normalization and Montgomery domain maintenance
  // 最终归一化并保持蒙哥马利域
  for i in 0..