c语⾔数学库math源码,math.h函数源代码hypot ( )/* -- C语⾔库函数源代码 - */
/*
hypot函数对于给定的直⾓三⾓形的两个直⾓边，
求其斜边的长度。
*/
//⼀般的常规算法：
double my_hypot01(double x, double y)
{
double hypotenuse;
x = fabs(x);
y = fabs(y);
IF (x < y)
{
double temp = x;
x = y;
y = temp;
}
if (x == 0.)
return 0.;
else
{
hypotenuse = y/x;
return x*sqrt(1.+hypotenuse*hypotenuse);
}
}
#define __SQRT_DBL_MAX 1.3e+154
#define __SQRT_DBL_MIN 2.3e-162
double my_hypot02(double x, double y)
{
double ratio;
double r, t, s, p, q;
x = fabs(x), y = fabs(y);
if (x < y)
{
double temp = x;
x = y;
y = temp;
}//保持x 是两个直⾓边中较长的边，y是较短的边。
if (y == 0.)
return x;
/*
主要考虑的是当x很⼤⽽y很⼩，那么它们的平⽅和将会造成
丢失⼩数的现象。⾸先要判断y是否是太⼩，x是不是太⼤。
如果出现这种情况则⽤，第⼀个公式来处理。其他的则⽤
这样可以让求出的斜边更加精确⼀些。
*/
if ((ratio = y / x) > __SQRT_DBL_MIN && x < __SQRT_DBL_MAX) return x * sqrt(1. + ratio*ratio);
else
{//使⽤3次迭代是增加精确度。
r = ratio*ratio, p = x, q = y;
do
{
t = 4.+ r;
if (t == 4.)
break;
s = r / t;
p += 2. * s * p;
q *= s;
r = (q / p) * (q / p);
} while (1);
return p;
}
}
struct complex
{
double x;
double y;
}
double cabs(struct complex x)
{
return hypot(z.x,z.y);
}//hypot 函数的封装(这⾥不再举调⽤的例⼦了。)
# define DBL_MAX 1.79769313486231e+308
# define DBL_MIN 2.22507385850721e-308
int main(void)
{
printf("hypot(3, 4)                    =%25.17e\n",hypot(3., 4.));
printf("hypot(3*10^150,4*10^150)        =%25.17g\n",hypot(3.e+150, 4.e+150));
printf("hypot(3*10^306,4*10^306)        =%25.17g\n",hypot(3.e+306, 4.e+306));
printf("hypot(3*10^-320,4*10^-320)      =%25.17g\n",hypot(3.e-320, 4.e-320));
printf("hypot(0.7*DBL_MAX,0.7*DBL_MAX) =%25.17g\n",hypot(0.7*DBL_MAX,0.7*DBL_MAX));
printf("hypot(DBL_MAX, 1.0)            =%25.17g\n",hypot(DBL_MAX, 1.0));
printf("hypot(1.0, DBL_MAX)            =%25.17g\n",hypot(1.0, DBL_MAX));
printf("hypot(0.0, DBL_MAX)            =%25.17g\n",hypot(0.0, DBL_MAX));
printf("\n************************************************************\n");
printf("hypot(3, 4)                    =%25.17e\n",my_hypot01(3., 4.));
printf("hypot(3*10^150,4*10^150)        =%25.17g\n",my_hypot01(3.e+150, 4.e+150));
printf("hypot(3*10^306,4*10^306)        =%25.17g\n",my_hypot01(3.e+306, 4.e+306));
printf("hypot(3*10^-320,4*10^-320)      =%25.17g\n",my_hypot01(3.e-320, 4.e-320));
printf("hypot(0.7*DBL_MAX,0.7*DBL_MAX) =%25.17g\n",my_hypot01(0.7*DBL_MAX,0.7*DBL_MAX)); printf("hypot(DBL_MAX, 1.0)            =%25.17g\n",my_hypot01(DBL_MAX, 1.0));
printf("hypot(1.0, DBL_MAX)            =%25.17g\n",my_hypot01(1.0, DBL_MAX));
printf("hypot(0.0, DBL_MAX)            =%25.17g\n",my_hypot01(0.0, DBL_MAX));
printf("\n************************************************************\n");
printf("hypot(3, 4)                    =%25.17e\n",my_hypot02(3., 4.));
printf("hypot(3*10^150,4*10^150)        =%25.17g\n",my_hypot02(3.e+150, 4.e+150));
printf("hypot(3*10^306,4*10^306)        =%25.17g\n",my_hypot02(3.e+306, 4.e+306));
printf("hypot(3*10^-320,4*10^-320)      =%25.17g\n",my_hypot02(3.e-320, 4.e-320));
printf("hypot(0.7*DBL_MAX,0.7*DBL_MAX) =%25.17g\n",my_hypot02(0.7*DBL_MAX,0.7*DBL_MAX)
); printf("hypot(DBL_MAX, 1.0)            =%25.17g\n",my_hypot02(DBL_MAX, 1.0));
printf("hypot(1.0, DBL_MAX)            =%25.17g\n",my_hypot02(1.0, DBL_MAX)); printf("hypot(0.0, DBL_MAX)            =%25.17g\n",my_hypot02(0.0, DBL_MAX)); return 0;
}
modf ( )/* -- C语⾔库函数源代码 - */
/*
将浮点数x分解成整数部分和⼩数部分。
返回⼩数部分，将整数部分存⼊* iptr所指内存中。
*/
double my_modf01(double x, double *iptr)
{
double ret = fmod(x,1.0);
*iptr = x - ret;
return ret;
}//这个函数算法⽐较简单，也容易让⼈理解。
//下⾯的这个函数理解起来就有点困难了。
typedef struct
{
unsigned int mantissal:32;
unsigned int mantissah:20;
unsigned int exponent:11;
unsigned int sign:1;
}double_t;//这个结构体在IEEE.h定义。
double my_modf02(double x, double *y)
{
double_t * z = (double_t *)&x;
double_t * iptr = (double_t *)y;
int j0;
unsigned int i;
j0 = z->exponent - 0x3ff;  /* exponent of x */
if(j0<20)
{/* integer part in high x */
if(j0<0)
{                  /* |x|<1 */
*y = 0.0;
iptr->sign = z->sign;
return x;
}
else
{
if ( z->mantissah == 0 && z->mantissal == 0 ) {
*y = x;
return 0.0;
}
i = (0x000fffff)>>j0;
iptr->sign = z->sign;
iptr->exponent = z->exponent;
iptr->mantissah = z->mantissah&(~i);
iptr->mantissal = 0;
if ( x == *y )
{
x = 0.0;
z->sign = iptr->sign;
return x;
源代码1080p在线}
return x - *y;
}
}
else if (j0>51)
{              /* no fraction part */
*y = x;
if ( isnan(x) || isinf(x) )
return x;
x = 0.0;
z->sign = iptr->sign;
return x;
}

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。