1/*  mpfr_ui_pow_ui -- compute the power between two machine integers
2
3Copyright 1999-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#include "mpfr-impl.h"
24
25int
26mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int k, unsigned long int n,
27                mpfr_rnd_t rnd)
28{
29  mpfr_exp_t err;
30  unsigned long m;
31  mpfr_t res;
32  mpfr_prec_t prec;
33  int size_n;
34  int inexact;
35  MPFR_ZIV_DECL (loop);
36  MPFR_SAVE_EXPO_DECL (expo);
37
38  MPFR_LOG_FUNC
39    (("k=%lu n=%lu rnd=%d", k, n, rnd),
40     ("y[%Pd]=%.*Rg inexact=%d",
41      mpfr_get_prec (x), mpfr_log_prec, x, inexact));
42
43  if (MPFR_UNLIKELY (n <= 1))
44    {
45      if (n == 1)
46        return mpfr_set_ui (x, k, rnd);     /* k^1 = k */
47      else
48        return mpfr_set_ui (x, 1, rnd);     /* k^0 = 1 for any k */
49    }
50  else if (MPFR_UNLIKELY (k <= 1))
51    {
52      if (k == 1)
53        return mpfr_set_ui (x, 1, rnd);     /* 1^n = 1 for any n > 0 */
54      else
55        return mpfr_set_ui (x, 0, rnd);     /* 0^n = 0 for any n > 0 */
56    }
57
58  for (size_n = 0, m = n; m != 0; size_n++, m >>= 1)
59    ;
60
61  MPFR_SAVE_EXPO_MARK (expo);
62  prec = MPFR_PREC (x) + 3 + size_n;
63  mpfr_init2 (res, prec);
64
65  MPFR_ZIV_INIT (loop, prec);
66  for (;;)
67    {
68      int i = size_n;
69      unsigned int inex_res;
70
71      inex_res = mpfr_set_ui (res, k, MPFR_RNDU);
72      err = 1;
73      /* now 2^(i-1) <= n < 2^i: i=1+floor(log2(n)) */
74      for (i -= 2; i >= 0; i--)
75        {
76          inex_res |= mpfr_sqr (res, res, MPFR_RNDU);
77          err++;
78          if (n & (1UL << i))
79            inex_res |= mpfr_mul_ui (res, res, k, MPFR_RNDU);
80        }
81
82      if (MPFR_UNLIKELY (MPFR_IS_INF (res)))
83        {
84          mpfr_t kf;
85          mpz_t z;
86          int size_k;
87          MPFR_BLOCK_DECL (flags);
88
89          /* Let's handle the overflow by calling mpfr_pow_z.
90             Alternatively, we could call mpfr_pow_ui; this would
91             need a bit shorter code below, but mpfr_pow_ui handles
92             the overflow by calling mpfr_pow_z, so that calling
93             mpfr_pow_z directly should be a bit more efficient. */
94
95          MPFR_ZIV_FREE (loop);
96          mpfr_clear (res);
97          for (size_k = 0, m = k; m != 0; size_k++, m >>= 1)
98            ;
99          mpfr_init2 (kf, size_k);
100          inexact = mpfr_set_ui (kf, k, MPFR_RNDN);
101          MPFR_ASSERTD (inexact == 0);
102          mpz_init (z);
103          mpz_set_ui (z, n);
104          MPFR_BLOCK (flags, inexact = mpfr_pow_z (x, kf, z, rnd););
105          mpz_clear (z);
106          mpfr_clear (kf);
107          MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, flags);
108          goto end;
109        }
110
111      /* since the loop is executed floor(log2(n)) times,
112         we have err = 1+floor(log2(n)).
113         Since prec >= MPFR_PREC(x) + 4 + floor(log2(n)), prec > err */
114      err = prec - err;
115
116      MPFR_LOG_VAR (res);
117      if (MPFR_LIKELY (!inex_res
118                       || MPFR_CAN_ROUND (res, err, MPFR_PREC (x), rnd)))
119        break;
120
121      /* Actualisation of the precision */
122      MPFR_ZIV_NEXT (loop, prec);
123      mpfr_set_prec (res, prec);
124    }
125  MPFR_ZIV_FREE (loop);
126
127  inexact = mpfr_set (x, res, rnd);
128
129  mpfr_clear (res);
130
131 end:
132  MPFR_SAVE_EXPO_FREE (expo);
133  return mpfr_check_range (x, inexact, rnd);
134}
135