Saturday, August 2, 2014

Fixed-point square-root

Simple implementation ...
Continuing from the previously posted fixed-point type ... , (UPDATE: My apology, because there was an error ... -- It's been fixed .. > <)

Output (reformatted, partial):
tRX4::sqrt(0.0000)=    0.0000   < 0          >   0          =sqrt(0)
tRX4::sqrt(0.1250)=    0.3515   <-0.00199088 >   0.353553   =sqrt(0.125)
tRX4::sqrt(0.2500)=    0.5000   < 0          >   0.5        =sqrt(0.25)
tRX4::sqrt(0.3750)=    0.6093   <-0.00299746 >   0.612372   =sqrt(0.375)
tRX4::sqrt(0.5000)=    0.7070   <-7.55191e-05>   0.707107   =sqrt(0.5)
tRX4::sqrt(0.6250)=    0.7890   <-0.00150692 >   0.790569   =sqrt(0.625)
tRX4::sqrt(0.7500)=    0.8632   <-0.00274414 >   0.866025   =sqrt(0.75)
tRX4::sqrt(0.8750)=    0.9335   <-0.00182062 >   0.935414   =sqrt(0.875)
tRX4::sqrt(1.0000)=    1.0000   < 0          >   1          =sqrt(1)
tRX4::sqrt(1.1250)=    1.0585   <-0.00206637 >   1.06066    =sqrt(1.125)
...

tRX4::sqrt(2.0000)=    1.4140   <-0.000151038>   1.41421    =sqrt(2)
tRX4::sqrt(2.1250)=    1.4570   <-0.000706673>   1.45774    =sqrt(2.125)
tRX4::sqrt(2.2500)=    1.5000   < 0          >   1.5        =sqrt(2.25)
tRX4::sqrt(2.3750)=    1.5390   <-0.00204098 >   1.5411     =sqrt(2.375)
...

tRX4::sqrt(3.8750)=    1.9648   <-0.00365818 >   1.9685     =sqrt(3.875)
tRX4::sqrt(4.0000)=    2.0000   < 0          >   2          =sqrt(4)
tRX4::sqrt(4.1250)=    2.0273   <-0.00366592 >   2.03101    =sqrt(4.125)
...

tRX4::sqrt(6.1250)=    2.4726   <-0.00221753 >   2.47487    =sqrt(6.125)
tRX4::sqrt(6.2500)=    2.5000   < 0          >   2.5        =sqrt(6.25)
tRX4::sqrt(6.3750)=    2.5234   <-0.00143862 >   2.52488    =sqrt(6.375)
...

tRX4::sqrt(8.0000)=    2.8281   <-0.000302076>   2.82843    =sqrt(8)

Code:
typedef uint32_t t4u;
typedef int32_t t4s;
typedef float t4f;

struct tRX4;/*The fixed-point type is a struct; Refer to the previous tutorial*/
typedef t4u tWU4;/*whole number, unsigned, 4 bytes*/
typedef t4s tII4;/*Integer, signed, 4 bytes*/

/*Would it ever be approriate, to implement it as a struct?
struct tWU4;
struct tII4;*/


static t4u wsLog2(t4u p)noexcept{static t4u const n[0x20]={
0x00,0x09,0x01,0x0A,0x0D,0x15,0x02,0x1D,0x0B,0x0E,0x10,0x12,0x16,0x19,0x03,0x1E,
0x08,0x0C,0x14,0x1C,0x0F,0x11,0x18,0x07,0x13,0x1B,0x17,0x06,0x1A,0x05,0x04,0x1F};
p|=p>>0x01,p|=p>>0x02,p|=p>>0x04,p|=p>>0x08,p|=p>>0x10;return n[(p*0x07C4ACDDU)>>0x1B];}

static tRX4 sqrt(tRX4 const&p)noexcept{
    tRX4 l=(tRX4)tWU4(1<<(wsLog2(t4u(p.g())>>tRX4::cFc)>>1));
    /*With mantissa of 8 bits (see previous tutorial page for fixed-point type),
    according to tests, the maximum loop-sum is 9.
    If you can guarantee that the parameter will never be zero,
    then you can slightly increase this loop-sum value for accuracy*/
    t4u i=0x09;

    while(i--){
        l=(p/l+l)/tWU4(2);}

    return l;}

#include<cmath>

static t4f fsf4(tRX4 const&p)noexcept{return(t4f)p.g()/(t4f)tRX4::c1.g();}

static v f()noexcept{
tRX4 l1(0),l2;t4f m1=0.F,m2;

for(;m1<=8.F;m1+=0.125F,l1+=tRX4(0.125F)){

    l2=sqrt(l1),m2=std::sqrt(m1);
    std::cout<<"tRX4::sqrt("<<l1<< ")=    "<<l2<<"   <" 

    <<(fsf4(l2)-m2)<<">   "<<m2<<"    =sqrt("<<m1<<")"<<std::endl;}}

int main(){f();getchar();return 0;}

No comments:

Post a Comment