Subj : Re: sqrt algo
To   : comp.programming
From : pete
Date : Thu Oct 06 2005 01:59 am

pete wrote:
> 
> himanshu bafna wrote:
> >
> > Hi
> >
> > Can anyone shed light on an Algorithm for
> > finding out the square root of
> > a number? to make it simple lets only consider integers.
> 
> /* BEGIN isqrt.c */

... and now for a float sqrt algorithm:

/* BEGIN ma_th.h */

#ifndef H_MA_TH_
#define H_MA_TH_

double sq_rt(double x);
double p_ow(double x, double y);
double l_og(double x);
double e_xp(double x);

#endif

/* END ma_th.h */


/* BEGIN ma_th.c */

#include <float.h>
#include <errno.h>
#include <math.h>

#include "ma_th.h"

double sq_rt(double x)
{
    return p_ow(x, 0.5);
}

double p_ow(double x, double y)
{
    double p;

    if (0 > x && floor(y) == ceil(y)) {
        if (floor(y / 2) == ceil(y / 2) ) {
            p =  e_xp(l_og(-x) * y);
        } else {
            p = -e_xp(l_og(-x) * y);
        }
    } else {
        if (x != 0 || 0 >= y) {
            p =  e_xp(l_og( x) * y);
        } else {
            p =  0.0;
        }
    }
    return p;
}

double l_og(double x)
{
    int n;
    double a, b, c, epsilon;
    static double A, B, C;
    
    if (DBL_MAX >= x) {
        if (x > 0) {
            if (1 > A) {
                A = 1.5;
                do {
                    B = A;
                    A = 1 / B + B / 2;
                } while (B > A);    
                B /= 2;             
                C = l_og(A);          
            }
            for (n = 0; x > A; x /= 2) {
                ++n;
            }
            while (B > x) {
                --n;
                x *= 2;
            }
            a = (x - 1) / (x + 1);
            x = C * n + a;
            c = a * a;
            n = 1;
            epsilon = DBL_EPSILON * x;
            if (0 > a) {
                if (epsilon > 0) {
                    epsilon = -epsilon;
                }
                do {
                    n += 2;
                    a *= c;
                    b = a / n;
                    x += b;
                } while (epsilon > b);
            } else {
                if (0 > epsilon) {
                    epsilon = -epsilon;
                }
                do {
                    n += 2;
                    a *= c;
                    b = a / n;
                    x += b;
                } while (b > epsilon);
            }
            x *= 2;
        } else {
            errno = 0 > x ? EDOM : ERANGE;
            x = -HUGE_VAL;
        }
    } else {
        errno = ERANGE;
    }
    return x;
}

double e_xp(double x)
{
    unsigned n, negative, square; 
    double b, e;
    static double x_max;

    if (0 > x) {
        x = -x;
        negative = 1;
    } else {
        negative = 0;
    }
    if (1 > x_max) {
        x_max = l_og(DBL_MAX);
    }
    if (x_max >= x) {        
        for (square = 0; x > 0.25; x /= 2) {
            ++square;
        }
        e = b = n = 1;
        do {
            b /= n++;
            b *= x;
            e += b;
        } while (b > DBL_EPSILON);
        while (square-- != 0) {
            e *= e;
        }
        if (negative) {
            e = 1 / e;
        }
    } else {
        errno = ERANGE;        
        e = negative ? 0.0 : HUGE_VAL;
    }
    return e;
}

/* END ma_th.c */

-- 
pete

.