rootofunity: Allow powers of primitive roots of unity.

* src/mpc.h: Change the prototype of rootofunity.
* doc/mpc.texi: Document the change.
* src/rootofunity.c [gcd]: New function.
  [rootofunity]: Adapt the algorithm. The special cases of n are not yet
    handled.
    Correct a mistake in the error analysis, leading to wrong results for
    large n.
* tests/trootofunity.c, tests/rootofunity.dsc, src/rootofunity.c: Adapt the
    tests.

6 files changed, 57 insertions, 22 deletions

diff --git a/doc/mpc.texi b/doc/mpc.texi
index aa629a1..7e77223 100644
--- a/doc/mpc.texi
+++ b/doc/mpc.texi
@@ -981,9 +981,9 @@ so that the imaginary part of the result lies in
 @math{]-\pi , \pi]} and @math{]-\pi/log(10) , \pi/log(10)]} respectively.
 @end deftypefun
 
-@deftypefun int mpc_rootofunity (mpc_t @var{rop}, unsigned long int @var{n}, mpc_rnd_t @var{rnd})
-Set @var{rop} to the standard primitive @var{n}-th root of unity,
-@math{\exp (2 \pi i / n)},
+@deftypefun int mpc_rootofunity (mpc_t @var{rop}, unsigned long int @var{n}, unsigned long int @var{k}, mpc_rnd_t @var{rnd})
+Set @var{rop} to the standard primitive @var{n}-th root of unity raised to the power @var{k}, that is,
+@math{\exp (2 \pi i k / n)},
 rounded according to @var{rnd} with the precision of @var{rop}.
 @end deftypefun
 
diff --git a/src/mpc.h b/src/mpc.h
index ead33d0..e7a2142 100644
--- a/src/mpc.h
+++ b/src/mpc.h
@@ -200,7 +200,7 @@ __MPC_DECLSPEC int  mpc_atan        (mpc_ptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_asinh       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_acosh       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_atanh       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_rootofunity (mpc_ptr, unsigned long int, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_rootofunity (mpc_ptr, unsigned long int, unsigned long int, mpc_rnd_t);
 __MPC_DECLSPEC void mpc_clear       (mpc_ptr);
 __MPC_DECLSPEC int  mpc_urandom     (mpc_ptr, gmp_randstate_t);
 __MPC_DECLSPEC void mpc_init2       (mpc_ptr, mpfr_prec_t);
diff --git a/src/rootofunity.c b/src/rootofunity.c
index bb3039b..a1b73da 100644
--- a/src/rootofunity.c
+++ b/src/rootofunity.c
@@ -19,11 +19,22 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 */
 
 #include "mpc-impl.h"
+#include <assert.h>
 
-/* put in rop the value of exp(2*i*pi/n) rounded according to rnd */
+static unsigned long int
+gcd (unsigned long int a, unsigned long int b)
+{
+   if (b == 0)
+      return a;
+   else return gcd (b, a % b);
+}
+
+/* put in rop the value of exp(2*i*pi*k/n) rounded according to rnd */
 int
-mpc_rootofunity (mpc_ptr rop, unsigned long int n, mpc_rnd_t rnd)
+mpc_rootofunity (mpc_ptr rop, unsigned long int n, unsigned long int k, mpc_rnd_t rnd)
 {
+   unsigned long int g;
+   mpq_t kn;
    mpfr_t t, s, c;
    mpfr_prec_t prec;
    int inex_re, inex_im;
@@ -35,8 +46,15 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, mpc_rnd_t rnd)
       return MPC_INEX (0, 0);
    }
 
-   /* We assume that only n=1, 2, 3, 4, 6, 8, 12 yield exact results and
-      treat them separately.                                          */
+   /* Eliminate common denominator. */
+   k %= n;
+   g = gcd (k, n);
+   k /= g;
+   n /= g;
+
+   /* We assume that only n=1, 2, 3, 4, 6 and 12 may yield exact results
+      and treat them separately; n=8 is also treated here for efficiency
+      reasons. */
    if (n == 1)
       return mpc_set_ui_ui (rop, 1, 0, rnd);
    else if (n == 2)
@@ -65,9 +83,12 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, mpc_rnd_t rnd)
 
    prec = MPC_MAX_PREC(rop);
 
-   mpfr_init2 (t, 2);
-   mpfr_init2 (s, 2);
-   mpfr_init2 (c, 2);
+   mpfr_init2 (t, 67); /* see the argument at the end of the following loop */
+   mpfr_init2 (s, 67);
+   mpfr_init2 (c, 67);
+   mpq_init (kn);
+   mpq_set_ui (kn, k, n);
+   mpq_mul_2exp (kn, kn, 1);
 
    do {
       prec += mpc_ceil_log2 (prec) + 5;
@@ -77,15 +98,27 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, mpc_rnd_t rnd)
       mpfr_set_prec (c, prec);
 
       mpfr_const_pi (t, MPFR_RNDN); /* error 0.5 ulp */
