summaryrefslogtreecommitdiff
path: root/src/xcms/cmsTrig.c
blob: a1e9b50cee655c164814bc2bb7800b2672089b4e (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
/* $Xorg: cmsTrig.c,v 1.3 2000/08/17 19:45:10 cpqbld Exp $ */

/*
 * Code and supporting documentation (c) Copyright 1990 1991 Tektronix, Inc.
 * 	All Rights Reserved
 * 
 * This file is a component of an X Window System-specific implementation
 * of Xcms based on the TekColor Color Management System.  Permission is
 * hereby granted to use, copy, modify, sell, and otherwise distribute this
 * software and its documentation for any purpose and without fee, provided
 * that this copyright, permission, and disclaimer notice is reproduced in
 * all copies of this software and in supporting documentation.  TekColor
 * is a trademark of Tektronix, Inc.
 * 
 * Tektronix makes no representation about the suitability of this software
 * for any purpose.  It is provided "as is" and with all faults.
 * 
 * TEKTRONIX DISCLAIMS ALL WARRANTIES APPLICABLE TO THIS SOFTWARE,
 * INCLUDING THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 * PARTICULAR PURPOSE.  IN NO EVENT SHALL TEKTRONIX BE LIABLE FOR ANY
 * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER
 * RESULTING FROM LOSS OF USE, DATA, OR PROFITS, WHETHER IN AN ACTION OF
 * CONTRACT, NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 * CONNECTION WITH THE USE OR THE PERFORMANCE OF THIS SOFTWARE.
 */

/*
 *	It should be pointed out that for simplicity's sake, the
 *	environment parameters are defined as floating point constants,
 *	rather than octal or hexadecimal initializations of allocated
 *	storage areas.  This means that the range of allowed numbers
 *	may not exactly match the hardware's capabilities.  For example,
 *	if the maximum positive double precision floating point number
 *	is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is
 *	defined to be 1.11E100 then the numbers between 1.11E100 and
 *	1.11...E100 are considered to be undefined.  For most
 *	applications, this will cause no problems.
 *
 *	An alternate method is to allocate a global static "double" variable,
 *	say "maxdouble", and use a union declaration and initialization
 *	to initialize it with the proper bits for the EXACT maximum value.
 *	This was not done because the only compilers available to the
 *	author did not fully support union initialization features.
 *
 */

/*
 *	EXTERNS
 */
extern double _XcmsSquareRoot();

/*
 *	FORWARD DECLARATIONS
 */
double _XcmsCosine();
static double _XcmsModulo();
static double _XcmsModuloF();
static double _XcmsPolynomial();
double _XcmsSine();
double _XcmsArcTangent();

/*
 *	DEFINES
 */

#ifdef __STDC__
#define Const const
#else
#define Const /**/
#endif
#define XCMS_MAXERROR       	0.000001
#define XCMS_MAXITER       	10000
#define XCMS_PI       		3.14159265358979323846264338327950
#define XCMS_TWOPI 		6.28318530717958620
#define XCMS_HALFPI		1.57079632679489660
#define XCMS_FOURTHPI		0.785398163397448280
#define XCMS_SIXTHPI		0.523598775598298820
#define XCMS_RADIANS(d)		((d) * XCMS_PI / 180.0)
#define XCMS_DEGREES(r)		((r) * 180.0 / XCMS_PI)
#define XCMS_X6_UNDERFLOWS	(4.209340e-52)	/* X**6 almost underflows */
#define XCMS_X16_UNDERFLOWS	(5.421010e-20)	/* X**16 almost underflows*/
#define XCMS_CHAR_BIT		8
#define XCMS_LONG_MAX		0x7FFFFFFF
#define XCMS_DEXPLEN		11
#define XCMS_NBITS(type)	(XCMS_CHAR_BIT * (int)sizeof(type))
#define XCMS_FABS(x)		((x) < 0.0 ? -(x) : (x))

/* XCMS_DMAXPOWTWO - largest power of two exactly representable as a double */
#ifdef _CRAY
#define XCMS_DMAXPOWTWO	((double)(1 < 47))
#else
#define XCMS_DMAXPOWTWO	((double)(XCMS_LONG_MAX) * \
	    (1L << (XCMS_NBITS(double)-XCMS_DEXPLEN) - XCMS_NBITS(long) + 1))
#endif

/*
 *	LOCAL VARIABLES
 */

static double Const cos_pcoeffs[] = {
    0.12905394659037374438e7,
   -0.37456703915723204710e6,
    0.13432300986539084285e5,
   -0.11231450823340933092e3
};

static double Const cos_qcoeffs[] = {
    0.12905394659037373590e7,
    0.23467773107245835052e5,
    0.20969518196726306286e3,
    1.0
};

