.module/boot=3/boot=4   square_root; 

{    SQRT.DSP - Calculate the squareroot

     INPUT:
       mr1 = msw of x in 16.16 format
       mr0 = lsw of x
       m1  = 5

     OUTPUT:
       sr1 = y in 8.8 unsigned format

     CALLED FUNCTIONS:
       None
     
     DESCRIPTION:
       Approximates the squareroot of x by a Taylor series
   
                y = sqrt(x)
   
     COMPUTATION TIME:
       75 cycles (maximum)
}

.const  BASE=h#0d49, SQRT2=h#5a82;

.var/pm sqrt_coeff[5];

.init   sqrt_coeff : h#5d1d00, h#a9ed00, h#46d600, h#ddaa00, h#072d00;

.entry  sqrt;

sqrt:   i6=^sqrt_coeff;  l6 = 0;              {pointer to coeff. buffer}
        se=exp mr1 (hi);                      {check for redundant bits}
        se=exp mr0 (lo);
        ax0=se, sr=norm mr1 (hi);             {remove redundant bits}
        sr=sr or norm mr0 (lo);
        my0=sr1, ar=pass sr1;
        if eq rts;
        mr=0;
        mr1=BASE;                             {load constant value}
        mf=ar*my0 (rnd), mx0=pm(i6,m5);       {mf = x**2}
        mr=mr+mx0*my0 (ss), mx0=pm(i6,m5);    {mr = BASE + c1*x}
        cntr=4;
        do approx until ce;
          mr=mr+mx0*mf (ss), mx0=pm(i6,m5);
        approx: mf=ar*mf (rnd);
        ay0=15;
        my0=mr1, ar=ax0+ay0;                 {se + 15 = 0?}
        if ne jump scale;                    {no, compute sqrt(s)}
        sr=ashift mr1 by -6 (hi);
        rts;
scale:  mr=0;
        mr1=SQRT2;                           {load 1/sqrt(2)}
        my1=mr1, ar=abs ar;
        ay0=ar;
        ar=ay0-1;
        if eq jump pwr_ok;
        cntr=ar;                             {compute (1/sqrt(2))^(se+15)}
        do compute until ce;
compute:  mr=mr1*my1 (rnd);
pwr_ok: if neg jump frac;
        ay1=h#0080;                          {load a 1 in 9.23 format}
        ay0=0;                               {compute reciprocal of mr}
        divs ay1, mr1;
        divq mr1; divq mr1; divq mr1;
        divq mr1; divq mr1; divq mr1;
        divq mr1; divq mr1; divq mr1;
        divq mr1; divq mr1; divq mr1;
        divq mr1; divq mr1; divq mr1;
        mx0=ay0;
        mr=0;
        mr0=h#2000;
        mr=mr+mx0*my0 (us);
        sr=ashift mr1 by 2 (hi);
        sr=sr or lshift mr0 by 2 (lo);
        rts;
frac:   mr=mr1*my0 (rnd);
        sr=ashift mr1 by -6 (hi);
        rts;
.endmod;