xcmsmath.c

早期freebsd实现

/*
 * square and cube roots by Newton's method
 *
 * $XConsortium: XcmsMath.c,v 1.8 91/08/14 15:13:41 rws Exp $
 *
 * Copyright 1990 Massachusetts Institute of Technology
 *
 * Permission to use, copy, modify, distribute, and sell this software and its
 * documentation for any purpose is hereby granted without fee, provided that
 * the above copyright notice appear in all copies and that both that
 * copyright notice and this permission notice appear in supporting
 * documentation, and that the name of M.I.T. not be used in advertising or
 * publicity pertaining to distribution of the software without specific,
 * written prior permission.  M.I.T. makes no representations about the
 * suitability of this software for any purpose.  It is provided "as is"
 * without express or implied warranty.
 *
 * M.I.T. DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL M.I.T.
 * BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
 * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN 
 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 *
 * Stephen Gildea, MIT X Consortium, January 1991
 */

#include "Xlibint.h"
#include "Xcmsint.h"

#if !defined(X_NOT_STDC_ENV) && (__STDC__ || !(defined(sun) || (defined(sony) && !defined(SYSTYPE_SYSV))))
#include <float.h>
#endif
#ifndef DBL_EPSILON
#define DBL_EPSILON 1e-6
#endif

#ifdef _X_ROOT_STATS
int cbrt_loopcount;
int sqrt_loopcount;
#endif

/* Newton's Method:  x_n+1 = x_n - ( f(x_n) / f'(x_n) ) */


/* for cube roots, x^3 - a = 0,  x_new = x - 1/3 (x - a/x^2) */

double
_XcmsCubeRoot(a)
    double a;
{
    register double abs_a, cur_guess, delta;

#ifdef DEBUG
    printf("_XcmsCubeRoot passed in %g\n", a);
#endif
#ifdef _X_ROOT_STATS
    cbrt_loopcount = 0;
#endif
    if (a == 0.)
	return 0.;

    abs_a = a<0. ? -a : a;	/* convert to positive to speed loop tests */

    /* arbitrary first guess */
    if (abs_a > 1.)
	cur_guess = abs_a/8.;
    else
	cur_guess = abs_a*8.;

    do {
#ifdef _X_ROOT_STATS
	cbrt_loopcount++;
#endif
	delta = (cur_guess - abs_a/(cur_guess*cur_guess))/3.;
	cur_guess -= delta;
	if (delta < 0.) delta = -delta;
    } while (delta >= cur_guess*DBL_EPSILON);

    if (a < 0.)
	cur_guess = -cur_guess;

#ifdef DEBUG
    printf("_XcmsCubeRoot returning %g\n", cur_guess);
#endif
    return cur_guess;
}
	


/* for square roots, x^2 - a = 0,  x_new = x - 1/2 (x - a/x) */

double
_XcmsSquareRoot(a)
    double a;
{
    register double cur_guess, delta;

#ifdef DEBUG
    printf("_XcmsSquareRoot passed in %g\n", a);
#endif
#ifdef _X_ROOT_STATS
    sqrt_loopcount = 0;
#endif
    if (a == 0.)
	return 0.;

    if (a < 0.) {
	/* errno = EDOM; */
	return 0.;
    }

    /* arbitrary first guess */
    if (a > 1.)
	cur_guess = a/4.;
    else
	cur_guess = a*4.;

    do {
#ifdef _X_ROOT_STATS
	sqrt_loopcount++;
#endif
	delta = (cur_guess - a/cur_guess)/2.;
	cur_guess -= delta;
	if (delta < 0.) delta = -delta;
    } while (delta >= cur_guess*DBL_EPSILON);

#ifdef DEBUG
    printf("_XcmsSquareRoot returning %g\n", cur_guess);
#endif
    return cur_guess;
}