-      mpfr_mul_2ui (t, t, 1u, MPFR_RNDN);
-      mpfr_div_ui (t, t, n, MPFR_RNDN); /* error 2*0.5+0.5=1.5 ulp */
+      mpfr_mul_q (t, t, kn, MPFR_RNDN); /* error 2*1.5+0.5=3.5 ulp */
       mpfr_sin_cos (s, c, t, MPFR_RNDN);
-         /* error (1.5*2^{Exp (t) - Exp (s resp.c)} + 0.5) ulp
-            <= 6.5 ulp for n>=3                             */
+         /* error (1.5*2^{Exp (t) - Exp (s resp. c)} + 0.5) ulp
+            We have 0<t<2*pi, so Exp (t) <= 3.
+            Unfortunately s or c can be close to 0.
+            Where the minimum of s and c over all primitive n-th roots of
+            unity is reached depends on n mod 4. 
+            To simplify the argument, we will consider the 4*n-th roots of
+            unity, which contain the n-th roots of unity and which are
+            symmmetrically distributed with respect to the real and imaginary
+            axes, so that it becomes enough to consider only s for k=1.
+            With n<2^64, the absolute values of all s or c are at least
+            sin (2 * pi * 2^(-64) / 4) > 2^(-64) of exponent at least -63.
+            So the error is bounded above by
+            (1.5*2^66+0.5) ulp < 2^67 ulp.
+            (This could be made dependent on n, which would be useful for
+            small n at small precision.) */
    }
-   while (   !mpfr_can_round (c, prec - 3, MPFR_RNDN, MPFR_RNDZ,
+   while (   !mpfr_can_round (c, prec - 67, MPFR_RNDN, MPFR_RNDZ,
                  MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == MPFR_RNDN))
-          || !mpfr_can_round (s, prec - 3, MPFR_RNDN, GMP_RNDZ,
+          || !mpfr_can_round (s, prec - 67, MPFR_RNDN, MPFR_RNDZ,
                  MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == MPFR_RNDN)));
 
    inex_re = mpfr_set (mpc_realref(rop), c, MPC_RND_RE(rnd));
@@ -94,6 +127,7 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, mpc_rnd_t rnd)
    mpfr_clear (t);
    mpfr_clear (s);
    mpfr_clear (c);
+   mpq_clear (kn);
 
    return MPC_INEX(inex_re, inex_im);
 }
diff --git a/tests/rootofunity.dat b/tests/rootofunity.dat
index f61848f..3742dd4 100644
--- a/tests/rootofunity.dat
+++ b/tests/rootofunity.dat
@@ -20,7 +20,8 @@
 # The line format respects the parameter order in function prototype as
 # follow:
 #
-# INEX_RE  INEX_IM  PREC_ROP_RE  ROP_RE  PREC_ROP_IM  ROP_IM  OP  RND_RE  RND_IM
+# INEX_RE  INEX_IM  PREC_ROP_RE  ROP_RE  PREC_ROP_IM  ROP_IM  OP1 OP2  RND_RE  RND_IM
 
 # special values
-0 0 53 nan 53 nan 0 N N
+0 0 53 nan 53 nan 0 0 N N
+0 0 53 nan 53 nan 0 1 N N
diff --git a/tests/rootofunity.dsc b/tests/rootofunity.dsc
index 613cdd9..87c0ddc 100644
--- a/tests/rootofunity.dsc
+++ b/tests/rootofunity.dsc
@@ -25,4 +25,5 @@ OUTPUT:
     mpc_ptr
 INPUT:
     unsigned long int
+    unsigned long int
     mpc_rnd_t
diff --git a/tests/trootofunity.c b/tests/trootofunity.c
index 7bb77d6..9624873 100644
--- a/tests/trootofunity.c
+++ b/tests/trootofunity.c
@@ -21,7 +21,7 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 #include "mpc-tests.h"
 
 #define MPC_FUNCTION_CALL                                               \
-  P[0].mpc_inex = mpc_rootofunity (P[1].mpc, P[2].ui, P[3].mpc_rnd)
+  P[0].mpc_inex = mpc_rootofunity (P[1].mpc, P[2].ui, P[3].ui, P[4].mpc_rnd)
 
 #include "data_check.tpl"
 #include "tgeneric.tpl"
@@ -41,7 +41,7 @@ check (unsigned long int n)
       mpc_set_prec (z, prec);
       mpc_set_prec (zero, prec);
 
-      mpc_rootofunity (z, n, MPC_RNDNN);
+      mpc_rootofunity (z, n, 1, MPC_RNDNN);
       mpc_pow_ui (zero, z, n, MPC_RNDNN);
       mpc_sub_ui (zero, zero, 1u, MPC_RNDNN);
       if (MPC_MAX (mpfr_get_exp (mpc_realref (zero)), mpfr_get_exp (mpc_imagref (zero)))
@@ -70,7 +70,6 @@ main (void)
 
    data_check_template ("rootofunity.dsc", "rootofunity.dat");
 
-
    /* Avoid checking roots of unity of high order at very low precision,
       so start only at 20. */
    tgeneric_template ("rootofunity.dsc", 20, 512, 7, 256);