static double Const sin_pcoeffs[] = {
    0.20664343336995858240e7,
   -0.18160398797407332550e6,
    0.35999306949636188317e4,
   -0.20107483294588615719e2
};

static double Const sin_qcoeffs[] = {
    0.26310659102647698963e7,
    0.39270242774649000308e5,
    0.27811919481083844087e3,
    1.0
};

/*
 *
 *  FUNCTION
 *
 *	_XcmsCosine   double precision cosine
 *
 *  KEY WORDS
 *
 *	cos
 *	machine independent routines
 *	trigonometric functions
 *	math libraries
 *
 *  DESCRIPTION
 *
 *	Returns double precision cosine of double precision
 *	floating point argument.
 *
 *  USAGE
 *
 *	double _XcmsCosine (x)
 *	double x;
 *
 *  REFERENCES
 *
 *	Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *	1968, pp. 112-120.
 *
 *  RESTRICTIONS
 *
 *	The sin and cos routines are interactive in the sense that
 *	in the process of reducing the argument to the range -PI/4
 *	to PI/4, each may call the other.  Ultimately one or the
 *	other uses a polynomial approximation on the reduced
 *	argument.  The sin approximation has a maximum relative error
 *	of 10**(-17.59) and the cos approximation has a maximum
 *	relative error of 10**(-16.18).
 *
 *	These error bounds assume exact arithmetic
 *	in the polynomial evaluation.  Additional rounding and
 *	truncation errors may occur as the argument is reduced
 *	to the range over which the polynomial approximation
 *	is valid, and as the polynomial is evaluated using
 *	finite-precision arithmetic.
 *	
 *  PROGRAMMER
 *
 *	Fred Fish
 *
 *  INTERNALS
 *
 *	Computes cos(x) from:
 *
 *		(1)	Reduce argument x to range -PI to PI.
 *
 *		(2)	If x > PI/2 then call cos recursively
 *			using relation cos(x) = -cos(x - PI).
 *
 *		(3)	If x < -PI/2 then call cos recursively
 *			using relation cos(x) = -cos(x + PI).
 *
 *		(4)	If x > PI/4 then call sin using
 *			relation cos(x) = sin(PI/2 - x).
 *
 *		(5)	If x < -PI/4 then call cos using
 *			relation cos(x) = sin(PI/2 + x).
 *
 *		(6)	If x would cause underflow in approx
 *			evaluation arithmetic then return
 *			sqrt(1.0 - x**2).
 *
 *		(7)	By now x has been reduced to range
 *			-PI/4 to PI/4 and the approximation
 *			from HART pg. 119 can be used:
 *
 *			cos(x) = ( p(y) / q(y) )
 *			Where:
 *
 *			y    =	x * (4/PI)
 *
 *			p(y) =  SUM [ Pj * (y**(2*j)) ]
 *			over j = {0,1,2,3}
 *
 *			q(y) =  SUM [ Qj * (y**(2*j)) ]
 *			over j = {0,1,2,3}
 *
 *			P0   =	0.12905394659037374438571854e+7
 *			P1   =	-0.3745670391572320471032359e+6
 *			P2   =	0.134323009865390842853673e+5
 *			P3   =	-0.112314508233409330923e+3
 *			Q0   =	0.12905394659037373590295914e+7
 *			Q1   =	0.234677731072458350524124e+5
 *			Q2   =	0.2096951819672630628621e+3
 *			Q3   =	1.0000...
 *			(coefficients from HART table #3843 pg 244)
 *
 *
 *	**** NOTE ****    The range reduction relations used in
 *	this routine depend on the final approximation being valid
 *	over the negative argument range in addition to the positive
 *	argument range.  The particular approximation chosen from
 *	HART satisfies this requirement, although not explicitly
 *	stated in the text.  This may not be true of other
 *	approximations given in the reference.
 *			
 */

double _XcmsCosine (x)
double x;
{
    auto double y;
    auto double yt2;
    double retval;

    if (x < -XCMS_PI || x > XCMS_PI) {
	x = _XcmsModulo (x, XCMS_TWOPI);
        if (x > XCMS_PI) {
	    x = x - XCMS_TWOPI;
        } else if (x < -XCMS_PI) {
	    x = x + XCMS_TWOPI;
        }
    }
    if (x > XCMS_HALFPI) {
	retval = -(_XcmsCosine (x - XCMS_PI));
    } else if (x < -XCMS_HALFPI) {
	retval = -(_XcmsCosine (x + XCMS_PI));
    } else if (x > XCMS_FOURTHPI) {
	retval = _XcmsSine (XCMS_HALFPI - x);
    } else if (x < -XCMS_FOURTHPI) {
	retval = _XcmsSine (XCMS_HALFPI + x);
    } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
	retval = _XcmsSquareRoot (1.0 - (x * x));
    } else {
	y = x / XCMS_FOURTHPI;
	yt2 = y * y;
	retval = _XcmsPolynomial (3, cos_pcoeffs, yt2) / _XcmsPolynomial (3, cos_qcoeffs, yt2);
    }
    return (retval);
}


