1
2 /*
3
4 Copyright 1990, 1998 The Open Group
5
6 Permission to use, copy, modify, distribute, and sell this software and its
7 documentation for any purpose is hereby granted without fee, provided that
8 the above copyright notice appear in all copies and that both that
9 copyright notice and this permission notice appear in supporting
10 documentation.
11
12 The above copyright notice and this permission notice shall be included in
13 all copies or substantial portions of the Software.
14
15 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18 OPEN GROUP BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
19 AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
20 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21
22 Except as contained in this notice, the name of The Open Group shall not be
23 used in advertising or otherwise to promote the sale, use or other dealings
24 in 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
48 int cbrt_loopcount;
49 int 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
57 double
_XcmsCubeRoot(double a)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
101 double
_XcmsSquareRoot(double a)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