cmsTrig.c revision c6e579a2 
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 */
52static double _XcmsModulo(double value, double base);
53static double _XcmsPolynomial(
54    register int order,
55    double const *coeffs,
56    double x);
57static 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#define XCMS_X6_UNDERFLOWS	(4.209340e-52)	/* X**6 almost underflows */
75#define XCMS_X16_UNDERFLOWS	(5.421010e-20)	/* X**16 almost underflows*/
76#define XCMS_CHAR_BIT		8
77#define XCMS_LONG_MAX		0x7FFFFFFF
78#define XCMS_DEXPLEN		11
79#define XCMS_NBITS(type)	(XCMS_CHAR_BIT * (int)sizeof(type))
80#define XCMS_FABS(x)		((x) < 0.0 ? -(x) : (x))
81
82/* XCMS_DMAXPOWTWO - largest power of two exactly representable as a double */
83#define XCMS_DMAXPOWTWO	((double)(XCMS_LONG_MAX) * \
84	    (1L << ((XCMS_NBITS(double)-XCMS_DEXPLEN) - XCMS_NBITS(int) + 1)))
85
86/*
87 *	LOCAL VARIABLES
88 */
89
90static double const cos_pcoeffs[] = {
91    0.12905394659037374438e7,
92   -0.37456703915723204710e6,
93    0.13432300986539084285e5,
94   -0.11231450823340933092e3
95};
96
97static double const cos_qcoeffs[] = {
98    0.12905394659037373590e7,
99    0.23467773107245835052e5,
100    0.20969518196726306286e3,
101    1.0
102};
103
104static double const sin_pcoeffs[] = {
105    0.20664343336995858240e7,
106   -0.18160398797407332550e6,
107    0.35999306949636188317e4,
108   -0.20107483294588615719e2
109};
110
111static double const sin_qcoeffs[] = {
112    0.26310659102647698963e7,
113    0.39270242774649000308e5,
114    0.27811919481083844087e3,
115    1.0
116};
117
118/*
119 *
120 *  FUNCTION
121 *
122 *	_XcmsCosine   double precision cosine
123 *
124 *  KEY WORDS
125 *
126 *	cos
127 *	machine independent routines
128 *	trigonometric functions
129 *	math libraries
130 *
131 *  DESCRIPTION
132 *
133 *	Returns double precision cosine of double precision
134 *	floating point argument.
135 *
136 *  USAGE
137 *
138 *	double _XcmsCosine (x)
139 *	double x;
140 *
141 *  REFERENCES
142 *
143 *	Computer Approximations, J.F. Hart et al, John Wiley & Sons,
144 *	1968, pp. 112-120.
145 *
146 *  RESTRICTIONS
147 *
148 *	The sin and cos routines are interactive in the sense that
149 *	in the process of reducing the argument to the range -PI/4
150 *	to PI/4, each may call the other.  Ultimately one or the
151 *	other uses a polynomial approximation on the reduced
152 *	argument.  The sin approximation has a maximum relative error
153 *	of 10**(-17.59) and the cos approximation has a maximum
154 *	relative error of 10**(-16.18).
155 *
156 *	These error bounds assume exact arithmetic
157 *	in the polynomial evaluation.  Additional rounding and
158 *	truncation errors may occur as the argument is reduced
159 *	to the range over which the polynomial approximation
160 *	is valid, and as the polynomial is evaluated using
161 *	finite-precision arithmetic.
162 *
163 *  PROGRAMMER
164 *
165 *	Fred Fish
166 *
167 *  INTERNALS
168 *
169 *	Computes cos(x) from:
170 *
171 *		(1)	Reduce argument x to range -PI to PI.
172 *
173 *		(2)	If x > PI/2 then call cos recursively
174 *			using relation cos(x) = -cos(x - PI).
175 *
176 *		(3)	If x < -PI/2 then call cos recursively
177 *			using relation cos(x) = -cos(x + PI).
178 *
179 *		(4)	If x > PI/4 then call sin using
180 *			relation cos(x) = sin(PI/2 - x).
181 *
182 *		(5)	If x < -PI/4 then call cos using
183 *			relation cos(x) = sin(PI/2 + x).
184 *
185 *		(6)	If x would cause underflow in approx
186 *			evaluation arithmetic then return
187 *			sqrt(1.0 - x**2).
188 *
189 *		(7)	By now x has been reduced to range
190 *			-PI/4 to PI/4 and the approximation
191 *			from HART pg. 119 can be used:
192 *
193 *			cos(x) = ( p(y) / q(y) )
194 *			Where:
195 *
196 *			y    =	x * (4/PI)
197 *
198 *			p(y) =  SUM [ Pj * (y**(2*j)) ]
199 *			over j = {0,1,2,3}
200 *
201 *			q(y) =  SUM [ Qj * (y**(2*j)) ]
202 *			over j = {0,1,2,3}
203 *
204 *			P0   =	0.12905394659037374438571854e+7
205 *			P1   =	-0.3745670391572320471032359e+6
206 *			P2   =	0.134323009865390842853673e+5
207 *			P3   =	-0.112314508233409330923e+3
208 *			Q0   =	0.12905394659037373590295914e+7
209 *			Q1   =	0.234677731072458350524124e+5
210 *			Q2   =	0.2096951819672630628621e+3
211 *			Q3   =	1.0000...
212 *			(coefficients from HART table #3843 pg 244)
213 *
214 *
215 *	**** NOTE ****    The range reduction relations used in
216 *	this routine depend on the final approximation being valid
217 *	over the negative argument range in addition to the positive
218 *	argument range.  The particular approximation chosen from
219 *	HART satisfies this requirement, although not explicitly
220 *	stated in the text.  This may not be true of other
221 *	approximations given in the reference.
222 *
223 */
224
225double _XcmsCosine(double x)
226{
227    auto double y;
228    auto double yt2;
229    double retval;
230
231    if (x < -XCMS_PI || x > XCMS_PI) {
232	x = _XcmsModulo (x, XCMS_TWOPI);
233        if (x > XCMS_PI) {
234	    x = x - XCMS_TWOPI;
235        } else if (x < -XCMS_PI) {
236	    x = x + XCMS_TWOPI;
237        }
238    }
239    if (x > XCMS_HALFPI) {
240	retval = -(_XcmsCosine (x - XCMS_PI));
241    } else if (x < -XCMS_HALFPI) {
242	retval = -(_XcmsCosine (x + XCMS_PI));
243    } else if (x > XCMS_FOURTHPI) {
244	retval = _XcmsSine (XCMS_HALFPI - x);
245    } else if (x < -XCMS_FOURTHPI) {
246	retval = _XcmsSine (XCMS_HALFPI + x);
247    } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
248	retval = _XcmsSquareRoot (1.0 - (x * x));
249    } else {
250	y = x / XCMS_FOURTHPI;
251	yt2 = y * y;
252	retval = _XcmsPolynomial (3, cos_pcoeffs, yt2) / _XcmsPolynomial (3, cos_qcoeffs, yt2);
253    }
254    return (retval);
255}
256
257
258/*
259 *  FUNCTION
260 *
261 *	_XcmsModulo   double precision modulo
262 *
263 *  KEY WORDS
264 *
265 *	_XcmsModulo
266 *	machine independent routines
267 *	math libraries
268 *
269 *  DESCRIPTION
270 *
271 *	Returns double precision modulo of two double
272 *	precision arguments.
273 *
274 *  USAGE
275 *
276 *	double _XcmsModulo (value, base)
277 *	double value;
278 *	double base;
279 *
280 *  PROGRAMMER
281 *
282 *	Fred Fish
283 *
284 */
285static double _XcmsModulo(double value, double base)
286{
287    auto double intpart;
288
289    value /= base;
290    value = _XcmsModuloF (value, &intpart);
291    value *= base;
292    return(value);
293}
294
295
296/*
297 * frac = (double) _XcmsModuloF(double val, double *dp)
298 *	return fractional part of 'val'
299 *	set *dp to integer part of 'val'
300 *
301 * Note  -> only compiled for the CA or KA.  For the KB/MC,
302 * "math.c" instantiates a copy of the inline function
303 * defined in "math.h".
304 */
305static double
306_XcmsModuloF(
307    double val,
308    register double *dp)
309{
310	register double abs;
311	/*
312	 * Don't use a register for this.  The extra precision this results
313	 * in on some systems causes problems.
314	 */
315	double ip;
316
317	/* should check for illegal values here - nan, inf, etc */
318	abs = XCMS_FABS(val);
319	if (abs >= XCMS_DMAXPOWTWO) {
320		ip = val;
321	} else {
322		ip = abs + XCMS_DMAXPOWTWO;	/* dump fraction */
323		ip -= XCMS_DMAXPOWTWO;	/* restore w/o frac */
324		if (ip > abs)		/* if it rounds up */
325			ip -= 1.0;	/* fix it */
326		ip = XCMS_FABS(ip);
327	}
328	*dp = ip;
329	return (val - ip); /* signed fractional part */
330}
331
332
333/*
334 *  FUNCTION
335 *
336 *	_XcmsPolynomial   double precision polynomial evaluation
337 *
338 *  KEY WORDS
339 *
340 *	poly
341 *	machine independent routines
342 *	math libraries
343 *
344 *  DESCRIPTION
345 *
346 *	Evaluates a polynomial and returns double precision
347 *	result.  Is passed a the order of the polynomial,
348 *	a pointer to an array of double precision polynomial
349 *	coefficients (in ascending order), and the independent
350 *	variable.
351 *
352 *  USAGE
353 *
354 *	double _XcmsPolynomial (order, coeffs, x)
355 *	int order;
356 *	double *coeffs;
357 *	double x;
358 *
359 *  PROGRAMMER
360 *
361 *	Fred Fish
362 *
363 *  INTERNALS
364 *
365 *	Evalates the polynomial using recursion and the form:
366 *
367 *		P(x) = P0 + x(P1 + x(P2 +...x(Pn)))
368 *
369 */
370
371static double _XcmsPolynomial(
372    register int order,
373    double const *coeffs,
374    double x)
375{
376    auto double rtn_value;
377
378    coeffs += order;
379    rtn_value = *coeffs--;
380    while(order-- > 0)
381	rtn_value = *coeffs-- + (x * rtn_value);
382
383    return(rtn_value);
384}
385
386
387/*
388 *  FUNCTION
389 *
390 *	_XcmsSine	double precision sine
391 *
392 *  KEY WORDS
393 *
394 *	sin
395 *	machine independent routines
396 *	trigonometric functions
397 *	math libraries
398 *
399 *  DESCRIPTION
400 *
401 *	Returns double precision sine of double precision
402 *	floating point argument.
403 *
404 *  USAGE
405 *
406 *	double _XcmsSine (x)
407 *	double x;
408 *
409 *  REFERENCES
410 *
411 *	Computer Approximations, J.F. Hart et al, John Wiley & Sons,
412 *	1968, pp. 112-120.
413 *
414 *  RESTRICTIONS
415 *
416 *	The sin and cos routines are interactive in the sense that
417 *	in the process of reducing the argument to the range -PI/4
418 *	to PI/4, each may call the other.  Ultimately one or the
419 *	other uses a polynomial approximation on the reduced
420 *	argument.  The sin approximation has a maximum relative error
421 *	of 10**(-17.59) and the cos approximation has a maximum
422 *	relative error of 10**(-16.18).
423 *
424 *	These error bounds assume exact arithmetic
425 *	in the polynomial evaluation.  Additional rounding and
426 *	truncation errors may occur as the argument is reduced
427 *	to the range over which the polynomial approximation
428 *	is valid, and as the polynomial is evaluated using
429 *	finite-precision arithmetic.
430 *
431 *  PROGRAMMER
432 *
433 *	Fred Fish
434 *
435 *  INTERNALS
436 *
437 *	Computes sin(x) from:
438 *
439 *	  (1)	Reduce argument x to range -PI to PI.
440 *
441 *	  (2)	If x > PI/2 then call sin recursively
442 *		using relation sin(x) = -sin(x - PI).
443 *
444 *	  (3)	If x < -PI/2 then call sin recursively
445 *		using relation sin(x) = -sin(x + PI).
446 *
447 *	  (4)	If x > PI/4 then call cos using
448 *		relation sin(x) = cos(PI/2 - x).
449 *
450 *	  (5)	If x < -PI/4 then call cos using
451 *		relation sin(x) = -cos(PI/2 + x).
452 *
453 *	  (6)	If x is small enough that polynomial
454 *		evaluation would cause underflow
455 *		then return x, since sin(x)
456 *		approaches x as x approaches zero.
457 *
458 *	  (7)	By now x has been reduced to range
459 *		-PI/4 to PI/4 and the approximation
460 *		from HART pg. 118 can be used:
461 *
462 *		sin(x) = y * ( p(y) / q(y) )
463 *		Where:
464 *
465 *		y    =  x * (4/PI)
466 *
467 *		p(y) =  SUM [ Pj * (y**(2*j)) ]
468 *		over j = {0,1,2,3}
469 *
470 *		q(y) =  SUM [ Qj * (y**(2*j)) ]
471 *		over j = {0,1,2,3}
472 *
473 *		P0   =  0.206643433369958582409167054e+7
474 *		P1   =  -0.18160398797407332550219213e+6
475 *		P2   =  0.359993069496361883172836e+4
476 *		P3   =  -0.2010748329458861571949e+2
477 *		Q0   =  0.263106591026476989637710307e+7
478 *		Q1   =  0.3927024277464900030883986e+5
479 *		Q2   =  0.27811919481083844087953e+3
480 *		Q3   =  1.0000...
481 *		(coefficients from HART table #3063 pg 234)
482 *
483 *
484 *	**** NOTE ****	  The range reduction relations used in
485 *	this routine depend on the final approximation being valid
486 *	over the negative argument range in addition to the positive
487 *	argument range.  The particular approximation chosen from
488 *	HART satisfies this requirement, although not explicitly
489 *	stated in the text.  This may not be true of other
490 *	approximations given in the reference.
491 *
492 */
493
494double
495_XcmsSine (double x)
496{
497    double y;
498    double yt2;
499    double retval;
500
501    if (x < -XCMS_PI || x > XCMS_PI) {
502	x = _XcmsModulo (x, XCMS_TWOPI);
503	if (x > XCMS_PI) {
504	    x = x - XCMS_TWOPI;
505	} else if (x < -XCMS_PI) {
506	    x = x + XCMS_TWOPI;
507	}
508    }
509    if (x > XCMS_HALFPI) {
510	retval = -(_XcmsSine (x - XCMS_PI));
511    } else if (x < -XCMS_HALFPI) {
512	retval = -(_XcmsSine (x + XCMS_PI));
513    } else if (x > XCMS_FOURTHPI) {
514	retval = _XcmsCosine (XCMS_HALFPI - x);
515    } else if (x < -XCMS_FOURTHPI) {
516	retval = -(_XcmsCosine (XCMS_HALFPI + x));
517    } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
518	retval = x;
519    } else {
520	y = x / XCMS_FOURTHPI;
521	yt2 = y * y;
522	retval = y * (_XcmsPolynomial (3, sin_pcoeffs, yt2) / _XcmsPolynomial(3, sin_qcoeffs, yt2));
523    }
524    return(retval);
525}
526
527
528/*
529 *	NAME
530 *		_XcmsArcTangent
531 *
532 *	SYNOPSIS
533 */
534double
535_XcmsArcTangent(double x)
536/*
537 *	DESCRIPTION
538 *		Computes the arctangent.
539 *		This is an implementation of the Gauss algorithm as
540 *		described in:
541 *		    Forman S. Acton, Numerical Methods That Work,
542 *			New York, NY, Harper & Row, 1970.
543 *
544 *	RETURNS
545 *		Returns the arctangent
546 */
547{
548    double ai, a1 = 0.0, bi, b1 = 0.0, l, d;
549    double maxerror;
550    int i;
551
552    if (x == 0.0)  {
553	return (0.0);
554    }
555    if (x < 1.0) {
556	maxerror = x * XCMS_MAXERROR;
557    } else {
558	maxerror = XCMS_MAXERROR;
559    }
560    ai = _XcmsSquareRoot( 1.0 / (1.0 + (x * x)) );
561    bi = 1.0;
562    for (i = 0; i < XCMS_MAXITER; i++) {
563	a1 = (ai + bi) / 2.0;
564	b1 = _XcmsSquareRoot((a1 * bi));
565	if (a1 == b1)
566	    break;
567	d = XCMS_FABS(a1 - b1);
568	if (d < maxerror)
569	    break;
570	ai = a1;
571	bi = b1;
572    }
573
574    l = ((a1 > b1) ? b1 : a1);
575
576    a1 = _XcmsSquareRoot(1 + (x * x));
577    return (x / (a1 * l));
578}
579