/*
 *  FUNCTION
 *
 *	_XcmsModulo   double precision modulo
 *
 *  KEY WORDS
 *
 *	_XcmsModulo
 *	machine independent routines
 *	math libraries
 *
 *  DESCRIPTION
 *
 *	Returns double precision modulo of two double
 *	precision arguments.
 *
 *  USAGE
 *
 *	double _XcmsModulo (value, base)
 *	double value;
 *	double base;
 *
 *  PROGRAMMER
 *
 *	Fred Fish
 *
 */
static double _XcmsModulo (value, base)
double value;
double base;
{
    auto double intpart;

    value /= base;
    value = _XcmsModuloF (value, &intpart);
    value *= base;
    return(value);
}


/*
 * frac = (double) _XcmsModuloF(double val, double *dp)
 *	return fractional part of 'val'
 *	set *dp to integer part of 'val'
 *
 * Note  -> only compiled for the CA or KA.  For the KB/MC,
 * "math.c" instantiates a copy of the inline function
 * defined in "math.h".
 */
static double
_XcmsModuloF(val, dp)
double val;
register double *dp;
{
	register double abs;
	register double ip;

	/* should check for illegal values here - nan, inf, etc */
	abs = XCMS_FABS(val);
	if (abs >= XCMS_DMAXPOWTWO) {
		ip = val;
	} else {
		ip = abs + XCMS_DMAXPOWTWO;	/* dump fraction */
		ip -= XCMS_DMAXPOWTWO;	/* restore w/o frac */
		if (ip > abs)		/* if it rounds up */
			ip -= 1.0;	/* fix it */
		ip = XCMS_FABS(ip);
	}
	*dp = ip;
	return (val - ip); /* signed fractional part */
}


/*
 *  FUNCTION
 *
 *	_XcmsPolynomial   double precision polynomial evaluation
 *
 *  KEY WORDS
 *
 *	poly
 *	machine independent routines
 *	math libraries
 *
 *  DESCRIPTION
 *
 *	Evaluates a polynomial and returns double precision
 *	result.  Is passed a the order of the polynomial,
 *	a pointer to an array of double precision polynomial
 *	coefficients (in ascending order), and the independent
 *	variable.
 *
 *  USAGE
 *
 *	double _XcmsPolynomial (order, coeffs, x)
 *	int order;
 *	double *coeffs;
 *	double x;
 *
 *  PROGRAMMER
 *
 *	Fred Fish
 *
 *  INTERNALS
 *
 *	Evalates the polynomial using recursion and the form:
 *
 *		P(x) = P0 + x(P1 + x(P2 +...x(Pn)))
 *
 */

static double _XcmsPolynomial (order, coeffs, x)
register int order;
double Const *coeffs;
double x;
{
    auto double rtn_value;

#if 0
    auto double curr_coeff;
    if (order <= 0) {
	rtn_value = *coeffs;
    } else {
	curr_coeff = *coeffs;	/* Bug in Unisoft's compiler.  Does not */
	coeffs++;		/* generate good code for *coeffs++ */
	rtn_value = curr_coeff + x * _XcmsPolynomial (--order, coeffs, x);
    }
#else /* ++jrb -- removed tail recursion */
    coeffs += order;
    rtn_value = *coeffs--;
    while(order-- > 0)
	rtn_value = *coeffs-- + (x * rtn_value);
#endif

    return(rtn_value);
}


