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