版權申明：本文為博主窗戶(Colin Cai)原創，歡迎轉帖。如要轉貼，必須注明原文網址 http://www.cnblogs.com/Colin-Cai/p/7223254.html 作者：窗戶 QQ：6679072 E-mail：6679072@qq.com
了解了浮點數的存儲以及手算平方根的原理，我們可以考慮程序實現了。
先實現一個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; }
然后再測，准確無誤。於是我們可以完工了。