1 /* Functions for evaluating Gamma(1/3) and Gamma(2/3). Used by mpfr_ai.
3 Copyright 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 #define MPFR_ACC_OR_MUL(v) \
29 if (v <= ULONG_MAX / acc) \
33 mpfr_mul_ui (y, y, acc, mode); acc = v; \
38 #define MPFR_ACC_OR_DIV(v) \
41 if (v <= ULONG_MAX / acc) \
45 mpfr_div_ui (y, y, acc, mode); acc = v; \
51 mpfr_mul_ui5 (mpfr_ptr y
, mpfr_srcptr x
,
52 unsigned long int v1
, unsigned long int v2
,
53 unsigned long int v3
, unsigned long int v4
,
54 unsigned long int v5
, mpfr_rnd_t mode
)
56 unsigned long int acc
= v1
;
57 mpfr_set (y
, x
, mode
);
62 mpfr_mul_ui (y
, y
, acc
, mode
);
66 mpfr_div_ui2 (mpfr_ptr y
, mpfr_srcptr x
,
67 unsigned long int v1
, unsigned long int v2
, mpfr_rnd_t mode
)
69 unsigned long int acc
= v1
;
70 mpfr_set (y
, x
, mode
);
72 mpfr_div_ui (y
, y
, acc
, mode
);
76 mpfr_div_ui8 (mpfr_ptr y
, mpfr_srcptr x
,
77 unsigned long int v1
, unsigned long int v2
,
78 unsigned long int v3
, unsigned long int v4
,
79 unsigned long int v5
, unsigned long int v6
,
80 unsigned long int v7
, unsigned long int v8
, mpfr_rnd_t mode
)
82 unsigned long int acc
= v1
;
83 mpfr_set (y
, x
, mode
);
91 mpfr_div_ui (y
, y
, acc
, mode
);
95 /* Gives an approximation of omega = Gamma(1/3)^6 * sqrt(10) / (12pi^4) */
96 /* using C. H. Brown's formula. */
97 /* The computed value s satisfies |s-omega| <= 2^{1-prec}*omega */
98 /* As usual, the variable s is supposed to be initialized. */
100 mpfr_Browns_const (mpfr_ptr s
, mpfr_prec_t prec
)
105 mpfr_prec_t working_prec
= prec
+ 10 + MPFR_INT_CEIL_LOG2 (2 + prec
/ 10);
107 mpfr_init2 (uk
, working_prec
);
108 mpfr_set_prec (s
, working_prec
);
110 mpfr_set_ui (uk
, 1, MPFR_RNDN
);
111 mpfr_set (s
, uk
, MPFR_RNDN
);
114 /* Invariants: uk ~ u(k-1) and s ~ sum(i=0..k-1, u(i)) */
117 mpfr_mul_ui5 (uk
, uk
, 6 * k
- 5, 6 * k
- 4, 6 * k
- 3, 6 * k
- 2,
118 6 * k
- 1, MPFR_RNDN
);
119 mpfr_div_ui8 (uk
, uk
, k
, k
, 3 * k
- 2, 3 * k
- 1, 3 * k
, 80, 160, 160,
121 MPFR_CHANGE_SIGN (uk
);
123 mpfr_add (s
, s
, uk
, MPFR_RNDN
);
125 if (MPFR_GET_EXP (uk
) + prec
<= MPFR_GET_EXP (s
) + 7)
133 /* Returns y such that |Gamma(1/3)-y| <= 2^{1-prec}*Gamma(1/3) */
135 mpfr_gamma_one_third (mpfr_ptr y
, mpfr_prec_t prec
)
137 mpfr_t tmp
, tmp2
, tmp3
;
139 mpfr_init2 (tmp
, prec
+ 9);
140 mpfr_init2 (tmp2
, prec
+ 9);
141 mpfr_init2 (tmp3
, prec
+ 4);
142 mpfr_set_prec (y
, prec
+ 2);
144 mpfr_const_pi (tmp
, MPFR_RNDN
);
145 mpfr_sqr (tmp
, tmp
, MPFR_RNDN
);
146 mpfr_sqr (tmp
, tmp
, MPFR_RNDN
);
147 mpfr_mul_ui (tmp
, tmp
, 12, MPFR_RNDN
);
149 mpfr_Browns_const (tmp2
, prec
+ 9);
150 mpfr_mul (tmp
, tmp
, tmp2
, MPFR_RNDN
);
152 mpfr_set_ui (tmp2
, 10, MPFR_RNDN
);
153 mpfr_sqrt (tmp2
, tmp2
, MPFR_RNDN
);
154 mpfr_div (tmp
, tmp
, tmp2
, MPFR_RNDN
);
156 mpfr_sqrt (tmp3
, tmp
, MPFR_RNDN
);
157 mpfr_cbrt (y
, tmp3
, MPFR_RNDN
);
165 /* Computes y1 and y2 such that: */
166 /* |y1-Gamma(1/3)| <= 2^{1-prec}Gamma(1/3) */
167 /* and |y2-Gamma(2/3)| <= 2^{1-prec}Gamma(2/3) */
169 /* Uses the formula Gamma(z)Gamma(1-z) = pi / sin(pi*z) */
170 /* to compute Gamma(2/3) from Gamma(1/3). */
172 mpfr_gamma_one_and_two_third (mpfr_ptr y1
, mpfr_ptr y2
, mpfr_prec_t prec
)
176 mpfr_init2 (temp
, prec
+ 4);
177 mpfr_set_prec (y2
, prec
+ 4);
179 mpfr_gamma_one_third (y1
, prec
+ 4);
181 mpfr_set_ui (temp
, 3, MPFR_RNDN
);
182 mpfr_sqrt (temp
, temp
, MPFR_RNDN
);
183 mpfr_mul (temp
, y1
, temp
, MPFR_RNDN
);
185 mpfr_const_pi (y2
, MPFR_RNDN
);
186 mpfr_mul_2ui (y2
, y2
, 1, MPFR_RNDN
);
188 mpfr_div (y2
, y2
, temp
, MPFR_RNDN
);