1 /* mpc_pow_ui -- Raise a complex number to an integer power. 2 3 Copyright (C) 2009, 2010, 2011, 2012 INRIA 4 5 This file is part of GNU MPC. 6 7 GNU MPC is free software; you can redistribute it and/or modify it under 8 the terms of the GNU Lesser General Public License as published by the 9 Free Software Foundation; either version 3 of the License, or (at your 10 option) any later version. 11 12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15 more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along with this program. If not, see http://www.gnu.org/licenses/ . 19 */ 20 21 #include <limits.h> /* for CHAR_BIT */ 22 #include "mpc-impl.h" 23 24 static int 25 mpc_pow_usi_naive (mpc_ptr z, mpc_srcptr x, unsigned long y, int sign, 26 mpc_rnd_t rnd) 27 { 28 int inex; 29 mpc_t t; 30 31 mpc_init3 (t, sizeof (unsigned long) * CHAR_BIT, MPFR_PREC_MIN); 32 if (sign > 0) 33 mpc_set_ui (t, y, MPC_RNDNN); /* exact */ 34 else 35 mpc_set_si (t, - (signed long) y, MPC_RNDNN); 36 inex = mpc_pow (z, x, t, rnd); 37 mpc_clear (t); 38 39 return inex; 40 } 41 42 43 int 44 mpc_pow_usi (mpc_ptr z, mpc_srcptr x, unsigned long y, int sign, 45 mpc_rnd_t rnd) 46 /* computes z = x^(sign*y) */ 47 { 48 int inex; 49 mpc_t t, x3; 50 mpfr_prec_t p, l, l0; 51 long unsigned int u; 52 int has3; /* non-zero if y has '11' in its binary representation */ 53 int loop, done; 54 55 /* let mpc_pow deal with special values */ 56 if (!mpc_fin_p (x) || mpfr_zero_p (mpc_realref (x)) || mpfr_zero_p (mpc_imagref(x)) 57 || y == 0) 58 return mpc_pow_usi_naive (z, x, y, sign, rnd); 59 /* easy special cases */ 60 else if (y == 1) { 61 if (sign > 0) 62 return mpc_set (z, x, rnd); 63 else 64 return mpc_ui_div (z, 1ul, x, rnd); 65 } 66 else if (y == 2 && sign > 0) 67 return mpc_sqr (z, x, rnd); 68 /* let mpc_pow treat potential over- and underflows */ 69 else { 70 mpfr_exp_t exp_r = mpfr_get_exp (mpc_realref (x)), 71 exp_i = mpfr_get_exp (mpc_imagref (x)); 72 if ( MPC_MAX (exp_r, exp_i) > mpfr_get_emax () / (mpfr_exp_t) y 73 /* heuristic for overflow */ 74 || MPC_MAX (-exp_r, -exp_i) > (-mpfr_get_emin ()) / (mpfr_exp_t) y 75 /* heuristic for underflow */ 76 ) 77 return mpc_pow_usi_naive (z, x, y, sign, rnd); 78 } 79 80 has3 = (y & (y >> 1)) != 0; 81 for (l = 0, u = y; u > 3; l ++, u >>= 1); 82 /* l>0 is the number of bits of y, minus 2, thus y has bits: 83 y_{l+1} y_l y_{l-1} ... y_1 y_0 */ 84 l0 = l + 2; 85 p = MPC_MAX_PREC(z) + l0 + 32; /* l0 ensures that y*2^{-p} <= 1 below */ 86 mpc_init2 (t, p); 87 if (has3) 88 mpc_init2 (x3, p); 89 90 loop = 0; 91 done = 0; 92 while (!done) { 93 loop++; 94 95 mpc_sqr (t, x, MPC_RNDNN); 96 if (has3) { 97 mpc_mul (x3, t, x, MPC_RNDNN); 98 if ((y >> l) & 1) /* y starts with 11... */ 99 mpc_set (t, x3, MPC_RNDNN); 100 } 101 while (l-- > 0) { 102 mpc_sqr (t, t, MPC_RNDNN); 103 if ((y >> l) & 1) { 104 if ((l > 0) && ((y >> (l-1)) & 1)) /* implies has3 <> 0 */ { 105 l--; 106 mpc_sqr (t, t, MPC_RNDNN); 107 mpc_mul (t, t, x3, MPC_RNDNN); 108 } 109 else 110 mpc_mul (t, t, x, MPC_RNDNN); 111 } 112 } 113 if (sign < 0) 114 mpc_ui_div (t, 1ul, t, MPC_RNDNN); 115 116 if (mpfr_zero_p (mpc_realref(t)) || mpfr_zero_p (mpc_imagref(t))) { 117 inex = mpc_pow_usi_naive (z, x, y, sign, rnd); 118 /* since mpfr_get_exp() is not defined for zero */ 119 done = 1; 120 } 121 else { 122 /* see error bound in algorithms.tex; we use y<2^l0 instead of y-1 123 also when sign>0 */ 124 mpfr_exp_t diff; 125 mpfr_prec_t er, ei; 126 127 diff = mpfr_get_exp (mpc_realref(t)) - mpfr_get_exp (mpc_imagref(t)); 128 /* the factor on the real part is 2+2^(-diff+2) <= 4 for diff >= 1 129 and < 2^(-diff+3) for diff <= 0 */ 130 er = (diff >= 1) ? l0 + 3 : l0 + (-diff) + 3; 131 /* the factor on the imaginary part is 2+2^(diff+2) <= 4 for diff <= -1 132 and < 2^(diff+3) for diff >= 0 */ 133 ei = (diff <= -1) ? l0 + 3 : l0 + diff + 3; 134 if (mpfr_can_round (mpc_realref(t), p - er, GMP_RNDN, GMP_RNDZ, 135 MPC_PREC_RE(z) + (MPC_RND_RE(rnd) == GMP_RNDN)) 136 && mpfr_can_round (mpc_imagref(t), p - ei, GMP_RNDN, GMP_RNDZ, 137 MPC_PREC_IM(z) + (MPC_RND_IM(rnd) == GMP_RNDN))) { 138 inex = mpc_set (z, t, rnd); 139 done = 1; 140 } 141 else if (loop == 1 && SAFE_ABS(mpfr_prec_t, diff) < MPC_MAX_PREC(z)) { 142 /* common case, make a second trial at higher precision */ 143 p += MPC_MAX_PREC(x); 144 mpc_set_prec (t, p); 145 if (has3) 146 mpc_set_prec (x3, p); 147 l = l0 - 2; 148 } 149 else { 150 /* stop the loop and use mpc_pow */ 151 inex = mpc_pow_usi_naive (z, x, y, sign, rnd); 152 done = 1; 153 } 154 } 155 } 156 157 mpc_clear (t); 158 if (has3) 159 mpc_clear (x3); 160 161 return inex; 162 } 163 164 165 int 166 mpc_pow_ui (mpc_ptr z, mpc_srcptr x, unsigned long y, mpc_rnd_t rnd) 167 { 168 return mpc_pow_usi (z, x, y, 1, rnd); 169 } 170