diff options
author | Paul Zimmermann <Paul.Zimmermann@inria.fr> | 2017-08-22 13:31:10 +0200 |
---|---|---|
committer | Paul Zimmermann <Paul.Zimmermann@inria.fr> | 2017-08-22 13:31:10 +0200 |
commit | 0c6c1efe96d177368df130556686e65d0e6bca54 (patch) | |
tree | 85195951c99fdb147d64ea39160dfb9036c43b39 | |
parent | 4c4778c54b389c7a38cbeb271afceebf8ed492e9 (diff) | |
download | mpc-git-0c6c1efe96d177368df130556686e65d0e6bca54.tar.gz |
started review of rootofunity.c
-rw-r--r-- | src/rootofunity.c | 169 |
1 files changed, 110 insertions, 59 deletions
diff --git a/src/rootofunity.c b/src/rootofunity.c index e4f533e..45d7cab 100644 --- a/src/rootofunity.c +++ b/src/rootofunity.c @@ -21,8 +21,8 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include "mpc-impl.h" #include <assert.h> -static unsigned long int -gcd (unsigned long int a, unsigned long int b) +static unsigned long +gcd (unsigned long a, unsigned long b) { if (b == 0) return a; @@ -31,9 +31,9 @@ gcd (unsigned long int a, unsigned long int 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, unsigned long int k, mpc_rnd_t rnd) +mpc_rootofunity (mpc_ptr rop, unsigned long n, unsigned long k, mpc_rnd_t rnd) { - unsigned long int g; + unsigned long g; mpq_t kn; mpfr_t t, s, c; mpfr_prec_t prec; @@ -53,68 +53,106 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, unsigned long int k, mpc_rnd_ k /= g; n /= g; + /* Now 0 <= k < n and gcd(k,n)=1. */ + /* 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); + { + /* necessarily k=0 thus we want exp(0)=1 */ + MPC_ASSERT (k == 0); + return mpc_set_ui_ui (rop, 1, 0, rnd); + } else if (n == 2) - return mpc_set_si_si (rop, -1, 0, rnd); + { + /* since gcd(k,n)=1, necessarily k=1, thus we want exp(i*pi)=-1 */ + MPC_ASSERT (k == 1); + return mpc_set_si_si (rop, -1, 0, rnd); + } else if (n == 4) - if (k == 1) + { + /* since gcd(k,n)=1, necessarily k=1 or k=3, thus we want + exp(2*i*pi/4)=i or exp(2*i*pi*3/4)=-i */ + MPC_ASSERT (k == 1 || k == 3); + if (k == 1) return mpc_set_ui_ui (rop, 0, 1, rnd); - else + 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)); - /* 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) + } + else if (n == 3 || n == 6) + { + MPC_ASSERT ((n == 3 && (k == 1 || k == 2)) || + (n == 6 && (k == 1 || k == 5))); + /* for n=3, necessarily k=1 or k=2: -1/2+/-1/2*sqrt(3)*i; + for n=6, necessarily k=1 or k=5: 1/2+/-1/2*sqrt(3)*i */ + inex_re = mpfr_set_si (mpc_realref (rop), (n == 3 ? -1 : 1), + MPC_RND_RE (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 == 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) + 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 == 12) + { + /* necessarily k=1, 5, 7, 11: + k=1: 1/2*sqrt(3) + 1/2*I + k=5: -1/2*sqrt(3) + 1/2*I + k=7: -1/2*sqrt(3) - 1/2*I + k=11: 1/2*sqrt(3) - 1/2*I */ + MPC_ASSERT (k == 1 || k == 5 || k == 7 || k == 11); + /* 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 == 8) { - rnd_re = MPC_RND_RE (rnd); - if (k == 3 || k == 5) + inex_re = mpfr_sqrt_ui (mpc_realref (rop), 3, rnd_re); + inex_im = mpfr_set_si (mpc_imagref (rop), k < 6 ? 1 : -1, + MPC_RND_IM (rnd)); + mpc_div_2ui (rop, rop, 1, 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 == 8) + { + /* k=1, 3, 5 or 7: + k=1: (1/2*I + 1/2)*sqrt(2) + k=3: (1/2*I - 1/2)*sqrt(2) + k=5: -(1/2*I + 1/2)*sqrt(2) + k=7: -(1/2*I - 1/2)*sqrt(2) */ + MPC_ASSERT (k == 1 || k == 3 || k == 5 || k == 7); + 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 = MPC_RND_IM (rnd); + if (k > 4) 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); - } + 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, 1, MPC_RNDNN); + if (k == 3 || k == 5) + { + mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN); + inex_re = -inex_re; + } + if (k > 4) + { + mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN); + inex_im = -inex_im; + } + return MPC_INEX (inex_re, inex_im); + } prec = MPC_MAX_PREC(rop); @@ -123,7 +161,7 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, unsigned long int k, mpc_rnd_ mpfr_init2 (c, 67); mpq_init (kn); mpq_set_ui (kn, k, n); - mpq_mul_2exp (kn, kn, 1); + mpq_mul_2exp (kn, kn, 1); /* kn=2*k/n < 2 */ do { prec += mpc_ceil_log2 (prec) + 5; @@ -132,18 +170,31 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, unsigned long int k, mpc_rnd_ mpfr_set_prec (s, prec); mpfr_set_prec (c, prec); - mpfr_const_pi (t, MPFR_RNDN); /* error 0.5 ulp */ - mpfr_mul_q (t, t, kn, MPFR_RNDN); /* error 2*1.5+0.5=3.5 ulp */ + mpfr_const_pi (t, MPFR_RNDN); /* error <= 0.5 ulp but since + ulp(t)=2^(2-prec), the absolute error + is bounded by 2^(1-prec), and the + relative error is bounded by + 2^(1-prec)/pi <= 0.64*2^(-prec) */ + mpfr_mul_q (t, t, kn, MPFR_RNDN); /* error <= 1.15 ulp(t) */ + /* Indeed, the error is bounded by 0.64*2^(-prec)*pi*kn + 0.5 ulp [1]. + Applying 2^(-prec)*|x| <= ulp(x) to x=pi*kn, we get a bound of: + 0.64*ulp(pi*kn)+0.5 ulp. + Now since ulp(pi*kn) <= 2*ulp(t), we get: |t-pi*kn| <= 1.78*ulp(t). + If we plug this into [1] we get: + 0.64*2^(-prec)*t + 0.64*2^(-prec)*1.78*ulp(t) + 0.5 ulp(t) + <= 0.64*ulp(t) + 0.64*2^(-prec)*1.78*ulp(t) + 0.5 ulp(t) + <= 1.15*ulp(t) for prec >= 7. */ + mpfr_sin_cos (s, c, t, MPFR_RNDN); /* 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, but with n<2^64, we lose at most about 64 bits: Where the minimum of s and c over all primitive n-th roots of - unity is reached depends on n mod 4. + 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 + symmetrically 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. |