20 #ifndef __MATHEMATICS_H_
21 #define __MATHEMATICS_H_
29 #ifndef _USE_MATH_DEFINES
30 #define _USE_MATH_DEFINES
36 #include <sys/types.h>
53 #define cygwin_log2 log2
58 #define M_PI 3.14159265358979323846
67 #define isfinite _isfinite
73 #define NAN (strtod("NAN",NULL))
77 #define RNG_SEED_SIZE 256
80 #define RADIX_STACK_SIZE 512
83 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
84 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
86 #ifndef DOXYGEN_SHOULD_SKIP_THIS
88 template <
class T>
struct radix_stack_t
101 template <
class T1,
class T2>
struct thread_qsort
111 int32_t* qsort_threads;
117 #endif // DOXYGEN_SHOULD_SKIP_THIS
119 #define COMPLEX128_ERROR_ONEARG(function) \
120 static inline complex128_t function(complex128_t a) \
122 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
124 return complex128_t(0.0, 0.0); \
127 #define COMPLEX128_STDMATH(function) \
128 static inline complex128_t function(complex128_t a) \
130 return std::function(a); \
160 static inline T
min(T a, T b)
162 return (a<=b) ? a : b;
167 static inline T
max(T a, T b)
169 return (a>=b) ? a : b;
174 static inline T
clamp(T value, T lb, T ub)
231 else return (a<0) ? (-1) : (+1);
236 static inline void swap(T &a,T &b)
246 static inline T
sq(T x)
293 tmp.i = 0x5f3759d5 - (tmp.i >> 1);
295 x = x*(1.5f - xhalf*x*x);
311 static inline int32_t
pow(
bool x, int32_t n)
316 static inline int32_t
pow(int32_t x, int32_t n)
438 for (i = 0; n != 0; i++)
500 for (
int i=1; i<len; i++)
501 area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
505 for (
int i=1; i<len; i++)
506 area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
515 for (
int i=2; i<=n; i++)
525 gettimeofday(&tv, NULL);
526 seed=(uint32_t) (4223517*getpid()*tv.tv_sec*tv.tv_usec);
539 static inline uint64_t
random(uint64_t min_value, uint64_t max_value)
544 static inline int64_t
random(int64_t min_value, int64_t max_value)
549 static inline uint32_t
random(uint32_t min_value, uint32_t max_value)
554 static inline int32_t
random(int32_t min_value, int32_t max_value)
588 rand_s = rand_u*rand_u + rand_v*rand_v;
589 }
while ((rand_s == 0) || (rand_s >= 1));
592 ret = rand_u*
sqrt(-2.0*
log(rand_s)/rand_s);
593 ret = std_dev*ret + mean;
610 rand_s = rand_u*rand_u + rand_v*rand_v;
611 }
while ((rand_s == 0) || (rand_s >= 1));
613 ret = rand_u*
sqrt(-2.0*
log(rand_s)/rand_s);
614 ret = std_dev*ret + mean;
650 static inline int64_t
nchoosek(int32_t n, int32_t k)
654 for (int32_t i=n-k+1; i<=n; i++)
699 for (from_idx=0; from_idx<values.
vlen; ++from_idx)
701 if (from_idx!=min_index)
703 values_without_X0[to_idx]=
exp(values[from_idx]-X0);
726 static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
732 static void qsort(T* output, int32_t size)
739 if (output[0] > output [1])
744 T
split=output[size/2];
747 int32_t right=size-1;
751 while (output[left] < split)
753 while (output[right] > split)
765 qsort(output,right+1);
768 qsort(&output[left],size-left);
776 for (int32_t i=0; i<size-1; i++)
780 while (j >= 0 && output[j] > value)
782 output[j+1] = output[j];
802 static inline uint8_t
byte(T word, uint16_t p)
804 return (word >> (
sizeof(T)-p-1) * 8) & 0xff;
810 SG_SERROR(
"CMath::byte():: Not supported for complex128_t\n");
817 static size_t count[256], nc, cmin;
821 T *an, *aj, *pile[256];
836 for (ak = array; ak < an; ak++) {
864 for (cp = count + cmin; nc > 0; cp++) {
879 pile[cp - count] = ak += *cp;
903 for (r = *aj; aj < (ak = --pile[c =
byte(r, i)]);)
916 SG_SERROR(
"CMath::radix_sort_helper():: Not supported for complex128_t\n");
936 if (*vector[0]>*vector[1])
937 swap(vector[0],vector[1]);
940 T*
split=vector[length/2];
943 int32_t right=length-1;
947 while (*vector[left]<*split)
949 while (*vector[right]>*split)
954 swap(vector[left],vector[right]);
961 qsort(vector,right+1);
964 qsort(&vector[left],length-left);
970 SG_SERROR(
"CMath::qsort():: Not supported for complex128_t\n");
974 template <
class T>
static void display_bits(T word, int32_t width=8*
sizeof(T))
977 for (
int i=0; i<width; i++)
979 T mask = ((T) 1)<<(
sizeof(T)*8-1);
996 SG_SERROR(
"CMath::display_bits():: Not supported for complex128_t\n");
1004 template <
class T1,
class T2>
1005 static void qsort_index(T1* output, T2* index, uint32_t size);
1011 SG_SERROR(
"CMath::qsort_index():: Not supported for complex128_t\n");
1019 template <
class T1,
class T2>
1021 T1* output, T2* index, int32_t size);
1028 SG_SERROR(
"CMath::qsort_backword_index():: \
1029 Not supported for complex128_t\n");
1039 template <
class T1,
class T2>
1040 inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1043 thread_qsort<T1,T2> t;
1049 t.num_threads=n_threads;
1050 parallel_qsort_index<T1,T2>(&t);
1056 uint32_t size, int32_t n_threads, int32_t limit=0)
1058 SG_SERROR(
"CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1062 template <
class T1,
class T2>
1069 static void min(
float64_t* output, T* index, int32_t size);
1074 SG_SERROR(
"CMath::min():: Not supported for complex128_t\n");
1081 float64_t* output, T* index, int32_t size, int32_t n);
1085 int32_t size, int32_t n)
1087 SG_SERROR(
"CMath::nmin():: Not supported for complex128_t\n");
1105 int32_t middle=(start+end)/2;
1107 if (output[middle]>elem)
1109 else if (output[middle]<elem)
1121 SG_SERROR(
"CMath::binary_search_helper():: Not supported for complex128_t\n");
1131 if (ind >= 0 && output[ind] == elem)
1139 SG_SERROR(
"CMath::binary_search():: Not supported for complex128_t\n");
1156 int32_t end=length-1;
1163 int32_t middle=(start+end)/2;
1165 if (*vector[middle]>*elem)
1167 else if (*vector[middle]<*elem)
1176 if (start>=0&&*vector[start]==*elem)
1185 SG_SERROR(
"CMath::binary_search():: Not supported for complex128_t\n");
1191 T* output, int32_t size, T elem)
1195 if (output[ind]<=elem)
1197 if (ind>0 && output[ind-1] <= elem)
1206 SG_SERROR(
"CMath::binary_search_max_lower_equal():: \
1207 Not supported for complex128_t\n");
1214 char * seq1,
char* seq2, int32_t l1, int32_t l2,
float64_t gapCost);
1245 inline static uint32_t get_log_accuracy()
1247 return CMath::LOGACCURACY;
1254 #if defined(isfinite) && !defined(SUNOS)
1265 static int is_nan(
double f);
1287 SG_SWARNING(
"INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1292 return diff >
LOGRANGE? p : p + logtable[(
int)(diff * LOGACCURACY)];
1293 return -diff >
LOGRANGE? q : q + logtable[(
int)(-diff * LOGACCURACY)];
1297 static
void init_log_table();
1300 static int32_t determine_logrange();
1303 static int32_t determine_logaccuracy(int32_t range);
1320 #ifdef USE_LOGSUMARRAY
1326 static inline float64_t logarithmic_sum_array(
1340 if (len%2==1) pp++ ;
1341 for (
register int32_t j=0; j < len>>1; j++)
1344 return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1346 #endif //USE_LOGSUMARRAY
1350 virtual const char*
get_name()
const {
return "Math"; }
1381 static int32_t LOGACCURACY;
1389 template <
class T1,
class T2>
1392 struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1393 T1* output=ps->output;
1394 T2* index=ps->index;
1395 int32_t size=ps->size;
1396 int32_t* qsort_threads=ps->qsort_threads;
1397 int32_t sort_limit=ps->sort_limit;
1398 int32_t num_threads=ps->num_threads;
1407 if (output[0] > output [1])
1409 swap(output[0], output[1]);
1410 swap(index[0], index[1]);
1415 T1
split=output[size/2];
1418 int32_t right=size-1;
1422 while (output[left] < split)
1424 while (output[right] > split)
1429 swap(output[left], output[right]);
1430 swap(index[left], index[right]);
1435 bool lthread_start=
false;
1436 bool rthread_start=
false;
1439 struct thread_qsort<T1,T2> t1;
1440 struct thread_qsort<T1,T2> t2;
1442 if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1444 else if (right+1> 1)
1451 t1.qsort_threads=qsort_threads;
1452 t1.sort_limit=sort_limit;
1453 t1.num_threads=num_threads;
1454 if (pthread_create(<hread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1456 lthread_start=
false;
1463 if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1464 qsort_index(&output[left],&index[left], size-left);
1465 else if (size-left> 1)
1469 t2.output=&output[left];
1470 t2.index=&index[left];
1472 t2.qsort_threads=qsort_threads;
1473 t2.sort_limit=sort_limit;
1474 t2.num_threads=num_threads;
1475 if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1477 rthread_start=
false;
1479 qsort_index(&output[left],&index[left], size-left);
1485 pthread_join(lthread, NULL);
1491 pthread_join(rthread, NULL);
1498 template <
class T1,
class T2>
1506 if (output[0] > output [1])
1508 swap(output[0],output[1]);
1509 swap(index[0],index[1]);
1514 T1
split=output[size/2];
1517 int32_t right=size-1;
1521 while (output[left] < split)
1523 while (output[right] > split)
1528 swap(output[left],output[right]);
1529 swap(index[left],index[right]);
1539 qsort_index(&output[left],&index[left], size-left);
1542 template <
class T1,
class T2>
1550 if (output[0] < output [1])
1552 swap(output[0],output[1]);
1553 swap(index[0],index[1]);
1559 T1
split=output[size/2];
1562 int32_t right=size-1;
1566 while (output[left] > split)
1568 while (output[right] < split)
1573 swap(output[left],output[right]);
1574 swap(index[left],index[right]);
1591 for (int32_t i=0; i<n; i++)
1592 min(&output[i], &index[i], size-i) ;
1604 int32_t min_index=0;
1605 for (int32_t i=1; i<size; i++)
1607 if (output[i]<min_elem)
1613 swap(output[0], output[min_index]);
1614 swap(index[0], index[min_index]);
1617 #define COMPLEX128_ERROR_ONEARG_T(function) \
1619 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
1621 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1623 return complex128_t(0.0, 0.0); \
1626 #define COMPLEX128_ERROR_TWOARGS_T(function) \
1628 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
1630 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1632 return complex128_t(0.0, 0.0); \
1635 #define COMPLEX128_ERROR_THREEARGS_T(function) \
1637 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
1639 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1641 return complex128_t(0.0, 0.0); \
1644 #define COMPLEX128_ERROR_SORT_T(function) \
1646 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
1648 SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1674 #undef COMPLEX128_ERROR_ONEARG
1675 #undef COMPLEX128_ERROR_ONEARG_T
1676 #undef COMPLEX128_ERROR_TWOARGS_T
1677 #undef COMPLEX128_ERROR_THREEARGS_T
1678 #undef COMPLEX128_STDMATH
1679 #undef COMPLEX128_ERROR_SORT_T