rootofunity: Handle powers of roots of order n = 3, 4, 6, 8, 12.

* src/rootofunity.c [mpc_rootofunity]: Handle powers for special n.
* tests/rootofunity.dat: Add corresponding tests, in particular with
  directed rounding.

2 files changed, 78 insertions, 10 deletions

diff --git a/src/rootofunity.c b/src/rootofunity.c
index a1b73da..84697c7 100644
--- a/src/rootofunity.c
+++ b/src/rootofunity.c
@@ -38,6 +38,7 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, unsigned long int k, mpc_rnd_
    mpfr_t t, s, c;
    mpfr_prec_t prec;
    int inex_re, inex_im;
+   mpfr_rnd_t rnd_re, rnd_im;
 
    if (n == 0) {
       /* Compute exp (0 + i*inf). */
@@ -60,24 +61,58 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, unsigned long int k, mpc_rnd_
    else if (n == 2)
       return mpc_set_si_si (rop, -1, 0, rnd);
    else if (n == 4)
-      return mpc_set_ui_ui (rop, 0, 1, rnd);
+      if (k == 1)
+         return mpc_set_ui_ui (rop, 0, 1, rnd);
+      else
+         return mpc_set_si_si (rop, 0, -1, rnd);
    else if (n == 3 || n == 6) {
       inex_re = mpfr_set_si (mpc_realref (rop), (n == 3 ? -1 : 1),
                              MPC_RND_RE (rnd));
-      inex_im = mpfr_sqrt_ui (mpc_imagref (rop), 3, MPC_RND_IM (rnd));
+      /* inverse the rounding mode for -sqrt(3)/2 for zeta_3^2 and zeta_6^5 */
+      rnd_im = MPC_RND_IM (rnd);
+      if (k != 1)
+         rnd_im = INV_RND (rnd_im);
+      inex_im = mpfr_sqrt_ui (mpc_imagref (rop), 3, rnd_im);
       mpc_div_2ui (rop, rop, 1, MPC_RNDNN);
+      if (k != 1) {
+         mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN);
+         inex_im = -inex_im;
+      }
       return MPC_INEX (inex_re, inex_im);
    }
-   else if (n == 8) {
-      inex_re = mpfr_sqrt_ui (mpc_realref (rop), 2, MPC_RND_RE (rnd));
-      inex_im = mpfr_sqrt_ui (mpc_imagref (rop), 2, MPC_RND_IM (rnd));
+   else if (n == 12) {
+      /* inverse the rounding mode for -sqrt(3)/2 for zeta_12^5 and zeta_12^7 */
+      rnd_re = MPC_RND_RE (rnd);
+      if (k == 5 || k == 7)
+         rnd_re = INV_RND (rnd_re);
+      inex_re = mpfr_sqrt_ui (mpc_realref (rop), 3, rnd_re);
+      inex_im = mpfr_set_si (mpc_imagref (rop), (k == 1 || k == 5 ? 1 : -1),
+         MPC_RND_IM (rnd));
       mpc_div_2ui (rop, rop, 1u, MPC_RNDNN);
+      if (k == 5 || k == 7) {
+         mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN);
+         inex_re = -inex_re;
+      }
       return MPC_INEX (inex_re, inex_im);
    }
-   else if (n == 12) {
-      inex_re = mpfr_sqrt_ui (mpc_realref (rop), 3, MPC_RND_RE (rnd));
-      inex_im = mpfr_set_ui (mpc_imagref (rop), 1, MPC_RND_IM (rnd));
+   else if (n == 8) {
+      rnd_re = MPC_RND_RE (rnd);
+      if (k == 3 || k == 5)
+         rnd_re = INV_RND (rnd_re);
+      rnd_im = MPC_RND_IM (rnd);
+      if (k == 5 || k == 7)
+         rnd_im = INV_RND (rnd_im);
+      inex_re = mpfr_sqrt_ui (mpc_realref (rop), 2, rnd_re);
+      inex_im = mpfr_sqrt_ui (mpc_imagref (rop), 2, rnd_im);
       mpc_div_2ui (rop, rop, 1u, MPC_RNDNN);
+      if (k == 3 || k == 5) {
+         mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN);
+         inex_re = -inex_re;
+      }
+      if (k == 5 || k == 7) {
+         mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN);
+         inex_im = -inex_im;
+      }
       return MPC_INEX (inex_re, inex_im);
    }
 
