xref: /netbsd-src/external/lgpl3/mpfr/dist/src/cache.c (revision 33881f779a77dce6440bdc44610d94de75bebefe)
1 /* mpfr_cache -- cache interface for multiple-precision constants in MPFR.
2 
3 Copyright 2004-2018 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 http://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 prec = MPFR_PREC (dest);
60   mpfr_prec_t pold;
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   pold = MPFR_PREC (cache->x);
73   if (MPFR_UNLIKELY (prec > pold))
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.
80          If there is no lock, there is no harm in this code */
81       pold = MPFR_PREC (cache->x);
82       if (MPFR_LIKELY (prec > pold))
83         {
84           /* No previous result in the cache or the precision of the previous
85              result is not sufficient. We increase the cache size by at least
86              10% to avoid invalidating the cache many times if one performs
87              several computations with small increase of precision. */
88           if (MPFR_UNLIKELY (pold == 0))  /* No previous result. */
89             mpfr_init2 (cache->x, prec);  /* as pold = prec below */
90           else
91             pold += pold / 10;
92 
93           /* Update the cache. */
94           if (pold < prec)
95             pold = prec;
96 
97           /* no need to keep the previous value */
98           mpfr_set_prec (cache->x, pold);
99           cache->inexact = (*cache->func) (cache->x, MPFR_RNDN);
100         }
101 
102       /* Free the cache in read-write mode */
103       /* Get the cache in read-only mode */
104       MPFR_LOCK_WRITE2READ(cache->lock);
105     }
106 
107   /* now pold >= prec is the precision of cache->x */
108   MPFR_ASSERTD (pold >= prec);
109   MPFR_ASSERTD (MPFR_PREC (cache->x) == pold);
110 
111   /* First, check if the cache has the exact value (unlikely).
112      Else the exact value is between (assuming x=cache->x > 0):
113        x and x+ulp(x) if cache->inexact < 0,
114        x-ulp(x) and x if cache->inexact > 0,
115      and abs(x-exact) <= ulp(x)/2. */
116 
117   /* we assume all cached constants are positive */
118   MPFR_ASSERTN (MPFR_IS_POS (cache->x)); /* TODO... */
119   sign = MPFR_SIGN (cache->x);
120   MPFR_EXP (dest) = MPFR_GET_EXP (cache->x);
121   MPFR_SET_SIGN (dest, sign);
122 
123   /* round cache->x from precision pold down to precision prec */
124   MPFR_RNDRAW_GEN (inexact, dest,
125                    MPFR_MANT (cache->x), pold, rnd, sign,
126                    if (MPFR_UNLIKELY (cache->inexact == 0))
127                      {
128                        if ((_sp[0] & _ulp) == 0)
129                          {
130                            inexact = -sign;
131                            goto trunc_doit;
132                          }
133                        else
134                          goto addoneulp;
135                      }
136                    else if (cache->inexact < 0)
137                      goto addoneulp;
138                    else /* cache->inexact > 0 */
139                      {
140                        inexact = -sign;
141                        goto trunc_doit;
142                      },
143                    if (MPFR_UNLIKELY (++MPFR_EXP (dest) > __gmpfr_emax))
144                      mpfr_overflow (dest, rnd, sign);
145                   );
146 
147   /* Rather a likely, this is a 100% succes rate for
148      all constants of MPFR */
149   if (MPFR_LIKELY (cache->inexact != 0))
150     {
151       switch (rnd)
152         {
153         case MPFR_RNDZ:
154         case MPFR_RNDD:
155           if (MPFR_UNLIKELY (inexact == 0))
156             {
157               inexact = cache->inexact;
158               if (inexact > 0)
159                 {
160                   mpfr_nextbelow (dest);
161                   inexact = -inexact;
162                 }
163             }
164           break;
165         case MPFR_RNDU:
166         case MPFR_RNDA:
167           if (MPFR_UNLIKELY (inexact == 0))
168             {
169               inexact = cache->inexact;
170               if (inexact < 0)
171                 {
172                   mpfr_nextabove (dest);
173                   inexact = -inexact;
174                 }
175             }
176           break;
177         default: /* MPFR_RNDN */
178           if (MPFR_UNLIKELY(inexact == 0))
179             inexact = cache->inexact;
180           break;
181         }
182     }
183 
184   MPFR_SAVE_EXPO_FREE (expo);
185 
186   /* Free the cache in read-only mode */
187   MPFR_UNLOCK_READ(cache->lock);
188 
189   return mpfr_check_range (dest, inexact, rnd);
190 }
191