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