淡泊明志 宁静致远
分类： C/C++
2006-12-08 15:13:49
【C语言库函数源代码】
【本程序在Dev C++ 4.9.9.2 下编译通过】
/*
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));
system("pause");
return 0;
}
