1#include "commons.h"
  2
  3#include <math.h>
  4#include <stdio.h>
  5#include <stdlib.h>
  6
  7EXPORT void* UNDEFINED_p_ii(int32_t n, int32_t m) { UNDEFINED(); }
  8EXPORT void* UNDEFINED_p_uu(uint32_t n, uint32_t m) { UNDEFINED(); }
  9EXPORT double* UNDEFINED_dp_pi(const void* p, int32_t n) { UNDEFINED(); }
 10EXPORT void* UNDEFINED_vp_pi(const void* p, int32_t n) { UNDEFINED(); }
 11EXPORT void* UNDEFINED_vp_pu(const void* p, uint32_t n) { UNDEFINED(); }
 12EXPORT void UNDEFINED_v_vpdp(const void* p, double* a) { UNDEFINED(); }
 13EXPORT void UNDEFINED_v_vpvp(const void* p, void* a) { UNDEFINED(); }
 14EXPORT double* NOT_IMPLEMENTED_dp_i(int32_t n) { NOT_IMPLEMENTED(); }
 15EXPORT void* NOT_IMPLEMENTED_vp_i(int32_t n) { NOT_IMPLEMENTED(); }
 16EXPORT void* NOT_IMPLEMENTED_vp_u(uint32_t n) { NOT_IMPLEMENTED(); }
 17EXPORT void NOT_IMPLEMENTED_v_dp(double* a) { NOT_IMPLEMENTED(); }
 18EXPORT void NOT_IMPLEMENTED_v_vp(void* p) { NOT_IMPLEMENTED(); }
 19EXPORT void NOT_IMPLEMENTED_v_idpdpdp(int32_t n, double* a, const double* b, const double* c) { NOT_IMPLEMENTED(); }
 20EXPORT void NOT_IMPLEMENTED_v_uvpcvpcvp(uint32_t n, void* r, const void* a, const void* b) { NOT_IMPLEMENTED(); }
 21EXPORT void NOT_IMPLEMENTED_v_uvpvpcvp(uint32_t n, void* a, void* b, const void* o) { NOT_IMPLEMENTED(); }
 22
 23#ifdef _WIN32
 24#define __always_inline inline __attribute((always_inline))
 25#endif
 26
 27void internal_accurate_sincos(double* rcos, double* rsin, double x) {
 28  double _4_x_over_pi = 4 * x / M_PI;
 29  int64_t int_part = ((int64_t)rint(_4_x_over_pi)) & 7;
 30  double frac_part = _4_x_over_pi - (double)(int_part);
 31  double frac_x = M_PI * frac_part / 4.;
 32  // compute the taylor series
 33  double cosp = 1.;
 34  double sinp = 0.;
 35  double powx = 1.;
 36  int64_t nn = 0;
 37  while (fabs(powx) > 1e-20) {
 38    ++nn;
 39    powx = powx * frac_x / (double)(nn);  // x^n/n!
 40    switch (nn & 3) {
 41      case 0:
 42        cosp += powx;
 43        break;
 44      case 1:
 45        sinp += powx;
 46        break;
 47      case 2:
 48        cosp -= powx;
 49        break;
 50      case 3:
 51        sinp -= powx;
 52        break;
 53      default:
 54        abort();  // impossible
 55    }
 56  }
 57  // final multiplication
 58  switch (int_part) {
 59    case 0:
 60      *rcos = cosp;
 61      *rsin = sinp;
 62      break;
 63    case 1:
 64      *rcos = M_SQRT1_2 * (cosp - sinp);
 65      *rsin = M_SQRT1_2 * (cosp + sinp);
 66      break;
 67    case 2:
 68      *rcos = -sinp;
 69      *rsin = cosp;
 70      break;
 71    case 3:
 72      *rcos = -M_SQRT1_2 * (cosp + sinp);
 73      *rsin = M_SQRT1_2 * (cosp - sinp);
 74      break;
 75    case 4:
 76      *rcos = -cosp;
 77      *rsin = -sinp;
 78      break;
 79    case 5:
 80      *rcos = -M_SQRT1_2 * (cosp - sinp);
 81      *rsin = -M_SQRT1_2 * (cosp + sinp);
 82      break;
 83    case 6:
 84      *rcos = sinp;
 85      *rsin = -cosp;
 86      break;
 87    case 7:
 88      *rcos = M_SQRT1_2 * (cosp + sinp);
 89      *rsin = -M_SQRT1_2 * (cosp - sinp);
 90      break;
 91    default:
 92      abort();  // impossible
 93  }
 94  if (fabs(cos(x) - *rcos) > 1e-10 || fabs(sin(x) - *rsin) > 1e-10) {
 95    printf("cos(%.17lf) =? %.17lf instead of %.17lf\n", x, *rcos, cos(x));
 96    printf("sin(%.17lf) =? %.17lf instead of %.17lf\n", x, *rsin, sin(x));
 97    printf("fracx = %.17lf\n", frac_x);
 98    printf("cosp = %.17lf\n", cosp);
 99    printf("sinp = %.17lf\n", sinp);
100    printf("nn = %d\n", (int)(nn));
101  }
102}
103
104double internal_accurate_cos(double x) {
105  double rcos, rsin;
106  internal_accurate_sincos(&rcos, &rsin, x);
107  return rcos;
108}
109double internal_accurate_sin(double x) {
110  double rcos, rsin;
111  internal_accurate_sincos(&rcos, &rsin, x);
112  return rsin;
113}
114
115EXPORT void spqlios_debug_free(void* addr) { free((uint8_t*)addr - 64); }
116
117EXPORT void* spqlios_debug_alloc(uint64_t size) { return (uint8_t*)malloc(size + 64) + 64; }
118
119EXPORT void spqlios_free(void* addr) {
120#ifndef NDEBUG
121  // in debug mode, we deallocated with spqlios_debug_free()
122  spqlios_debug_free(addr);
123#else
124  // in release mode, the function will free aligned memory
125#ifdef _WIN32
126  _aligned_free(addr);
127#else
128  free(addr);
129#endif
130#endif
131}
132
133EXPORT void* spqlios_alloc(uint64_t size) {
134#ifndef NDEBUG
135  // in debug mode, the function will not necessarily have any particular alignment
136  // it will also ensure that memory can only be deallocated with spqlios_free()
137  return spqlios_debug_alloc(size);
138#else
139  // in release mode, the function will return 64-bytes aligned memory
140#ifdef _WIN32
141  void* reps = _aligned_malloc((size + 63) & (UINT64_C(-64)), 64);
142#else
143  void* reps = aligned_alloc(64, (size + 63) & (UINT64_C(-64)));
144#endif
145  if (reps == 0) FATAL_ERROR("Out of memory");
146  return reps;
147#endif
148}
149
150EXPORT void* spqlios_alloc_custom_align(uint64_t align, uint64_t size) {
151#ifndef NDEBUG
152  // in debug mode, the function will not necessarily have any particular alignment
153  // it will also ensure that memory can only be deallocated with spqlios_free()
154  return spqlios_debug_alloc(size);
155#else
156  // in release mode, the function will return aligned memory
157#ifdef _WIN32
158  void* reps = _aligned_malloc(size, align);
159#else
160  void* reps = aligned_alloc(align, size);
161#endif
162  if (reps == 0) FATAL_ERROR("Out of memory");
163  return reps;
164#endif
165}