mingwrt/mingwex/math/hypot_generic.c

   1 /*
   2  * hypot_generic.c
   3  *
   4  * Compute the length of the hypotenuse of a right triangle, given the
   5  * lengths of the two perpendicular sides.
   6  *
   7  * $Id$
   8  *
   9  * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
  10  * Copyright 2015, MinGW.org Project
  11  *
  12  *
  13  * Permission is hereby granted, free of charge, to any person obtaining a
  14  * copy of this software and associated documentation files (the "Software"),
  15  * to deal in the Software without restriction, including without limitation
  16  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  17  * and/or sell copies of the Software, and to permit persons to whom the
  18  * Software is furnished to do so, subject to the following conditions:
  19  *
  20  * The above copyright notice and this permission notice (including the next
  21  * paragraph) shall be included in all copies or substantial portions of the
  22  * Software.
  23  *
  24  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  25  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  26  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  27  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  28  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  29  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  30  * DEALINGS IN THE SOFTWARE.
  31  *
  32  *
  33  * This module should be compiled separately for each supported function:
  34  *
  35  *   gcc -D FUNCTION=hypot  -o hypot.o  hypot_generic.c
  36  *   gcc -D FUNCTION=hypotl -o hypotl.o hypot_generic.c
  37  *   gcc -D FUNCTION=hypotf -o hypotf.o hypot_generic.c
  38  *
  39  */
  40 #include <math.h>
  41 #include <float.h>
  42 #include <errno.h>
  43
  44 #ifndef FUNCTION
  45 /* If user neglected to specify it, the default compilation is for
  46  * the hypot() function.
  47  */
  48 # define FUNCTION hypot
  49 #endif
  50
  51 #define argtype_hypot  double
  52 #define argtype_hypotl long double
  53 #define argtype_hypotf float
  54
  55 #define PASTE(PREFIX,SUFFIX)   PREFIX##SUFFIX
  56 #define mapname(PREFIX,SUFFIX) PASTE(PREFIX,SUFFIX)
  57
  58 #define ARGTYPE(FUNCTION) PASTE(argtype_,FUNCTION)
  59 #define ARGCAST(VALUE)    mapname(argcast_,FUNCTION)(VALUE)
  60
  61 #define argcast_hypot(VALUE)  VALUE
  62 #define argcast_hypotl(VALUE) PASTE(VALUE,L)
  63 #define argcast_hypotf(VALUE) PASTE(VALUE,F)
  64
  65 #define mapfunc(NAME) mapname(mapfunc_,FUNCTION)(NAME)
  66
  67 #define mapfunc_hypot(NAME)  NAME
  68 #define mapfunc_hypotl(NAME) PASTE(NAME,l)
  69 #define mapfunc_hypotf(NAME) PASTE(NAME,f)
  70
  71 #define hypot_argval_max  DBL_MAX
  72 #define hypotf_argval_max FLT_MAX
  73 #define hypotl_argval_max LDBL_MAX
  74
  75 #define hypot_result_overflow  HUGE_VAL
  76 #define hypotf_result_overflow HUGE_VALF
  77 #define hypotl_result_overflow HUGE_VALL
  78
  79 /* Prefer fast version of the sqrt() function.
  80  */
  81 #include "fastmath.h"
  82 #define sqrt __fast_sqrt
  83
  84 /* Define the generic function implementation.
  85  */
  86 #define METHODL 1
  87 ARGTYPE(FUNCTION) FUNCTION( ARGTYPE(FUNCTION) x, ARGTYPE(FUNCTION) y )
  88 {
  89 # if ARGCAST(METHOD)
  90   /*
  91    * By default, we implement the actual computation for the
  92    * hypotl() function only, (since only METHODL is explicitly
  93    * defined as non-zero); to do likewise for the hypot() and
  94    * hypotf() functions, compile with -D METHOD and -D METHODF
  95    * set respectively.
  96    *
  97    */
  98   if( isfinite(x) && isfinite(y) )
  99   {
 100     /* For all finite values of x and y, we may compute hypot(x,y)
 101      * in terms of their respective magnitudes; we apply Pythagoras'
 102      * theorum, but to avoid a possible overflow when computing the
 103      * sum of the squares of x and y, we use a transformed variant
 104      * of the conventional formula, based on the assumption that
 105      * the arguments are ordered such that |x| >= |y|...
 106      */
 107     ARGTYPE(FUNCTION) h;
 108     if( (h = x = mapfunc(fabs)( x )) < (y = mapfunc(fabs)( y )) )
 109     {
 110       /* ...otherwise, we transpose them.
 111        */
 112       x = y; y = h;
 113     }
 114     /* After ensuring that |x| >= |y|, we must perform one more
 115      * sanity check to ensure that we don't have |x| == 0.0, (in
 116      * which case |y| must also be 0.0, and the quotient |y|/|x|
 117      * becomes indeterminate).
 118      */
 119     if( x == ARGCAST(0.0) )
 120       /*
 121        * For this specific case, we may avoid the need to compute
 122        * the indeterminate quotient, by simple evaluation of the
 123        * normal expression of Pythagoras' theorum:
 124        *
 125        *   hypot = sqrt( x * x + y * y )
 126        *
 127        * directly, to obtain the expected result of zero.
 128        */
 129       return ARGCAST(0.0);
 130
 131     /* We can now be confident that the ratio...
 132      */
 133     h = y / x;
 134     /*
 135      * ...is determinate, and must lie in the range 0.0 .. 1.0, so
 136      * we may safely compute the overflow protected transformation
 137      * of the expression of Pythagoras' theorum:
 138      *
 139      *   hypot = x * sqrt( 1.0 + (y / x) * (y / x))
 140      *
 141      * without any possibility of overflow, when squaring the y / x
 142      * ratio term; however, there is still a slight possibility that
 143      * x itself is so large that multiplication by the square root
 144      * term, (which must lie in the range 1.0 .. sqrt(2.0)), may
 145      * cause overflow in the final result...
 146      */
 147     h = mapfunc(sqrt)( ARGCAST(1.0) + h * h );
 148     if( x > (mapname(FUNCTION,_argval_max) / h) )
 149     {
 150       /* ...in which case, we set errno, and return the appropriate
 151        * representation of HUGE_VAL...
 152        */
 153       errno = ERANGE;
 154       return mapname(FUNCTION,_result_overflow);
 155     }
 156     /* ...otherwise, we may compute and return the requisite result;
 157      * we may be certain that this is unaffected by OVERFLOW.  On the
 158      * other hand, UNDERFLOW may have occurred in the intermediate
 159      * computation of y/x; however, the effect is insignificant.
 160      */
 161     return x * h;
 162   }
 163   /* If we get to here, one or both of x and y was indefinite; POSIX
 164    * requires us to return +inf, if either x or y is infinite, (even
 165    * if the other is nan)...
 166    */
 167   if( isinf( x ) ) return mapfunc(copysign)( x, ARGCAST(+1.0) );
 168   if( isinf( y ) ) return mapfunc(copysign)( y, ARGCAST(+1.0) );
 169   /*
 170    * ...otherwise, combine them to force an indefinite result.
 171    */
 172   return x + y;
 173
 174 #else
 175   /* For the hypot() and hypotf() implementations, we prefer to
 176    * compute an intermediate result in the long double domain...
 177    */
 178   long double intermediate = hypotl( (long double)(x), (long double)(y) );
 179
 180   /* ...then check that this may be represented within the
 181    * effective domain of the return data type...
 182    */
 183   if(  isfinite( intermediate )
 184   &&  (intermediate > (long double)(mapname(FUNCTION,_argval_max))) )
 185   {
 186     /* ...otherwise flag a RANGE error, and return the
 187      * appropriate representation of HUGE_VAL...
 188      */
 189     errno = ERANGE;
 190     return mapname(FUNCTION,_result_overflow);
 191   }
 192   /* ...or the valid imtermediate result, cast to the
 193    * appropriate data type for the function.
 194    */
 195   return (ARGTYPE(FUNCTION))(intermediate);
 196
 197 #endif
 198 }
 199
 200 /* $RCSfile$: end of file */