lib/hypot.c

   1 /* Hypotenuse of a right-angled triangle.
   2    Copyright (C) 2012-2024 Free Software Foundation, Inc.
   3
   4    This file is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU Lesser General Public License as
   6    published by the Free Software Foundation, either version 3 of the
   7    License, or (at your option) any later version.
   8
   9    This file is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU Lesser General Public License for more details.
  13
  14    You should have received a copy of the GNU Lesser General Public License
  15    along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
  16
  17 /* Written by Bruno Haible <bruno@clisp.org>, 2012.  */
  18
  19 #include <config.h>
  20
  21 /* Specification.  */
  22 #include <math.h>
  23
  24 double
  25 hypot (double x, double y)
  26 {
  27   if (isfinite (x) && isfinite (y))
  28     {
  29       /* Determine absolute values.  */
  30       x = fabs (x);
  31       y = fabs (y);
  32
  33       {
  34         /* Find the bigger and the smaller one.  */
  35         double a;
  36         double b;
  37
  38         if (x >= y)
  39           {
  40             a = x;
  41             b = y;
  42           }
  43         else
  44           {
  45             a = y;
  46             b = x;
  47           }
  48         /* Now 0 <= b <= a.  */
  49
  50         {
  51           int e;
  52           double an;
  53           double bn;
  54
  55           /* Write a = an * 2^e, b = bn * 2^e with 0 <= bn <= an < 1.  */
  56           an = frexp (a, &e);
  57           bn = ldexp (b, - e);
  58
  59           {
  60             double cn;
  61
  62             /* Through the normalization, no unneeded overflow or underflow
  63                will occur here.  */
  64             cn = sqrt (an * an + bn * bn);
  65             return ldexp (cn, e);
  66           }
  67         }
  68       }
  69     }
  70   else
  71     {
  72       if (isinf (x) || isinf (y))
  73         /* x or y is infinite.  Return +Infinity.  */
  74         return HUGE_VAL;
  75       else
  76         /* x or y is NaN.  Return NaN.  */
  77         return x + y;
  78     }
  79 }