1 /* mpfr_cache -- cache interface for multiple-precision constants in MPFR.
2
3 Copyright 2004-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 #include "mpfr-impl.h"
24
25 #if 0 /* this function is not used/documented/tested so far, it could be
26 useful if some user wants to add a new constant to mpfr, and
27 implement a cache mechanism for that constant */
28 void
29 mpfr_init_cache (mpfr_cache_t cache, int (*func)(mpfr_ptr, mpfr_rnd_t))
30 {
31 MPFR_PREC (cache->x) = 0; /* Invalid prec to detect that the cache is not
32 valid. Maybe add a flag? */
33 cache->func = func;
34 }
35 #endif
36
37 void
38 mpfr_clear_cache (mpfr_cache_t cache)
39 {
40 if (MPFR_UNLIKELY (MPFR_PREC (cache->x) != 0))
41 {
42 /* Get the cache in read-write mode */
43 MPFR_LOCK_WRITE(cache->lock);
44
45 if (MPFR_LIKELY (MPFR_PREC (cache->x) != 0))
46 {
47 mpfr_clear (cache->x);
48 MPFR_PREC (cache->x) = 0;
49 }
50
51 /* Free the cache in read-write mode */
52 MPFR_UNLOCK_WRITE(cache->lock);
53 }
54 }
55
56 int
57 mpfr_cache (mpfr_ptr dest, mpfr_cache_t cache, mpfr_rnd_t rnd)
58 {
59 mpfr_prec_t dprec = MPFR_PREC (dest);
60 mpfr_prec_t cprec; /* precision of the cache */
61 int inexact, sign;
62 MPFR_SAVE_EXPO_DECL (expo);
63
64 /* Call the initialisation function of the cache if it's needed */
65 MPFR_DEFERRED_INIT_CALL(cache);
66
67 MPFR_SAVE_EXPO_MARK (expo);
68
69 /* Get the cache in read-only mode */
70 MPFR_LOCK_READ(cache->lock);
71 /* Read the precision within the cache */
72 cprec = MPFR_PREC (cache->x);
73 if (MPFR_UNLIKELY (dprec > cprec))
74 {
75 /* Free the cache in read-only mode */
76 /* And get the cache in read-write mode */
77 MPFR_LOCK_READ2WRITE(cache->lock);
78
79 /* Retest the precision once we get the lock (since it might have
80 changed). If there is no lock, there is no harm in this code. */
81 cprec = MPFR_PREC (cache->x);
82 if (MPFR_LIKELY (dprec > cprec))
83 {
84 /* No previous result in the cache or the precision of the
85 previous result is not sufficient. */
86 if (MPFR_UNLIKELY (cprec == 0)) /* No previous result. */
87 {
88 cprec = dprec;
89 mpfr_init2 (cache->x, cprec);
90 }
91 else
92 {
93 /* We increase the cache size by at least 10% to avoid
94 invalidating the cache many times if one performs
95 several computations with small increase of precision. */
96 cprec += cprec / 10;
97 if (cprec < dprec)
98 cprec = dprec;
99 /* no need to keep the previous value */
100 mpfr_set_prec (cache->x, cprec);
101 }
102
103 cache->inexact = (*cache->func) (cache->x, MPFR_RNDN);
104 }
105
106 /* Free the cache in read-write mode */
107 /* Get the cache in read-only mode */
108 MPFR_LOCK_WRITE2READ(cache->lock);
109 }
110
111 /* now cprec >= dprec is the precision of cache->x */
112 MPFR_ASSERTD (cprec >= dprec);
113 MPFR_ASSERTD (MPFR_PREC (cache->x) == cprec);
114
115 /* First, check if the cache has the exact value (unlikely).
116 Else the exact value is between (assuming x=cache->x > 0):
117 x and x+ulp(x) if cache->inexact < 0,
118 x-ulp(x) and x if cache->inexact > 0,
119 and abs(x-exact) <= ulp(x)/2. */
120
121 /* we assume all cached constants are positive */
122 MPFR_ASSERTN (MPFR_IS_POS (cache->x)); /* TODO... */
123 sign = MPFR_SIGN (cache->x);
124 MPFR_EXP (dest) = MPFR_GET_EXP (cache->x);
125 MPFR_SET_SIGN (dest, sign);
126
127 /* round cache->x from precision cprec down to precision dprec;
128 since we are in extended exponent range, for the values considered
129 here, an overflow is not possible (and wouldn't make much sense). */
130 MPFR_RNDRAW_GEN (inexact, dest,
131 MPFR_MANT (cache->x), cprec, rnd, sign,
132 if (MPFR_UNLIKELY (cache->inexact == 0))
133 {
134 if ((_sp[0] & _ulp) == 0)
135 {
136 inexact = -sign;
137 goto trunc_doit;
138 }
139 else
140 goto addoneulp;
141 }
142 else if (cache->inexact < 0)
143 goto addoneulp;
144 else /* cache->inexact > 0 */
145 {
146 inexact = -sign;
147 goto trunc_doit;
148 },
149 MPFR_EXP (dest) ++;
150 MPFR_ASSERTD (MPFR_EXP (dest) <= __gmpfr_emax);
151 );
152
153 /* Rather a likely, this is a 100% success rate for
154 all constants of MPFR */
155 if (MPFR_LIKELY (cache->inexact != 0))
156 {
157 switch (rnd)
158 {
159 case MPFR_RNDZ:
160 case MPFR_RNDD:
161 if (MPFR_UNLIKELY (inexact == 0))
162 {
163 inexact = cache->inexact;
164 if (inexact > 0)
165 {
166 mpfr_nextbelow (dest);
167 inexact = -inexact;
168 }
169 }
170 break;
171 case MPFR_RNDU:
172 case MPFR_RNDA:
173 if (MPFR_UNLIKELY (inexact == 0))
174 {
175 inexact = cache->inexact;
176 if (inexact < 0)
177 {
178 mpfr_nextabove (dest);
179 inexact = -inexact;
180 }
181 }
182 break;
183 default: /* MPFR_RNDN */
184 if (MPFR_UNLIKELY(inexact == 0))
185 inexact = cache->inexact;
186 break;
187 }
188 }
189
190 MPFR_SAVE_EXPO_FREE (expo);
191
192 /* Free the cache in read-only mode */
193 MPFR_UNLOCK_READ(cache->lock);
194
195 return mpfr_check_range (dest, inexact, rnd);
196 }