1//===-- Utils which wrap MPC ----------------------------------------------===//
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 "MPCUtils.h"
10
11#include "src/__support/CPP/array.h"
12#include "src/__support/CPP/stringstream.h"
13#include "utils/MPCWrapper/mpc_inc.h"
14#include "utils/MPFRWrapper/MPCommon.h"
15
16#include <stdint.h>
17
18template <typename T> using FPBits = LIBC_NAMESPACE::fputil::FPBits<T>;
19
20namespace LIBC_NAMESPACE_DECL {
21namespace testing {
22namespace mpc {
23
24static inline cpp::string str(RoundingMode mode) {
25 switch (mode) {
26 case RoundingMode::Upward:
27 return "MPFR_RNDU";
28 case RoundingMode::Downward:
29 return "MPFR_RNDD";
30 case RoundingMode::TowardZero:
31 return "MPFR_RNDZ";
32 case RoundingMode::Nearest:
33 return "MPFR_RNDN";
34 }
35}
36
37class MPCNumber {
38private:
39 unsigned int precision;
40 mpc_t value;
41 mpc_rnd_t mpc_rounding;
42
43public:
44 explicit MPCNumber(unsigned int p) : precision(p), mpc_rounding(MPC_RNDNN) {
45 mpc_init2(value, precision);
46 }
47
48 MPCNumber() : precision(256), mpc_rounding(MPC_RNDNN) {
49 mpc_init2(value, 256);
50 }
51
52 MPCNumber(unsigned int p, mpc_rnd_t rnd) : precision(p), mpc_rounding(rnd) {
53 mpc_init2(value, precision);
54 }
55
56 template <typename XType,
57 cpp::enable_if_t<cpp::is_same_v<_Complex float, XType>, bool> = 0>
58 MPCNumber(XType x,
59 unsigned int precision = mpfr::ExtraPrecision<float>::VALUE,
60 RoundingMode rnd = RoundingMode::Nearest)
61 : precision(precision),
62 mpc_rounding(MPC_RND(mpfr::get_mpfr_rounding_mode(rnd),
63 mpfr::get_mpfr_rounding_mode(rnd))) {
64 mpc_init2(value, precision);
65 Complex<float> x_c = cpp::bit_cast<Complex<float>>(x);
66 mpfr_t real, imag;
67 mpfr_init2(real, precision);
68 mpfr_init2(imag, precision);
69 mpfr_set_flt(real, x_c.real, mpfr::get_mpfr_rounding_mode(rnd));
70 mpfr_set_flt(imag, x_c.imag, mpfr::get_mpfr_rounding_mode(rnd));
71 mpc_set_fr_fr(value, real, imag, mpc_rounding);
72 mpfr_clear(real);
73 mpfr_clear(imag);
74 }
75
76 template <typename XType,
77 cpp::enable_if_t<cpp::is_same_v<_Complex double, XType>, bool> = 0>
78 MPCNumber(XType x,
79 unsigned int precision = mpfr::ExtraPrecision<double>::VALUE,
80 RoundingMode rnd = RoundingMode::Nearest)
81 : precision(precision),
82 mpc_rounding(MPC_RND(mpfr::get_mpfr_rounding_mode(rnd),
83 mpfr::get_mpfr_rounding_mode(rnd))) {
84 mpc_init2(value, precision);
85 Complex<double> x_c = cpp::bit_cast<Complex<double>>(x);
86 mpc_set_d_d(value, x_c.real, x_c.imag, mpc_rounding);
87 }
88
89 MPCNumber(const MPCNumber &other)
90 : precision(other.precision), mpc_rounding(other.mpc_rounding) {
91 mpc_init2(value, precision);
92 mpc_set(value, other.value, mpc_rounding);
93 }
94
95 ~MPCNumber() { mpc_clear(value); }
96
97 MPCNumber &operator=(const MPCNumber &rhs) {
98 precision = rhs.precision;
99 mpc_rounding = rhs.mpc_rounding;
100 mpc_init2(value, precision);
101 mpc_set(value, rhs.value, mpc_rounding);
102 return *this;
103 }
104
105 void setValue(mpc_t val) const { mpc_set(val, value, mpc_rounding); }
106
107 mpc_t &getValue() { return value; }
108
109 MPCNumber carg() const {
110 mpfr_t res;
111 MPCNumber result(precision, mpc_rounding);
112
113 mpfr_init2(res, precision);
114
115 mpc_arg(res, value, MPC_RND_RE(mpc_rounding));
116 mpc_set_fr(result.value, res, mpc_rounding);
117
118 mpfr_clear(res);
119
120 return result;
121 }
122
123 MPCNumber cproj() const {
124 MPCNumber result(precision, mpc_rounding);
125 mpc_proj(result.value, value, mpc_rounding);
126 return result;
127 }
128};
129
130namespace internal {
131
132template <typename InputType>
133cpp::enable_if_t<cpp::is_complex_v<InputType>, MPCNumber>
134unary_operation(Operation op, InputType input, unsigned int precision,
135 RoundingMode rounding) {
136 MPCNumber mpcInput(input, precision, rounding);
137 switch (op) {
138 case Operation::Carg:
139 return mpcInput.carg();
140 case Operation::Cproj:
141 return mpcInput.cproj();
142 default:
143 __builtin_unreachable();
144 }
145}
146
147template <typename InputType, typename OutputType>
148bool compare_unary_operation_single_output_same_type(Operation op,
149 InputType input,
150 OutputType libc_result,
151 double ulp_tolerance,
152 RoundingMode rounding) {
153
154 unsigned int precision =
155 mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance);
156
157 MPCNumber mpc_result;
158 mpc_result = unary_operation(op, input, precision, rounding);
159
160 mpc_t mpc_result_val;
161 mpc_init2(mpc_result_val, precision);
162 mpc_result.setValue(mpc_result_val);
163
164 mpfr_t real, imag;
165 mpfr_init2(real, precision);
166 mpfr_init2(imag, precision);
167 mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));
168 mpc_imag(imag, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));
169
170 mpfr::MPFRNumber mpfr_real(real, precision, rounding);
171 mpfr::MPFRNumber mpfr_imag(imag, precision, rounding);
172
173 double ulp_real = mpfr_real.ulp(
174 (cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)).real);
175 double ulp_imag = mpfr_imag.ulp(
176 (cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)).imag);
177 mpc_clear(mpc_result_val);
178 mpfr_clear(real);
179 mpfr_clear(imag);
180 return (ulp_real <= ulp_tolerance) && (ulp_imag <= ulp_tolerance);
181}
182
183template bool compare_unary_operation_single_output_same_type(
184 Operation, _Complex float, _Complex float, double, RoundingMode);
185template bool compare_unary_operation_single_output_same_type(
186 Operation, _Complex double, _Complex double, double, RoundingMode);
187
188template <typename InputType, typename OutputType>
189bool compare_unary_operation_single_output_different_type(
190 Operation op, InputType input, OutputType libc_result, double ulp_tolerance,
191 RoundingMode rounding) {
192
193 unsigned int precision =
194 mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance);
195
196 MPCNumber mpc_result;
197 mpc_result = unary_operation(op, input, precision, rounding);
198
199 mpc_t mpc_result_val;
200 mpc_init2(mpc_result_val, precision);
201 mpc_result.setValue(mpc_result_val);
202
203 mpfr_t real;
204 mpfr_init2(real, precision);
205 mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));
206
207 mpfr::MPFRNumber mpfr_real(real, precision, rounding);
208
209 double ulp_real = mpfr_real.ulp(libc_result);
210 mpc_clear(mpc_result_val);
211 mpfr_clear(real);
212 return (ulp_real <= ulp_tolerance);
213}
214
215template bool compare_unary_operation_single_output_different_type(
216 Operation, _Complex float, float, double, RoundingMode);
217template bool compare_unary_operation_single_output_different_type(
218 Operation, _Complex double, double, double, RoundingMode);
219
220template <typename InputType, typename OutputType>
221void explain_unary_operation_single_output_different_type_error(
222 Operation op, InputType input, OutputType libc_result, double ulp_tolerance,
223 RoundingMode rounding) {
224
225 unsigned int precision =
226 mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance);
227
228 MPCNumber mpc_result;
229 mpc_result = unary_operation(op, input, precision, rounding);
230
231 mpc_t mpc_result_val;
232 mpc_init2(mpc_result_val, precision);
233 mpc_result.setValue(mpc_result_val);
234
235 mpfr_t real;
236 mpfr_init2(real, precision);
237 mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));
238
239 mpfr::MPFRNumber mpfr_result(real, precision, rounding);
240 mpfr::MPFRNumber mpfrLibcResult(libc_result, precision, rounding);
241 mpfr::MPFRNumber mpfrInputReal(
242 cpp::bit_cast<Complex<make_real_t<InputType>>>(input).real, precision,
243 rounding);
244 mpfr::MPFRNumber mpfrInputImag(
245 cpp::bit_cast<Complex<make_real_t<InputType>>>(input).imag, precision,
246 rounding);
247
248 cpp::array<char, 2048> msg_buf;
249 cpp::StringStream msg(msg_buf);
250 msg << "Match value not within tolerance value of MPFR result:\n"
251 << " Input: " << mpfrInputReal.str() << " + " << mpfrInputImag.str()
252 << "i\n"
253 << " Rounding mode: " << str(rounding) << '\n'
254 << " Libc: " << mpfrLibcResult.str() << '\n'
255 << " MPC: " << mpfr_result.str() << '\n'
256 << '\n'
257 << " ULP error: " << mpfr_result.ulp_as_mpfr_number(libc_result).str()
258 << '\n';
259 tlog << msg.str();
260 mpc_clear(mpc_result_val);
261 mpfr_clear(real);
262}
263
264template void explain_unary_operation_single_output_different_type_error(
265 Operation, _Complex float, float, double, RoundingMode);
266template void explain_unary_operation_single_output_different_type_error(
267 Operation, _Complex double, double, double, RoundingMode);
268
269template <typename InputType, typename OutputType>
270void explain_unary_operation_single_output_same_type_error(
271 Operation op, InputType input, OutputType libc_result, double ulp_tolerance,
272 RoundingMode rounding) {
273
274 unsigned int precision =
275 mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance);
276
277 MPCNumber mpc_result;
278 mpc_result = unary_operation(op, input, precision, rounding);
279
280 mpc_t mpc_result_val;
281 mpc_init2(mpc_result_val, precision);
282 mpc_result.setValue(mpc_result_val);
283
284 mpfr_t real, imag;
285 mpfr_init2(real, precision);
286 mpfr_init2(imag, precision);
287 mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));
288 mpc_imag(imag, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));
289
290 mpfr::MPFRNumber mpfr_real(real, precision, rounding);
291 mpfr::MPFRNumber mpfr_imag(imag, precision, rounding);
292 mpfr::MPFRNumber mpfrLibcResultReal(
293 cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result).real,
294 precision, rounding);
295 mpfr::MPFRNumber mpfrLibcResultImag(
296 cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result).imag,
297 precision, rounding);
298 mpfr::MPFRNumber mpfrInputReal(
299 cpp::bit_cast<Complex<make_real_t<InputType>>>(input).real, precision,
300 rounding);
301 mpfr::MPFRNumber mpfrInputImag(
302 cpp::bit_cast<Complex<make_real_t<InputType>>>(input).imag, precision,
303 rounding);
304
305 cpp::array<char, 2048> msg_buf;
306 cpp::StringStream msg(msg_buf);
307 msg << "Match value not within tolerance value of MPFR result:\n"
308 << " Input: " << mpfrInputReal.str() << " + " << mpfrInputImag.str()
309 << "i\n"
310 << " Rounding mode: " << str(rounding) << " , " << str(rounding) << '\n'
311 << " Libc: " << mpfrLibcResultReal.str() << " + "
312 << mpfrLibcResultImag.str() << "i\n"
313 << " MPC: " << mpfr_real.str() << " + " << mpfr_imag.str() << "i\n"
314 << '\n'
315 << " ULP error: "
316 << mpfr_real
317 .ulp_as_mpfr_number(
318 cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)
319 .real)
320 .str()
321 << " , "
322 << mpfr_imag
323 .ulp_as_mpfr_number(
324 cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)
325 .imag)
326 .str()
327 << '\n';
328 tlog << msg.str();
329 mpc_clear(mpc_result_val);
330 mpfr_clear(real);
331 mpfr_clear(imag);
332}
333
334template void explain_unary_operation_single_output_same_type_error(
335 Operation, _Complex float, _Complex float, double, RoundingMode);
336template void explain_unary_operation_single_output_same_type_error(
337 Operation, _Complex double, _Complex double, double, RoundingMode);
338
339} // namespace internal
340
341} // namespace mpc
342} // namespace testing
343} // namespace LIBC_NAMESPACE_DECL
344

source code of libc/utils/MPCWrapper/MPCUtils.cpp