implemented rootofunity

git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/branches/rootsunity@1194 211d60ee-9f03-0410-a15a-8952a2c7a4e4

5 files changed, 193 insertions, 40 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 5cea5c2..a9bfa9c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -25,7 +25,8 @@ libmpc_la_SOURCES = mpc-impl.h abs.c acos.c acosh.c add.c add_fr.c         \
   fma.c fr_div.c fr_sub.c get_prec2.c get_prec.c get_version.c get_x.c     \
   imag.c init2.c init3.c inp_str.c log.c log10.c mem.c mul_2exp.c mul.c    \
   mul_fr.c mul_i.c mul_si.c mul_ui.c neg.c norm.c out_str.c pow.c pow_fr.c \
-  pow_ld.c pow_d.c pow_si.c pow_ui.c pow_z.c proj.c real.c urandom.c set.c \
+  pow_ld.c pow_d.c pow_si.c pow_ui.c pow_z.c proj.c real.c rootofunity.c   \
+  urandom.c set.c \
   set_prec.c set_str.c set_x.c set_x_x.c sin.c sin_cos.c sinh.c sqr.c      \
   sqrt.c strtoc.c sub.c sub_fr.c sub_ui.c swap.c tan.c tanh.c uceil_log2.c \
   ui_div.c ui_ui_sub.c
