|
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; } 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: fprem\n\t" "fstsw %%ax\n\t" "sahf\n\t" "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) = %f\n",x,y,fmod(x,y)); printf("my_fmod01(%f,%f) = %f\n",x,y,my_fmod01(x,y)); printf("my_fmod02(%f,%f) = %f\n",x,y,my_fmod01(x,y)); printf("\n******************************************\n"); x = -55.968,y = 8.8; printf("fmod(%f,%f) = %f\n",x,y,fmod(x,y)); printf("my_fmod01(%f,%f) = %f\n",x,y,my_fmod01(x,y)); printf("my_fmod02(%f,%f) = %f\n",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=%.2f\n",x); x=my_ldexp01(-y,6); // 1.0*2^6 printf("2^6=%.2f\n",x); x=my_ldexp02(y,6); // 1.0*2^6 printf("2^6=%.2f\n",x); x=my_ldexp02(-y,6); // 1.0*2^6 printf("2^6=%.2f\n",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^%d\n",y,x,exp); x = my_frexp01(-y,&exp); printf("%f=%.2f*2^%d\n",y,x,exp); printf("\n************************\n"); x = my_frexp02(y,&exp); printf("%f=%.2f*2^%d\n",y,x,exp); x = my_frexp02(-y,&exp); printf("%f=%.2f*2^%d\n",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) = %f\n",a,tanh(a)); printf("my_tanh(%f) = %f\n",a,my_tanh(a)); a = -0.5; printf("tanh(%f) = %f\n",a,tanh(a)); printf("my_tanh(%f) = %f\n",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) = %f\n",a,sinh(a)); printf("my_sinh(%f) = %f\n",a,my_sinh(a)); a = -0.5; printf("sinh(%f) = %f\n",a,sinh(a)); printf("my_sinh(%f) = %f\n",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) = %f\n",a,cosh(a)); printf("my_cosh(%f) = %f\n",a,my_cosh(a)); a = -0.5; printf("cosh(%f) = %f\n",a,cosh(a)); printf("my_cosh(%f) = %f\n",a,my_cosh(a)); return 0; } acos( )/* -- C语言库函数源代码 - */ double atan2 (double x, double y) { register double ret; __asm__( "fpatan\n\t" "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) = %f\n",a,acos(a)); printf("my_acos(%f) = %f\n",a,my_acos(a)); a = -0.5; printf("acos(%f) = %f\n",a,acos(a)); printf("my_acos(%f) = %f\n",a,my_acos(a)); return 0; } |
|
|