diff --git a/tests/rootofunity.dat b/tests/rootofunity.dat
index 3742dd4..3ee9e86 100644
--- a/tests/rootofunity.dat
+++ b/tests/rootofunity.dat
@@ -23,5 +23,38 @@
 # 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 0 N N
-0 0 53 nan 53 nan 0 1 N N
+0 0 53 nan 53 nan  0 0 N N
+0 0 53 nan 53 nan  0 1 N N
+
+# roots of unity of small order, some of them with exact results
+0 0 53 1   53 0   1 1 Z U
+0 0 53 -1  53 0   2 1 N N
+0 0 53 0   53 1   4 1 U Z
+0 0 53 -1  53 0   4 2 U Z
+0 0 53 0   53 -1  4 3 U Z
+0 0 53 1   53 0   4 4 U Z
+0 - 53 -0.5 53  0b1.1011101101100111101011101000010110000100110010101010e-1  3 1 N D
+0 + 53 -0.5 53  0b1.1011101101100111101011101000010110000100110010101011e-1  3 1 N U
+0 - 53 -0.5 53 -0b1.1011101101100111101011101000010110000100110010101011e-1  3 2 N D
+0 + 53 -0.5 53 -0b1.1011101101100111101011101000010110000100110010101010e-1  3 2 N U
+0 - 53  0.5 53  0b1.1011101101100111101011101000010110000100110010101010e-1  6 1 D Z
+0 - 53  0.5 53 -0b1.1011101101100111101011101000010110000100110010101011e-1  6 5 D D
+0 + 53  0.5 53 -0b1.1011101101100111101011101000010110000100110010101010e-1  6 5 D U
+- 0 53  0b1.1011101101100111101011101000010110000100110010101010e-1 53  0.5  12  1 N Z
+- 0 53  0b1.1011101101100111101011101000010110000100110010101010e-1 53  0.5  12  1 D Z
++ 0 53  0b1.1011101101100111101011101000010110000100110010101011e-1 53  0.5  12  1 U Z
+- 0 53 -0b1.1011101101100111101011101000010110000100110010101011e-1 53  0.5  12  5 D Z
++ 0 53 -0b1.1011101101100111101011101000010110000100110010101010e-1 53  0.5  12  5 U Z
+- 0 53 -0b1.1011101101100111101011101000010110000100110010101011e-1 53 -0.5  12  7 D Z
++ 0 53 -0b1.1011101101100111101011101000010110000100110010101010e-1 53 -0.5  12  7 U Z
+- 0 53  0b1.1011101101100111101011101000010110000100110010101010e-1 53 -0.5  12 11 D Z
++ 0 53  0b1.1011101101100111101011101000010110000100110010101011e-1 53 -0.5  12 11 U Z
++ + 53  0b1.0110101000001001111001100110011111110011101111001101e-1 53  0b1.0110101000001001111001100110011111110011101111001101e-1  8 1 N N
++ - 53  0b1.0110101000001001111001100110011111110011101111001101e-1 53  0b1.0110101000001001111001100110011111110011101111001100e-1  8 1 U D
+- + 53  0b1.0110101000001001111001100110011111110011101111001100e-1 53  0b1.0110101000001001111001100110011111110011101111001101e-1  8 1 D U
++ - 53 -0b1.0110101000001001111001100110011111110011101111001100e-1 53  0b1.0110101000001001111001100110011111110011101111001100e-1  8 3 U D
+- + 53 -0b1.0110101000001001111001100110011111110011101111001101e-1 53  0b1.0110101000001001111001100110011111110011101111001101e-1  8 3 D U
++ - 53 -0b1.0110101000001001111001100110011111110011101111001100e-1 53 -0b1.0110101000001001111001100110011111110011101111001101e-1  8 5 U D
+- + 53 -0b1.0110101000001001111001100110011111110011101111001101e-1 53 -0b1.0110101000001001111001100110011111110011101111001100e-1  8 5 D U
++ - 53  0b1.0110101000001001111001100110011111110011101111001101e-1 53 -0b1.0110101000001001111001100110011111110011101111001101e-1  8 7 U D
+- + 53  0b1.0110101000001001111001100110011111110011101111001100e-1 53 -0b1.0110101000001001111001100110011111110011101111001100e-1  8 7 D U