- 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.17en",hypot(3., 4.));
- printf("hypot(3*10^150,4*10^150) =%25.17gn",hypot(3.e+150, 4.e+150));
- printf("hypot(3*10^306,4*10^306) =%25.17gn",hypot(3.e+306, 4.e+306));
- printf("hypot(3*10^-320,4*10^-320) =%25.17gn",hypot(3.e-320, 4.e-320));
- printf("hypot(0.7*DBL_MAX,0.7*DBL_MAX) =%25.17gn",hypot(0.7*DBL_MAX,0.7*DBL_MAX));
- printf("hypot(DBL_MAX, 1.0) =%25.17gn",hypot(DBL_MAX, 1.0));
- printf("hypot(1.0, DBL_MAX) =%25.17gn",hypot(1.0, DBL_MAX));
- printf("hypot(0.0, DBL_MAX) =%25.17gn",hypot(0.0, DBL_MAX));
- printf("n************************************************************n");
- printf("hypot(3, 4) =%25.17en",my_hypot01(3., 4.));
- printf("hypot(3*10^150,4*10^150) =%25.17gn",my_hypot01(3.e+150, 4.e+150));
- printf("hypot(3*10^306,4*10^306) =%25.17gn",my_hypot01(3.e+306, 4.e+306));
- printf("hypot(3*10^-320,4*10^-320) =%25.17gn",my_hypot01(3.e-320, 4.e-320));
- printf("hypot(0.7*DBL_MAX,0.7*DBL_MAX) =%25.17gn",my_hypot01(0.7*DBL_MAX,0.7*DBL_MAX));
- printf("hypot(DBL_MAX, 1.0) =%25.17gn",my_hypot01(DBL_MAX, 1.0));
- printf("hypot(1.0, DBL_MAX) =%25.17gn",my_hypot01(1.0, DBL_MAX));
- printf("hypot(0.0, DBL_MAX) =%25.17gn",my_hypot01(0.0, DBL_MAX));
- printf("n************************************************************n");
- printf("hypot(3, 4) =%25.17en",my_hypot02(3., 4.));
- printf("hypot(3*10^150,4*10^150) =%25.17gn",my_hypot02(3.e+150, 4.e+150));
- printf("hypot(3*10^306,4*10^306) =%25.17gn",my_hypot02(3.e+306, 4.e+306));
- printf("hypot(3*10^-320,4*10^-320) =%25.17gn",my_hypot02(3.e-320, 4.e-320));
- printf("hypot(0.7*DBL_MAX,0.7*DBL_MAX) =%25.17gn",my_hypot02(0.7*DBL_MAX,0.7*DBL_MAX));
- printf("hypot(DBL_MAX, 1.0) =%25.17gn",my_hypot02(DBL_MAX, 1.0));
- printf("hypot(1.0, DBL_MAX) =%25.17gn",my_hypot02(1.0, DBL_MAX));
- printf("hypot(0.0, DBL_MAX) =%25.17gn",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;
- }
- 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;
- }
- else
- { /* fraction part in low x */
- i = ((unsigned)(0xffffffff))>>(j0-20);
- iptr->sign = z->sign;
- iptr->exponent = z->exponent;
- iptr->mantissah = z->mantissah;
- iptr->mantissal = z->mantissal&(~i);
- if ( x == *y )
- {
- x = 0.0;
- z->sign = iptr->sign;
- return x;
- }
- return x - *y;
- }
- }
- //下面是两个要用到的函数
- int isnan(double d)
- {
- union
- {
- unsigned long long l;
- double d;
- } u;
- u.d=d;
- return (u.l==0x7FF8000000000000ll || u.l==0x7FF0000000000000ll || u.l==0xfff8000000000000ll);
- }
- int isinf(double d)
- {
- union
- {
- unsigned long long l;
- double d;
- } u;
- u.d=d;
- return (u.l==0x7FF0000000000000ll?1:u.l==0xFFF0000000000000ll?-1:0);
- }
- int main()
- {
- double x,y;
-
- x = 12345678901.1234567;
- printf("%f = (%f) + (%f) n",x,y,modf(x,&y));
- printf("%f = (%f) + (%f) n",x,y,my_modf01(x,&y));
- printf("%f = (%f) + (%f) n",x,y,my_modf02(x,&y));
-
- printf("n******************************************n");
-
- printf("%f = (%f) + (%f) n",-x,y,modf(-x,&y));
- printf("%f = (%f) + (%f) n",-x,y,my_modf01(-x,&y));
- printf("%f = (%f) + (%f) n",-x,y,my_modf02(-x,&y));
-
- return 0;
- }
- fmod ( )/* -- C语言库函数源代码 - */
- /*
- 计算x/y的余数。返回x-n*y,符号同y。
- n=[x/y](向离开零的方向取整)
- */
- double my_fmod01(double x, double y)
- {
- register double ret;
- __asm__(
- "1: fpremnt"
- "fstsw %%axnt"
- "sahfnt"
- "jp 1b"
- : "=t" (ret)
- : "0" (x), "u" (y)
- : "ax", "cc"
- );
- return ret;
- }
- double my_fmod02(double x, double y)
- {
- double temp, ret;
-
- if (y == 0.0)
- return 0.0;
- temp = floor(x/y);
- ret = x - temp * y;
- if ((x < 0.0) != (y < 0.0))
- ret = ret - y;
- return ret;
- }
- int main()
- {
- double x,y;
-
- x = 80.8,y = 3.0;
- printf("fmod(%f,%f) = %fn",x,y,fmod(x,y));
- printf("my_fmod01(%f,%f) = %fn",x,y,my_fmod01(x,y));
- printf("my_fmod02(%f,%f) = %fn",x,y,my_fmod01(x,y));
-
- printf("n******************************************n");
-
- x = -55.968,y = 8.8;
- printf("fmod(%f,%f) = %fn",x,y,fmod(x,y));
- printf("my_fmod01(%f,%f) = %fn",x,y,my_fmod01(x,y));
- printf("my_fmod02(%f,%f) = %fn",x,y,my_fmod01(x,y));
-
- return 0;
- }
- ldexp ( )/* -- C语言库函数源代码 - */
- /*
- 装载浮点数,v是尾数,e为指数。
- 如:x=ldexp(1.0,6);则表示要转载的浮点数是1.0*2^6
- */
- double my_ldexp01(double v, int e)
- {
- double two = 2.0;
- if (e < 0)
- {
- e = -e; /*这句话和后面的if语句是用来对两位溢出码的机器*/
- if (e < 0) return 0.0;
- while (e > 0)
- {
- if (e & 1) v /= two;
- two *= two;
- e >>= 1;
- }
- }
- else if (e > 0)
- {
- while (e > 0)
- {
- if (e & 1) v *= two;
- two *= two;
- e >>= 1;
- }
- }
- return v;
- }
- double my_ldexp02(double v, int e)
- {
- double temp1, texp, temp2;
- texp = e;
- __asm__(
- "fscale "
- : "=u" (temp2), "=t" (temp1)
- : "0" (texp), "1" (v)
- );
- return (temp1);
- }
- main()
- {
- float x,y;
-
- y = 1.0;
- x=my_ldexp01(y,6); // 1.0*2^6
- printf("2^6=%.2fn",x);
- x=my_ldexp01(-y,6); // 1.0*2^6
- printf("2^6=%.2fn",x);
- x=my_ldexp02(y,6); // 1.0*2^6
- printf("2^6=%.2fn",x);
- x=my_ldexp02(-y,6); // 1.0*2^6
- printf("2^6=%.2fn",x);
-
- return 0;
- }
- frexp ( )/* -- C语言库函数源代码 - */
- /*
- 把浮点数x分解成尾数和指数。x=m*2^exptr,m为规格化小数。
- 返回尾数m,并将指数存入exptr中。
- */
- double my_frexp01(double x, int *exptr)
- {
- union
- {
- double d;
- unsigned char c[8];
- } u;
- u.d = x;
- //得到移码,并减去1022得到指数值。
- *exptr = (int)(((u.c[7] & 0x7f) << 4) | (u.c[6] >> 4)) - 1022;
- //把指数部分置为0x03FE
- u.c[7] &= 0x80;
- u.c[7] |= 0x3f;
- u.c[6] &= 0x0f;
- u.c[6] |= 0xe0;
- return u.d;
- }
- double my_frexp02(double x, int *eptr)
- {
- union
- {
- double v;
- struct
- {
- unsigned mantissa2 : 32;
- unsigned mantissa1 : 20;
- unsigned exponent : 11;
- unsigned sign : 1;
- } s;
- } u;
- if (x)
- {
- u.v = x;
- //得到移码,并减去1022得到指数值。
- *eptr = u.s.exponent - 1022;
- //把指数部分置为0x03FE
- u.s.exponent = 1022;
- return(u.v);
- }
- else
- {
- *eptr = 0;
- return((double)0);
- }
- }
- main()
- {
- float x,y;
- int exp;
-
- y = 64.0;
- x = my_frexp01(y,&exp);
- printf("%f=%.2f*2^%dn",y,x,exp);
- x = my_frexp01(-y,&exp);
- printf("%f=%.2f*2^%dn",y,x,exp);
-
- printf("n************************n");
-
- x = my_frexp02(y,&exp);
- printf("%f=%.2f*2^%dn",y,x,exp);
- x = my_frexp02(-y,&exp);
- printf("%f=%.2f*2^%dn",y,x,exp);
- return 0;
- }
- tanh ( )/* -- C语言库函数源代码 - */
- double my_tanh(double x)
- {
- double ret,temp;
- if (x > 50)
- return 1;
- else if (x < -50)
- return -1;
- else
- {
- ret = exp(x);
- temp = 1.0 / ret;
- return ( (ret - temp) / (ret + temp));
- }
- }//计算x的双曲正切值。
- int main()
- {
- double a = 0.5;
- printf("tanh(%f) = %fn",a,tanh(a));
- printf("my_tanh(%f) = %fn",a,my_tanh(a));
- a = -0.5;
- printf("tanh(%f) = %fn",a,tanh(a));
- printf("my_tanh(%f) = %fn",a,my_tanh(a));
-
- return 0;
- }
- sinh ( )/* -- C语言库函数源代码 - */
- double my_sinh(double x)
- {
- double ret;
- if(x >= 0.0)
- {
- ret = exp(x);
- return (ret - 1.0/ret) / 2.0;
- }
- else
- {
- ret = exp(-x);
- return (1.0/ret - ret) / 2.0;
- }
- }//计算x的双曲正弦值。
- int main()
- {
- double a = 0.5;
- printf("sinh(%f) = %fn",a,sinh(a));
- printf("my_sinh(%f) = %fn",a,my_sinh(a));
- a = -0.5;
- printf("sinh(%f) = %fn",a,sinh(a));
- printf("my_sinh(%f) = %fn",a,my_sinh(a));
-
- return 0;
- }
- cosh ( )/* -- C语言库函数源代码 - */
- double my_cosh(double x)
- {
- double ret;
- ret = exp(fabs(x));
- return (ret + 1.0/ret) / 2.0;
- }//计算x的双曲余弦值。
- int main()
- {
- double a = 0.5;
- printf("cosh(%f) = %fn",a,cosh(a));
- printf("my_cosh(%f) = %fn",a,my_cosh(a));
- a = -0.5;
- printf("cosh(%f) = %fn",a,cosh(a));
- printf("my_cosh(%f) = %fn",a,my_cosh(a));
-
- return 0;
- }
- acos( )/* -- C语言库函数源代码 - */
- double atan2 (double x, double y)
- {
- register double ret;
- __asm__(
- "fpatannt"
- "fld %%st(0)"
- : "=t" (ret)
- : "0" (y), "u" (x)
- );
-
- return ret;
- }//求x / y的反正切值。
- double my_acos(double x)
- {
- return atan2 (sqrt (1.0 - x * x), x);
- }//求x的反余弦值。
- int main()
- {
- double a = 0.5;
- printf("acos(%f) = %fn",a,acos(a));
- printf("my_acos(%f) = %fn",a,my_acos(a));
- a = -0.5;
- printf("acos(%f) = %fn",a,acos(a));
- printf("my_acos(%f) = %fn",a,my_acos(a));
-
- return 0;
- }
复制代码
0
评分
-
参与人数 1 | 威望 +5 |
+5 |
积分 +5 |
+5 |
收起
理由
|
何生
| + 5 |
+ 5 |
+ 5 |
+ 5 |
谢谢分享^-^ |
查看全部评分
|
|
|
|