1
2 /*
3 * Code and supporting documentation (c) Copyright 1990 1991 Tektronix, Inc.
4 * All Rights Reserved
5 *
6 * This file is a component of an X Window System-specific implementation
7 * of Xcms based on the TekColor Color Management System. Permission is
8 * hereby granted to use, copy, modify, sell, and otherwise distribute this
9 * software and its documentation for any purpose and without fee, provided
10 * that this copyright, permission, and disclaimer notice is reproduced in
11 * all copies of this software and in supporting documentation. TekColor
12 * is a trademark of Tektronix, Inc.
13 *
14 * Tektronix makes no representation about the suitability of this software
15 * for any purpose. It is provided "as is" and with all faults.
16 *
17 * TEKTRONIX DISCLAIMS ALL WARRANTIES APPLICABLE TO THIS SOFTWARE,
18 * INCLUDING THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
19 * PARTICULAR PURPOSE. IN NO EVENT SHALL TEKTRONIX BE LIABLE FOR ANY
20 * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER
21 * RESULTING FROM LOSS OF USE, DATA, OR PROFITS, WHETHER IN AN ACTION OF
22 * CONTRACT, NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
23 * CONNECTION WITH THE USE OR THE PERFORMANCE OF THIS SOFTWARE.
24 */
25
26 /*
27 * It should be pointed out that for simplicity's sake, the
28 * environment parameters are defined as floating point constants,
29 * rather than octal or hexadecimal initializations of allocated
30 * storage areas. This means that the range of allowed numbers
31 * may not exactly match the hardware's capabilities. For example,
32 * if the maximum positive double precision floating point number
33 * is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is
34 * defined to be 1.11E100 then the numbers between 1.11E100 and
35 * 1.11...E100 are considered to be undefined. For most
36 * applications, this will cause no problems.
37 *
38 * An alternate method is to allocate a global static "double" variable,
39 * say "maxdouble", and use a union declaration and initialization
40 * to initialize it with the proper bits for the EXACT maximum value.
41 * This was not done because the only compilers available to the
42 * author did not fully support union initialization features.
43 *
44 */
45
46 #ifdef HAVE_CONFIG_H
47 #include <config.h>
48 #endif
49 #include "Xcmsint.h"
50
51 /* forward/static */
52 static double _XcmsModulo(double value, double base);
53 static double _XcmsPolynomial(
54 register int order,
55 double const *coeffs,
56 double x);
57 static double
58 _XcmsModuloF(
59 double val,
60 register double *dp);
61
62 /*
63 * DEFINES
64 */
65 #define XCMS_MAXERROR 0.000001
66 #define XCMS_MAXITER 10000
67 #define XCMS_PI 3.14159265358979323846264338327950
68 #define XCMS_TWOPI 6.28318530717958620
69 #define XCMS_HALFPI 1.57079632679489660
70 #define XCMS_FOURTHPI 0.785398163397448280
71 #define XCMS_SIXTHPI 0.523598775598298820
72 #define XCMS_RADIANS(d) ((d) * XCMS_PI / 180.0)
73 #define XCMS_DEGREES(r) ((r) * 180.0 / XCMS_PI)
74 #ifdef __vax__
75 #define XCMS_X6_UNDERFLOWS (3.784659e-07) /* X**6 almost underflows*/
76 #else
77 #define XCMS_X6_UNDERFLOWS (4.209340e-52) /* X**6 almost underflows */
78 #endif
79 #define XCMS_X16_UNDERFLOWS (5.421010e-20) /* X**16 almost underflows*/
80 #define XCMS_CHAR_BIT 8
81 #define XCMS_LONG_MAX 0x7FFFFFFF
82 #define XCMS_DEXPLEN 11
83 #define XCMS_NBITS(type) (XCMS_CHAR_BIT * (int)sizeof(type))
84 #define XCMS_FABS(x) ((x) < 0.0 ? -(x) : (x))
85
86 /* XCMS_DMAXPOWTWO - largest power of two exactly representable as a double */
87 #define XCMS_DMAXPOWTWO ((double)(XCMS_LONG_MAX) * \
88 (1L << ((XCMS_NBITS(double)-XCMS_DEXPLEN) - XCMS_NBITS(int) + 1)))
89
90 /*
91 * LOCAL VARIABLES
92 */
93
94 static double const cos_pcoeffs[] = {
95 0.12905394659037374438e7,
96 -0.37456703915723204710e6,
97 0.13432300986539084285e5,
98 -0.11231450823340933092e3
99 };
100
101 static double const cos_qcoeffs[] = {
102 0.12905394659037373590e7,
103 0.23467773107245835052e5,
104 0.20969518196726306286e3,
105 1.0
106 };
107
108 static double const sin_pcoeffs[] = {
109 0.20664343336995858240e7,
110 -0.18160398797407332550e6,
111 0.35999306949636188317e4,
112 -0.20107483294588615719e2
113 };
114
115 static double const sin_qcoeffs[] = {
116 0.26310659102647698963e7,
117 0.39270242774649000308e5,
118 0.27811919481083844087e3,
119 1.0
120 };
121
122 /*
123 *
124 * FUNCTION
125 *
126 * _XcmsCosine double precision cosine
127 *
128 * KEY WORDS
129 *
130 * cos
131 * machine independent routines
132 * trigonometric functions
133 * math libraries
134 *
135 * DESCRIPTION
136 *
137 * Returns double precision cosine of double precision
138 * floating point argument.
139 *
140 * USAGE
141 *
142 * double _XcmsCosine (x)
143 * double x;
144 *
145 * REFERENCES
146 *
147 * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
148 * 1968, pp. 112-120.
149 *
150 * RESTRICTIONS
151 *
152 * The sin and cos routines are interactive in the sense that
153 * in the process of reducing the argument to the range -PI/4
154 * to PI/4, each may call the other. Ultimately one or the
155 * other uses a polynomial approximation on the reduced
156 * argument. The sin approximation has a maximum relative error
157 * of 10**(-17.59) and the cos approximation has a maximum
158 * relative error of 10**(-16.18).
159 *
160 * These error bounds assume exact arithmetic
161 * in the polynomial evaluation. Additional rounding and
162 * truncation errors may occur as the argument is reduced
163 * to the range over which the polynomial approximation
164 * is valid, and as the polynomial is evaluated using
165 * finite-precision arithmetic.
166 *
167 * PROGRAMMER
168 *
169 * Fred Fish
170 *
171 * INTERNALS
172 *
173 * Computes cos(x) from:
174 *
175 * (1) Reduce argument x to range -PI to PI.
176 *
177 * (2) If x > PI/2 then call cos recursively
178 * using relation cos(x) = -cos(x - PI).
179 *
180 * (3) If x < -PI/2 then call cos recursively
181 * using relation cos(x) = -cos(x + PI).
182 *
183 * (4) If x > PI/4 then call sin using
184 * relation cos(x) = sin(PI/2 - x).
185 *
186 * (5) If x < -PI/4 then call cos using
187 * relation cos(x) = sin(PI/2 + x).
188 *
189 * (6) If x would cause underflow in approx
190 * evaluation arithmetic then return
191 * sqrt(1.0 - x**2).
192 *
193 * (7) By now x has been reduced to range
194 * -PI/4 to PI/4 and the approximation
195 * from HART pg. 119 can be used:
196 *
197 * cos(x) = ( p(y) / q(y) )
198 * Where:
199 *
200 * y = x * (4/PI)
201 *
202 * p(y) = SUM [ Pj * (y**(2*j)) ]
203 * over j = {0,1,2,3}
204 *
205 * q(y) = SUM [ Qj * (y**(2*j)) ]
206 * over j = {0,1,2,3}
207 *
208 * P0 = 0.12905394659037374438571854e+7
209 * P1 = -0.3745670391572320471032359e+6
210 * P2 = 0.134323009865390842853673e+5
211 * P3 = -0.112314508233409330923e+3
212 * Q0 = 0.12905394659037373590295914e+7
213 * Q1 = 0.234677731072458350524124e+5
214 * Q2 = 0.2096951819672630628621e+3
215 * Q3 = 1.0000...
216 * (coefficients from HART table #3843 pg 244)
217 *
218 *
219 * **** NOTE **** The range reduction relations used in
220 * this routine depend on the final approximation being valid
221 * over the negative argument range in addition to the positive
222 * argument range. The particular approximation chosen from
223 * HART satisfies this requirement, although not explicitly
224 * stated in the text. This may not be true of other
225 * approximations given in the reference.
226 *
227 */
228
_XcmsCosine(double x)229 double _XcmsCosine(double x)
230 {
231 auto double y;
232 auto double yt2;
233 double retval;
234
235 if (x < -XCMS_PI || x > XCMS_PI) {
236 x = _XcmsModulo (x, XCMS_TWOPI);
237 if (x > XCMS_PI) {
238 x = x - XCMS_TWOPI;
239 } else if (x < -XCMS_PI) {
240 x = x + XCMS_TWOPI;
241 }
242 }
243 if (x > XCMS_HALFPI) {
244 retval = -(_XcmsCosine (x - XCMS_PI));
245 } else if (x < -XCMS_HALFPI) {
246 retval = -(_XcmsCosine (x + XCMS_PI));
247 } else if (x > XCMS_FOURTHPI) {
248 retval = _XcmsSine (XCMS_HALFPI - x);
249 } else if (x < -XCMS_FOURTHPI) {
250 retval = _XcmsSine (XCMS_HALFPI + x);
251 } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
252 retval = _XcmsSquareRoot (1.0 - (x * x));
253 } else {
254 y = x / XCMS_FOURTHPI;
255 yt2 = y * y;
256 retval = _XcmsPolynomial (3, cos_pcoeffs, yt2) / _XcmsPolynomial (3, cos_qcoeffs, yt2);
257 }
258 return (retval);
259 }
260
261
262 /*
263 * FUNCTION
264 *
265 * _XcmsModulo double precision modulo
266 *
267 * KEY WORDS
268 *
269 * _XcmsModulo
270 * machine independent routines
271 * math libraries
272 *
273 * DESCRIPTION
274 *
275 * Returns double precision modulo of two double
276 * precision arguments.
277 *
278 * USAGE
279 *
280 * double _XcmsModulo (value, base)
281 * double value;
282 * double base;
283 *
284 * PROGRAMMER
285 *
286 * Fred Fish
287 *
288 */
_XcmsModulo(double value,double base)289 static double _XcmsModulo(double value, double base)
290 {
291 auto double intpart;
292
293 value /= base;
294 value = _XcmsModuloF (value, &intpart);
295 value *= base;
296 return(value);
297 }
298
299
300 /*
301 * frac = (double) _XcmsModuloF(double val, double *dp)
302 * return fractional part of 'val'
303 * set *dp to integer part of 'val'
304 *
305 * Note -> only compiled for the CA or KA. For the KB/MC,
306 * "math.c" instantiates a copy of the inline function
307 * defined in "math.h".
308 */
309 static double
_XcmsModuloF(double val,register double * dp)310 _XcmsModuloF(
311 double val,
312 register double *dp)
313 {
314 register double abs;
315 /*
316 * Don't use a register for this. The extra precision this results
317 * in on some systems causes problems.
318 */
319 double ip;
320
321 /* should check for illegal values here - nan, inf, etc */
322 abs = XCMS_FABS(val);
323 if (abs >= XCMS_DMAXPOWTWO) {
324 ip = val;
325 } else {
326 ip = abs + XCMS_DMAXPOWTWO; /* dump fraction */
327 ip -= XCMS_DMAXPOWTWO; /* restore w/o frac */
328 if (ip > abs) /* if it rounds up */
329 ip -= 1.0; /* fix it */
330 ip = XCMS_FABS(ip);
331 }
332 *dp = ip;
333 return (val - ip); /* signed fractional part */
334 }
335
336
337 /*
338 * FUNCTION
339 *
340 * _XcmsPolynomial double precision polynomial evaluation
341 *
342 * KEY WORDS
343 *
344 * poly
345 * machine independent routines
346 * math libraries
347 *
348 * DESCRIPTION
349 *
350 * Evaluates a polynomial and returns double precision
351 * result. Is passed a the order of the polynomial,
352 * a pointer to an array of double precision polynomial
353 * coefficients (in ascending order), and the independent
354 * variable.
355 *
356 * USAGE
357 *
358 * double _XcmsPolynomial (order, coeffs, x)
359 * int order;
360 * double *coeffs;
361 * double x;
362 *
363 * PROGRAMMER
364 *
365 * Fred Fish
366 *
367 * INTERNALS
368 *
369 * Evalates the polynomial using recursion and the form:
370 *
371 * P(x) = P0 + x(P1 + x(P2 +...x(Pn)))
372 *
373 */
374
_XcmsPolynomial(register int order,double const * coeffs,double x)375 static double _XcmsPolynomial(
376 register int order,
377 double const *coeffs,
378 double x)
379 {
380 auto double rtn_value;
381
382 coeffs += order;
383 rtn_value = *coeffs--;
384 while(order-- > 0)
385 rtn_value = *coeffs-- + (x * rtn_value);
386
387 return(rtn_value);
388 }
389
390
391 /*
392 * FUNCTION
393 *
394 * _XcmsSine double precision sine
395 *
396 * KEY WORDS
397 *
398 * sin
399 * machine independent routines
400 * trigonometric functions
401 * math libraries
402 *
403 * DESCRIPTION
404 *
405 * Returns double precision sine of double precision
406 * floating point argument.
407 *
408 * USAGE
409 *
410 * double _XcmsSine (x)
411 * double x;
412 *
413 * REFERENCES
414 *
415 * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
416 * 1968, pp. 112-120.
417 *
418 * RESTRICTIONS
419 *
420 * The sin and cos routines are interactive in the sense that
421 * in the process of reducing the argument to the range -PI/4
422 * to PI/4, each may call the other. Ultimately one or the
423 * other uses a polynomial approximation on the reduced
424 * argument. The sin approximation has a maximum relative error
425 * of 10**(-17.59) and the cos approximation has a maximum
426 * relative error of 10**(-16.18).
427 *
428 * These error bounds assume exact arithmetic
429 * in the polynomial evaluation. Additional rounding and
430 * truncation errors may occur as the argument is reduced
431 * to the range over which the polynomial approximation
432 * is valid, and as the polynomial is evaluated using
433 * finite-precision arithmetic.
434 *
435 * PROGRAMMER
436 *
437 * Fred Fish
438 *
439 * INTERNALS
440 *
441 * Computes sin(x) from:
442 *
443 * (1) Reduce argument x to range -PI to PI.
444 *
445 * (2) If x > PI/2 then call sin recursively
446 * using relation sin(x) = -sin(x - PI).
447 *
448 * (3) If x < -PI/2 then call sin recursively
449 * using relation sin(x) = -sin(x + PI).
450 *
451 * (4) If x > PI/4 then call cos using
452 * relation sin(x) = cos(PI/2 - x).
453 *
454 * (5) If x < -PI/4 then call cos using
455 * relation sin(x) = -cos(PI/2 + x).
456 *
457 * (6) If x is small enough that polynomial
458 * evaluation would cause underflow
459 * then return x, since sin(x)
460 * approaches x as x approaches zero.
461 *
462 * (7) By now x has been reduced to range
463 * -PI/4 to PI/4 and the approximation
464 * from HART pg. 118 can be used:
465 *
466 * sin(x) = y * ( p(y) / q(y) )
467 * Where:
468 *
469 * y = x * (4/PI)
470 *
471 * p(y) = SUM [ Pj * (y**(2*j)) ]
472 * over j = {0,1,2,3}
473 *
474 * q(y) = SUM [ Qj * (y**(2*j)) ]
475 * over j = {0,1,2,3}
476 *
477 * P0 = 0.206643433369958582409167054e+7
478 * P1 = -0.18160398797407332550219213e+6
479 * P2 = 0.359993069496361883172836e+4
480 * P3 = -0.2010748329458861571949e+2
481 * Q0 = 0.263106591026476989637710307e+7
482 * Q1 = 0.3927024277464900030883986e+5
483 * Q2 = 0.27811919481083844087953e+3
484 * Q3 = 1.0000...
485 * (coefficients from HART table #3063 pg 234)
486 *
487 *
488 * **** NOTE **** The range reduction relations used in
489 * this routine depend on the final approximation being valid
490 * over the negative argument range in addition to the positive
491 * argument range. The particular approximation chosen from
492 * HART satisfies this requirement, although not explicitly
493 * stated in the text. This may not be true of other
494 * approximations given in the reference.
495 *
496 */
497
498 double
_XcmsSine(double x)499 _XcmsSine (double x)
500 {
501 double y;
502 double yt2;
503 double retval;
504
505 if (x < -XCMS_PI || x > XCMS_PI) {
506 x = _XcmsModulo (x, XCMS_TWOPI);
507 if (x > XCMS_PI) {
508 x = x - XCMS_TWOPI;
509 } else if (x < -XCMS_PI) {
510 x = x + XCMS_TWOPI;
511 }
512 }
513 if (x > XCMS_HALFPI) {
514 retval = -(_XcmsSine (x - XCMS_PI));
515 } else if (x < -XCMS_HALFPI) {
516 retval = -(_XcmsSine (x + XCMS_PI));
517 } else if (x > XCMS_FOURTHPI) {
518 retval = _XcmsCosine (XCMS_HALFPI - x);
519 } else if (x < -XCMS_FOURTHPI) {
520 retval = -(_XcmsCosine (XCMS_HALFPI + x));
521 } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
522 retval = x;
523 } else {
524 y = x / XCMS_FOURTHPI;
525 yt2 = y * y;
526 retval = y * (_XcmsPolynomial (3, sin_pcoeffs, yt2) / _XcmsPolynomial(3, sin_qcoeffs, yt2));
527 }
528 return(retval);
529 }
530
531
532 /*
533 * NAME
534 * _XcmsArcTangent
535 *
536 * SYNOPSIS
537 */
538 double
_XcmsArcTangent(double x)539 _XcmsArcTangent(double x)
540 /*
541 * DESCRIPTION
542 * Computes the arctangent.
543 * This is an implementation of the Gauss algorithm as
544 * described in:
545 * Forman S. Acton, Numerical Methods That Work,
546 * New York, NY, Harper & Row, 1970.
547 *
548 * RETURNS
549 * Returns the arctangent
550 */
551 {
552 double ai, a1 = 0.0, bi, b1 = 0.0, l, d;
553 double maxerror;
554 int i;
555
556 if (x == 0.0) {
557 return (0.0);
558 }
559 if (x < 1.0) {
560 maxerror = x * XCMS_MAXERROR;
561 } else {
562 maxerror = XCMS_MAXERROR;
563 }
564 ai = _XcmsSquareRoot( 1.0 / (1.0 + (x * x)) );
565 bi = 1.0;
566 for (i = 0; i < XCMS_MAXITER; i++) {
567 a1 = (ai + bi) / 2.0;
568 b1 = _XcmsSquareRoot((a1 * bi));
569 if (a1 == b1)
570 break;
571 d = XCMS_FABS(a1 - b1);
572 if (d < maxerror)
573 break;
574 ai = a1;
575 bi = b1;
576 }
577
578 l = ((a1 > b1) ? b1 : a1);
579
580 a1 = _XcmsSquareRoot(1 + (x * x));
581 return (x / (a1 * l));
582 }
583