source: trunk/minix/lib/ack/libp/sqt.c@ 10

/*
 * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 * See the copyright notice in the ACK home directory, in the file "Copyright".
 *
 * Author: Ceriel J.H. Jacobs
 */

/* $Header: /cvsup/minix/src/lib/ack/libp/sqt.c,v 1.1 2005/10/10 15:27:47 beng Exp $ */
#define __NO_DEFS
#include <math.h>
#include <pc_err.h>
extern  _trp();

#define NITER   5

static double
Ldexp(fl,exp)
        double fl;
        int exp;
{
        extern double _fef();
        int sign = 1;
        int currexp;

        if (fl<0) {
                fl = -fl;
                sign = -1;
        }
        fl = _fef(fl,&currexp);
        exp += currexp;
        if (exp > 0) {
                while (exp>30) {
                        fl *= (double) (1L << 30);
                        exp -= 30;
                }
                fl *= (double) (1L << exp);
        }
        else    {
                while (exp<-30) {
                        fl /= (double) (1L << 30);
                        exp += 30;
                }
                fl /= (double) (1L << -exp);
        }
        return sign * fl;
}

double
_sqt(x)
        double x;
{
        extern double _fef();
        int exponent;
        double val;

        if (x <= 0) {
                if (x < 0) _trp(ESQT);
                return 0;
        }

        val = _fef(x, &exponent);
        if (exponent & 1) {
                exponent--;
                val *= 2;
        }
        val = Ldexp(val + 1.0, exponent/2 - 1);
        /* was: val = (val + 1.0)/2.0; val = Ldexp(val, exponent/2); */
        for (exponent = NITER - 1; exponent >= 0; exponent--) {
                val = (val + x / val) / 2.0;
        }
        return val;
}