1 //===-- Utils which wrap MPFR ---------------------------------------------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 
9 #include "MPFRUtils.h"
10 
11 #include "utils/FPUtil/FPBits.h"
12 #include "utils/FPUtil/TestHelpers.h"
13 
14 #include "llvm/ADT/StringExtras.h"
15 #include "llvm/ADT/StringRef.h"
16 
17 #include <memory>
18 #include <mpfr.h>
19 #include <stdint.h>
20 #include <string>
21 
22 template <typename T> using FPBits = __llvm_libc::fputil::FPBits<T>;
23 
24 namespace __llvm_libc {
25 namespace testing {
26 namespace mpfr {
27 
28 class MPFRNumber {
29   // A precision value which allows sufficiently large additional
30   // precision even compared to quad-precision floating point values.
31   static constexpr unsigned int mpfrPrecision = 128;
32 
33   mpfr_t value;
34 
35 public:
36   MPFRNumber() { mpfr_init2(value, mpfrPrecision); }
37 
38   // We use explicit EnableIf specializations to disallow implicit
39   // conversions. Implicit conversions can potentially lead to loss of
40   // precision.
41   template <typename XType,
42             cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0>
43   explicit MPFRNumber(XType x) {
44     mpfr_init2(value, mpfrPrecision);
45     mpfr_set_flt(value, x, MPFR_RNDN);
46   }
47 
48   template <typename XType,
49             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
50   explicit MPFRNumber(XType x) {
51     mpfr_init2(value, mpfrPrecision);
52     mpfr_set_d(value, x, MPFR_RNDN);
53   }
54 
55   template <typename XType,
56             cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0>
57   explicit MPFRNumber(XType x) {
58     mpfr_init2(value, mpfrPrecision);
59     mpfr_set_ld(value, x, MPFR_RNDN);
60   }
61 
62   template <typename XType,
63             cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
64   explicit MPFRNumber(XType x) {
65     mpfr_init2(value, mpfrPrecision);
66     mpfr_set_sj(value, x, MPFR_RNDN);
67   }
68 
69   MPFRNumber(const MPFRNumber &other) {
70     mpfr_set(value, other.value, MPFR_RNDN);
71   }
72 
73   ~MPFRNumber() {
74     mpfr_clear(value);
75   }
76 
77   MPFRNumber &operator=(const MPFRNumber &rhs) {
78     mpfr_set(value, rhs.value, MPFR_RNDN);
79     return *this;
80   }
81 
82   MPFRNumber abs() const {
83     MPFRNumber result;
84     mpfr_abs(result.value, value, MPFR_RNDN);
85     return result;
86   }
87 
88   MPFRNumber ceil() const {
89     MPFRNumber result;
90     mpfr_ceil(result.value, value);
91     return result;
92   }
93 
94   MPFRNumber cos() const {
95     MPFRNumber result;
96     mpfr_cos(result.value, value, MPFR_RNDN);
97     return result;
98   }
99 
100   MPFRNumber exp() const {
101     MPFRNumber result;
102     mpfr_exp(result.value, value, MPFR_RNDN);
103     return result;
104   }
105 
106   MPFRNumber exp2() const {
107     MPFRNumber result;
108     mpfr_exp2(result.value, value, MPFR_RNDN);
109     return result;
110   }
111 
112   MPFRNumber floor() const {
113     MPFRNumber result;
114     mpfr_floor(result.value, value);
115     return result;
116   }
117 
118   MPFRNumber frexp(int &exp) {
119     MPFRNumber result;
120     mpfr_exp_t resultExp;
121     mpfr_frexp(&resultExp, result.value, value, MPFR_RNDN);
122     exp = resultExp;
123     return result;
124   }
125 
126   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
127     MPFRNumber remainder;
128     long q;
129     mpfr_remquo(remainder.value, &q, value, divisor.value, MPFR_RNDN);
130     quotient = q;
131     return remainder;
132   }
133 
134   MPFRNumber round() const {
135     MPFRNumber result;
136     mpfr_round(result.value, value);
137     return result;
138   }
139 
140   MPFRNumber sin() const {
141     MPFRNumber result;
142     mpfr_sin(result.value, value, MPFR_RNDN);
143     return result;
144   }
145 
146   MPFRNumber sqrt() const {
147     MPFRNumber result;
148     mpfr_sqrt(result.value, value, MPFR_RNDN);
149     return result;
150   }
151 
152   MPFRNumber trunc() const {
153     MPFRNumber result;
154     mpfr_trunc(result.value, value);
155     return result;
156   }
157 
158   std::string str() const {
159     // 200 bytes should be more than sufficient to hold a 100-digit number
160     // plus additional bytes for the decimal point, '-' sign etc.
161     constexpr size_t printBufSize = 200;
162     char buffer[printBufSize];
163     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
164     llvm::StringRef ref(buffer);
165     ref = ref.trim();
166     return ref.str();
167   }
168 
169   // These functions are useful for debugging.
170   template <typename T> T as() const;
171 
172   template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
173   template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
174   template <> long double as<long double>() const {
175     return mpfr_get_ld(value, MPFR_RNDN);
176   }
177 
178   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
179 
180   // Return the ULP (units-in-the-last-place) difference between the
181   // stored MPFR and a floating point number.
182   //
183   // We define:
184   //   ULP(mpfr_value, value) = abs(mpfr_value - value) / eps(value)
185   //
186   // Remarks:
187   // 1. ULP < 0.5 will imply that the value is correctly rounded.
188   // 2. We expect that this value and the value to be compared (the [input]
189   //    argument) are reasonable close, and we will provide an upper bound
190   //    of ULP value for testing.  Morever, most of the fractional parts of
191   //    ULP value do not matter much, so using double as the return type
192   //    should be good enough.
193   template <typename T>
194   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
195     fputil::FPBits<T> bits(input);
196     MPFRNumber mpfrInput(input);
197 
198     // abs(value - input)
199     mpfr_sub(mpfrInput.value, value, mpfrInput.value, MPFR_RNDN);
200     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
201 
202     // get eps(input)
203     int epsExponent = bits.exponent - fputil::FPBits<T>::exponentBias -
204                       fputil::MantissaWidth<T>::value;
205     if (bits.exponent == 0) {
206       // correcting denormal exponent
207       ++epsExponent;
208     } else if ((bits.mantissa == 0) && (bits.exponent > 1) &&
209                mpfr_less_p(value, mpfrInput.value)) {
210       // when the input is exactly 2^n, distance (epsilon) between the input
211       // and the next floating point number is different from the distance to
212       // the previous floating point number.  So in that case, if the correct
213       // value from MPFR is smaller than the input, we use the smaller epsilon
214       --epsExponent;
215     }
216 
217     // Since eps(value) is of the form 2^e, instead of dividing such number,
218     // we multiply by its inverse 2^{-e}.
219     mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
220 
221     return mpfrInput.as<double>();
222   }
223 };
224 
225 namespace internal {
226 
227 template <typename InputType>
228 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
229 unaryOperation(Operation op, InputType input) {
230   MPFRNumber mpfrInput(input);
231   switch (op) {
232   case Operation::Abs:
233     return mpfrInput.abs();
234   case Operation::Ceil:
235     return mpfrInput.ceil();
236   case Operation::Cos:
237     return mpfrInput.cos();
238   case Operation::Exp:
239     return mpfrInput.exp();
240   case Operation::Exp2:
241     return mpfrInput.exp2();
242   case Operation::Floor:
243     return mpfrInput.floor();
244   case Operation::Round:
245     return mpfrInput.round();
246   case Operation::Sin:
247     return mpfrInput.sin();
248   case Operation::Sqrt:
249     return mpfrInput.sqrt();
250   case Operation::Trunc:
251     return mpfrInput.trunc();
252   default:
253     __builtin_unreachable();
254   }
255 }
256 
257 template <typename InputType>
258 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
259 unaryOperationTwoOutputs(Operation op, InputType input, int &output) {
260   MPFRNumber mpfrInput(input);
261   switch (op) {
262   case Operation::Frexp:
263     return mpfrInput.frexp(output);
264   default:
265     __builtin_unreachable();
266   }
267 }
268 
269 template <typename InputType>
270 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
271 binaryOperationTwoOutputs(Operation op, InputType x, InputType y, int &output) {
272   MPFRNumber inputX(x), inputY(y);
273   switch (op) {
274   case Operation::RemQuo:
275     return inputX.remquo(inputY, output);
276   default:
277     __builtin_unreachable();
278   }
279 }
280 
281 template <typename T>
282 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
283                                             testutils::StreamWrapper &OS) {
284   MPFRNumber mpfrInput(input);
285   MPFRNumber mpfrResult = unaryOperation(op, input);
286   MPFRNumber mpfrMatchValue(matchValue);
287   FPBits<T> inputBits(input);
288   FPBits<T> matchBits(matchValue);
289   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
290   OS << "Match value not within tolerance value of MPFR result:\n"
291      << "  Input decimal: " << mpfrInput.str() << '\n';
292   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
293   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
294   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
295                                               OS);
296   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
297   __llvm_libc::fputil::testing::describeValue(
298       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
299   OS << '\n';
300   OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
301      << '\n';
302 }
303 
304 template void
305 explainUnaryOperationSingleOutputError<float>(Operation op, float, float,
306                                               testutils::StreamWrapper &);
307 template void
308 explainUnaryOperationSingleOutputError<double>(Operation op, double, double,
309                                                testutils::StreamWrapper &);
310 template void explainUnaryOperationSingleOutputError<long double>(
311     Operation op, long double, long double, testutils::StreamWrapper &);
312 
313 template <typename T>
314 void explainUnaryOperationTwoOutputsError(Operation op, T input,
315                                           const BinaryOutput<T> &libcResult,
316                                           testutils::StreamWrapper &OS) {
317   MPFRNumber mpfrInput(input);
318   FPBits<T> inputBits(input);
319   int mpfrIntResult;
320   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
321 
322   if (mpfrIntResult != libcResult.i) {
323     OS << "MPFR integral result: " << mpfrIntResult << '\n'
324        << "Libc integral result: " << libcResult.i << '\n';
325   } else {
326     OS << "Integral result from libc matches integral result from MPFR.\n";
327   }
328 
329   MPFRNumber mpfrMatchValue(libcResult.f);
330   OS << "Libc floating point result is not within tolerance value of the MPFR "
331      << "result.\n\n";
332 
333   OS << "            Input decimal: " << mpfrInput.str() << "\n\n";
334 
335   OS << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
336   __llvm_libc::fputil::testing::describeValue(
337       " Libc floating point bits: ", libcResult.f, OS);
338   OS << "\n\n";
339 
340   OS << "              MPFR result: " << mpfrResult.str() << '\n';
341   __llvm_libc::fputil::testing::describeValue(
342       "             MPFR rounded: ", mpfrResult.as<T>(), OS);
343   OS << '\n'
344      << "                ULP error: "
345      << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
346 }
347 
348 template void explainUnaryOperationTwoOutputsError<float>(
349     Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &);
350 template void
351 explainUnaryOperationTwoOutputsError<double>(Operation, double,
352                                              const BinaryOutput<double> &,
353                                              testutils::StreamWrapper &);
354 template void explainUnaryOperationTwoOutputsError<long double>(
355     Operation, long double, const BinaryOutput<long double> &,
356     testutils::StreamWrapper &);
357 
358 template <typename T>
359 void explainBinaryOperationTwoOutputsError(Operation op,
360                                            const BinaryInput<T> &input,
361                                            const BinaryOutput<T> &libcResult,
362                                            testutils::StreamWrapper &OS) {
363   MPFRNumber mpfrX(input.x);
364   MPFRNumber mpfrY(input.y);
365   FPBits<T> xbits(input.x);
366   FPBits<T> ybits(input.y);
367   int mpfrIntResult;
368   MPFRNumber mpfrResult =
369       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
370   MPFRNumber mpfrMatchValue(libcResult.f);
371 
372   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
373      << "MPFR integral result: " << mpfrIntResult << '\n'
374      << "Libc integral result: " << libcResult.i << '\n'
375      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
376      << "               MPFR result: " << mpfrResult.str() << '\n';
377   __llvm_libc::fputil::testing::describeValue(
378       "Libc floating point result bits: ", libcResult.f, OS);
379   __llvm_libc::fputil::testing::describeValue(
380       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
381   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
382 }
383 
384 template void explainBinaryOperationTwoOutputsError<float>(
385     Operation, const BinaryInput<float> &, const BinaryOutput<float> &,
386     testutils::StreamWrapper &);
387 template void explainBinaryOperationTwoOutputsError<double>(
388     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
389     testutils::StreamWrapper &);
390 template void explainBinaryOperationTwoOutputsError<long double>(
391     Operation, const BinaryInput<long double> &,
392     const BinaryOutput<long double> &, testutils::StreamWrapper &);
393 
394 template <typename T>
395 bool compareUnaryOperationSingleOutput(Operation op, T input, T libcResult,
396                                        double ulpError) {
397   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
398   // is rounded to the nearest even.
399   MPFRNumber mpfrResult = unaryOperation(op, input);
400   double ulp = mpfrResult.ulp(libcResult);
401   bool bitsAreEven = ((FPBits<T>(libcResult).bitsAsUInt() & 1) == 0);
402   return (ulp < ulpError) ||
403          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
404 }
405 
406 template bool compareUnaryOperationSingleOutput<float>(Operation, float, float,
407                                                        double);
408 template bool compareUnaryOperationSingleOutput<double>(Operation, double,
409                                                         double, double);
410 template bool compareUnaryOperationSingleOutput<long double>(Operation,
411                                                              long double,
412                                                              long double,
413                                                              double);
414 
415 template <typename T>
416 bool compareUnaryOperationTwoOutputs(Operation op, T input,
417                                      const BinaryOutput<T> &libcResult,
418                                      double ulpError) {
419   int mpfrIntResult;
420   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
421   double ulp = mpfrResult.ulp(libcResult.f);
422 
423   if (mpfrIntResult != libcResult.i)
424     return false;
425 
426   bool bitsAreEven = ((FPBits<T>(libcResult.f).bitsAsUInt() & 1) == 0);
427   return (ulp < ulpError) ||
428          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
429 }
430 
431 template bool
432 compareUnaryOperationTwoOutputs<float>(Operation, float,
433                                        const BinaryOutput<float> &, double);
434 template bool
435 compareUnaryOperationTwoOutputs<double>(Operation, double,
436                                         const BinaryOutput<double> &, double);
437 template bool compareUnaryOperationTwoOutputs<long double>(
438     Operation, long double, const BinaryOutput<long double> &, double);
439 
440 template <typename T>
441 bool compareBinaryOperationTwoOutputs(Operation op, const BinaryInput<T> &input,
442                                       const BinaryOutput<T> &libcResult,
443                                       double ulpError) {
444   int mpfrIntResult;
445   MPFRNumber mpfrResult =
446       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
447   double ulp = mpfrResult.ulp(libcResult.f);
448 
449   if (mpfrIntResult != libcResult.i) {
450     if (op == Operation::RemQuo) {
451       if ((0x7 & mpfrIntResult) != (0x7 & libcResult.i))
452         return false;
453     } else {
454       return false;
455     }
456   }
457 
458   bool bitsAreEven = ((FPBits<T>(libcResult.f).bitsAsUInt() & 1) == 0);
459   return (ulp < ulpError) ||
460          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
461 }
462 
463 template bool
464 compareBinaryOperationTwoOutputs<float>(Operation, const BinaryInput<float> &,
465                                         const BinaryOutput<float> &, double);
466 template bool
467 compareBinaryOperationTwoOutputs<double>(Operation, const BinaryInput<double> &,
468                                          const BinaryOutput<double> &, double);
469 template bool compareBinaryOperationTwoOutputs<long double>(
470     Operation, const BinaryInput<long double> &,
471     const BinaryOutput<long double> &, double);
472 
473 } // namespace internal
474 
475 } // namespace mpfr
476 } // namespace testing
477 } // namespace __llvm_libc
478