1 /*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12 #include <math.h>
13 #include <inttypes.h>
14 #include <errno.h>
15
16 /* NAN builtins for gcc, as they are not part of math.h */
17 #ifndef NANF
18 #define NANF __builtin_nanf ("")
19 #endif
20 #ifndef NANL
21 #define NANL __builtin_nanl ("")
22 #endif
23
24 extern double bsd__ieee754_fmod(double, double);
25 extern float bsd__ieee754_fmodf(float, float);
26 extern double bsd__ieee754_pow(double, double);
27 extern float bsd__ieee754_powf(float, float);
28 extern double bsd__ieee754_remainder(double, double);
29 extern float bsd__ieee754_remainderf(float, float);
30
softmath_fact(double number)31 static inline double softmath_fact(double number)
32 {
33 if (number <= 0)
34 return 1;
35
36 return number * softmath_fact(number - 1);
37 }
38
softmath_expf(float x)39 static inline float softmath_expf(float x)
40 {
41 float result = 0.0;
42 int n;
43
44 for(n = 0; n < 64; n++)
45 {
46 result += (bsd__ieee754_powf(x, n) / softmath_fact(n));
47 if (isnan(result) || isinf(result))
48 break;
49 }
50
51 return result;
52 }
53
softmath_log(double x)54 static inline double softmath_log(double x)
55 {
56 int n, aprox = 8 / (x / 2);
57 double result = 0.0;
58
59 if (x == 0.0)
60 return -HUGE_VALF;
61 else if (x < 0.0001)
62 aprox = 32768;
63
64 if (aprox < 30)
65 aprox = 30;
66
67 for(n = 0; n < aprox; n++)
68 {
69 result += bsd__ieee754_pow((x - 1.0) / (x + 1.0), 2 * n + 1) * (1.0 / (2.0 * n + 1.0));
70 if (isinf(result))
71 break;
72 }
73 result *= 2;
74
75 return result;
76 }
77
softmath_logf(float x)78 static inline float softmath_logf(float x)
79 {
80 int n, aprox = 8 / (x / 2);
81 float result = 0.0;
82
83 if (x == 0.0)
84 return -HUGE_VALF;
85 else if (x < 0.0001)
86 aprox = 32768;
87
88 if (aprox < 30)
89 aprox = 30;
90
91 for(n = 0; n < aprox; n++)
92 {
93 result += bsd__ieee754_powf((x - 1.0) / (x + 1.0), 2 * n + 1) * (1.0 / (2.0 * n + 1.0));
94 if (isinf(result))
95 break;
96 }
97 result *= 2;
98
99 return result;
100 }
101