124 lines
3.1 KiB
C
124 lines
3.1 KiB
C
|
/***
|
||
|
|
*hypot.c - hypotenuse and complex absolute value
|
||
|
|
*
|
||
|
|
* Copyright (c) 1991-2001, Microsoft Corporation. All rights reserved.
|
||
|
|
*
|
||
|
|
*Purpose:
|
||
|
|
*
|
||
|
|
*Revision History:
|
||
|
|
* 08-15-91 GDP written
|
||
|
|
* 10-20-91 GDP removed inline assembly for calling sqrt
|
||
|
|
*
|
||
|
|
*******************************************************************************/
|
||
|
|
#include <math.h>
|
||
|
|
#include <trans.h>
|
||
|
|
|
||
|
|
static double _hypothlp(double x, double y, int who);
|
||
|
|
|
||
|
|
/*
|
||
|
|
* Function name: hypot
|
||
|
|
*
|
||
|
|
* Arguments: x, y - double
|
||
|
|
*
|
||
|
|
* Description: hypot returns sqrt(x*x + y*y), taking precautions against
|
||
|
|
* unwarrented over and underflows.
|
||
|
|
*
|
||
|
|
* Side Effects: no global data is used or affected.
|
||
|
|
*
|
||
|
|
* Author: written R.K. Wyss, Microsoft, Sept. 9, 1983
|
||
|
|
*
|
||
|
|
* History:
|
||
|
|
* 03/13/89 WAJ Minor changes to source.
|
||
|
|
* 04/13/89 WAJ Now uses _cdecl for _CDECL
|
||
|
|
* 06/07/91 JCR ANSI naming (_hypot)
|
||
|
|
* 08/26/91 GDP NaN support, error handling
|
||
|
|
* 01/13/91 GDP IEEE exceptions support
|
||
|
|
*/
|
||
|
|
|
||
|
|
double _hypot(double x, double y)
|
||
|
|
{
|
||
|
|
return _hypothlp(x,y,OP_HYPOT);
|
||
|
|
}
|
||
|
|
|
||
|
|
/***
|
||
|
|
*double _cabs(struct _complex z) - absolute value of a complex number
|
||
|
|
*
|
||
|
|
*Purpose:
|
||
|
|
*
|
||
|
|
*Entry:
|
||
|
|
*
|
||
|
|
*Exit:
|
||
|
|
*
|
||
|
|
*Exceptions:
|
||
|
|
*******************************************************************************/
|
||
|
|
double _cabs(struct _complex z)
|
||
|
|
{
|
||
|
|
return( _hypothlp(z.x, z.y, OP_CABS ) );
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
static double _hypothlp(double x, double y, int who)
|
||
|
|
{
|
||
|
|
double max;
|
||
|
|
double result, sum;
|
||
|
|
uintptr_t savedcw;
|
||
|
|
int exp1, exp2, newexp;
|
||
|
|
|
||
|
|
/* save user fp control word */
|
||
|
|
savedcw = _maskfp();
|
||
|
|
|
||
|
|
/* check for infinity or NAN */
|
||
|
|
if (IS_D_SPECIAL(x) || IS_D_SPECIAL(y)){
|
||
|
|
if (IS_D_SNAN(x) || IS_D_SNAN(y)){
|
||
|
|
return _except2(FP_I,who,x,y,_d_snan2(x,y),savedcw);
|
||
|
|
}
|
||
|
|
if (IS_D_QNAN(x) || IS_D_QNAN(y)){
|
||
|
|
return _handle_qnan2(who,x,y,savedcw);
|
||
|
|
}
|
||
|
|
/* there is at least one infinite argument ... */
|
||
|
|
RETURN(savedcw,D_INF);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
/* take the absolute values of x and y, compute the max, and then scale by
|
||
|
|
max to prevent over or underflowing */
|
||
|
|
|
||
|
|
if ( x < 0.0 )
|
||
|
|
x = - x;
|
||
|
|
|
||
|
|
if ( y < 0.0 )
|
||
|
|
y = - y;
|
||
|
|
|
||
|
|
max = ( ( y > x ) ? y : x );
|
||
|
|
|
||
|
|
if ( max == 0.0 )
|
||
|
|
RETURN(savedcw, 0.0 );
|
||
|
|
|
||
|
|
x /= max; //this may pollute the fp status word (underflow flag)
|
||
|
|
y /= max;
|
||
|
|
|
||
|
|
sum = x*x + y*y;
|
||
|
|
|
||
|
|
result = _decomp(sqrt(sum),&exp1) * _decomp(max,&exp2);
|
||
|
|
newexp = exp1 + exp2 + _get_exp(result);
|
||
|
|
|
||
|
|
// in case of overflow or underflow
|
||
|
|
// adjusting exp by IEEE_ADJUST will certainly
|
||
|
|
// bring the result in the representable range
|
||
|
|
|
||
|
|
if (newexp > MAXEXP) {
|
||
|
|
result = _set_exp(result, newexp - IEEE_ADJUST);
|
||
|
|
return _except2(FP_O | FP_P, who, x, y, result, savedcw);
|
||
|
|
}
|
||
|
|
if (newexp < MINEXP) {
|
||
|
|
result = _set_exp(result, newexp + IEEE_ADJUST);
|
||
|
|
return _except2(FP_U | FP_P, who, x, y, result, savedcw);
|
||
|
|
}
|
||
|
|
|
||
|
|
result = _set_exp(result, newexp);
|
||
|
|
// fix needed: P exception is raised even if the result is exact
|
||
|
|
|
||
|
|
RETURN_INEXACT2(who, x, y, result, savedcw);
|
||
|
|
}
|