标签：运算 四舍五入 回忆一下 bsp 测试结果 理论 a.out 进制 简单
了解了浮点数的存储以及手算平方根的原理，我们可以考虑程序实现了。
先实现一个64位整数的平方根，根据之前的手算平方根，程序也不是那么难写了。
#include <stdint.h>
uint64_t _sqrt_u64(uint64_t a) { int i; uint64_t res; uint64_t remain; //0的平方根是0,特殊处理一下 if(a == 0ull) return 0ull; //找到最高位的1，并且产生平方根结果最高位的1 for(i=62;;i-=2) if(a&(3ull<<i)) { res = 1ull; remain = ((a&(3ull<<i))>>i) - 1ull; i -= 2; break; } //根据手算平方根的原理，依次产生各位结果 for(;i>=0;i-=2) { //右移动两位，并把a接着的两位并入remain remain = (remain<<2)|((a&(3ull<<i))>>i); if(((res<<2)|1ull) <= remain) { //产生新一位的1 remain = remain - ((res<<2)|1ull); res = (res<<1)|1ull; } else { //产生新一位的0 res <<= 1; } } return res; }
其实，可以合在一起写，代码会短一些，但效率会低那么一点点，而且编译器应该不太容易优化。
#include <stdint.h>
uint64_t _sqrt_u64(uint64_t a) { int i; uint64_t res; uint64_t remain; res = remain = 0ull; for(i=62;i>=0;i-=2) { remain = (remain<<2)|((a&(3ull<<i))>>i); if(((res<<2)|1ull) <= remain) { remain = remain - ((res<<2)|1ull); res = (res<<1)|1ull; } else { res <<= 1; } } return res; }
不过，我们不需要这个结果。
为了验证其正确性，我们来写个C语言的main
#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>
uint64_t _sqrt_u64(uint64_t a);
int main()
{
uint64_t a, b;
scanf("%" PRIu64, &a);
b = _sqrt_u64(a);
printf("%" PRIu64 "\n",b);
return 0;
}
我们shell程序测试一下，我们当然不可能测试过每一个64bits的数，这个运算量太大，不现实。我们可以用随机取一部分来测试。
#!/bin/bash #编译 gcc -O2 -s sqrt_u64.c main_sqrt_u64.c -o a.out #随机测试10000个数 for((i=0;i<10000;i++));do #随机产生bits 0~64，如果是0，代表测试的数就是0 #如果不是0，则代表要产生的数二进制可以有多少位 let bits=RANDOM%65 if [ $bits -eq 0 ]; then x=0 y=0 else #产生一个bits位的二进制数x x=$({ #最高位1 echo -n 1 #之后每位随机产生 for((j=1;j<bits;j++));do let x=RANDOM%2 echo -n $x done }) #用bc将x转换成十进制 x=$(echo ‘obase=10;ibase=2;‘"$x" | bc) #用bc计算x的平方根取整，理论上和我们的C语言计算一致 y=$(echo ‘sqrt(‘"$x"‘)‘ | bc) fi #z是我们的C语言计算结果 z=$(echo $x | ./a.out) #比较，如果不一致，就报错 if [ $y -ne $z ];then echo $x $y $z error exit 1 fi done echo OK
测试结果表明，我们的C语言还是可以得到正确的结果的。
再来回忆下第一节里讲过的浮点数结构，
S(1bits) | N(8bits) | A(23bits)
对于浮点数a*2n,
1<=a<2，n为整数，
如果n是偶数，
那么a*2n的平方根是sqrt(a)*2n/2，也满足1<=sqrt(a)<2，n/2是整数；
如果n为奇数，
那么a*2n的平方根是sqrt(2*a)*2(n-1)/2，也满足1<=sqrt(2*a)<2，(n-1)/2是整数。
所以此处要用a或者2*a来开平方根，
回忆一下浮点数的结构，单精度浮点数的精度是23位。
表示的是科学计数法a*2n的a减去1的部分，那么加上整数1可以用二进制24位表示。
于是，我们就想，一个二进制48位或47位长的数，平方根是二进制24位。那么，我们就可以用一个48位或47位的二进制整数的平方根计算结果的小数部分。
nan/inf/-inf以及负数的平方根都是nan，
0.0的平方根是0.0,
-0.0的平方根是-0.0（可能只是某些库里是这样的），
以上都可以在计算的时候特殊化一下。
规格数（就是用科学计数法表示的浮点数）的平方根也是规格数，
S=0，N=0，A>0代表的是A*2-149，也就是(A*2)*2-150，
我们稍微计算一下，可以明白，所有的此类数的平方根都在规格数表示的范围内。
于是，有了以下的代码。
#include <stdint.h>
static uint32_t _sqrt_(uint64_t a) { int i; uint64_t res; uint64_t remain; res = remain = 0ull; //之前整数平方根被直接优化，我们只需要求47位或者48位整数的平方根 for(i=46;i>=0;i-=2) { remain = (remain<<2)|((a&(3ull<<i))>>i); if(((res<<2)|1ull) <= remain) { remain = remain - ((res<<2)|1ull); res = (res<<1)|1ull; } else { res <<= 1; } } return (uint32_t)res; } float mysqrtf(float f) { union { float f; uint32_t u; } n; uint32_t N,A; int _N, i; uint64_t _A; n.f = f; if(n.u == 0x80000000 || n.u == 0x00000000) /* 0.0/-0.0 */ return n.f; N = (n.u&(0xff<<23))>>23; if(N==0xff||(n.u&0x80000000)) { /* inf/-inf/nan/ f < 0.0*/ n.u = 0x7fc00000; /* nan */ return n.f; } if(N!=0x0) { /* 用科学计数法表示的规格数 */ A = (n.u&0x7fffff)|0x800000; _N = (int)N - 127; if(N&0x1) { _A = (uint64_t)A<<23; } else { _A = (uint64_t)A<<24; _N--; } } else { //A*2^(-149)这种表示方式的浮点数 //还是需要找最高位 for(i=22;;i--) if(n.u&((0x1)<<i)) break; //然后需要移位，要区分奇数和偶数 if(i&0x1) { _N = i-149; _A = (uint64_t)n.u << (46-i); } else { _N = i-150; _A = (uint64_t)n.u << (47-i); } } //小数部分 A = _sqrt_(_A); //指数部分 N = (uint32_t)(_N/2+127); //得到结果 n.u = (A&0x7fffff)|(N<<23); return n.f; }
同样，也写个测试用的程序，对inf/-inf/nan/0.0/-0.0以及负数不测了，这些很简单。
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <inttypes.h>
int main(int argc, char **argv)
{
union {
float f;
uint32_t u;
} n;
uint32_t A,N;
float f,f2;
int i;
srand((unsigned)time(NULL));
//随机10000个数据
for(i=0;i<10000;i++) {
N = rand()%256;
if(N==255)
N=254;
A = 0x0;
A |= rand()%256;
A |= (rand()%256)>>8;
A |= (rand()%256)>>16;
n.u = (A&0x7fffff)|(N<<23);
f = sqrtf(n.f);
f2 = mysqrtf(n.f);
printf("%.60f %.60f\n",f,f2);
}
return 0;
}
结果发现，我们的程序和数学库里的sqrtf结果有细微差别。
于是，我们决定再加个小东西，就是四舍五入。之前我们用的是47位或者48位数开平方，为了四舍五入，我们需要多一位，于是就用49位或者50位数开平方。
修改一下mysqrtf，增加两位拿去开平方，_sqrt_也动一下。
#include <stdint.h>
static uint32_t _sqrt_(uint64_t a) { int i; uint64_t res; uint64_t remain; res = remain = 0ull; //之前整数平方根被直接优化，我们只需要求49位或者50位整数的平方根 for(i=48;i>=0;i-=2) {//这里之前是46,改成48 remain = (remain<<2)|((a&(3ull<<i))>>i); if(((res<<2)|1ull) <= remain) { remain = remain - ((res<<2)|1ull); res = (res<<1)|1ull; } else { res <<= 1; } } return (uint32_t)res; }
float mysqrtf(float f) { union { float f; uint32_t u; } n; uint32_t N,A; int _N, i; uint64_t _A; n.f = f; if(n.u == 0x80000000 || n.u == 0x00000000) /* 0.0/-0.0 */ return n.f; N = (n.u&(0xff<<23))>>23; if(N==0xff||(n.u&0x80000000)) { /* inf/-inf/nan/ f < 0.0*/ n.u = 0x7fc00000; /* nan */ return n.f; } if(N!=0x0) { /* 用科学计数法表示的规格数 */ A = (n.u&0x7fffff)|0x800000; _N = (int)N - 127; if(N&0x1) { _A = (uint64_t)A<<25; } else { _A = (uint64_t)A<<26; _N--; } } else { //A*2^(-149)这种表示方式的浮点数 //还是需要找最高位 for(i=22;;i--) if(n.u&((0x1)<<i)) break; //然后需要移位，要区分奇数和偶数 if(i&0x1) { _N = i-149; _A = (uint64_t)n.u << (48-i); } else { _N = i-150; _A = (uint64_t)n.u << (49-i); } } //小数部分 A = _sqrt_(_A); //四舍五入 A = (A+(A&0x1))>>1; //指数部分 N = (uint32_t)(_N/2+127); //得到结果 n.u = (A&0x7fffff)|(N<<23); return n.f; }
然后再测，准确无误。于是我们可以完工了。
标签：运算 四舍五入 回忆一下 bsp 测试结果 理论 a.out 进制 简单
原文地址：http://www.cnblogs.com/Colin-Cai/p/7223254.html
