C language implementation of square root (iii.)--Final program implementation

Knowing the storage of floating-point numbers and the principle of square root of hand, we can consider the implementation of the program.

The square root of a 64-bit integer is implemented first, and the program is not so difficult to write according to the square root of the previous hand.

#include <stdint.h>
uint64_t _sqrt_u64 (uint64_t a) {int i;        uint64_t Res;        Uint64_t remain;        The square root of 0 is 0, special handling if (a = = 0ull) return 0ull; Find the highest 1 and produce a square root result with a maximum bit of 1 for (i=62;;                        i-=2) if (a& (3ull<<i)) {res = 1ull;                        remain = ((a& (3ull<<i)) >>i)-1ull;                        I-= 2;                Break  }//According to the principle of the square root of the hand, generate the results for (; i>=0;i-=2) {//Right move two-bit, and the two bits of a followed by the remain remain = (REMAIN&LT;&LT;2) |                ((a& (3ull<<i)) >>i); if ((res<<2) |1ull) <= remain) {//Generate a new bit of 1 remain = remain-((RES&L                        T;&LT;2) |1ull);                res = (res<<1) |1ull;                } else {//generates a new bit of 0 res <<= 1; }} return res;}

In fact, you can write together, the code will be shorter, but the efficiency will be low a little bit, and the compiler should not be easy to optimize.

#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;}

However, we do not need this result.

To verify its correctness, let's write a C-language 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;}

Our shell program tests, we certainly can not test every 64bits number, this calculation is too large, unrealistic. We can use a random part to test.

#!/bin/bash# CompilingGCC-o2-s sqrt_u64.c main_sqrt_u64.c-o a.out# random test 10,000 numbers for((i=0;i<10000; i++)); Do#随机产生bits0~ -, if it is 0, the number that represents the test is 0 #如果不是0, which represents how many bits of the binary to generate=random% $        if[$bits-eq0]; Thenx=0y=0        Else#产生一个bits位的二进制数x x=$ ({#最高位1Echo-N1#之后每位随机产生 for((j=1; j<bits;j++)); DoLet x=random%2                                Echo-N $x Done}) #用bc将x转换成十进制 x=$(Echo 'obase=10;ibase=2;' "$x"|BC) #用bc计算x的平方根取整, theoretically consistent with our C language calculation y=$(Echo 'sqrt (' "$x"')'|BC)fi#z是我们的C语言计算结果 z=$(Echo$x |./a.out) #比较, if not consistent, the errorif[$y-ne $z]; Then                Echo$x $y $z Error Exit1        fi DoneEchoOk

The test results show that our C language can still get the correct results.

Recall the floating-point structure mentioned in the first section,

　S (1bits) |  N (8bits) | A (23bits)

For floating-point a*2n,

1<=a<2,n is an integer,

If n is an even number,

Then the square root of a*2n is sqrt (a) *2N/2, also satisfies 1<=sqrt (a) <2,N/2 is an integer;

If n is an odd number,

Then the square root of a*2n is sqrt (2*a) (n-1)/2, also satisfies 1<=sqrt (2*a) <2, (n-1)/2 is an integer.

So here we use a or 2*a to open the square root,

Recall the structure of floating-point numbers, the precision of single-precision floating-point numbers is 23 bits.

Represents a portion of the scientific notation a*2n minus 1, then an integer 1 can be represented by a binary 24 bit.

So, we thought, a binary 48-bit or 47-bit long, square root is a binary 24-bit. Then we can calculate the fractional part of the result with a square root of a 48-bit or 47-bit binary integer.

The square root of nan/inf/-inf and negative numbers is Nan,

The square root of 0.0 is 0.0,

The square root of 0.0 is-0.0 (which is probably just the case in some libraries),

All of these can be specialized at the time of calculation.

The square root of the number of specifications (which is the floating-point number represented by scientific notation) is also the size,

S=0,n=0,a>0 stands for a*2-149, which is (a*2) *2-150,

Let's take a little bit of calculation to see that the square root of all such numbers is within the range indicated by the number of specifications.

As a result, the following code is available.

　　

#include <stdint.h>
Static uint32_t _sqrt_ (uint64_t a) {int i;        uint64_t Res;        Uint64_t remain;        res = remain = 0ull; Before the square root of an integer is directly optimized, we simply require a 47-bit or 48-bit integer for the square root for (i=46;i>=0;i-=2) {remain = (REMAIN&LT;&LT;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) {/* The number of specifications represented by scientific notation */A = (N.U&AMP;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) This representation of the floating-point number//or need to find the highest bit for (i=22;;                i--) if (n.u& (0x1) <<i) break;                        The shift is then required to differentiate between odd and even if (i&0x1) {_n = i-149;                _a = (uint64_t) n.u << (46-i);                        } else {_n = i-150;                _a = (uint64_t) n.u << (47-i);        }}//Decimal part A = _sqrt_ (_a);        Exponent part N = (uint32_t) (_n/2+127); Get Results n.u = (A&AMP;0X7FFFFF) |        (n<<23); RetUrn N.f;} 

Also, write a test program, for inf/-inf/nan/0.0/-0.0 and negative, these are very simple.

#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));        Random 10,000 data 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;}

It turns out that the SQRTF results in our program and in the math library are slightly different.

So, we decided to add a small thing, that is, rounding. Before we used a 47-bit or 48-digit open square, in order to round, we need one more, so we use 49 or 50 digits square.

Modify the Mysqrtf, add two bit take to open square, _sqrt_ also move a bit.

#include <stdint.h>
Static uint32_t _sqrt_ (uint64_t a) {        int i;        uint64_t Res;        uint64_t remain;        res = remain = 0ull;        Before the square root of an integer is directly optimized, we simply require a 49-bit or 50-bit integer for the square root for        (i=48;i>=0;i-=2) {//Here is 46, changed to                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) {/* The number of specifications represented by scientific notation */A = (N.U&AMP;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) This representation of the floating-point number//or need to find the highest bit for (i=22;;                i--) if (n.u& (0x1) <<i) break; Then you need toTo shift, distinguish between odd and even if (i&0x1) {_n = i-149;                _a = (uint64_t) n.u << (48-i);                        } else {_n = i-150;                _a = (uint64_t) n.u << (49-i);        }}//Decimal part A = _sqrt_ (_a);        Rounding A = (A + (a&0x1)) >>1;        Exponent part N = (uint32_t) (_n/2+127); Get Results n.u = (A&AMP;0X7FFFFF) |        (n<<23); return N.F;}

Then test again, accurate. So we can finish it.

C language implementation of square root (iii.)--Final program implementation