jcs.org
Headers and library sources from all versions of Lightspeed C and THINK C
1
2/*
3
4 Math library for LightspeedC
5
6 (C) Copyright 1986 THINK Technologies. All rights reserved.
7
8 For details, refer to Harbison & Steele's "C: A Reference Manual",
9 Chapter 11.
10
11 Only one version of each function is defined by this library.
12 Whether error checking is done is determined by how this file
13 is compiled. Error checking is ENABLED when the following
14 #define is COMMENTED OUT.
15
16*/
17
18/*#define _NOERRORCHECK_*/
19#include "math.h"
20#include "sane.h"
21
22/* useful constants */
23static double Zero = 0.0;
24static double One = 1.0;
25static double Two = 2.0;
26static double MinusOne = -1.0;
27static double MinusTwo = -2.0;
28static double PointFive = 0.5;
29static double PointTwoFive = 0.25;
30static double Pi = PI;
31static double Pi2 = PI2;
32static double Log2Ten = 3.321928094887362348;
33
34static short _Max[] = { 0x7FFE, 0x7FFF, 0xFFFF, 0xFFFF, 0xFFFF };
35static short _MinusMax[] = { 0xFFFE, 0x7FFF, 0xFFFF, 0xFFFF, 0xFFFF };
36#define Max (* (double *) _Max)
37#define MinusMax (* (double *) _MinusMax)
38
39/* seed for pseudo-random number generator */
40static unsigned long seed = 1;
41
42
43/* environment word masks */
44#define ROUND 0x6000
45#define ROUND_UP 0x2000
46#define ROUND_DOWN 0x4000
47#define ERROR 0x0F00
48
49
50/* ---------- error checking ---------- */
51
52#ifdef _ERRORCHECK_
53
54#define DomainCheck(test, result) if (test) { \
55 errno = EDOM; \
56 return(result); \
57 }
58#define RangeCheck(target) if (GetState() & ERROR) { \
59 errno = ERANGE; \
60 target = Max; \
61 }
62
63#else _ERRORCHECK_
64
65#define ClearExceptions()
66#define DomainCheck(test, result)
67#define RangeCheck(target)
68
69#endif _ERRORCHECK_
70
71/* ---------- environment ---------- */
72
73
74/*
75 * GetState - query environment word
76 *
77 */
78
79static
80GetState()
81{
82 int state;
83
84 GetEnvironment(&state);
85 return(state);
86}
87
88
89/*
90 * SetState - define environment word
91 *
92 */
93
94static
95SetState(state)
96{
97 SetEnvironment(&state);
98}
99
100
101/*
102 * ClearExceptions - clear error conditions
103 *
104 * This must be called if RangeCheck is to be used.
105 *
106 */
107
108#ifdef _ERRORCHECK_
109static
110ClearExceptions()
111{
112 int state;
113
114 GetEnvironment(&state);
115 state &= ~ERROR;
116 SetEnvironment(&state);
117}
118#endif _ERRORCHECK_
119
120
121/*
122 * SetRound - define rounding direction
123 *
124 * The previous environment word is returned. The rounding direction can
125 * be restored by passing this value to SetState.
126 *
127 */
128
129static
130SetRound(direction)
131{
132 int state;
133
134 GetEnvironment(&state);
135 SetState((state & ~ROUND) | direction);
136 return(state);
137}
138
139
140/*
141 * xfersign - transfer sign from one floating number to another
142 *
143 */
144
145static
146xfersign(x, yp)
147double x, *yp;
148{
149 asm {
150 movea.l yp,a0
151 bclr #7,(a0)
152 tst.w x
153 bpl.s @1
154 bset #7,(a0)
155@1 }
156}
157
158
159/* ---------- math functions (alphabetically) ---------- */
160
161
162/*
163 * abs - absolute value of an integer
164 *
165 */
166
167int abs(x)
168int x;
169{
170 return(x < 0 ? -x : x);
171}
172
173
174/*
175 * acos - inverse circular cosine
176 *
177 */
178
179double acos(x)
180double x;
181{
182 DomainCheck(x > One || x < MinusOne, Zero);
183 if (x == MinusOne)
184 return(Pi);
185 return(Two * atan(sqrt((One - x) / (One + x))));
186}
187
188
189/*
190 * asin - inverse circular sine
191 *
192 */
193
194double asin(x)
195double x;
196{
197 double y = fabs(x);
198
199 DomainCheck(y > One, Zero);
200 if (y == One) {
201 y = Pi2;
202 xfersign(x, &y);
203 return(y);
204 }
205 if (y > PointFive) {
206 y = One - y;
207 y = Two * y - y * y;
208 }
209 else
210 y = One - y * y;
211 return(atan(x / sqrt(y)));
212}
213
214
215/*
216 * atan - inverse circular tangent
217 *
218 */
219
220double atan(x)
221double x;
222{
223 elems68k(&x, FOATANX);
224 return(x);
225}
226
227
228/*
229 * atan2 - inverse circular tangent (y/x)
230 *
231 */
232
233double atan2(y, x)
234double y, x;
235{
236 double z;
237
238 if (x == 0) {
239 DomainCheck(y == 0, Zero);
240 z = Pi2;
241 }
242 else {
243 z = atan(fabs(y / x));
244 if (x < 0)
245 z = Pi - z;
246 }
247 xfersign(y, &z);
248 return(z);
249}
250
251
252/*
253 * ceil - round up to an integer
254 *
255 */
256
257double ceil(x)
258double x;
259{
260 int state = SetRound(ROUND_UP);
261 fp68k(&x, FORTI);
262 SetState(state);
263 return(x);
264}
265
266
267/*
268 * cos - circular cosine
269 *
270 */
271
272double cos(x)
273double x;
274{
275 elems68k(&x, FOCOSX);
276 return(x);
277}
278
279
280/*
281 * cosh - hyperbolic cosine
282 *
283 */
284
285double cosh(x)
286double x;
287{
288 ClearExceptions();
289 x = PointFive * exp(fabs(x));
290 x += PointTwoFive / x;
291 RangeCheck(x);
292 return(x);
293}
294
295
296/*
297 * exp - exponential function
298 *
299 */
300
301double exp(x)
302double x;
303{
304 ClearExceptions();
305 elems68k(&x, FOEXPX);
306 RangeCheck(x);
307 return(x);
308}
309
310
311/*
312 * fabs - absolute value of a floating number
313 *
314 */
315
316double fabs(x)
317double x;
318{
319 fp68k(&x, FOABS);
320 return(x);
321}
322
323
324/*
325 * floor - round down to an integer
326 *
327 */
328
329double floor(x)
330double x;
331{
332 int state = SetRound(ROUND_DOWN);
333
334 fp68k(&x, FORTI);
335 SetState(state);
336 return(x);
337}
338
339
340/*
341 * fmod - remainder function
342 *
343 * This computes a value z, with the same sign as x, such that for some
344 * integer k, k*y + z == x.
345 *
346 */
347
348double fmod(x, y)
349double x, y;
350{
351 double z = x;
352
353 fp68k(&y, FOABS);
354 fp68k(&y, &z, FOREM);
355
356 if (x > 0 && z < 0)
357 z += y;
358 else if (x < 0 && z > 0)
359 z -= y;
360 return(z);
361}
362
363
364/*
365 * frexp - split floating number into fraction/exponent
366 *
367 * This computes a value z, where 0.5 <= fabs(z) < 1.0, and an integer n such
368 * that z*(2^n) == x.
369 *
370 */
371
372double frexp(x, nptr)
373double x;
374register int *nptr;
375{
376 double y = fabs(x), z = Two;
377
378 if (y == 0) {
379 *nptr = 0;
380 return(Zero);
381 }
382
383 elems68k(&y, FOLOG2X);
384
385 y -= *nptr = y;
386
387 elems68k(&y, &z, FOPWRY);
388
389 if (z >= One) {
390 z *= PointFive;
391 ++*nptr;
392 }
393 else if (z < PointFive) {
394 z += z;
395 --*nptr;
396 }
397 xfersign(x, &z);
398 return(z);
399}
400
401
402/*
403 * labs - absolute value of a long integer
404 *
405 */
406
407long labs(x)
408long x;
409{
410 return(x < 0 ? -x : x);
411}
412
413
414/*
415 * ldexp - combine fraction/exponent into a floating number
416 *
417 */
418
419double ldexp(x, n)
420double x;
421int n;
422{
423 fp68k(&n, &x, FOSCALB);
424 return(x);
425}
426
427
428/*
429 * log - natural logarithm
430 *
431 */
432
433double log(x)
434double x;
435{
436 DomainCheck(x <= 0, MinusMax);
437 elems68k(&x, FOLNX);
438 return(x);
439}
440
441
442/*
443 * log10 - logarithm base 10
444 *
445 */
446
447double log10(x)
448double x;
449{
450 DomainCheck(x <= 0, MinusMax);
451 elems68k(&x, FOLOG2X); /* LOG2 is much faster than LN */
452 return(x / Log2Ten);
453}
454
455
456/*
457 * modf - split a floating number into fraction/integer
458 *
459 */
460
461double modf(x, nptr)
462double x;
463register double *nptr;
464{
465 *nptr = x;
466 fp68k(nptr, FOTTI);
467 return(x - *nptr);
468}
469
470
471/*
472 * pow - power function (exponentiation)
473 *
474 */
475
476double pow(x, y)
477double x, y;
478{
479 int odd = 0;
480 double z;
481
482 ClearExceptions();
483 if (x == 0) {
484 if (y <= 0) {
485#ifdef _ERRORCHECK_
486 errno = EDOM;
487#endif
488 return (MinusMax);
489 }
490 return(Zero);
491 }
492 if (y == 0)
493 return(One);
494 if (x < 0) {
495 if (modf(y, &y)) {
496#ifdef _ERRORCHECK_
497 errno = EDOM;
498#endif
499 return (MinusMax);
500 }
501 x = -x;
502 odd = fmod(y, Two);
503 }
504 elems68k(&y, &x, FOPWRY);
505
506 RangeCheck(x);
507 return(odd ? -x : x);
508}
509
510/*
511 * rand - pseudo-random number generator (ANSI C standard)
512 *
513 */
514
515int rand()
516{
517 seed = seed * 1103515245 + 12345;
518 asm {
519 move.w seed,d0 ; high word of long
520 andi.w #0x7FFF,d0 ; remove high bit
521 }
522}
523
524
525/*
526 * sin - circular sine
527 *
528 */
529
530double sin(x)
531double x;
532{
533 elems68k(&x, FOSINX);
534 return(x);
535}
536
537
538/*
539 * sinh - hyperbolic sine
540 *
541 */
542
543double sinh(x)
544double x;
545{
546 double y = fabs(x);
547
548 ClearExceptions();
549 elems68k(&y, FOEXP1X);
550 y += y / (y + One);
551 y *= PointFive;
552 RangeCheck(y);
553 xfersign(x, &y);
554 return(y);
555}
556
557
558/*
559 * sqrt - square root
560 *
561 */
562
563double sqrt(x)
564double x;
565{
566 DomainCheck(x < 0, Zero);
567 fp68k(&x, FOSQRT);
568 return(x);
569}
570
571
572/*
573 * srand - seed pseudo-random number generator
574 *
575 */
576
577void srand(n)
578unsigned n;
579{
580 seed = n;
581}
582
583
584/*
585 * tan - circular tangent
586 *
587 */
588
589double tan(x)
590double x;
591{
592 ClearExceptions();
593 elems68k(&x, FOTANX);
594 RangeCheck(x);
595 return(x);
596}
597
598
599/*
600 * tanh - hyperbolic tangent
601 *
602 */
603
604double tanh(x)
605double x;
606{
607 double y = MinusTwo * fabs(x);
608 elems68k(&y, FOEXP1X);
609 y = -y / (y + Two);
610 xfersign(x, &y);
611 return(y);
612}