1 /*
2 Copyright (C) 2013 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "arb.h"
13
14 static void
bsplit(arb_t y,const fmpz_t p,const fmpz_t q,ulong a,ulong b,slong prec)15 bsplit(arb_t y, const fmpz_t p, const fmpz_t q, ulong a, ulong b, slong prec)
16 {
17 if (b - a <= 8)
18 {
19 fmpz_t t, u;
20 ulong c;
21
22 fmpz_init(t);
23 fmpz_init(u);
24
25 fmpz_mul_ui(t, q, a);
26 fmpz_add(t, t, p);
27 fmpz_set(u, t);
28
29 for (c = a + 1; c < b; c++)
30 {
31 fmpz_add(u, u, q);
32 fmpz_mul(t, t, u);
33 }
34
35 arb_set_round_fmpz(y, t, prec);
36
37 fmpz_clear(t);
38 fmpz_clear(u);
39 }
40 else
41 {
42 arb_t w;
43 ulong m = a + (b - a) / 2;
44 arb_init(w);
45
46 bsplit(y, p, q, a, m, prec);
47 bsplit(w, p, q, m, b, prec);
48
49 arb_mul(y, y, w, prec);
50 arb_clear(w);
51 }
52 }
53
54 void
arb_rising_fmpq_ui(arb_t y,const fmpq_t x,ulong n,slong prec)55 arb_rising_fmpq_ui(arb_t y, const fmpq_t x, ulong n, slong prec)
56 {
57 if (n == 0)
58 {
59 arb_one(y);
60 }
61 else if (n == 1)
62 {
63 arb_set_fmpq(y, x, prec);
64 }
65 else
66 {
67 slong wp;
68
69 wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
70
71 bsplit(y, fmpq_numref(x), fmpq_denref(x), 0, n, wp);
72
73 if (fmpz_is_one(fmpq_denref(x)))
74 {
75 arb_set_round(y, y, prec);
76 }
77 else
78 {
79 arb_t t;
80 arb_init(t);
81 arb_set_fmpz(t, fmpq_denref(x));
82 arb_pow_ui(t, t, n, wp);
83 arb_div(y, y, t, prec);
84 arb_clear(t);
85 }
86 }
87 }
88
89