Add positive root function of two algebraic numbers

include/algebraic_number.h

@@ -48,7 +48,7 @@ void lp_algebraic_number_construct(lp_algebraic_number_t* a, lp_upolynomial_t* f
 /** Construct a zero algebraic number */
 void lp_algebraic_number_construct_zero(lp_algebraic_number_t* a);
 
-/** Construct a zero algebraic number */
+/** Construct a one algebraic number */
 void lp_algebraic_number_construct_one(lp_algebraic_number_t* a);
 
 /** Construct a copy of the algebraic number. */
@@ -142,6 +142,11 @@ void lp_algebraic_number_div(lp_algebraic_number_t* div, const lp_algebraic_numb
 /** Exponentiation */
 void lp_algebraic_number_pow(lp_algebraic_number_t* pow, const lp_algebraic_number_t* a, unsigned n);
 
+/** Square, cubic, ... roots positive
+    The function suppose that the algebraic number a is positive.
+ */
+void lp_algebraic_number_positive_root(lp_algebraic_number_t* root, const lp_algebraic_number_t* a, unsigned n);
+
 /** Returns true if a rational number (not complete) */
 int lp_algebraic_number_is_rational(const lp_algebraic_number_t* a);
 
@@ -154,7 +159,7 @@ void lp_algebraic_number_ceiling(const lp_algebraic_number_t* a, lp_integer_t* a
 /** Returns the floor of the number */
 void lp_algebraic_number_floor(const lp_algebraic_number_t* a, lp_integer_t* a_floor);
 
-/** Returns a hash of the of the dyadic approximation of a */
+/** Returns a hash of the dyadic approximation of a */
 size_t lp_algebraic_number_hash_approx(const lp_algebraic_number_t* a, unsigned precision);
 
 #ifdef __cplusplus

include/upolynomial.h

@@ -173,6 +173,12 @@ int lp_upolynomial_cmp(const lp_upolynomial_t* p, const lp_upolynomial_t* q);
  */
 lp_upolynomial_t* lp_upolynomial_subst_x_neg(const lp_upolynomial_t* f);
 
+/**
+ * Replace the polynomial with f(x**n).
+ * The function suppose that the maximum degree * n doesn't overflow.
+ */
+void lp_upolynomial_subst_x_pow_in_place(lp_upolynomial_t* f, size_t n);
+
 /**
  * Returns the polynomial -f.
  */

python/polypyAlgebraicNumber.h

@@ -29,10 +29,10 @@ typedef struct {
   lp_algebraic_number_t a;
 } AlgebraicNumber;
 
-/** Methods on coefficient rings */
+/** Methods on algebraic numbers */
 extern PyMethodDef AlgebraicNumber_methods[];
 
-/** Definition of the CoefficientRing type */
+/** Definition of the AlgebraicNumber type */
 extern PyTypeObject AlgebraicNumberType;
 
 /** Create an algebraic object (makes a copy of a) */

python/polypyAlgebraicNumber2.c

@@ -65,9 +65,13 @@ AlgebraicNumber_div(PyObject* self, PyObject* args);
 static PyObject*
 AlgebraicNumber_pow(PyObject* self, PyObject* args);
 
+static PyObject*
+AlgebraicNumber_positive_root(PyObject* self, PyObject* args);
+
 PyMethodDef AlgebraicNumber_methods[] = {
     {"refine", (PyCFunction)AlgebraicNumber_refine, METH_NOARGS, "Refines the number to half the interval"},
     {"to_double", (PyCFunction)AlgebraicNumber_to_double, METH_NOARGS, "Returns the approximation of the algebraic number"},
+    {"positive_root", (PyCFunction)AlgebraicNumber_positive_root, METH_VARARGS, "Returns the positive root of the number is positive"},
     {NULL}  /* Sentinel */
 };
 
@@ -446,3 +450,29 @@ AlgebraicNumber_pow(PyObject* self, PyObject* other) {
 
   return result;
 }
+
+static PyObject*
+AlgebraicNumber_positive_root(PyObject* self, PyObject* args) {
+  UPolynomialObject* p = (UPolynomialObject*) self;
+  if (p) {
+    // Get the argument
+    if (PyTuple_Size(args) == 1) {
+      // Get n
+      PyObject* other = PyTuple_GetItem(args, 0);
+      AlgebraicNumber* a1 = (AlgebraicNumber*) self;
+      long n = PyLong_AsLong(other);
+
+      lp_algebraic_number_t root;
+      lp_algebraic_number_construct_zero(&root);
+      lp_algebraic_number_positive_root(&root, &a1->a, n);
+      PyObject* result = PyAlgebraicNumber_create(&root);
+      lp_algebraic_number_destruct(&root);
+
+      return result;
+    } else {
+      Py_RETURN_NONE;
+    }
+  } else {
+    Py_RETURN_NONE;
+  }
+}

python/polypyAlgebraicNumber3.c

@@ -62,9 +62,13 @@ AlgebraicNumber_div(PyObject* self, PyObject* args);
 static PyObject*
 AlgebraicNumber_pow(PyObject* self, PyObject* args);
 
+static PyObject*
+AlgebraicNumber_positive_root(PyObject* self, PyObject* args);
+
 PyMethodDef AlgebraicNumber_methods[] = {
     {"refine", (PyCFunction)AlgebraicNumber_refine, METH_NOARGS, "Refines the number to half the interval"},
     {"to_double", (PyCFunction)AlgebraicNumber_to_double, METH_NOARGS, "Returns the approximation of the algebraic number"},
+    {"positive_root", (PyCFunction)AlgebraicNumber_positive_root, METH_VARARGS, "Returns the positive root of the number is positive"},
     {NULL}  /* Sentinel */
 };
 
@@ -429,3 +433,29 @@ AlgebraicNumber_pow(PyObject* self, PyObject* other) {
 
   return result;
 }
+
+static PyObject*
+AlgebraicNumber_positive_root(PyObject* self, PyObject* args) {
+  UPolynomialObject* p = (UPolynomialObject*) self;
+  if (p) {
+    // Get the argument
+    if (PyTuple_Size(args) == 1) {
+      // Get n
+      PyObject* other = PyTuple_GetItem(args, 0);
+      AlgebraicNumber* a1 = (AlgebraicNumber*) self;
+      long n = PyLong_AsLong(other);
+
+      lp_algebraic_number_t root;
+      lp_algebraic_number_construct_zero(&root);
+      lp_algebraic_number_positive_root(&root, &a1->a, n);
+      PyObject* result = PyAlgebraicNumber_create(&root);
+      lp_algebraic_number_destruct(&root);
+
+      return result;
+    } else {
+      Py_RETURN_NONE;
+    }
+  } else {
+    Py_RETURN_NONE;
+  }
+}

src/interval/arithmetic.c

@@ -1254,3 +1254,34 @@ void dyadic_interval_pow(lp_dyadic_interval_t* P, const lp_dyadic_interval_t* I,
     }
   }
 }
+
+void dyadic_interval_root_overapprox(lp_dyadic_interval_t* P, const lp_dyadic_interval_t* I, unsigned n, unsigned prec) {
+  assert ( n > 0 );
+  if (n == 1) {
+    // I^1 = I
+    lp_dyadic_interval_assign(P,I);
+  } else {
+    if(dyadic_rational_root_approx(&P->a,&I->a,n,prec,0/*floor*/)) {
+       // root of lower bound is exact
+       if (I->is_point) {
+         // Plain root
+         if (!P->is_point) {
+           dyadic_rational_destruct(&P->b);
+           P->is_point = 1;
+           P->a_open = P->b_open = 0;
+         }
+         return;
+       }
+    }
+    if (P->is_point) {
+      P->is_point = 0;
+      dyadic_rational_construct(&P->b);
+    }
+    if (I->is_point) {
+      dyadic_rational_root_approx(&P->b,&I->a,n,prec,1/*ceil*/);
+    } else {
+      dyadic_rational_root_approx(&P->b,&I->b,n,prec,1);
+    }
+    P->a_open = P->b_open = 0;
+  }
+}

src/interval/arithmetic.h

@@ -39,3 +39,9 @@ void dyadic_interval_sub(lp_dyadic_interval_t* S, const lp_dyadic_interval_t* I1
 void dyadic_interval_mul(lp_dyadic_interval_t* P, const lp_dyadic_interval_t* I1, const lp_dyadic_interval_t* I2);
 
 void dyadic_interval_pow(lp_dyadic_interval_t* P, const lp_dyadic_interval_t* I, unsigned n);
+
+/** Compute an overapproximation of the n th positive root of the reals in
+    I. The precision of the approximation is at least prec (roughly the
+    overapproximation of the result is at least 2^(-prec/n) ) .
+ */
+void dyadic_interval_root_overapprox(lp_dyadic_interval_t* P, const lp_dyadic_interval_t* I, unsigned n, unsigned prec);

src/number/algebraic_number.c

@@ -571,6 +571,19 @@ void lp_algebraic_number_to_rational(const lp_algebraic_number_t* a_const, lp_ra
     return;
   }
 
+  if (lp_upolynomial_degree(a_const->f) == 1) {
+    // If degree 1, we're directly rational
+    if(a_const->f->size==2){
+      mpz_neg(&q->_mp_num,&a_const->f->monomials[0].coefficient);
+      mpz_set(&q->_mp_den,&a_const->f->monomials[1].coefficient);
+    } else {
+      /* a_const == 0, is it possible without a_const->f==0? */
+      assert (a_const->f->size<=1);
+      mpq_set_ui(q,0,1);
+    }
+    return;
+  }
+
   // We do the necessary refinement on a copy
   lp_algebraic_number_t a;
   lp_algebraic_number_construct_copy(&a, a_const);
@@ -1144,6 +1157,116 @@ void lp_algebraic_number_pow(lp_algebraic_number_t* pow, const lp_algebraic_numb
   }
 }
 
+/**
+ * Compute the positive root of an algebraic number.
+ *
+ * f(x) = a_p x^p + ... + a_0
+ *
+ * with f(a) = 0, a only root in (l, u), (l, u) positive
+ *
+ * NOTE: l, u could be 0
+ *
+ * g(x) = f(x^n)
+ *
+ * with g(root(a,n)) = 0, only root in 0 < root(l,n) < root(a,n) < root(u,n)
+ *
+ * but, root(l,n) and root(u,n) are not dyadic rationals. So we overapproximate
+ * the resulting interval and refine the overapproximation until there is a
+ * unique root in g in the resulting interval.
+ */
+void lp_algebraic_number_positive_root(
+    lp_algebraic_number_t* op, const lp_algebraic_number_t* a, unsigned n)
+{
+  assert(0 < n);
+  assert(lp_algebraic_number_sgn(a) >= 0);
+
+  if (trace_is_enabled("algebraic_number")) {
+    tracef("a = "); lp_algebraic_number_print(a, trace_out); tracef("\n");
+    tracef("root = %d",n); tracef("\n");
+  }
+
+  lp_upolynomial_t* f;
+
+  if (a->f) {
+    f = lp_upolynomial_construct_copy(a->f);
+  } else {
+    assert(a->I.is_point);
+    // x = p/q -> q*x - p = 0
+    lp_integer_t coefs[2];
+    integer_construct(&coefs[0]);
+    integer_construct(&coefs[1]);
+    integer_neg(lp_Z, &coefs[0], &a->I.a.a);
+    dyadic_rational_get_den(&a->I.a, &coefs[1]);
+    f = lp_upolynomial_construct(lp_Z,1,coefs);
+    lp_integer_destruct(&coefs[0]);
+    lp_integer_destruct(&coefs[1]);
+  }
+
+  lp_upolynomial_subst_x_pow_in_place(f,n);
+
+  // Get the roots of f
+  size_t f_roots_size = 0;
+  lp_algebraic_number_t* f_roots = malloc(sizeof(lp_algebraic_number_t)*lp_upolynomial_degree(f));
+  lp_upolynomial_roots_isolate(f, f_roots, &f_roots_size);
+  lp_upolynomial_delete(f);
+
+  // Interval for the result
+  lp_dyadic_interval_t I;
+  lp_dyadic_interval_construct_zero(&I);
+
+  unsigned long prec = 0;
+
+  // Approximate the result
+  while (f_roots_size > 1) {
+
+    // Root the interval
+    dyadic_interval_root_overapprox(&I, &a->I, n, prec);
+
+    if (trace_is_enabled("algebraic_number")) {
+      tracef("a = "); lp_algebraic_number_print(a, trace_out); tracef("\n");
+      tracef("I = "); lp_dyadic_interval_print(&I, trace_out); tracef("\n");
+      size_t i;
+      for (i = 0; i < f_roots_size; ++ i) {
+        tracef("f[%zu] = ", i); lp_algebraic_number_print(f_roots + i, trace_out); tracef("\n");
+      }
+    }
+
+    // Filter the roots over I
+    filter_roots(f_roots, &f_roots_size, &I);
+
+    // If more than one root, we need to refine a and b
+    if (f_roots_size > 1) {
+      if (a->f) {
+        lp_algebraic_number_refine_const_internal(a);
+      }
+      size_t i;
+      for (i = 0; i < f_roots_size; ++ i) {
+        if ((f_roots + i)->f) {
+          // not a point, refine it
+          lp_algebraic_number_refine_const_internal(f_roots + i);
+        }
+      }
+    }
+
+    prec += 1;
+  }
+
+  assert(f_roots_size == 1);
+
+  // Our root is the only one left
+  lp_algebraic_number_destruct(op);
+  *op = *f_roots;
+
+  if (trace_is_enabled("algebraic_number")) {
+    tracef("op = "); lp_algebraic_number_print(op, trace_out); tracef("\n");
+  }
+
+  // Remove temps
+  lp_dyadic_interval_destruct(&I);
+  free(f_roots);
+}
+
+
 void lp_algebraic_number_get_dyadic_midpoint(const lp_algebraic_number_t* a, lp_dyadic_rational_t* q) {
   if (a->I.is_point) {
     lp_dyadic_rational_assign(q, &a->I.a);

src/number/dyadic_rational.h

@@ -302,6 +302,36 @@ void dyadic_rational_pow(lp_dyadic_rational_t* pow, const lp_dyadic_rational_t*
   pow->n = a->n * n;
 }
 
+/* Approximate the positive root of a. The precision of the approximation is at
+   least prec (roughly the overapproximation of the result is at least
+   2^(-prec/n) ). The rounding is done toward plus infinity if ceil is true,
+   otherwise toward minus infinity. If the result is non-null the computation is exact */
+static inline
+int dyadic_rational_root_approx(lp_dyadic_rational_t* pow, const lp_dyadic_rational_t* a, unsigned long n, unsigned long prec, int ceil) {
+  assert(dyadic_rational_is_normalized(a));
+  assert(mpz_sgn(&a->a) >= 0);
+
+  if (mpz_sgn(&a->a) == 0){
+    /* pow = 0 = a */
+    mpz_set(&pow->a,&a->a);
+    pow->n=a->n;
+    return 1;
+  }
+
+  /* ensure denominator is a perfect root and bigger than prec */
+  unsigned long k = (a->n < prec? prec : a->n);
+  k += k%n;
+  pow->n = k/n;
+  mpz_mul_2exp(&pow->a,&a->a, k - a->n);
+
+  int exact = mpz_root(&pow->a,&pow->a,n);
+
+  if(ceil && !exact) mpz_add_ui(&pow->a,&pow->a,1);
+
+  dyadic_rational_normalize(pow);
+  return exact;
+}
+
 static inline
 void dyadic_rational_div_2exp(lp_dyadic_rational_t* div, const lp_dyadic_rational_t* a, unsigned long n) {
   assert(dyadic_rational_is_normalized(a));

src/upolynomial/upolynomial.c

@@ -1186,6 +1186,13 @@ void lp_upolynomial_reverse_in_place(lp_upolynomial_t* p) {
   }
 }
 
+void lp_upolynomial_subst_x_pow_in_place(lp_upolynomial_t* p, size_t n) {
+  size_t i;
+  for (i = 0; i < p->size; ++ i) {
+    p->monomials[i].degree *= n;
+  }
+}
+
 lp_upolynomial_t* lp_upolynomial_subst_x_neg(const lp_upolynomial_t* f) {
 
   lp_upolynomial_t* neg = lp_upolynomial_construct_copy(f);