From b17945e0ad92be088e36bad22fe41438776b6f35 Mon Sep 17 00:00:00 2001 From: enge Date: Wed, 27 Jun 2012 14:39:37 +0000 Subject: implemented rootofunity git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/branches/rootsunity@1194 211d60ee-9f03-0410-a15a-8952a2c7a4e4 --- src/Makefile.am | 3 +- src/mpc.h | 77 +++++++++++++++++++++++----------------------- src/rootofunity.c | 86 ++++++++++++++++++++++++++++++++++++++++++++++++++++ tests/Makefile.am | 3 +- tests/trootofunity.c | 64 ++++++++++++++++++++++++++++++++++++++ 5 files changed, 193 insertions(+), 40 deletions(-) create mode 100644 src/rootofunity.c create mode 100644 tests/trootofunity.c 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; +} -- cgit v1.2.1