ren1244的部落格

流程圖：

程式碼：

1. #include <stdint.h>
3. double msqrt(double x)
4. {//直式開方，適用 IEEE754 double
5.     int i;
6.     uint64_t r=0,k;
7.     
8.     //取得浮點數的s,e,f
9.     uint32_t s=(*(uint64_t*)&x)>>63;
10.     int32_t e=(*(uint64_t*)&x)>>52&0x7FF;
11.     uint64_t f=(*(uint64_t*)&x)&0xFFFFFFFFFFFFFllu;
12.     e-=1023;
13.     
14.     
15.     if(e==1024)
16.         return x;
17.     if(e==-1023)
18.     {
19.         if(f==0)
20.             return 0.0;
21.         else
22.         {
23.             //正規化
24.             while(!(f&0x10000000000000llu))
25.             {
26.                 f<<=1;
27.                 --e;
28.             }
29.         }
30.     }
31.     if(s==1)
32.     {
33.         k=0xFFFFFFFFFFFFFFFF;
34.         return *(double*)&k; //return NaN;
35.     }
36.     
37.     //位值調整，並把e開根號
38.     f=(e&1)?(f|0x10000000000000llu)<<1:(f|0x10000000000000llu);
39.     e=(e>>1)|(e&0x80000000);
40.     
41.     //對f開根號
42.     for(i=52;i>=0;i-=2)
43.     {
44.         k=(r<<2|1)<<i;
45.         if(f>=k)
46.         {
47.             f-=k;
48.             r=r<<1|1;
49.         }
50.         else
51.             r<<=1;
52.     }
53.     for(i=0;i<26;++i)
54.     {
55.         f<<=2;
56.         k=r<<2|1;
57.         if(f>=k)
58.         {
59.             f-=k;
60.             r=r<<1|1;
61.         }
62.         else
63.             r<<=1;
64.     }
65.     //四捨五入
66.     if( (f<<2) >= (r<<2|1) )
67.         ++r;
68.     
69.     //e,f合併為double回傳
70.     k=e+1023;
71.     k=(r&0xFFFFFFFFFFFFFllu)|k<<52;
72.     return *(double*)&k;
73. }