gcd_coefficients()

TowerPolynomial DPoly::gcd_coefficients	(	int	level,
		const TowerPolynomial	f,
		const TowerPolynomial	g,
		TowerPolynomial &	result_u,
		TowerPolynomial &	result_v )

Definition at line 1159 of file dpoly.cpp.

{
  // Assumption:
  //  f and g are non-zero
  TowerPolynomial v1, v2, v3;
  TowerPolynomial u1, u2, u3;
  TowerPolynomial q = nullptr;
 
  v1 = nullptr;
  v2 = from_long(level, 1);
  v3 = copy(level, g);
 
  u1 = from_long(level, 1);
  u2 = nullptr;
  u3 = copy(level, f);
 
  if (v3 == nullptr || (u3 != nullptr && v3->deg > u3->deg))
    {
      swap_poly(u1, v1);
      swap_poly(u2, v2);
      swap_poly(u3, v3);
    }
 
// At the end of the loop:
// u1 * f + u2 * g == u3
// v1 * f + v2 * g == v3
// (and (v1,v2,v3) is close to the gcd
 
#ifdef DEBUGGCD
  if (level == 1)
    {
      printf("u1 = %s\n", to_string(level, u1));
      printf("u2 = %s\n", to_string(level, u2));
      printf("u3 = %s\n", to_string(level, u3));
 
      printf("v1 = %s\n", to_string(level, v1));
      printf("v2 = %s\n", to_string(level, v2));
      printf("v3 = %s\n", to_string(level, v3));
    }
#endif
 
  while (v3 != nullptr)
    {
      if (!make_monic3(level, v1, v2, v3))
        {
          // deallocate some polynomials, then return 0.  No monic gcd
          dealloc_poly(q);
          dealloc_poly(u1);
          dealloc_poly(u2);
          dealloc_poly(u3);
          dealloc_poly(v1);
          dealloc_poly(v2);
          dealloc_poly(v3);
          result_u = nullptr;
          result_v = nullptr;
          return nullptr;
        }
      q = division_in_place_monic(
          level,
          u3,
          v3);  // u3 := u3 - q*v3, as v3 is monic, this is always possible
 
#ifdef DEBUGGCD
      if (level == 1)
        {
          printf("q = %s\n", to_string(level, q));
          printf("u3 = %s\n", to_string(level, u3));
        }
#endif
 
      negate_in_place(level, q);
      TowerPolynomial a = mult(level, q, v1, false);
      TowerPolynomial b = mult(level, q, v2, false);
      add_in_place(level, u1, a);
      add_in_place(level, u2, b);
 
#ifdef DEBUGGCD
      if (level == 1)
        {
          printf("u1 = %s\n", to_string(level, u1));
          printf("u2 = %s\n", to_string(level, u2));
          printf("u3 = %s\n", to_string(level, u3));
        }
#endif
      dealloc_poly(a);  // MES: totally wipeout polys a, b here!
      dealloc_poly(b);
      swap_poly(u1, v1);
      swap_poly(u2, v2);
      swap_poly(u3, v3);
    }
 
#ifdef DEBUGGCD
  if (level == 1)
    {
      printf("u1 = %s\n", to_string(level, u1));
      printf("u2 = %s\n", to_string(level, u2));
      printf("u3 = %s\n", to_string(level, u3));
    }
#endif
 
  // At this point u3 is monic, and is the monic gcd
 
  // v3 is already 0 here.
  dealloc_poly(v1);
  dealloc_poly(v2);
 
  result_u = u1;
  result_v = u2;
  return u3;
}

References add_in_place(), copy(), dealloc_poly(), division_in_place_monic(), from_long(), make_monic3(), mult(), negate_in_place(), swap_poly(), and to_string().

Referenced by invert().