c_math.c 4.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236
  1. #include "c_math.h"
  2. #include "c_types.h"
  3. #include "user_config.h"
  4. double floor(double x)
  5. {
  6. return (double) (x < 0.f ? ((int) x == x ? x : (((int) x) - 1)) : ((int) x));
  7. }
  8. #define MAXEXP 2031 /* (MAX_EXP * 16) - 1 */
  9. #define MINEXP -2047 /* (MIN_EXP * 16) - 1 */
  10. #define HUGE MAXFLOAT
  11. double a1[] ICACHE_STORE_ATTR ICACHE_RODATA_ATTR =
  12. {
  13. 1.0,
  14. 0.95760328069857365,
  15. 0.91700404320467123,
  16. 0.87812608018664974,
  17. 0.84089641525371454,
  18. 0.80524516597462716,
  19. 0.77110541270397041,
  20. 0.73841307296974966,
  21. 0.70710678118654752,
  22. 0.67712777346844637,
  23. 0.64841977732550483,
  24. 0.62092890603674203,
  25. 0.59460355750136054,
  26. 0.56939431737834583,
  27. 0.54525386633262883,
  28. 0.52213689121370692,
  29. 0.50000000000000000
  30. };
  31. double a2[] ICACHE_STORE_ATTR ICACHE_RODATA_ATTR =
  32. {
  33. 0.24114209503420288E-17,
  34. 0.92291566937243079E-18,
  35. -0.15241915231122319E-17,
  36. -0.35421849765286817E-17,
  37. -0.31286215245415074E-17,
  38. -0.44654376565694490E-17,
  39. 0.29306999570789681E-17,
  40. 0.11260851040933474E-17
  41. };
  42. double p1 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.833333333333332114e-1;
  43. double p2 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.125000000005037992e-1;
  44. double p3 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.223214212859242590e-2;
  45. double p4 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.434457756721631196e-3;
  46. double q1 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.693147180559945296e0;
  47. double q2 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.240226506959095371e0;
  48. double q3 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.555041086640855953e-1;
  49. double q4 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.961812905951724170e-2;
  50. double q5 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.133335413135857847e-2;
  51. double q6 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.154002904409897646e-3;
  52. double q7 ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.149288526805956082e-4;
  53. double k ICACHE_STORE_ATTR ICACHE_RODATA_ATTR = 0.442695040888963407;
  54. double pow(double x, double y)
  55. {
  56. double frexp(), g, ldexp(), r, u1, u2, v, w, w1, w2, y1, y2, z;
  57. int iw1, m, p;
  58. bool flipsignal = false;
  59. if (y == 0.0)
  60. return (1.0);
  61. if (x <= 0.0)
  62. {
  63. if (x == 0.0)
  64. {
  65. if (y > 0.0)
  66. return 0.0;
  67. //cmemsg(FP_POWO, &y);
  68. //return(HUGE);
  69. }
  70. else
  71. {
  72. //cmemsg(FP_POWN, &x);
  73. x = -x;
  74. if (y != (int) y) { // if y is fractional, then this woud result in a complex number
  75. return NAN;
  76. }
  77. flipsignal = ((int) y) & 1;
  78. }
  79. }
  80. g = frexp(x, &m);
  81. p = 0;
  82. if (g <= a1[8])
  83. p = 8;
  84. if (g <= a1[p + 4])
  85. p += 4;
  86. if (g <= a1[p + 2])
  87. p += 2;
  88. p++;
  89. z = ((g - a1[p]) - a2[p / 2]) / (g + a1[p]);
  90. z += z;
  91. v = z * z;
  92. r = (((p4 * v + p3) * v + p2) * v + p1) * v * z;
  93. r += k * r;
  94. u2 = (r + z * k) + z;
  95. u1 = 0.0625 * (double)(16 * m - p);
  96. y1 = 0.0625 * (double)((int)(16.0 * y));
  97. y2 = y - y1;
  98. w = u2 * y + u1 * y2;
  99. w1 = 0.0625 * (double)((int)(16.0 * w));
  100. w2 = w - w1;
  101. w = w1 + u1 * y1;
  102. w1 = 0.0625 * (double)((int)(16.0 * w));
  103. w2 += (w - w1);
  104. w = 0.0625 * (double)((int)(16.0 * w2));
  105. iw1 = 16.0 * (w1 + w);
  106. w2 -= w;
  107. while (w2 > 0.0)
  108. {
  109. iw1++;
  110. w2 -= 0.0625;
  111. }
  112. if (iw1 > MAXEXP)
  113. {
  114. //cmemsg(FP_POWO, &y);
  115. return (HUGE);
  116. }
  117. if (iw1 < MINEXP)
  118. {
  119. //cmemsg(FP_POWU, &y);
  120. return (0.0);
  121. }
  122. m = iw1 / 16;
  123. if (iw1 >= 0)
  124. m++;
  125. p = 16 * m - iw1;
  126. z = ((((((q7 * w2 + q6) * w2 + q5) * w2 + q4) * w2 + q3) * w2 + q2) * w2 + q1) * w2;
  127. z = a1[p] + a1[p] * z;
  128. double res = ldexp(z, m);
  129. return flipsignal ? -res : res;
  130. }
  131. #if 0
  132. #ifndef __math_68881
  133. double atan(double x)
  134. {
  135. return x;
  136. }
  137. double cos(double x)
  138. {
  139. return x;
  140. }
  141. double sin(double x)
  142. {
  143. return x;
  144. }
  145. double tan(double x)
  146. {
  147. return x;
  148. }
  149. double tanh(double x)
  150. {
  151. return x;
  152. }
  153. double frexp(double x, int *y)
  154. {
  155. return x;
  156. }
  157. double modf(double x, double *y)
  158. {
  159. return x;
  160. }
  161. double ceil(double x)
  162. {
  163. return x;
  164. }
  165. double fabs(double x)
  166. {
  167. return x;
  168. }
  169. double floor(double x)
  170. {
  171. return x;
  172. }
  173. #endif /* ! defined (__math_68881) */
  174. /* Non reentrant ANSI C functions. */
  175. #ifndef _REENT_ONLY
  176. #ifndef __math_68881
  177. double acos(double x)
  178. {
  179. return x;
  180. }
  181. double asin(double x)
  182. {
  183. return x;
  184. }
  185. double atan2(double x, double y)
  186. {
  187. return x;
  188. }
  189. double cosh(double x)
  190. {
  191. return x;
  192. }
  193. double sinh(double x)
  194. {
  195. return x;
  196. }
  197. double exp(double x)
  198. {
  199. return x;
  200. }
  201. double ldexp(double x, int y)
  202. {
  203. return x;
  204. }
  205. double log(double x)
  206. {
  207. return x;
  208. }
  209. double log10(double x)
  210. {
  211. return x;
  212. }
  213. double pow(double x, double y)
  214. {
  215. return x;
  216. }
  217. double sqrt(double x)
  218. {
  219. return x;
  220. }
  221. double fmod(double x, double y)
  222. {
  223. return x;
  224. }
  225. #endif /* ! defined (__math_68881) */
  226. #endif /* ! defined (_REENT_ONLY) */
  227. #endif