1 //===-- Unittests for sqrtf128---------------------------------------------===//
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 "SqrtTest.h"
10
11 #include "src/__support/uint128.h"
12 #include "src/math/sqrtf128.h"
13
14 LIST_SQRT_TESTS(float128, LIBC_NAMESPACE::sqrtf128);
15
TEST_F(LlvmLibcSqrtTest,HardToRound)16 TEST_F(LlvmLibcSqrtTest, HardToRound) {
17 using LIBC_NAMESPACE::fputil::testing::RoundingMode;
18 using FPBits = LIBC_NAMESPACE::fputil::FPBits<float128>;
19
20 // Since there is no exact half cases for square root I encode the
21 // round direction in the sign of the result. E.g. if the number is
22 // negative it means that the exact root is below the rounded value
23 // (the absolute value). Thus I can test not only hard to round
24 // cases for the round to nearest mode but also the directional
25 // modes.
26 float128 HARD_TO_ROUND[][2] = {
27 {0x0.000000dee2f5b6a26c8f07f05442p-16382q,
28 -0x1.ddbd8763a617cff753e2a31083p-8204q},
29 {0x0.000000c86d174c5ad8ae54a548e7p-16382q,
30 0x1.c507bb538940719890851ec1ca88p-8204q},
31 {0x0.000020ab15cfe0b8e488e128f535p-16382q,
32 -0x1.6dccb402560213bc0d62d62e910bp-8201q},
33 {0x0.0000219e97732a9970f2511989bap-16382q,
34 0x1.73163d28be706f4b5052791e28a5p-8201q},
35 {0x0.000026e477546ae99ef57066f9fdp-16382q,
36 -0x1.8f20dd0d0c570a23ea59bc2bf009p-8201q},
37 {0x0.00002d0f88d27a496b3e533f5067p-16382q,
38 0x1.ad9d4abe9f047225a7352bcc52c1p-8201q},
39 {0x1.0000000000000000000000000001p+0q, 0x1p+0q},
40 {0x1.0000000000000000000000000002p+0q,
41 -0x1.0000000000000000000000000001p+0q},
42 {0x1.0000000000000000000000000003p+0q,
43 0x1.0000000000000000000000000001p+0q},
44 {0x1.0000000000000000000000000005p+0q,
45 0x1.0000000000000000000000000002p+0q},
46 {0x1.0000000000000000000000000006p+0q,
47 -0x1.0000000000000000000000000003p+0q},
48 {0x1.1d4c381cbf3a0aa15b9aee344892p+0q,
49 0x1.0e408c3fadc5e64b449c63673f4bp+0q},
50 {0x1.2af17a4ae6f93d11310c49c11b59p+0q,
51 -0x1.14a3bdf0ea5231f12d421a5dbe33p+0q},
52 {0x1.96f893bf29fb91e0fbe19a46d0c8p+0q,
53 0x1.42c6bf6202e66f2295807dee44d9p+0q},
54 {0x1.97fb3839925b66804c429289cce8p+0q,
55 -0x1.432d4049ac1c85a241f333d326e9p+0q},
56 {0x1.be1d900eaeb1533f0f19cc15c7e6p+0q,
57 0x1.51f1715154da44f3bf11f3d96c2dp+0q},
58 {0x1.c4f5074269525063a26051a0ad27p+0q,
59 0x1.54864e9b1daa4d9135ff00663366p+0q},
60 {0x1.035cb5f298a801dc4be9b1f8cd97p+1q,
61 -0x1.6c688775bffcb3f507ba11d0abb9p+0q},
62 {0x1.274be02380427e709beab4dedeb4p+1q,
63 -0x1.84d5763281f2318422392e506b1cp+0q},
64 {0x1.64e797cfdbaa3f7e2f33279dbc6p+1q,
65 0x1.ab79b164e255b26eca00ff99cc99p+0q},
66 {0x1.693a741358c9dac44a570a7e9f6cp+1q,
67 0x1.ae0e8eaeab25bb0c40ee0c2693d3p+0q},
68 {0x1.8275db3fc4d822596047adcb71b9p+1q,
69 -0x1.bcd2bfb653e37a5dbe0ccc2cd917p+0q},
70 {0x1.83280bb98c4a7b88bd6f535899d9p+1q,
71 0x1.bd39409dfd1990dd6a7f8211bb27p+0q},
72 {0x1.d78d8352b48608b510bfd5c75315p+1q,
73 -0x1.eb5c420f15adce0ed2bde5a241cep+0q},
74 {0x1.e3e4774f564b526edff84ce46668p+1q,
75 0x1.f1bf73c0523a19b4bb639c98c0b5p+0q},
76 {0x1.fffffffffffffffffffffffffffap+1q,
77 -0x1.fffffffffffffffffffffffffffdp+0q},
78 {0x1.fffffffffffffffffffffffffffbp+1q,
79 0x1.fffffffffffffffffffffffffffdp+0q},
80 {0x1.fffffffffffffffffffffffffffdp+1q,
81 0x1.fffffffffffffffffffffffffffep+0q},
82 {0x1.fffffffffffffffffffffffffffep+1q,
83 -0x1.ffffffffffffffffffffffffffffp+0q},
84 {0x1.ffffffffffffffffffffffffffffp+1q,
85 0x1.ffffffffffffffffffffffffffffp+0q},
86 };
87
88 auto rnd = [](float128 x, RoundingMode rm) -> float128 {
89 bool is_neg = x < 0;
90 float128 y = is_neg ? -x : x;
91 FPBits ybits(y);
92
93 if (is_neg &&
94 (rm == RoundingMode::Downward || rm == RoundingMode::TowardZero))
95 return FPBits(ybits.uintval() - 1).get_val();
96 if (!is_neg && (rm == RoundingMode::Upward))
97 return FPBits(ybits.uintval() + 1).get_val();
98
99 return y;
100 };
101
102 for (auto &t : HARD_TO_ROUND) {
103 EXPECT_FP_EQ_ALL_ROUNDING(
104 rnd(t[1], RoundingMode::Nearest), rnd(t[1], RoundingMode::Upward),
105 rnd(t[1], RoundingMode::Downward), rnd(t[1], RoundingMode::TowardZero),
106 LIBC_NAMESPACE::sqrtf128(t[0]));
107 }
108
109 // Exact results for subnormal arguments
110 float128 EXACT_SUBNORMAL[][2] = {
111 {0x0.0000000000000000000000000001p-16382q, 0x1p-8247q},
112 {0x0.0000000000000000000000000004p-16382q, 0x1p-8246q},
113 {0x0.0000000000001000000000000000p-16382q, 0x1p-8217q},
114 {0x0.0000000000010000000000000000p-16382q, 0x1p-8215q},
115 {0x0.0000000000100000000000000000p-16382q, 0x1p-8213q},
116 };
117
118 for (auto t : EXACT_SUBNORMAL)
119 EXPECT_FP_EQ_ALL_ROUNDING(t[1], LIBC_NAMESPACE::sqrtf128(t[0]));
120
121 // Check exact cases starting from small numbers
122 for (unsigned k = 1; k < 100 * 100; ++k) {
123 unsigned k2 = k * k;
124 float128 x = static_cast<float128>(k2);
125 float128 y = static_cast<float128>(k);
126 EXPECT_FP_EQ_ALL_ROUNDING(y, LIBC_NAMESPACE::sqrtf128(x));
127 };
128
129 // Then from the largest number.
130 uint64_t k0 = 101904826760412362ULL;
131 for (uint64_t k = k0; k > k0 - 10000; --k) {
132 float128 k_f128 = static_cast<float128>(k);
133 float128 x = k_f128 * k_f128;
134 float128 y = static_cast<float128>(k);
135 EXPECT_FP_EQ_ALL_ROUNDING(y, LIBC_NAMESPACE::sqrtf128(x));
136 }
137 }
138