gcd_poly.b

fun is_zero(A: [f64; D], D: u64) -> bool {
  for i: u64 = 0 .. D {
    if A[i] != 0. { return false; }
  }
  return true;
}

fun deg(A: [f64; DA], DA: u64) -> i64 {
  for decr i: i64 = ((i64) DA - 1) .. -1i64 {
    if A[i] != 0. { return i; }
  }
  return -1i64;
}

fun d_comb_lin(A: mut [f64; DA], B: [f64; DB], x: f64, d: i64, DA: u64, DB: u64) {
  let db: i64 = deg(B, DB);
  if db >= 0 {
    for i: u64 = 0 .. ((u64) db + 1) {
      A[d + (i64) i] = A[d + (i64) i] + x * B[i];
    }
  }
}

fun modulo(A: mut [f64; DA], B: [f64; DB], DA: u64, DB: u64) {
  let da: i64 = deg(A, DA);
  let db: i64 = deg(B, DB);
  if da >= 0 && db >= 0 {
    let k: f64 = 0.;
    for decr i: i64 = (da - db) .. -1i64 {
      if A[i + db] == 0. { continue; }
      k = - (A[i + db] / B[db]);
      d_comb_lin(A, B, k, i, DA, DB);
    }
  }
}

fun gcd(A: mut [f64; DA], B: mut [f64; DB], DA: u64, DB: u64) {
  if !is_zero(B, DB) {
    modulo(A, B, DA, DB);
    gcd(B, A, DB, DA);
  }
}