cmsMath.c revision 1ab64890 
1/* $Xorg: cmsMath.c,v 1.4 2001/02/09 02:03:39 xorgcvs Exp $ */
2
3/*
4
5Copyright 1990, 1998  The Open Group
6
7Permission to use, copy, modify, distribute, and sell this software and its
8documentation for any purpose is hereby granted without fee, provided that
9the above copyright notice appear in all copies and that both that
10copyright notice and this permission notice appear in supporting
11documentation.
12
13The above copyright notice and this permission notice shall be included in
14all copies or substantial portions of the Software.
15
16THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
19OPEN GROUP BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
20AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
21CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22
23Except as contained in this notice, the name of The Open Group shall not be
24used in advertising or otherwise to promote the sale, use or other dealings
25in this Software without prior written authorization from The Open Group.
26
27*/
28/* $XFree86: xc/lib/X11/cmsMath.c,v 3.4 2001/07/29 05:01:11 tsi Exp $ */
29
30/*
31 * Stephen Gildea, MIT X Consortium, January 1991
32 */
33
34#ifdef HAVE_CONFIG_H
35#include <config.h>
36#endif
37#include "Xlibint.h"
38#include "Xcmsint.h"
39
40#include <float.h>
41#ifndef DBL_EPSILON
42#define DBL_EPSILON 1e-6
43#endif
44
45#ifdef _X_ROOT_STATS
46int cbrt_loopcount;
47int sqrt_loopcount;
48#endif
49
50/* Newton's Method:  x_n+1 = x_n - ( f(x_n) / f'(x_n) ) */
51
52
53/* for cube roots, x^3 - a = 0,  x_new = x - 1/3 (x - a/x^2) */
54
55double
56_XcmsCubeRoot(double a)
57{
58    register double abs_a, cur_guess, delta;
59
60#ifdef DEBUG
61    printf("_XcmsCubeRoot passed in %g\n", a);
62#endif
63#ifdef _X_ROOT_STATS
64    cbrt_loopcount = 0;
65#endif
66    if (a == 0.)
67	return 0.;
68
69    abs_a = a<0. ? -a : a;	/* convert to positive to speed loop tests */
70
71    /* arbitrary first guess */
72    if (abs_a > 1.)
73	cur_guess = abs_a/8.;
74    else
75	cur_guess = abs_a*8.;
76
77    do {
78#ifdef _X_ROOT_STATS
79	cbrt_loopcount++;
80#endif
81	delta = (cur_guess - abs_a/(cur_guess*cur_guess))/3.;
82	cur_guess -= delta;
83	if (delta < 0.) delta = -delta;
84    } while (delta >= cur_guess*DBL_EPSILON);
85
86    if (a < 0.)
87	cur_guess = -cur_guess;
88
89#ifdef DEBUG
90    printf("_XcmsCubeRoot returning %g\n", cur_guess);
91#endif
92    return cur_guess;
93}
94
95
96
97/* for square roots, x^2 - a = 0,  x_new = x - 1/2 (x - a/x) */
98
99double
100_XcmsSquareRoot(double a)
101{
102    register double cur_guess, delta;
103
104#ifdef DEBUG
105    printf("_XcmsSquareRoot passed in %g\n", a);
106#endif
107#ifdef _X_ROOT_STATS
108    sqrt_loopcount = 0;
109#endif
110    if (a == 0.)
111	return 0.;
112
113    if (a < 0.) {
114	/* errno = EDOM; */
115	return 0.;
116    }
117
118    /* arbitrary first guess */
119    if (a > 1.)
120	cur_guess = a/4.;
121    else
122	cur_guess = a*4.;
123
124    do {
125#ifdef _X_ROOT_STATS
126	sqrt_loopcount++;
127#endif
128	delta = (cur_guess - a/cur_guess)/2.;
129	cur_guess -= delta;
130	if (delta < 0.) delta = -delta;
131    } while (delta >= cur_guess*DBL_EPSILON);
132
133#ifdef DEBUG
134    printf("_XcmsSquareRoot returning %g\n", cur_guess);
135#endif
136    return cur_guess;
137}
138
139