1 //===-- lib/Evaluate/real.cpp ---------------------------------------------===//
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 "flang/Evaluate/real.h"
10 #include "int-power.h"
11 #include "flang/Common/idioms.h"
12 #include "flang/Decimal/decimal.h"
13 #include "flang/Parser/characters.h"
14 #include "llvm/Support/raw_ostream.h"
15 #include <limits>
16
17 namespace Fortran::evaluate::value {
18
Compare(const Real & y) const19 template <typename W, int P> Relation Real<W, P>::Compare(const Real &y) const {
20 if (IsNotANumber() || y.IsNotANumber()) { // NaN vs x, x vs NaN
21 return Relation::Unordered;
22 } else if (IsInfinite()) {
23 if (y.IsInfinite()) {
24 if (IsNegative()) { // -Inf vs +/-Inf
25 return y.IsNegative() ? Relation::Equal : Relation::Less;
26 } else { // +Inf vs +/-Inf
27 return y.IsNegative() ? Relation::Greater : Relation::Equal;
28 }
29 } else { // +/-Inf vs finite
30 return IsNegative() ? Relation::Less : Relation::Greater;
31 }
32 } else if (y.IsInfinite()) { // finite vs +/-Inf
33 return y.IsNegative() ? Relation::Greater : Relation::Less;
34 } else { // two finite numbers
35 bool isNegative{IsNegative()};
36 if (isNegative != y.IsNegative()) {
37 if (word_.IOR(y.word_).IBCLR(bits - 1).IsZero()) {
38 return Relation::Equal; // +/-0.0 == -/+0.0
39 } else {
40 return isNegative ? Relation::Less : Relation::Greater;
41 }
42 } else {
43 // same sign
44 Ordering order{evaluate::Compare(Exponent(), y.Exponent())};
45 if (order == Ordering::Equal) {
46 order = GetSignificand().CompareUnsigned(y.GetSignificand());
47 }
48 if (isNegative) {
49 order = Reverse(order);
50 }
51 return RelationFromOrdering(order);
52 }
53 }
54 }
55
56 template <typename W, int P>
Add(const Real & y,Rounding rounding) const57 ValueWithRealFlags<Real<W, P>> Real<W, P>::Add(
58 const Real &y, Rounding rounding) const {
59 ValueWithRealFlags<Real> result;
60 if (IsNotANumber() || y.IsNotANumber()) {
61 result.value = NotANumber(); // NaN + x -> NaN
62 if (IsSignalingNaN() || y.IsSignalingNaN()) {
63 result.flags.set(RealFlag::InvalidArgument);
64 }
65 return result;
66 }
67 bool isNegative{IsNegative()};
68 bool yIsNegative{y.IsNegative()};
69 if (IsInfinite()) {
70 if (y.IsInfinite()) {
71 if (isNegative == yIsNegative) {
72 result.value = *this; // +/-Inf + +/-Inf -> +/-Inf
73 } else {
74 result.value = NotANumber(); // +/-Inf + -/+Inf -> NaN
75 result.flags.set(RealFlag::InvalidArgument);
76 }
77 } else {
78 result.value = *this; // +/-Inf + x -> +/-Inf
79 }
80 return result;
81 }
82 if (y.IsInfinite()) {
83 result.value = y; // x + +/-Inf -> +/-Inf
84 return result;
85 }
86 int exponent{Exponent()};
87 int yExponent{y.Exponent()};
88 if (exponent < yExponent) {
89 // y is larger in magnitude; simplify by reversing operands
90 return y.Add(*this, rounding);
91 }
92 if (exponent == yExponent && isNegative != yIsNegative) {
93 Ordering order{GetSignificand().CompareUnsigned(y.GetSignificand())};
94 if (order == Ordering::Less) {
95 // Same exponent, opposite signs, and y is larger in magnitude
96 return y.Add(*this, rounding);
97 }
98 if (order == Ordering::Equal) {
99 // x + (-x) -> +0.0 unless rounding is directed downwards
100 if (rounding.mode == common::RoundingMode::Down) {
101 result.value.word_ = result.value.word_.IBSET(bits - 1); // -0.0
102 }
103 return result;
104 }
105 }
106 // Our exponent is greater than y's, or the exponents match and y is not
107 // of the opposite sign and greater magnitude. So (x+y) will have the
108 // same sign as x.
109 Fraction fraction{GetFraction()};
110 Fraction yFraction{y.GetFraction()};
111 int rshift = exponent - yExponent;
112 if (exponent > 0 && yExponent == 0) {
113 --rshift; // correct overshift when only y is subnormal
114 }
115 RoundingBits roundingBits{yFraction, rshift};
116 yFraction = yFraction.SHIFTR(rshift);
117 bool carry{false};
118 if (isNegative != yIsNegative) {
119 // Opposite signs: subtract via addition of two's complement of y and
120 // the rounding bits.
121 yFraction = yFraction.NOT();
122 carry = roundingBits.Negate();
123 }
124 auto sum{fraction.AddUnsigned(yFraction, carry)};
125 fraction = sum.value;
126 if (isNegative == yIsNegative && sum.carry) {
127 roundingBits.ShiftRight(sum.value.BTEST(0));
128 fraction = fraction.SHIFTR(1).IBSET(fraction.bits - 1);
129 ++exponent;
130 }
131 NormalizeAndRound(
132 result, isNegative, exponent, fraction, rounding, roundingBits);
133 return result;
134 }
135
136 template <typename W, int P>
Multiply(const Real & y,Rounding rounding) const137 ValueWithRealFlags<Real<W, P>> Real<W, P>::Multiply(
138 const Real &y, Rounding rounding) const {
139 ValueWithRealFlags<Real> result;
140 if (IsNotANumber() || y.IsNotANumber()) {
141 result.value = NotANumber(); // NaN * x -> NaN
142 if (IsSignalingNaN() || y.IsSignalingNaN()) {
143 result.flags.set(RealFlag::InvalidArgument);
144 }
145 } else {
146 bool isNegative{IsNegative() != y.IsNegative()};
147 if (IsInfinite() || y.IsInfinite()) {
148 if (IsZero() || y.IsZero()) {
149 result.value = NotANumber(); // 0 * Inf -> NaN
150 result.flags.set(RealFlag::InvalidArgument);
151 } else {
152 result.value = Infinity(isNegative);
153 }
154 } else {
155 auto product{GetFraction().MultiplyUnsigned(y.GetFraction())};
156 std::int64_t exponent{CombineExponents(y, false)};
157 if (exponent < 1) {
158 int rshift = 1 - exponent;
159 exponent = 1;
160 bool sticky{false};
161 if (rshift >= product.upper.bits + product.lower.bits) {
162 sticky = !product.lower.IsZero() || !product.upper.IsZero();
163 } else if (rshift >= product.lower.bits) {
164 sticky = !product.lower.IsZero() ||
165 !product.upper
166 .IAND(product.upper.MASKR(rshift - product.lower.bits))
167 .IsZero();
168 } else {
169 sticky = !product.lower.IAND(product.lower.MASKR(rshift)).IsZero();
170 }
171 product.lower = product.lower.SHIFTRWithFill(product.upper, rshift);
172 product.upper = product.upper.SHIFTR(rshift);
173 if (sticky) {
174 product.lower = product.lower.IBSET(0);
175 }
176 }
177 int leadz{product.upper.LEADZ()};
178 if (leadz >= product.upper.bits) {
179 leadz += product.lower.LEADZ();
180 }
181 int lshift{leadz};
182 if (lshift > exponent - 1) {
183 lshift = exponent - 1;
184 }
185 exponent -= lshift;
186 product.upper = product.upper.SHIFTLWithFill(product.lower, lshift);
187 product.lower = product.lower.SHIFTL(lshift);
188 RoundingBits roundingBits{product.lower, product.lower.bits};
189 NormalizeAndRound(result, isNegative, exponent, product.upper, rounding,
190 roundingBits, true /*multiply*/);
191 }
192 }
193 return result;
194 }
195
196 template <typename W, int P>
Divide(const Real & y,Rounding rounding) const197 ValueWithRealFlags<Real<W, P>> Real<W, P>::Divide(
198 const Real &y, Rounding rounding) const {
199 ValueWithRealFlags<Real> result;
200 if (IsNotANumber() || y.IsNotANumber()) {
201 result.value = NotANumber(); // NaN / x -> NaN, x / NaN -> NaN
202 if (IsSignalingNaN() || y.IsSignalingNaN()) {
203 result.flags.set(RealFlag::InvalidArgument);
204 }
205 } else {
206 bool isNegative{IsNegative() != y.IsNegative()};
207 if (IsInfinite()) {
208 if (y.IsInfinite()) {
209 result.value = NotANumber(); // Inf/Inf -> NaN
210 result.flags.set(RealFlag::InvalidArgument);
211 } else { // Inf/x -> Inf, Inf/0 -> Inf
212 result.value = Infinity(isNegative);
213 }
214 } else if (y.IsZero()) {
215 if (IsZero()) { // 0/0 -> NaN
216 result.value = NotANumber();
217 result.flags.set(RealFlag::InvalidArgument);
218 } else { // x/0 -> Inf, Inf/0 -> Inf
219 result.value = Infinity(isNegative);
220 result.flags.set(RealFlag::DivideByZero);
221 }
222 } else if (IsZero() || y.IsInfinite()) { // 0/x, x/Inf -> 0
223 if (isNegative) {
224 result.value.word_ = result.value.word_.IBSET(bits - 1);
225 }
226 } else {
227 // dividend and divisor are both finite and nonzero numbers
228 Fraction top{GetFraction()}, divisor{y.GetFraction()};
229 std::int64_t exponent{CombineExponents(y, true)};
230 Fraction quotient;
231 bool msb{false};
232 if (!top.BTEST(top.bits - 1) || !divisor.BTEST(divisor.bits - 1)) {
233 // One or two subnormals
234 int topLshift{top.LEADZ()};
235 top = top.SHIFTL(topLshift);
236 int divisorLshift{divisor.LEADZ()};
237 divisor = divisor.SHIFTL(divisorLshift);
238 exponent += divisorLshift - topLshift;
239 }
240 for (int j{1}; j <= quotient.bits; ++j) {
241 if (NextQuotientBit(top, msb, divisor)) {
242 quotient = quotient.IBSET(quotient.bits - j);
243 }
244 }
245 bool guard{NextQuotientBit(top, msb, divisor)};
246 bool round{NextQuotientBit(top, msb, divisor)};
247 bool sticky{msb || !top.IsZero()};
248 RoundingBits roundingBits{guard, round, sticky};
249 if (exponent < 1) {
250 std::int64_t rshift{1 - exponent};
251 for (; rshift > 0; --rshift) {
252 roundingBits.ShiftRight(quotient.BTEST(0));
253 quotient = quotient.SHIFTR(1);
254 }
255 exponent = 1;
256 }
257 NormalizeAndRound(
258 result, isNegative, exponent, quotient, rounding, roundingBits);
259 }
260 }
261 return result;
262 }
263
264 template <typename W, int P>
SQRT(Rounding rounding) const265 ValueWithRealFlags<Real<W, P>> Real<W, P>::SQRT(Rounding rounding) const {
266 ValueWithRealFlags<Real> result;
267 if (IsNotANumber()) {
268 result.value = NotANumber();
269 if (IsSignalingNaN()) {
270 result.flags.set(RealFlag::InvalidArgument);
271 }
272 } else if (IsNegative()) {
273 if (IsZero()) {
274 // SQRT(-0) == -0 in IEEE-754.
275 result.value.word_ = result.value.word_.IBSET(bits - 1);
276 } else {
277 result.value = NotANumber();
278 }
279 } else if (IsInfinite()) {
280 // SQRT(+Inf) == +Inf
281 result.value = Infinity(false);
282 } else {
283 int expo{UnbiasedExponent()};
284 if (expo < -1 || expo > 1) {
285 // Reduce the range to [0.5 .. 4.0) by dividing by an integral power
286 // of four to avoid trouble with very large and very small values
287 // (esp. truncation of subnormals).
288 // SQRT(2**(2a) * x) = SQRT(2**(2a)) * SQRT(x) = 2**a * SQRT(x)
289 Real scaled;
290 int adjust{expo / 2};
291 scaled.Normalize(false, expo - 2 * adjust + exponentBias, GetFraction());
292 result = scaled.SQRT(rounding);
293 result.value.Normalize(false,
294 result.value.UnbiasedExponent() + adjust + exponentBias,
295 result.value.GetFraction());
296 return result;
297 }
298 // Compute the square root of the reduced value with the slow but
299 // reliable bit-at-a-time method. Start with a clear significand and
300 // half of the unbiased exponent, and then try to set significand bits
301 // in descending order of magnitude without exceeding the exact result.
302 expo = expo / 2 + exponentBias;
303 result.value.Normalize(false, expo, Fraction::MASKL(1));
304 Real initialSq{result.value.Multiply(result.value).value};
305 if (Compare(initialSq) == Relation::Less) {
306 // Initial estimate is too large; this can happen for values just
307 // under 1.0.
308 --expo;
309 result.value.Normalize(false, expo, Fraction::MASKL(1));
310 }
311 for (int bit{significandBits - 1}; bit >= 0; --bit) {
312 Word word{result.value.word_};
313 result.value.word_ = word.IBSET(bit);
314 auto squared{result.value.Multiply(result.value, rounding)};
315 if (squared.flags.test(RealFlag::Overflow) ||
316 squared.flags.test(RealFlag::Underflow) ||
317 Compare(squared.value) == Relation::Less) {
318 result.value.word_ = word;
319 }
320 }
321 // The computed square root has a square that's not greater than the
322 // original argument. Check this square against the square of the next
323 // larger Real and return that one if its square is closer in magnitude to
324 // the original argument.
325 Real resultSq{result.value.Multiply(result.value).value};
326 Real diff{Subtract(resultSq).value.ABS()};
327 if (diff.IsZero()) {
328 return result; // exact
329 }
330 Real ulp;
331 ulp.Normalize(false, expo, Fraction::MASKR(1));
332 Real nextAfter{result.value.Add(ulp).value};
333 auto nextAfterSq{nextAfter.Multiply(nextAfter)};
334 if (!nextAfterSq.flags.test(RealFlag::Overflow) &&
335 !nextAfterSq.flags.test(RealFlag::Underflow)) {
336 Real nextAfterDiff{Subtract(nextAfterSq.value).value.ABS()};
337 if (nextAfterDiff.Compare(diff) == Relation::Less) {
338 result.value = nextAfter;
339 if (nextAfterDiff.IsZero()) {
340 return result; // exact
341 }
342 }
343 }
344 result.flags.set(RealFlag::Inexact);
345 }
346 return result;
347 }
348
349 // HYPOT(x,y) = SQRT(x**2 + y**2) by definition, but those squared intermediate
350 // values are susceptible to over/underflow when computed naively.
351 // Assuming that x>=y, calculate instead:
352 // HYPOT(x,y) = SQRT(x**2 * (1+(y/x)**2))
353 // = ABS(x) * SQRT(1+(y/x)**2)
354 template <typename W, int P>
HYPOT(const Real & y,Rounding rounding) const355 ValueWithRealFlags<Real<W, P>> Real<W, P>::HYPOT(
356 const Real &y, Rounding rounding) const {
357 ValueWithRealFlags<Real> result;
358 if (IsNotANumber() || y.IsNotANumber()) {
359 result.flags.set(RealFlag::InvalidArgument);
360 result.value = NotANumber();
361 } else if (ABS().Compare(y.ABS()) == Relation::Less) {
362 return y.HYPOT(*this);
363 } else if (IsZero()) {
364 return result; // x==y==0
365 } else {
366 auto yOverX{y.Divide(*this, rounding)}; // y/x
367 bool inexact{yOverX.flags.test(RealFlag::Inexact)};
368 auto squared{yOverX.value.Multiply(yOverX.value, rounding)}; // (y/x)**2
369 inexact |= squared.flags.test(RealFlag::Inexact);
370 Real one;
371 one.Normalize(false, exponentBias, Fraction::MASKL(1)); // 1.0
372 auto sum{squared.value.Add(one, rounding)}; // 1.0 + (y/x)**2
373 inexact |= sum.flags.test(RealFlag::Inexact);
374 auto sqrt{sum.value.SQRT()};
375 inexact |= sqrt.flags.test(RealFlag::Inexact);
376 result = sqrt.value.Multiply(ABS(), rounding);
377 if (inexact) {
378 result.flags.set(RealFlag::Inexact);
379 }
380 }
381 return result;
382 }
383
384 template <typename W, int P>
ToWholeNumber(common::RoundingMode mode) const385 ValueWithRealFlags<Real<W, P>> Real<W, P>::ToWholeNumber(
386 common::RoundingMode mode) const {
387 ValueWithRealFlags<Real> result{*this};
388 if (IsNotANumber()) {
389 result.flags.set(RealFlag::InvalidArgument);
390 result.value = NotANumber();
391 } else if (IsInfinite()) {
392 result.flags.set(RealFlag::Overflow);
393 } else {
394 constexpr int noClipExponent{exponentBias + binaryPrecision - 1};
395 if (Exponent() < noClipExponent) {
396 Real adjust; // ABS(EPSILON(adjust)) == 0.5
397 adjust.Normalize(IsSignBitSet(), noClipExponent, Fraction::MASKL(1));
398 // Compute ival=(*this + adjust), losing any fractional bits; keep flags
399 result = Add(adjust, Rounding{mode});
400 result.flags.reset(RealFlag::Inexact); // result *is* exact
401 // Return (ival-adjust) with original sign in case we've generated a zero.
402 result.value =
403 result.value.Subtract(adjust, Rounding{common::RoundingMode::ToZero})
404 .value.SIGN(*this);
405 }
406 }
407 return result;
408 }
409
410 template <typename W, int P>
Normalize(bool negative,int exponent,const Fraction & fraction,Rounding rounding,RoundingBits * roundingBits)411 RealFlags Real<W, P>::Normalize(bool negative, int exponent,
412 const Fraction &fraction, Rounding rounding, RoundingBits *roundingBits) {
413 int lshift{fraction.LEADZ()};
414 if (lshift == fraction.bits /* fraction is zero */ &&
415 (!roundingBits || roundingBits->empty())) {
416 // No fraction, no rounding bits -> +/-0.0
417 exponent = lshift = 0;
418 } else if (lshift < exponent) {
419 exponent -= lshift;
420 } else if (exponent > 0) {
421 lshift = exponent - 1;
422 exponent = 0;
423 } else if (lshift == 0) {
424 exponent = 1;
425 } else {
426 lshift = 0;
427 }
428 if (exponent >= maxExponent) {
429 // Infinity or overflow
430 if (rounding.mode == common::RoundingMode::TiesToEven ||
431 rounding.mode == common::RoundingMode::TiesAwayFromZero ||
432 (rounding.mode == common::RoundingMode::Up && !negative) ||
433 (rounding.mode == common::RoundingMode::Down && negative)) {
434 word_ = Word{maxExponent}.SHIFTL(significandBits); // Inf
435 } else {
436 // directed rounding: round to largest finite value rather than infinity
437 // (x86 does this, not sure whether it's standard behavior)
438 word_ = Word{word_.MASKR(word_.bits - 1)}.IBCLR(significandBits);
439 }
440 if (negative) {
441 word_ = word_.IBSET(bits - 1);
442 }
443 RealFlags flags{RealFlag::Overflow};
444 if (!fraction.IsZero()) {
445 flags.set(RealFlag::Inexact);
446 }
447 return flags;
448 }
449 word_ = Word::ConvertUnsigned(fraction).value;
450 if (lshift > 0) {
451 word_ = word_.SHIFTL(lshift);
452 if (roundingBits) {
453 for (; lshift > 0; --lshift) {
454 if (roundingBits->ShiftLeft()) {
455 word_ = word_.IBSET(lshift - 1);
456 }
457 }
458 }
459 }
460 if constexpr (isImplicitMSB) {
461 word_ = word_.IBCLR(significandBits);
462 }
463 word_ = word_.IOR(Word{exponent}.SHIFTL(significandBits));
464 if (negative) {
465 word_ = word_.IBSET(bits - 1);
466 }
467 return {};
468 }
469
470 template <typename W, int P>
Round(Rounding rounding,const RoundingBits & bits,bool multiply)471 RealFlags Real<W, P>::Round(
472 Rounding rounding, const RoundingBits &bits, bool multiply) {
473 int origExponent{Exponent()};
474 RealFlags flags;
475 bool inexact{!bits.empty()};
476 if (inexact) {
477 flags.set(RealFlag::Inexact);
478 }
479 if (origExponent < maxExponent &&
480 bits.MustRound(rounding, IsNegative(), word_.BTEST(0) /* is odd */)) {
481 typename Fraction::ValueWithCarry sum{
482 GetFraction().AddUnsigned(Fraction{}, true)};
483 int newExponent{origExponent};
484 if (sum.carry) {
485 // The fraction was all ones before rounding; sum.value is now zero
486 sum.value = sum.value.IBSET(binaryPrecision - 1);
487 if (++newExponent >= maxExponent) {
488 flags.set(RealFlag::Overflow); // rounded away to an infinity
489 }
490 }
491 flags |= Normalize(IsNegative(), newExponent, sum.value);
492 }
493 if (inexact && origExponent == 0) {
494 // inexact subnormal input: signal Underflow unless in an x86-specific
495 // edge case
496 if (rounding.x86CompatibleBehavior && Exponent() != 0 && multiply &&
497 bits.sticky() &&
498 (bits.guard() ||
499 (rounding.mode != common::RoundingMode::Up &&
500 rounding.mode != common::RoundingMode::Down))) {
501 // x86 edge case in which Underflow fails to signal when a subnormal
502 // inexact multiplication product rounds to a normal result when
503 // the guard bit is set or we're not using directed rounding
504 } else {
505 flags.set(RealFlag::Underflow);
506 }
507 }
508 return flags;
509 }
510
511 template <typename W, int P>
NormalizeAndRound(ValueWithRealFlags<Real> & result,bool isNegative,int exponent,const Fraction & fraction,Rounding rounding,RoundingBits roundingBits,bool multiply)512 void Real<W, P>::NormalizeAndRound(ValueWithRealFlags<Real> &result,
513 bool isNegative, int exponent, const Fraction &fraction, Rounding rounding,
514 RoundingBits roundingBits, bool multiply) {
515 result.flags |= result.value.Normalize(
516 isNegative, exponent, fraction, rounding, &roundingBits);
517 result.flags |= result.value.Round(rounding, roundingBits, multiply);
518 }
519
MapRoundingMode(common::RoundingMode rounding)520 inline enum decimal::FortranRounding MapRoundingMode(
521 common::RoundingMode rounding) {
522 switch (rounding) {
523 case common::RoundingMode::TiesToEven:
524 break;
525 case common::RoundingMode::ToZero:
526 return decimal::RoundToZero;
527 case common::RoundingMode::Down:
528 return decimal::RoundDown;
529 case common::RoundingMode::Up:
530 return decimal::RoundUp;
531 case common::RoundingMode::TiesAwayFromZero:
532 return decimal::RoundCompatible;
533 }
534 return decimal::RoundNearest; // dodge gcc warning about lack of result
535 }
536
MapFlags(decimal::ConversionResultFlags flags)537 inline RealFlags MapFlags(decimal::ConversionResultFlags flags) {
538 RealFlags result;
539 if (flags & decimal::Overflow) {
540 result.set(RealFlag::Overflow);
541 }
542 if (flags & decimal::Inexact) {
543 result.set(RealFlag::Inexact);
544 }
545 if (flags & decimal::Invalid) {
546 result.set(RealFlag::InvalidArgument);
547 }
548 return result;
549 }
550
551 template <typename W, int P>
Read(const char * & p,Rounding rounding)552 ValueWithRealFlags<Real<W, P>> Real<W, P>::Read(
553 const char *&p, Rounding rounding) {
554 auto converted{
555 decimal::ConvertToBinary<P>(p, MapRoundingMode(rounding.mode))};
556 const auto *value{reinterpret_cast<Real<W, P> *>(&converted.binary)};
557 return {*value, MapFlags(converted.flags)};
558 }
559
DumpHexadecimal() const560 template <typename W, int P> std::string Real<W, P>::DumpHexadecimal() const {
561 if (IsNotANumber()) {
562 return "NaN0x"s + word_.Hexadecimal();
563 } else if (IsNegative()) {
564 return "-"s + Negate().DumpHexadecimal();
565 } else if (IsInfinite()) {
566 return "Inf"s;
567 } else if (IsZero()) {
568 return "0.0"s;
569 } else {
570 Fraction frac{GetFraction()};
571 std::string result{"0x"};
572 char intPart = '0' + frac.BTEST(frac.bits - 1);
573 result += intPart;
574 result += '.';
575 int trailz{frac.TRAILZ()};
576 if (trailz >= frac.bits - 1) {
577 result += '0';
578 } else {
579 int remainingBits{frac.bits - 1 - trailz};
580 int wholeNybbles{remainingBits / 4};
581 int lostBits{remainingBits - 4 * wholeNybbles};
582 if (wholeNybbles > 0) {
583 std::string fracHex{frac.SHIFTR(trailz + lostBits)
584 .IAND(frac.MASKR(4 * wholeNybbles))
585 .Hexadecimal()};
586 std::size_t field = wholeNybbles;
587 if (fracHex.size() < field) {
588 result += std::string(field - fracHex.size(), '0');
589 }
590 result += fracHex;
591 }
592 if (lostBits > 0) {
593 result += frac.SHIFTR(trailz)
594 .IAND(frac.MASKR(lostBits))
595 .SHIFTL(4 - lostBits)
596 .Hexadecimal();
597 }
598 }
599 result += 'p';
600 int exponent = Exponent() - exponentBias;
601 result += Integer<32>{exponent}.SignedDecimal();
602 return result;
603 }
604 }
605
606 template <typename W, int P>
AsFortran(llvm::raw_ostream & o,int kind,bool minimal) const607 llvm::raw_ostream &Real<W, P>::AsFortran(
608 llvm::raw_ostream &o, int kind, bool minimal) const {
609 if (IsNotANumber()) {
610 o << "(0._" << kind << "/0.)";
611 } else if (IsInfinite()) {
612 if (IsNegative()) {
613 o << "(-1._" << kind << "/0.)";
614 } else {
615 o << "(1._" << kind << "/0.)";
616 }
617 } else {
618 using B = decimal::BinaryFloatingPointNumber<P>;
619 B value{word_.template ToUInt<typename B::RawType>()};
620 char buffer[common::MaxDecimalConversionDigits(P) +
621 EXTRA_DECIMAL_CONVERSION_SPACE];
622 decimal::DecimalConversionFlags flags{}; // default: exact representation
623 if (minimal) {
624 flags = decimal::Minimize;
625 }
626 auto result{decimal::ConvertToDecimal<P>(buffer, sizeof buffer, flags,
627 static_cast<int>(sizeof buffer), decimal::RoundNearest, value)};
628 const char *p{result.str};
629 if (DEREF(p) == '-' || *p == '+') {
630 o << *p++;
631 }
632 int expo{result.decimalExponent};
633 if (*p != '0') {
634 --expo;
635 }
636 o << *p << '.' << (p + 1);
637 if (expo != 0) {
638 o << 'e' << expo;
639 }
640 o << '_' << kind;
641 }
642 return o;
643 }
644
645 template class Real<Integer<16>, 11>;
646 template class Real<Integer<16>, 8>;
647 template class Real<Integer<32>, 24>;
648 template class Real<Integer<64>, 53>;
649 template class Real<Integer<80>, 64>;
650 template class Real<Integer<128>, 113>;
651 } // namespace Fortran::evaluate::value
652