cmsTrig.c revision b4ee4795 
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#ifdef _CRAY
84#define XCMS_DMAXPOWTWO	((double)(1 < 47))
85#else
86#define XCMS_DMAXPOWTWO	((double)(XCMS_LONG_MAX) * \
87	    (1L << ((XCMS_NBITS(double)-XCMS_DEXPLEN) - XCMS_NBITS(int) + 1)))
88#endif
89
90/*
91 *	LOCAL VARIABLES
92 */
93
94static double const cos_pcoeffs[] = {
95    0.12905394659037374438e7,
96   -0.37456703915723204710e6,
97    0.13432300986539084285e5,
98   -0.11231450823340933092e3
99};
100
101static double const cos_qcoeffs[] = {
102    0.12905394659037373590e7,
103    0.23467773107245835052e5,
104    0.20969518196726306286e3,
105    1.0
106};
107
108static double const sin_pcoeffs[] = {
109    0.20664343336995858240e7,
110   -0.18160398797407332550e6,
111    0.35999306949636188317e4,
112   -0.20107483294588615719e2
113};
114
115static 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
229double _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 */
289static 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 */
309static double
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
375static double _XcmsPolynomial(
376    register int order,
377    double const *coeffs,
378    double x)
379{
380    auto double rtn_value;
381
382#if 0
383    auto double curr_coeff;
384    if (order <= 0) {
385	rtn_value = *coeffs;
386    } else {
387	curr_coeff = *coeffs;	/* Bug in Unisoft's compiler.  Does not */
388	coeffs++;		/* generate good code for *coeffs++ */
389	rtn_value = curr_coeff + x * _XcmsPolynomial (--order, coeffs, x);
390    }
391#else /* ++jrb -- removed tail recursion */
392    coeffs += order;
393    rtn_value = *coeffs--;
394    while(order-- > 0)
395	rtn_value = *coeffs-- + (x * rtn_value);
396#endif
397
398    return(rtn_value);
399}
400
401
402/*
403 *  FUNCTION
404 *
405 *	_XcmsSine	double precision sine
406 *
407 *  KEY WORDS
408 *
409 *	sin
410 *	machine independent routines
411 *	trigonometric functions
412 *	math libraries
413 *
414 *  DESCRIPTION
415 *
416 *	Returns double precision sine of double precision
417 *	floating point argument.
418 *
419 *  USAGE
420 *
421 *	double _XcmsSine (x)
422 *	double x;
423 *
424 *  REFERENCES
425 *
426 *	Computer Approximations, J.F. Hart et al, John Wiley & Sons,
427 *	1968, pp. 112-120.
428 *
429 *  RESTRICTIONS
430 *
431 *	The sin and cos routines are interactive in the sense that
432 *	in the process of reducing the argument to the range -PI/4
433 *	to PI/4, each may call the other.  Ultimately one or the
434 *	other uses a polynomial approximation on the reduced
435 *	argument.  The sin approximation has a maximum relative error
436 *	of 10**(-17.59) and the cos approximation has a maximum
437 *	relative error of 10**(-16.18).
438 *
439 *	These error bounds assume exact arithmetic
440 *	in the polynomial evaluation.  Additional rounding and
441 *	truncation errors may occur as the argument is reduced
442 *	to the range over which the polynomial approximation
443 *	is valid, and as the polynomial is evaluated using
444 *	finite-precision arithmetic.
445 *
446 *  PROGRAMMER
447 *
448 *	Fred Fish
449 *
450 *  INTERNALS
451 *
452 *	Computes sin(x) from:
453 *
454 *	  (1)	Reduce argument x to range -PI to PI.
455 *
456 *	  (2)	If x > PI/2 then call sin recursively
457 *		using relation sin(x) = -sin(x - PI).
458 *
459 *	  (3)	If x < -PI/2 then call sin recursively
460 *		using relation sin(x) = -sin(x + PI).
461 *
462 *	  (4)	If x > PI/4 then call cos using
463 *		relation sin(x) = cos(PI/2 - x).
464 *
465 *	  (5)	If x < -PI/4 then call cos using
466 *		relation sin(x) = -cos(PI/2 + x).
467 *
468 *	  (6)	If x is small enough that polynomial
469 *		evaluation would cause underflow
470 *		then return x, since sin(x)
471 *		approaches x as x approaches zero.
472 *
473 *	  (7)	By now x has been reduced to range
474 *		-PI/4 to PI/4 and the approximation
475 *		from HART pg. 118 can be used:
476 *
477 *		sin(x) = y * ( p(y) / q(y) )
478 *		Where:
479 *
480 *		y    =  x * (4/PI)
481 *
482 *		p(y) =  SUM [ Pj * (y**(2*j)) ]
483 *		over j = {0,1,2,3}
484 *
485 *		q(y) =  SUM [ Qj * (y**(2*j)) ]
486 *		over j = {0,1,2,3}
487 *
488 *		P0   =  0.206643433369958582409167054e+7
489 *		P1   =  -0.18160398797407332550219213e+6
490 *		P2   =  0.359993069496361883172836e+4
491 *		P3   =  -0.2010748329458861571949e+2
492 *		Q0   =  0.263106591026476989637710307e+7
493 *		Q1   =  0.3927024277464900030883986e+5
494 *		Q2   =  0.27811919481083844087953e+3
495 *		Q3   =  1.0000...
496 *		(coefficients from HART table #3063 pg 234)
497 *
498 *
499 *	**** NOTE ****	  The range reduction relations used in
500 *	this routine depend on the final approximation being valid
501 *	over the negative argument range in addition to the positive
502 *	argument range.  The particular approximation chosen from
503 *	HART satisfies this requirement, although not explicitly
504 *	stated in the text.  This may not be true of other
505 *	approximations given in the reference.
506 *
507 */
508
509double
510_XcmsSine (double x)
511{
512    double y;
513    double yt2;
514    double retval;
515
516    if (x < -XCMS_PI || x > XCMS_PI) {
517	x = _XcmsModulo (x, XCMS_TWOPI);
518	if (x > XCMS_PI) {
519	    x = x - XCMS_TWOPI;
520	} else if (x < -XCMS_PI) {
521	    x = x + XCMS_TWOPI;
522	}
523    }
524    if (x > XCMS_HALFPI) {
525	retval = -(_XcmsSine (x - XCMS_PI));
526    } else if (x < -XCMS_HALFPI) {
527	retval = -(_XcmsSine (x + XCMS_PI));
528    } else if (x > XCMS_FOURTHPI) {
529	retval = _XcmsCosine (XCMS_HALFPI - x);
530    } else if (x < -XCMS_FOURTHPI) {
531	retval = -(_XcmsCosine (XCMS_HALFPI + x));
532    } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
533	retval = x;
534    } else {
535	y = x / XCMS_FOURTHPI;
536	yt2 = y * y;
537	retval = y * (_XcmsPolynomial (3, sin_pcoeffs, yt2) / _XcmsPolynomial(3, sin_qcoeffs, yt2));
538    }
539    return(retval);
540}
541
542
543/*
544 *	NAME
545 *		_XcmsArcTangent
546 *
547 *	SYNOPSIS
548 */
549double
550_XcmsArcTangent(double x)
551/*
552 *	DESCRIPTION
553 *		Computes the arctangent.
554 *		This is an implementation of the Gauss algorithm as
555 *		described in:
556 *		    Forman S. Acton, Numerical Methods That Work,
557 *			New York, NY, Harper & Row, 1970.
558 *
559 *	RETURNS
560 *		Returns the arctangent
561 */
562{
563    double ai, a1 = 0.0, bi, b1 = 0.0, l, d;
564    double maxerror;
565    int i;
566
567    if (x == 0.0)  {
568	return (0.0);
569    }
570    if (x < 1.0) {
571	maxerror = x * XCMS_MAXERROR;
572    } else {
573	maxerror = XCMS_MAXERROR;
574    }
575    ai = _XcmsSquareRoot( 1.0 / (1.0 + (x * x)) );
576    bi = 1.0;
577    for (i = 0; i < XCMS_MAXITER; i++) {
578	a1 = (ai + bi) / 2.0;
579	b1 = _XcmsSquareRoot((a1 * bi));
580	if (a1 == b1)
581	    break;
582	d = XCMS_FABS(a1 - b1);
583	if (d < maxerror)
584	    break;
585	ai = a1;
586	bi = b1;
587    }
588
589    l = ((a1 > b1) ? b1 : a1);
590
591    a1 = _XcmsSquareRoot(1 + (x * x));
592    return (x / (a1 * l));
593}
594