libc/src/math/generic/hypotf.cpp

   1 //===-- Implementation of hypotf function ---------------------------------===//
   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 #include "src/math/hypotf.h"
   9 #include "src/__support/FPUtil/BasicOperations.h"
  10 #include "src/__support/FPUtil/FPBits.h"
  11 #include "src/__support/FPUtil/sqrt.h"
  12 #include "src/__support/common.h"
  13
  14 namespace __llvm_libc {
  15
  16 LLVM_LIBC_FUNCTION(float, hypotf, (float x, float y)) {
  17   using DoubleBits = fputil::FPBits<double>;
  18   using FPBits = fputil::FPBits<float>;
  19
  20   FPBits x_bits(x), y_bits(y);
  21
  22   uint16_t x_exp = x_bits.get_unbiased_exponent();
  23   uint16_t y_exp = y_bits.get_unbiased_exponent();
  24   uint16_t exp_diff = (x_exp > y_exp) ? (x_exp - y_exp) : (y_exp - x_exp);
  25
  26   if (exp_diff >= fputil::MantissaWidth<float>::VALUE + 2) {
  27     return fputil::abs(x) + fputil::abs(y);
  28   }
  29
  30   double xd = static_cast<double>(x);
  31   double yd = static_cast<double>(y);
  32
  33   // These squares are exact.
  34   double x_sq = xd * xd;
  35   double y_sq = yd * yd;
  36
  37   // Compute the sum of squares.
  38   double sum_sq = x_sq + y_sq;
  39
  40   // Compute the rounding error with Fast2Sum algorithm:
  41   // x_sq + y_sq = sum_sq - err
  42   double err = (x_sq >= y_sq) ? (sum_sq - x_sq) - y_sq : (sum_sq - y_sq) - x_sq;
  43
  44   // Take sqrt in double precision.
  45   DoubleBits result(fputil::sqrt(sum_sq));
  46
  47   if (!DoubleBits(sum_sq).is_inf_or_nan()) {
  48     // Correct rounding.
  49     double r_sq = static_cast<double>(result) * static_cast<double>(result);
  50     double diff = sum_sq - r_sq;
  51     constexpr uint64_t mask = 0x0000'0000'3FFF'FFFFULL;
  52     uint64_t lrs = result.uintval() & mask;
  53
  54     if (lrs == 0x0000'0000'1000'0000ULL && err < diff) {
  55       result.bits |= 1ULL;
  56     } else if (lrs == 0x0000'0000'3000'0000ULL && err > diff) {
  57       result.bits -= 1ULL;
  58     }
  59   } else {
  60     FPBits bits_x(x), bits_y(y);
  61     if (bits_x.is_inf_or_nan() || bits_y.is_inf_or_nan()) {
  62       if (bits_x.is_inf() || bits_y.is_inf())
  63         return static_cast<float>(FPBits::inf());
  64       if (bits_x.is_nan())
  65         return x;
  66       return y;
  67     }
  68   }
  69
  70   return static_cast<float>(static_cast<double>(result));
  71 }
  72
  73 } // namespace __llvm_libc