1 /* Functions for evaluating Gamma(1/3) and Gamma(2/3). Used by mpfr_ai.
2
3 Copyright 2010-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
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.
12
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.
17
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 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 #define MPFR_ACC_OR_MUL(v) \
27 do \
28 { \
29 if (v <= ULONG_MAX / acc) \
30 acc *= v; \
31 else \
32 { \
33 mpfr_mul_ui (y, y, acc, mode); acc = v; \
34 } \
35 } \
36 while (0)
37
38 #define MPFR_ACC_OR_DIV(v) \
39 do \
40 { \
41 if (v <= ULONG_MAX / acc) \
42 acc *= v; \
43 else \
44 { \
45 mpfr_div_ui (y, y, acc, mode); acc = v; \
46 } \
47 } \
48 while (0)
49
50 static void
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)
55 {
56 unsigned long int acc = v1;
57 mpfr_set (y, x, mode);
58 MPFR_ACC_OR_MUL (v2);
59 MPFR_ACC_OR_MUL (v3);
60 MPFR_ACC_OR_MUL (v4);
61 MPFR_ACC_OR_MUL (v5);
62 mpfr_mul_ui (y, y, acc, mode);
63 }
64
65 void
66 mpfr_div_ui2 (mpfr_ptr y, mpfr_srcptr x,
67 unsigned long int v1, unsigned long int v2, mpfr_rnd_t mode)
68 {
69 unsigned long int acc = v1;
70 mpfr_set (y, x, mode);
71 MPFR_ACC_OR_DIV (v2);
72 mpfr_div_ui (y, y, acc, mode);
73 }
74
75 static void
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)
81 {
82 unsigned long int acc = v1;
83 mpfr_set (y, x, mode);
84 MPFR_ACC_OR_DIV (v2);
85 MPFR_ACC_OR_DIV (v3);
86 MPFR_ACC_OR_DIV (v4);
87 MPFR_ACC_OR_DIV (v5);
88 MPFR_ACC_OR_DIV (v6);
89 MPFR_ACC_OR_DIV (v7);
90 MPFR_ACC_OR_DIV (v8);
91 mpfr_div_ui (y, y, acc, mode);
92 }
93
94
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. */
99 static void
100 mpfr_Browns_const (mpfr_ptr s, mpfr_prec_t prec)
101 {
102 mpfr_t uk;
103 unsigned long int k;
104
105 mpfr_prec_t working_prec = prec + 10 + MPFR_INT_CEIL_LOG2 (2 + prec / 10);
106
107 mpfr_init2 (uk, working_prec);
108 mpfr_set_prec (s, working_prec);
109
110 mpfr_set_ui (uk, 1, MPFR_RNDN);
111 mpfr_set (s, uk, MPFR_RNDN);
112 k = 1;
113
114 /* Invariants: uk ~ u(k-1) and s ~ sum(i=0..k-1, u(i)) */
115 for (;;)
116 {
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,
120 MPFR_RNDN);
121 MPFR_CHANGE_SIGN (uk);
122
123 mpfr_add (s, s, uk, MPFR_RNDN);
124 k++;
125 if (MPFR_GET_EXP (uk) + prec <= MPFR_GET_EXP (s) + 7)
126 break;
127 }
128
129 mpfr_clear (uk);
130 return;
131 }
132
133 /* Returns y such that |Gamma(1/3)-y| <= 2^{1-prec}*Gamma(1/3) */
134 static void
135 mpfr_gamma_one_third (mpfr_ptr y, mpfr_prec_t prec)
136 {
137 mpfr_t tmp, tmp2, tmp3;
138
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);
143
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);
148
149 mpfr_Browns_const (tmp2, prec + 9);
150 mpfr_mul (tmp, tmp, tmp2, MPFR_RNDN);
151
152 mpfr_set_ui (tmp2, 10, MPFR_RNDN);
153 mpfr_sqrt (tmp2, tmp2, MPFR_RNDN);
154 mpfr_div (tmp, tmp, tmp2, MPFR_RNDN);
155
156 mpfr_sqrt (tmp3, tmp, MPFR_RNDN);
157 mpfr_cbrt (y, tmp3, MPFR_RNDN);
158
159 mpfr_clear (tmp);
160 mpfr_clear (tmp2);
161 mpfr_clear (tmp3);
162 return;
163 }
164
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) */
168 /* */
169 /* Uses the formula Gamma(z)Gamma(1-z) = pi / sin(pi*z) */
170 /* to compute Gamma(2/3) from Gamma(1/3). */
171 void
172 mpfr_gamma_one_and_two_third (mpfr_ptr y1, mpfr_ptr y2, mpfr_prec_t prec)
173 {
174 mpfr_t temp;
175
176 mpfr_init2 (temp, prec + 4);
177 mpfr_set_prec (y2, prec + 4);
178
179 mpfr_gamma_one_third (y1, prec + 4);
180
181 mpfr_set_ui (temp, 3, MPFR_RNDN);
182 mpfr_sqrt (temp, temp, MPFR_RNDN);
183 mpfr_mul (temp, y1, temp, MPFR_RNDN);
184
185 mpfr_const_pi (y2, MPFR_RNDN);
186 mpfr_mul_2ui (y2, y2, 1, MPFR_RNDN);
187
188 mpfr_div (y2, y2, temp, MPFR_RNDN);
189
190 mpfr_clear (temp);
191 }