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
mpc_pow_usi_naive(mpc_ptr z,mpc_srcptr x,unsigned long y,int sign,mpc_rnd_t rnd)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
mpc_pow_usi(mpc_ptr z,mpc_srcptr x,unsigned long y,int sign,mpc_rnd_t rnd)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
mpc_pow_ui(mpc_ptr z,mpc_srcptr x,unsigned long y,mpc_rnd_t rnd)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