diff --git a/src/mpc.h b/src/mpc.h
index c98c3fb..179ca29 100644
--- a/src/mpc.h
+++ b/src/mpc.h
@@ -180,53 +180,54 @@ __MPC_DECLSPEC int  mpc_fma       (mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_srcptr,
 
 __MPC_DECLSPEC void mpc_set_nan   (mpc_ptr);
 
-__MPC_DECLSPEC int  mpc_real      (mpfr_ptr, mpc_srcptr, mpfr_rnd_t);
-__MPC_DECLSPEC int  mpc_imag      (mpfr_ptr, mpc_srcptr, mpfr_rnd_t);
-__MPC_DECLSPEC int  mpc_arg       (mpfr_ptr, mpc_srcptr, mpfr_rnd_t);
-__MPC_DECLSPEC int  mpc_proj      (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_cmp       (mpc_srcptr, mpc_srcptr);
-__MPC_DECLSPEC int  mpc_cmp_si_si (mpc_srcptr, long int, long int);
-__MPC_DECLSPEC int  mpc_exp       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_log       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_log10     (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_sin       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_cos       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_sin_cos   (mpc_ptr, mpc_ptr, mpc_srcptr, mpc_rnd_t, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_tan       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_sinh      (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_cosh      (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_tanh      (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_asin      (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_acos      (mpc_ptr, mpc_srcptr, mpc_rnd_t);
-__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 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);
-__MPC_DECLSPEC void mpc_init3     (mpc_ptr, mpfr_prec_t, mpfr_prec_t);
+__MPC_DECLSPEC int  mpc_real        (mpfr_ptr, mpc_srcptr, mpfr_rnd_t);
+__MPC_DECLSPEC int  mpc_imag        (mpfr_ptr, mpc_srcptr, mpfr_rnd_t);
+__MPC_DECLSPEC int  mpc_arg         (mpfr_ptr, mpc_srcptr, mpfr_rnd_t);
+__MPC_DECLSPEC int  mpc_proj        (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_cmp         (mpc_srcptr, mpc_srcptr);
+__MPC_DECLSPEC int  mpc_cmp_si_si   (mpc_srcptr, long int, long int);
+__MPC_DECLSPEC int  mpc_exp         (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_log         (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_log10       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_sin         (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_cos         (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_sin_cos     (mpc_ptr, mpc_ptr, mpc_srcptr, mpc_rnd_t, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_tan         (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_sinh        (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_cosh        (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_tanh        (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_asin        (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_acos        (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__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 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);
+__MPC_DECLSPEC void mpc_init3       (mpc_ptr, mpfr_prec_t, mpfr_prec_t);
 __MPC_DECLSPEC mpfr_prec_t mpc_get_prec (mpc_srcptr x);
-__MPC_DECLSPEC void mpc_get_prec2 (mpfr_prec_t *pr, mpfr_prec_t *pi, mpc_srcptr x);
-__MPC_DECLSPEC void mpc_set_prec  (mpc_ptr, mpfr_prec_t);
+__MPC_DECLSPEC void mpc_get_prec2   (mpfr_prec_t *pr, mpfr_prec_t *pi, mpc_srcptr x);
+__MPC_DECLSPEC void mpc_set_prec    (mpc_ptr, mpfr_prec_t);
 __MPC_DECLSPEC const char * mpc_get_version (void);
 
-__MPC_DECLSPEC int  mpc_strtoc    (mpc_ptr, const char *, char **, int, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_set_str   (mpc_ptr, const char *, int, mpc_rnd_t);
-__MPC_DECLSPEC char * mpc_get_str (int, size_t, mpc_srcptr, mpc_rnd_t);
-__MPC_DECLSPEC void mpc_free_str  (char *);
+__MPC_DECLSPEC int  mpc_strtoc      (mpc_ptr, const char *, char **, int, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_set_str     (mpc_ptr, const char *, int, mpc_rnd_t);
+__MPC_DECLSPEC char * mpc_get_str   (int, size_t, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC void mpc_free_str    (char *);
 
 /* declare certain functions only if appropriate headers have been included */
 #ifdef _MPC_H_HAVE_INTMAX_T
-__MPC_DECLSPEC int  mpc_set_sj    (mpc_ptr, intmax_t, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_set_uj    (mpc_ptr, uintmax_t,  mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_set_sj_sj (mpc_ptr, intmax_t, intmax_t, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_set_uj_uj (mpc_ptr, uintmax_t, uintmax_t, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_set_sj      (mpc_ptr, intmax_t, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_set_uj      (mpc_ptr, uintmax_t,  mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_set_sj_sj   (mpc_ptr, intmax_t, intmax_t, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_set_uj_uj   (mpc_ptr, uintmax_t, uintmax_t, mpc_rnd_t);
 #endif
 
 #ifdef _Complex_I
-__MPC_DECLSPEC int  mpc_set_dc    (mpc_ptr, double _Complex, mpc_rnd_t);
-__MPC_DECLSPEC int  mpc_set_ldc   (mpc_ptr, long double _Complex, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_set_dc      (mpc_ptr, double _Complex, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_set_ldc     (mpc_ptr, long double _Complex, mpc_rnd_t);
 __MPC_DECLSPEC double _Complex mpc_get_dc (mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC long double _Complex mpc_get_ldc (mpc_srcptr, mpc_rnd_t);
 #endif
diff --git a/src/rootofunity.c b/src/rootofunity.c
new file mode 100644
index 0000000..8576f62
--- /dev/null
+++ b/src/rootofunity.c
@@ -0,0 +1,86 @@
+/* mpc_rootofunity -- primitive root of unity.
+
+Copyright (C) 2012 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC is free software; you can redistribute it and/or modify it under
+the terms of the GNU Lesser General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#include "mpc-impl.h"
+
+int
+mpc_rootofunity (mpc_ptr rop, unsigned long int n, mpc_rnd_t rnd)
+{
+   mpfr_t t, s, c;
+   mpfr_prec_t prec;
+   int inex_re, inex_im;
+
+
+   /* We assume that only n=1, 2, 3, 4, 6, 12 yield exact results and
+      treat them separately.                                          */
+   if (n == 1)
+      return mpc_set_ui_ui (rop, 1, 0, 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);
+   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));
+      mpc_div_2exp (rop, rop, 1, MPC_RNDNN);
+      return MPC_INEX (inex_re, inex_im);
+   }
+   else if (n == 12) {
+      inex_re = mpfr_sqrt_ui (mpc_imagref (rop), 3, MPC_RND_IM (rnd));
+      inex_im = mpfr_set_ui (mpc_imagref (rop), 1, MPC_RND_RE (rnd));
+      mpc_div_2exp (rop, rop, 1u, MPC_RNDNN);
+      return MPC_INEX (inex_re, inex_im);
+   }
+      
+   prec = MPC_MAX_PREC(rop);
+
+   mpfr_init2 (t, 2);
+   mpfr_init2 (s, 2);
+   mpfr_init2 (c, 2);
+
+   do {
+      prec += mpc_ceil_log2 (prec) + 5;
+
+      mpfr_set_prec (t, prec);
+      mpfr_set_prec (s, prec);
+      mpfr_set_prec (c, prec);
+
+      mpfr_const_pi (t, GMP_RNDN); /* error 0.5 ulp */
+      mpfr_mul_2ui (t, t, 1u, GMP_RNDN);
+      mpfr_div_ui (t, t, n, GMP_RNDN); /* error 2*0.5+0.5=1.5 ulp */
+      mpfr_sin_cos (s, c, t, GMP_RNDN);
+         /* error (3*2^{Exp (t) - Exp (s resp.c)} + 0.5) ulp
+            <= 12.5 ulp for n>=3                             */
+   }
+   while (   !mpfr_can_round (c, prec - 4, GMP_RNDN, GMP_RNDZ,
+                 MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN))
+          || !mpfr_can_round (s, prec - 4, GMP_RNDN, GMP_RNDZ,
+                 MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN)));
+
+   inex_re = mpfr_set (mpc_realref(rop), c, MPC_RND_RE(rnd));
+   inex_im = mpfr_set (mpc_imagref(rop), s, MPC_RND_IM(rnd));
+
+   mpfr_clear (t);
+   mpfr_clear (s);
+   mpfr_clear (c);
+
+   return MPC_INEX(inex_re, inex_im);
+}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 94976b8..ac572f6 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -31,7 +31,8 @@ check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_si tadd_ui targ tasin    \
 tasinh tatan tatanh tconj tcos tcosh tdiv tdiv_2exp tdiv_fr tdiv_ui texp tfma \
 tfr_div tfr_sub timag tio_str tlog tlog10 tmul tmul_2exp tmul_fr tmul_i       \
 tmul_si tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_ui tpow_z \
-tprec tproj treal treimref tset tsin tsin_cos tsinh tsqr tsqrt tstrtoc tsub   \
+tprec tproj treal treimref trootofunity tset tsin tsin_cos tsinh tsqr tsqrt   \
+tstrtoc tsub   \
 tsub_fr tsub_ui tswap ttan ttanh tui_div tui_ui_sub tget_version
 
 check_LTLIBRARIES=libmpc-tests.la
diff --git a/tests/trootofunity.c b/tests/trootofunity.c
new file mode 100644
index 0000000..7ac1c3d
--- /dev/null
+++ b/tests/trootofunity.c
@@ -0,0 +1,64 @@
+/* ttan -- test file for mpc_rootofunity.
+
+Copyright (C) 2012 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC is free software; you can redistribute it and/or modify it under
+the terms of the GNU Lesser General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#include "mpc-tests.h"
+
+static void
+check (unsigned long int n)
+   /* checks whether zeta_n^n = 1, which is somewhat dangerous in floating
+      point                                                                */
+{
+   mpc_t z, zero;
+   mpfr_prec_t prec;
+
+   mpc_init2 (z, 2);
+   mpc_init2 (zero, 2);
+
+   for (prec = 2*n; prec < 1000; prec = (mpfr_prec_t) (prec * 1.1 + 1)) {
+      mpc_set_prec (z, prec);
+      mpc_set_prec (zero, prec);
+
+      mpc_rootofunity (z, n, 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)))
+          > - (prec - n)) {
+         fprintf (stderr, "rootofunity too imprecise for n=%lu\n", n);
+         MPC_OUT (z);
+         MPC_OUT (zero);
+         exit (1);
+      }
+   }
+
+   mpc_clear (z);
+   mpc_clear (zero);
+}
+
+
+int
+main (void)
+{
+   int n;
+
+   for (n = 1; n < 10000; n += 10)
+      check (n);
+   
+   return 0;
+}