htswanalysis/MACS/lib/gsl/gsl-1.11/sys/hypot.c

   1 /* sys/hypot.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough, Patrick Alken
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <math.h>
  22 #include <gsl/gsl_math.h>
  23
  24 double gsl_hypot (const double x, const double y);
  25 double gsl_hypot3 (const double x, const double y, const double z);
  26
  27 double gsl_hypot (const double x, const double y)
  28 {
  29   double xabs = fabs(x) ;
  30   double yabs = fabs(y) ;
  31   double min, max;
  32
  33   if (xabs < yabs) {
  34     min = xabs ;
  35     max = yabs ;
  36   } else {
  37     min = yabs ;
  38     max = xabs ;
  39   }
  40
  41   if (min == 0)
  42     {
  43       return max ;
  44     }
  45
  46   {
  47     double u = min / max ;
  48     return max * sqrt (1 + u * u) ;
  49   }
  50 }
  51
  52 double
  53 gsl_hypot3(const double x, const double y, const double z)
  54 {
  55   double xabs = fabs(x);
  56   double yabs = fabs(y);
  57   double zabs = fabs(z);
  58   double w = GSL_MAX(xabs, GSL_MAX(yabs, zabs));
  59
  60   if (w == 0.0)
  61     {
  62       return (0.0);
  63     }
  64   else
  65     {
  66       double r = w * sqrt((xabs / w) * (xabs / w) +
  67                           (yabs / w) * (yabs / w) +
  68                           (zabs / w) * (zabs / w));
  69       return r;
  70     }
  71 }