/*
 *  FUNCTION
 *
 *	_XcmsSine	double precision sine
 *
 *  KEY WORDS
 *
 *	sin
 *	machine independent routines
 *	trigonometric functions
 *	math libraries
 *
 *  DESCRIPTION
 *
 *	Returns double precision sine of double precision
 *	floating point argument.
 *
 *  USAGE
 *
 *	double _XcmsSine (x)
 *	double x;
 *
 *  REFERENCES
 *
 *	Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *	1968, pp. 112-120.
 *
 *  RESTRICTIONS
 *
 *	The sin and cos routines are interactive in the sense that
 *	in the process of reducing the argument to the range -PI/4
 *	to PI/4, each may call the other.  Ultimately one or the
 *	other uses a polynomial approximation on the reduced
 *	argument.  The sin approximation has a maximum relative error
 *	of 10**(-17.59) and the cos approximation has a maximum
 *	relative error of 10**(-16.18).
 *
 *	These error bounds assume exact arithmetic
 *	in the polynomial evaluation.  Additional rounding and
 *	truncation errors may occur as the argument is reduced
 *	to the range over which the polynomial approximation
 *	is valid, and as the polynomial is evaluated using
 *	finite-precision arithmetic.
 *	
 *  PROGRAMMER
 *
 *	Fred Fish
 *
 *  INTERNALS
 *
 *	Computes sin(x) from:
 *
 *	  (1)	Reduce argument x to range -PI to PI.
 *
 *	  (2)	If x > PI/2 then call sin recursively
 *		using relation sin(x) = -sin(x - PI).
 *
 *	  (3)	If x < -PI/2 then call sin recursively
 *		using relation sin(x) = -sin(x + PI).
 *
 *	  (4)	If x > PI/4 then call cos using
 *		relation sin(x) = cos(PI/2 - x).
 *
 *	  (5)	If x < -PI/4 then call cos using
 *		relation sin(x) = -cos(PI/2 + x).
 *
 *	  (6)	If x is small enough that polynomial
 *		evaluation would cause underflow
 *		then return x, since sin(x)
 *		approaches x as x approaches zero.
 *
 *	  (7)	By now x has been reduced to range
 *		-PI/4 to PI/4 and the approximation
 *		from HART pg. 118 can be used:
 *
 *		sin(x) = y * ( p(y) / q(y) )
 *		Where:
 *
 *		y    =  x * (4/PI)
 *
 *		p(y) =  SUM [ Pj * (y**(2*j)) ]
 *		over j = {0,1,2,3}
 *
 *		q(y) =  SUM [ Qj * (y**(2*j)) ]
 *		over j = {0,1,2,3}
 *
 *		P0   =  0.206643433369958582409167054e+7
 *		P1   =  -0.18160398797407332550219213e+6
 *		P2   =  0.359993069496361883172836e+4
 *		P3   =  -0.2010748329458861571949e+2
 *		Q0   =  0.263106591026476989637710307e+7
 *		Q1   =  0.3927024277464900030883986e+5
 *		Q2   =  0.27811919481083844087953e+3
 *		Q3   =  1.0000...
 *		(coefficients from HART table #3063 pg 234)
 *
 *
 *	**** NOTE ****	  The range reduction relations used in
 *	this routine depend on the final approximation being valid
 *	over the negative argument range in addition to the positive
 *	argument range.  The particular approximation chosen from
 *	HART satisfies this requirement, although not explicitly
 *	stated in the text.  This may not be true of other
 *	approximations given in the reference.
 *			
 */

double
_XcmsSine (x)
double x;
{
    double y;
    double yt2;
    double retval;

    if (x < -XCMS_PI || x > XCMS_PI) {
	x = _XcmsModulo (x, XCMS_TWOPI);
	if (x > XCMS_PI) {
	    x = x - XCMS_TWOPI;
	} else if (x < -XCMS_PI) {
	    x = x + XCMS_TWOPI;
	}
    }
    if (x > XCMS_HALFPI) {
	retval = -(_XcmsSine (x - XCMS_PI));
    } else if (x < -XCMS_HALFPI) {
	retval = -(_XcmsSine (x + XCMS_PI));
    } else if (x > XCMS_FOURTHPI) {
	retval = _XcmsCosine (XCMS_HALFPI - x);
    } else if (x < -XCMS_FOURTHPI) {
	retval = -(_XcmsCosine (XCMS_HALFPI + x));
    } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
	retval = x;
    } else {
	y = x / XCMS_FOURTHPI;
	yt2 = y * y;
	retval = y * (_XcmsPolynomial (3, sin_pcoeffs, yt2) / _XcmsPolynomial(3, sin_qcoeffs, yt2));
    }
    return(retval);
}


/*
 *	NAME
 *		_XcmsArcTangent
 *
 *	SYNOPSIS
 */
double
_XcmsArcTangent(x)
    double x;
/*
 *	DESCRIPTION
 *		Computes the arctangent.
 *		This is an implementation of the Gauss algorithm as
 *		described in:
 *		    Forman S. Acton, Numerical Methods That Work,
 *			New York, NY, Harper & Row, 1970.
 *
 *	RETURNS
 *		Returns the arctangent 
 */
{
    double ai, a1, bi, b1, l, d;
    double maxerror;
    int i;

    if (x == 0.0)  {
	return (0.0);
    }
    if (x < 1.0) {
	maxerror = x * XCMS_MAXERROR;
    } else {
	maxerror = XCMS_MAXERROR;
    }
    ai = _XcmsSquareRoot( 1.0 / (1.0 + (x * x)) );
    bi = 1.0;
    for (i = 0; i < XCMS_MAXITER; i++) {
	a1 = (ai + bi) / 2.0;
	b1 = _XcmsSquareRoot((a1 * bi));
	if (a1 == b1)
	    break;
	d = XCMS_FABS(a1 - b1);
	if (d < maxerror)
	    break;
	ai = a1;
	bi = b1;
    }

    l = ((a1 > b1) ? b1 : a1);

    a1 = _XcmsSquareRoot(1 + (x * x));
    return (x / (a1 * l));
}