cmsMath.c revision 2f44462d 
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#ifdef DEBUG
41#include <stdio.h>
42#endif
43
44#include <float.h>
45#ifndef DBL_EPSILON
46#define DBL_EPSILON 1e-6
47#endif
48
49#ifdef _X_ROOT_STATS
50int cbrt_loopcount;
51int sqrt_loopcount;
52#endif
53
54/* Newton's Method:  x_n+1 = x_n - ( f(x_n) / f'(x_n) ) */
55
56
57/* for cube roots, x^3 - a = 0,  x_new = x - 1/3 (x - a/x^2) */
58
59double
60_XcmsCubeRoot(double a)
61{
62    register double abs_a, cur_guess, delta;
63
64#ifdef DEBUG
65    printf("_XcmsCubeRoot passed in %g\n", a);
66#endif
67#ifdef _X_ROOT_STATS
68    cbrt_loopcount = 0;
69#endif
70    if (a == 0.)
71	return 0.;
72
73    abs_a = a<0. ? -a : a;	/* convert to positive to speed loop tests */
74
75    /* arbitrary first guess */
76    if (abs_a > 1.)
77	cur_guess = abs_a/8.;
78    else
79	cur_guess = abs_a*8.;
80
81    do {
82#ifdef _X_ROOT_STATS
83	cbrt_loopcount++;
84#endif
85	delta = (cur_guess - abs_a/(cur_guess*cur_guess))/3.;
86	cur_guess -= delta;
87	if (delta < 0.) delta = -delta;
88    } while (delta >= cur_guess*DBL_EPSILON);
89
90    if (a < 0.)
91	cur_guess = -cur_guess;
92
93#ifdef DEBUG
94    printf("_XcmsCubeRoot returning %g\n", cur_guess);
95#endif
96    return cur_guess;
97}
98
99
100
101/* for square roots, x^2 - a = 0,  x_new = x - 1/2 (x - a/x) */
102
103double
104_XcmsSquareRoot(double a)
105{
106    register double cur_guess, delta;
107
108#ifdef DEBUG
109    printf("_XcmsSquareRoot passed in %g\n", a);
110#endif
111#ifdef _X_ROOT_STATS
112    sqrt_loopcount = 0;
113#endif
114    if (a == 0.)
115	return 0.;
116
117    if (a < 0.) {
118	/* errno = EDOM; */
119	return 0.;
120    }
121
122    /* arbitrary first guess */
123    if (a > 1.)
124	cur_guess = a/4.;
125    else
126	cur_guess = a*4.;
127
128    do {
129#ifdef _X_ROOT_STATS
130	sqrt_loopcount++;
131#endif
132	delta = (cur_guess - a/cur_guess)/2.;
133	cur_guess -= delta;
134	if (delta < 0.) delta = -delta;
135    } while (delta >= cur_guess*DBL_EPSILON);
136
137#ifdef DEBUG
138    printf("_XcmsSquareRoot returning %g\n", cur_guess);
139#endif
140    return cur_guess;
141}
142
143