SHOGUN  3.2.1
 全部  命名空间 文件 函数 变量 类型定义 枚举 枚举值 友元 宏定义  
Math.h
浏览该文件的文档.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Thoralf Klein
8  * Written (W) 2013 Soumyajit De
9  * Written (W) 2012 Fernando Jose Iglesias Garcia
10  * Written (W) 2011 Siddharth Kherada
11  * Written (W) 2011 Justin Patera
12  * Written (W) 2011 Alesis Novik
13  * Written (W) 2011-2013 Heiko Strathmann
14  * Written (W) 1999-2009 Soeren Sonnenburg
15  * Written (W) 1999-2008 Gunnar Raetsch
16  * Written (W) 2007 Konrad Rieck
17  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
18  */
19 
20 #ifndef __MATHEMATICS_H_
21 #define __MATHEMATICS_H_
22 
23 #include <shogun/base/SGObject.h>
24 #include <shogun/lib/common.h>
25 #include <shogun/io/SGIO.h>
26 #include <shogun/base/Parallel.h>
28 
29 #ifndef _USE_MATH_DEFINES
30 #define _USE_MATH_DEFINES
31 #endif
32 
33 #include <math.h>
34 #include <stdio.h>
35 #include <float.h>
36 #include <sys/types.h>
37 #ifndef _WIN32
38 #include <unistd.h>
39 #endif
40 
41 #ifdef HAVE_PTHREAD
42 #include <pthread.h>
43 #endif
44 
45 #ifdef SUNOS
46 #include <ieeefp.h>
47 #endif
48 
50 #ifdef log2
51 #define cygwin_log2 log2
52 #undef log2
53 #endif
54 
55 #ifndef M_PI
56 #define M_PI 3.14159265358979323846
57 #endif
58 
59 #ifdef _WIN32
60 #ifndef isnan
61 #define isnan _isnan
62 #endif
63 
64 #ifndef isfinite
65 #define isfinite _isfinite
66 #endif
67 #endif //_WIN32
68 
69 /* Size of RNG seed */
70 #define RNG_SEED_SIZE 256
71 
72 /* Maximum stack size */
73 #define RADIX_STACK_SIZE 512
74 
75 /* Stack macros */
76 #define radix_push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
77 #define radix_pop(a, n, i) a = (--sp)->sa, n = sp->sn, i = sp->si
78 
79 #ifndef DOXYGEN_SHOULD_SKIP_THIS
80 
81 template <class T> struct radix_stack_t
82 {
84  T *sa;
86  size_t sn;
88  uint16_t si;
89 };
90 
92 
94 template <class T1, class T2> struct thread_qsort
95 {
97  T1* output;
99  T2* index;
101  uint32_t size;
102 
104  int32_t* qsort_threads;
106  int32_t sort_limit;
108  int32_t num_threads;
109 };
110 #endif // DOXYGEN_SHOULD_SKIP_THIS
111 
112 #define COMPLEX128_ERROR_ONEARG(function) \
113 static inline complex128_t function(complex128_t a) \
114 { \
115  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
116  #function);\
117  return complex128_t(0.0, 0.0); \
118 }
119 
120 #define COMPLEX128_STDMATH(function) \
121 static inline complex128_t function(complex128_t a) \
122 { \
123  return std::function(a); \
124 }
125 
126 namespace shogun
127 {
129  extern CRandom* sg_rand;
130  class CSGObject;
133 class CMath : public CSGObject
134 {
135  public:
139  CMath();
141 
143  virtual ~CMath();
145 
149 
151  //
152  template <class T>
153  static inline T min(T a, T b)
154  {
155  return (a<=b) ? a : b;
156  }
157 
159  template <class T>
160  static inline T max(T a, T b)
161  {
162  return (a>=b) ? a : b;
163  }
164 
166  template <class T>
167  static inline T clamp(T value, T lb, T ub)
168  {
169  if (value<=lb)
170  return lb;
171  else if (value>=ub)
172  return ub;
173  else
174  return value;
175  }
176 
178  template <class T>
179  static inline T abs(T a)
180  {
181  // can't be a>=0?(a):(-a), because compiler complains about
182  // 'comparison always true' when T is unsigned
183  if (a==0)
184  return 0;
185  else if (a>0)
186  return a;
187  else
188  return -a;
189  }
190 
192  static inline float64_t abs(complex128_t a)
193  {
194  float64_t a_real=a.real();
195  float64_t a_imag=a.imag();
196  return (CMath::sqrt(a_real*a_real+a_imag*a_imag));
197  }
199 
202 
209  template <class T>
210  static inline bool fequals_abs(const T& a, const T& b,
211  const float64_t eps)
212  {
213  const T diff = CMath::abs<T>((a-b));
214  return (diff < eps);
215  }
216 
226  template <class T>
227  static inline bool fequals(const T& a, const T& b,
228  const float64_t eps, bool tolerant=false)
229  {
230  const T absA = CMath::abs<T>(a);
231  const T absB = CMath::abs<T>(b);
232  const T diff = CMath::abs<T>((a-b));
233  T comp;
234 
235  // Handle this separately since NAN is unordered
237  return true;
238 
239  // Required for JSON Serialization Tests
240  if (tolerant)
241  return CMath::fequals_abs<T>(a, b, eps);
242 
243  // handles float32_t and float64_t separately
244  if (sizeof(T) == 4)
246 
247  else
249 
250  if (a==b)
251  return true;
252 
253  // both a and b are 0 and relative error is less meaningful
254  else if ( (a==0) || (b==0) || (diff < comp) )
255  return (diff<(eps * comp));
256 
257  // use max(relative error, diff) to handle large eps
258  else
259  {
260  T check = ((diff/(absA + absB)) > diff)?
261  (diff/(absA + absB)):diff;
262  return (check < eps);
263  }
264  }
265 
266  static inline float64_t round(float64_t d)
267  {
268  return ::floor(d+0.5);
269  }
270 
271  static inline float64_t floor(float64_t d)
272  {
273  return ::floor(d);
274  }
275 
276  static inline float64_t ceil(float64_t d)
277  {
278  return ::ceil(d);
279  }
280 
282  template <class T>
283  static inline T sign(T a)
284  {
285  if (a==0)
286  return 0;
287  else return (a<0) ? (-1) : (+1);
288  }
289 
291  template <class T>
292  static inline void swap(T &a,T &b)
293  {
294  /* register for fast swaps */
295  register T c=a;
296  a=b;
297  b=c;
298  }
299 
301  template <class T>
302  static inline T sq(T x)
303  {
304  return x*x;
305  }
306 
308  static inline float32_t sqrt(float32_t x)
309  {
310  return ::sqrtf(x);
311  }
312 
314  static inline float64_t sqrt(float64_t x)
315  {
316  return ::sqrt(x);
317  }
318 
320  static inline floatmax_t sqrt(floatmax_t x)
321  {
322  //fall back to double precision sqrt if sqrtl is not
323  //available
324 #ifdef HAVE_SQRTL
325  return ::sqrtl(x);
326 #else
327  return ::sqrt(x);
328 #endif
329  }
330 
333 
334 
335  static inline float32_t invsqrt(float32_t x)
336  {
337  union float_to_int
338  {
339  float32_t f;
340  int32_t i;
341  };
342 
343  float_to_int tmp;
344  tmp.f=x;
345 
346  float32_t xhalf = 0.5f * x;
347  // store floating-point bits in integer tmp.i
348  // and do initial guess for Newton's method
349  tmp.i = 0x5f3759d5 - (tmp.i >> 1);
350  x = tmp.f; // convert new bits into float
351  x = x*(1.5f - xhalf*x*x); // One round of Newton's method
352  return x;
353  }
354 
356  static inline floatmax_t powl(floatmax_t x, floatmax_t n)
357  {
358  //fall back to double precision pow if powl is not
359  //available
360 #ifdef HAVE_POWL
361  return ::powl((long double) x, (long double) n);
362 #else
363  return ::pow((double) x, (double) n);
364 #endif
365  }
366 
367  static inline int32_t pow(bool x, int32_t n)
368  {
369  return 0;
370  }
371 
372  static inline int32_t pow(int32_t x, int32_t n)
373  {
374  ASSERT(n>=0)
375  int32_t result=1;
376  while (n--)
377  result*=x;
378 
379  return result;
380  }
381 
382  static inline float64_t pow(float64_t x, int32_t n)
383  {
384  if (n>=0)
385  {
386  float64_t result=1;
387  while (n--)
388  result*=x;
389 
390  return result;
391  }
392  else
393  return ::pow((double)x, (double)n);
394  }
395 
396  static inline float64_t pow(float64_t x, float64_t n)
397  {
398  return ::pow((double) x, (double) n);
399  }
400 
402  static inline complex128_t pow(complex128_t x, int32_t n)
403  {
404  return std::pow(x, n);
405  }
406 
408  {
409  return std::pow(x, n);
410  }
411 
413  {
414  return std::pow(x, n);
415  }
416 
418  {
419  return std::pow(x, n);
420  }
421 
422  static inline float64_t exp(float64_t x)
423  {
424  return ::exp((double) x);
425  }
426 
429 
430 
431  static inline float64_t tan(float64_t x)
432  {
433  return ::tan((double) x);
434  }
435 
438 
439 
440  static inline float64_t atan(float64_t x)
441  {
442  return ::atan((double) x);
443  }
444 
447 
448 
449  static inline float64_t atan2(float64_t x, float64_t y)
450  {
451  return ::atan2((double) x, (double) y);
452  }
453 
456 
457 
458  static inline float64_t tanh(float64_t x)
459  {
460  return ::tanh((double) x);
461  }
462 
465 
466  static inline float64_t log10(float64_t v)
467  {
468  return ::log(v)/::log(10.0);
469  }
470 
473 
474  static inline float64_t log2(float64_t v)
475  {
476 #ifdef HAVE_LOG2
477  return ::log2(v);
478 #else
479  return ::log(v)/::log(2.0);
480 #endif //HAVE_LOG2
481  }
482 
483  static inline float64_t log(float64_t v)
484  {
485  return ::log(v);
486  }
487 
490 
491  static inline index_t floor_log(index_t n)
492  {
493  index_t i;
494  for (i = 0; n != 0; i++)
495  n >>= 1;
496 
497  return i;
498  }
499 
500  static inline float64_t sin(float64_t x)
501  {
502  return ::sin(x);
503  }
504 
507 
508  static inline float64_t asin(float64_t x)
509  {
510  return ::asin(x);
511  }
512 
515 
516  static inline float64_t sinh(float64_t x)
517  {
518  return ::sinh(x);
519  }
520 
523 
524  static inline float64_t cos(float64_t x)
525  {
526  return ::cos(x);
527  }
528 
531 
532  static inline float64_t acos(float64_t x)
533  {
534  return ::acos(x);
535  }
536 
539 
540  static inline float64_t cosh(float64_t x)
541  {
542  return ::cosh(x);
543  }
544 
547 
548  static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
549  {
550  ASSERT(len>0 && xy)
551 
552  float64_t area = 0.0;
553 
554  if (!reversed)
555  {
556  for (int i=1; i<len; i++)
557  area += 0.5*(xy[2*i]-xy[2*(i-1)])*(xy[2*i+1]+xy[2*(i-1)+1]);
558  }
559  else
560  {
561  for (int i=1; i<len; i++)
562  area += 0.5*(xy[2*i+1]-xy[2*(i-1)+1])*(xy[2*i]+xy[2*(i-1)]);
563  }
564 
565  return area;
566  }
567 
568  static bool strtof(const char* str, float32_t* float_result);
569  static bool strtod(const char* str, float64_t* double_result);
570  static bool strtold(const char* str, floatmax_t* long_double_result);
571 
572  static inline int64_t factorial(int32_t n)
573  {
574  int64_t res=1;
575  for (int i=2; i<=n; i++)
576  res*=i ;
577  return res ;
578  }
579 
580  static void init_random(uint32_t initseed=0)
581  {
582  if (initseed==0)
584  else
585  seed=initseed;
586 
588  }
589 
590  static inline uint64_t random()
591  {
592  return sg_rand->random_64();
593  }
594 
595  static inline uint64_t random(uint64_t min_value, uint64_t max_value)
596  {
597  return sg_rand->random(min_value, max_value);
598  }
599 
600  static inline int64_t random(int64_t min_value, int64_t max_value)
601  {
602  return sg_rand->random(min_value, max_value);
603  }
604 
605  static inline uint32_t random(uint32_t min_value, uint32_t max_value)
606  {
607  return sg_rand->random(min_value, max_value);
608  }
609 
610  static inline int32_t random(int32_t min_value, int32_t max_value)
611  {
612  return sg_rand->random(min_value, max_value);
613  }
614 
615  static inline float32_t random(float32_t min_value, float32_t max_value)
616  {
617  return sg_rand->random(min_value, max_value);
618  }
619 
620  static inline float64_t random(float64_t min_value, float64_t max_value)
621  {
622  return sg_rand->random(min_value, max_value);
623  }
624 
625  static inline floatmax_t random(floatmax_t min_value, floatmax_t max_value)
626  {
627  return sg_rand->random(min_value, max_value);
628  }
629 
633  static inline float32_t normal_random(float32_t mean, float32_t std_dev)
634  {
635  // sets up variables & makes sure rand_s.range == (0,1)
636  float32_t ret;
637  float32_t rand_u;
638  float32_t rand_v;
639  float32_t rand_s;
640  do
641  {
642  rand_u = CMath::random(-1.0, 1.0);
643  rand_v = CMath::random(-1.0, 1.0);
644  rand_s = rand_u*rand_u + rand_v*rand_v;
645  } while ((rand_s == 0) || (rand_s >= 1));
646 
647  // the meat & potatos, and then the mean & standard deviation shifting...
648  ret = rand_u*sqrt(-2.0*log(rand_s)/rand_s);
649  ret = std_dev*ret + mean;
650  return ret;
651  }
652 
656  static inline float64_t normal_random(float64_t mean, float64_t std_dev)
657  {
658  return sg_rand->normal_distrib(mean, std_dev);
659  }
660 
663  static inline float32_t randn_float()
664  {
665  return normal_random(0.0, 1.0);
666  }
667 
670  static inline float64_t randn_double()
671  {
672  return sg_rand->std_normal_distrib();
673  }
674 
675  template <class T>
676  static int32_t get_num_nonzero(T* vec, int32_t len)
677  {
678  int32_t nnz = 0;
679  for (index_t i=0; i<len; ++i)
680  nnz += vec[i] != 0;
681 
682  return nnz;
683  }
684 
685  static int32_t get_num_nonzero(complex128_t* vec, int32_t len)
686  {
687  int32_t nnz=0;
688  for (index_t i=0; i<len; ++i)
689  nnz+=vec[i]!=0.0;
690  return nnz;
691  }
692 
693  static inline int64_t nchoosek(int32_t n, int32_t k)
694  {
695  int64_t res=1;
696 
697  for (int32_t i=n-k+1; i<=n; i++)
698  res*=i;
699 
700  return res/factorial(k);
701  }
702 
710  static void linspace(float64_t* output, float64_t start, float64_t end, int32_t n = 100);
711 
718  template <class T>
719  static T log_sum_exp(SGVector<T> values)
720  {
721  REQUIRE(values.vector, "Values are empty");
722 
723  /* find minimum element index */
724  index_t min_index=0;
725  T X0=values[0];
726  if (values.vlen>1)
727  {
728  for (index_t i=1; i<values.vlen; ++i)
729  {
730  if (values[i]<X0)
731  {
732  X0=values[i];
733  min_index=i;
734  }
735  }
736  }
737 
738  /* remove element from vector copy and compute log sum exp */
739  SGVector<T> values_without_X0(values.vlen-1);
740  index_t from_idx=0;
741  index_t to_idx=0;
742  for (from_idx=0; from_idx<values.vlen; ++from_idx)
743  {
744  if (from_idx!=min_index)
745  {
746  values_without_X0[to_idx]=exp(values[from_idx]-X0);
747  to_idx++;
748  }
749  }
750 
751  return X0+log(SGVector<T>::sum(values_without_X0)+1);
752  }
753 
760  template <class T>
761  static T log_mean_exp(SGVector<T> values)
762  {
763  return log_sum_exp(values) - log(values.vlen);
764  }
765 
769  static void sort(int32_t *a, int32_t cols, int32_t sort_col=0);
770  static void sort(float64_t *a, int32_t*idx, int32_t N);
771 
774  template <class T>
775  static void qsort(T* output, int32_t size)
776  {
777  if (size<=1)
778  return;
779 
780  if (size==2)
781  {
782  if (output[0] > output [1])
783  CMath::swap(output[0],output[1]);
784  return;
785  }
786  //T split=output[random(0,size-1)];
787  T split=output[size/2];
788 
789  int32_t left=0;
790  int32_t right=size-1;
791 
792  while (left<=right)
793  {
794  while (output[left] < split)
795  left++;
796  while (output[right] > split)
797  right--;
798 
799  if (left<=right)
800  {
801  CMath::swap(output[left],output[right]);
802  left++;
803  right--;
804  }
805  }
806 
807  if (right+1> 1)
808  qsort(output,right+1);
809 
810  if (size-left> 1)
811  qsort(&output[left],size-left);
812  }
813 
816  template <class T>
817  static void insertion_sort(T* output, int32_t size)
818  {
819  for (int32_t i=0; i<size-1; i++)
820  {
821  int32_t j=i-1;
822  T value=output[i];
823  while (j >= 0 && output[j] > value)
824  {
825  output[j+1] = output[j];
826  j--;
827  }
828  output[j+1]=value;
829  }
830  }
831 
833  template <class T>
834  inline static void radix_sort(T* array, int32_t size)
835  {
836  radix_sort_helper(array,size,0);
837  }
838 
839  /*
840  * Inline function to extract the byte at position p (from left)
841  * of an 64 bit integer. The function is somewhat identical to
842  * accessing an array of characters via [].
843  */
844  template <class T>
845  static inline uint8_t byte(T word, uint16_t p)
846  {
847  return (word >> (sizeof(T)-p-1) * 8) & 0xff;
848  }
849 
851  static inline uint8_t byte(complex128_t word, uint16_t p)
852  {
853  SG_SERROR("CMath::byte():: Not supported for complex128_t\n");
854  return uint8_t(0);
855  }
856 
857  template <class T>
858  static void radix_sort_helper(T* array, int32_t size, uint16_t i)
859  {
860  static size_t count[256], nc, cmin;
861  T *ak;
862  uint8_t c=0;
863  radix_stack_t<T> s[RADIX_STACK_SIZE], *sp, *olds, *bigs;
864  T *an, *aj, *pile[256];
865  size_t *cp, cmax;
866 
867  /* Push initial array to stack */
868  sp = s;
869  radix_push(array, size, i);
870 
871  /* Loop until all digits have been sorted */
872  while (sp>s) {
873  radix_pop(array, size, i);
874  an = array + size;
875 
876  /* Make character histogram */
877  if (nc == 0) {
878  cmin = 0xff;
879  for (ak = array; ak < an; ak++) {
880  c = byte(*ak, i);
881  count[c]++;
882  if (count[c] == 1) {
883  /* Determine smallest character */
884  if (c < cmin)
885  cmin = c;
886  nc++;
887  }
888  }
889 
890  /* Sort recursively if stack size too small */
891  if (sp + nc > s + RADIX_STACK_SIZE) {
892  radix_sort_helper(array, size, i);
893  continue;
894  }
895  }
896 
897  /*
898  * Set pile[]; push incompletely sorted bins onto stack.
899  * pile[] = pointers to last out-of-place element in bins.
900  * Before permuting: pile[c-1] + count[c] = pile[c];
901  * during deal: pile[c] counts down to pile[c-1].
902  */
903  olds = bigs = sp;
904  cmax = 2;
905  ak = array;
906  pile[0xff] = an;
907  for (cp = count + cmin; nc > 0; cp++) {
908  /* Find next non-empty pile */
909  while (*cp == 0)
910  cp++;
911  /* Pile with several entries */
912  if (*cp > 1) {
913  /* Determine biggest pile */
914  if (*cp > cmax) {
915  cmax = *cp;
916  bigs = sp;
917  }
918 
919  if (i < sizeof(T)-1)
920  radix_push(ak, *cp, (uint16_t) (i + 1));
921  }
922  pile[cp - count] = ak += *cp;
923  nc--;
924  }
925 
926  /* Play it safe -- biggest bin last. */
927  swap(*olds, *bigs);
928 
929  /*
930  * Permute misplacements home. Already home: everything
931  * before aj, and in pile[c], items from pile[c] on.
932  * Inner loop:
933  * r = next element to put in place;
934  * ak = pile[r[i]] = location to put the next element.
935  * aj = bottom of 1st disordered bin.
936  * Outer loop:
937  * Once the 1st disordered bin is done, ie. aj >= ak,
938  * aj<-aj + count[c] connects the bins in array linked list;
939  * reset count[c].
940  */
941  aj = array;
942  while (aj<an)
943  {
944  T r;
945 
946  for (r = *aj; aj < (ak = --pile[c = byte(r, i)]);)
947  swap(*ak, r);
948 
949  *aj = r;
950  aj += count[c];
951  count[c] = 0;
952  }
953  }
954  }
955 
957  static void radix_sort_helper(complex128_t* array, int32_t size, uint16_t i)
958  {
959  SG_SERROR("CMath::radix_sort_helper():: Not supported for complex128_t\n");
960  }
961 
971  template <class T>
972  static void qsort(T** vector, index_t length)
973  {
974  if (length<=1)
975  return;
976 
977  if (length==2)
978  {
979  if (*vector[0]>*vector[1])
980  swap(vector[0],vector[1]);
981  return;
982  }
983  T* split=vector[length/2];
984 
985  int32_t left=0;
986  int32_t right=length-1;
987 
988  while (left<=right)
989  {
990  while (*vector[left]<*split)
991  ++left;
992  while (*vector[right]>*split)
993  --right;
994 
995  if (left<=right)
996  {
997  swap(vector[left],vector[right]);
998  ++left;
999  --right;
1000  }
1001  }
1002 
1003  if (right+1>1)
1004  qsort(vector,right+1);
1005 
1006  if (length-left>1)
1007  qsort(&vector[left],length-left);
1008  }
1009 
1011  static void qsort(complex128_t** vector, index_t length)
1012  {
1013  SG_SERROR("CMath::qsort():: Not supported for complex128_t\n");
1014  }
1015 
1017  template <class T> static void display_bits(T word, int32_t width=8*sizeof(T))
1018  {
1019  ASSERT(width>=0)
1020  for (int i=0; i<width; i++)
1021  {
1022  T mask = ((T) 1)<<(sizeof(T)*8-1);
1023  while (mask)
1024  {
1025  if (mask & word)
1026  SG_SPRINT("1")
1027  else
1028  SG_SPRINT("0")
1029 
1030  mask>>=1;
1031  }
1032  }
1033  }
1034 
1036  static void display_bits(complex128_t word,
1037  int32_t width=8*sizeof(complex128_t))
1038  {
1039  SG_SERROR("CMath::display_bits():: Not supported for complex128_t\n");
1040  }
1041 
1047  template <class T1,class T2>
1048  static void qsort_index(T1* output, T2* index, uint32_t size);
1049 
1051  template <class T>
1052  static void qsort_index(complex128_t* output, T* index, uint32_t size)
1053  {
1054  SG_SERROR("CMath::qsort_index():: Not supported for complex128_t\n");
1055  }
1056 
1062  template <class T1,class T2>
1063  static void qsort_backward_index(
1064  T1* output, T2* index, int32_t size);
1065 
1067  template <class T>
1069  complex128_t* output, T* index, uint32_t size)
1070  {
1071  SG_SERROR("CMath::qsort_backword_index():: \
1072  Not supported for complex128_t\n");
1073  }
1074 
1082  template <class T1,class T2>
1083  inline static void parallel_qsort_index(T1* output, T2* index, uint32_t size, int32_t n_threads, int32_t limit=262144)
1084  {
1085  int32_t n=0;
1086  thread_qsort<T1,T2> t;
1087  t.output=output;
1088  t.index=index;
1089  t.size=size;
1090  t.qsort_threads=&n;
1091  t.sort_limit=limit;
1092  t.num_threads=n_threads;
1093  parallel_qsort_index<T1,T2>(&t);
1094  }
1095 
1097  template <class T>
1098  inline static void parallel_qsort_index(complex128_t* output, T* index,
1099  uint32_t size, int32_t n_threads, int32_t limit=0)
1100  {
1101  SG_SERROR("CMath::parallel_qsort_index():: Not supported for complex128_t\n");
1102  }
1103 
1104 
1105  template <class T1,class T2>
1106  static void* parallel_qsort_index(void* p);
1107 
1108 
1109  /* finds the smallest element in output and puts that element as the
1110  first element */
1111  template <class T>
1112  static void min(float64_t* output, T* index, int32_t size);
1113 
1115  static void min(float64_t* output, complex128_t* index, int32_t size)
1116  {
1117  SG_SERROR("CMath::min():: Not supported for complex128_t\n");
1118  }
1119 
1120  /* finds the n smallest elements in output and puts these elements as the
1121  first n elements */
1122  template <class T>
1123  static void nmin(
1124  float64_t* output, T* index, int32_t size, int32_t n);
1125 
1127  static void nmin(float64_t* output, complex128_t* index,
1128  int32_t size, int32_t n)
1129  {
1130  SG_SERROR("CMath::nmin():: Not supported for complex128_t\n");
1131  }
1132 
1133 
1134 
1135  /* finds an element in a sorted array via binary search
1136  * returns -1 if not found */
1137  template <class T>
1138  static int32_t binary_search_helper(T* output, int32_t size, T elem)
1139  {
1140  int32_t start=0;
1141  int32_t end=size-1;
1142 
1143  if (size<1)
1144  return 0;
1145 
1146  while (start<end)
1147  {
1148  int32_t middle=(start+end)/2;
1149 
1150  if (output[middle]>elem)
1151  end=middle-1;
1152  else if (output[middle]<elem)
1153  start=middle+1;
1154  else
1155  return middle;
1156  }
1157 
1158  return start;
1159  }
1160 
1162  static int32_t binary_search_helper(complex128_t* output, int32_t size, complex128_t elem)
1163  {
1164  SG_SERROR("CMath::binary_search_helper():: Not supported for complex128_t\n");
1165  return int32_t(0);
1166  }
1167 
1168  /* finds an element in a sorted array via binary search
1169  * * returns -1 if not found */
1170  template <class T>
1171  static inline int32_t binary_search(T* output, int32_t size, T elem)
1172  {
1173  int32_t ind = binary_search_helper(output, size, elem);
1174  if (ind >= 0 && output[ind] == elem)
1175  return ind;
1176  return -1;
1177  }
1178 
1180  static inline int32_t binary_search(complex128_t* output, int32_t size, complex128_t elem)
1181  {
1182  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1183  return int32_t(-1);
1184  }
1185 
1186 
1187  /* Finds an element in a sorted array of pointers via binary search
1188  * Every element is dereferenced once before being compared
1189  *
1190  * @param array array of pointers to search in (assumed being sorted)
1191  * @param length length of array
1192  * @param elem pointer to element to search for
1193  * @return index of elem, -1 if not found */
1194  template<class T>
1195  static inline int32_t binary_search(T** vector, index_t length,
1196  T* elem)
1197  {
1198  int32_t start=0;
1199  int32_t end=length-1;
1200 
1201  if (length<1)
1202  return -1;
1203 
1204  while (start<end)
1205  {
1206  int32_t middle=(start+end)/2;
1207 
1208  if (*vector[middle]>*elem)
1209  end=middle-1;
1210  else if (*vector[middle]<*elem)
1211  start=middle+1;
1212  else
1213  {
1214  start=middle;
1215  break;
1216  }
1217  }
1218 
1219  if (start>=0&&*vector[start]==*elem)
1220  return start;
1221 
1222  return -1;
1223  }
1224 
1226  static inline int32_t binary_search(complex128_t** vector, index_t length, complex128_t* elem)
1227  {
1228  SG_SERROR("CMath::binary_search():: Not supported for complex128_t\n");
1229  return int32_t(-1);
1230  }
1231 
1232  template <class T>
1234  T* output, int32_t size, T elem)
1235  {
1236  int32_t ind = binary_search_helper(output, size, elem);
1237 
1238  if (output[ind]<=elem)
1239  return ind;
1240  if (ind>0 && output[ind-1] <= elem)
1241  return ind-1;
1242  return -1;
1243  }
1244 
1247  int32_t size, complex128_t elem)
1248  {
1249  SG_SERROR("CMath::binary_search_max_lower_equal():: \
1250  Not supported for complex128_t\n");
1251  return int32_t(-1);
1252  }
1253 
1256  static float64_t Align(
1257  char * seq1, char* seq2, int32_t l1, int32_t l2, float64_t gapCost);
1258 
1259 
1261 
1264  {
1265  return c.real();
1266  }
1267 
1270  {
1271  return c.imag();
1272  }
1273 
1275  inline static uint32_t get_seed()
1276  {
1277  return CMath::seed;
1278  }
1279 
1281  inline static uint32_t get_log_range()
1282  {
1283  return CMath::LOGRANGE;
1284  }
1285 
1286 #ifdef USE_LOGCACHE
1287  inline static uint32_t get_log_accuracy()
1289  {
1290  return CMath::LOGACCURACY;
1291  }
1292 #endif
1293 
1295  static int is_finite(double f);
1296 
1298  static int is_infinity(double f);
1299 
1301  static int is_nan(double f);
1302 
1313 #ifdef USE_LOGCACHE
1314  static inline float64_t logarithmic_sum(float64_t p, float64_t q)
1315  {
1316  float64_t diff;
1317 
1318  if (!CMath::is_finite(p))
1319  return q;
1320 
1321  if (!CMath::is_finite(q))
1322  {
1323  SG_SWARNING("INVALID second operand to logsum(%f,%f) expect undefined results\n", p, q)
1324  return NOT_A_NUMBER;
1325  }
1326  diff = p - q;
1327  if (diff > 0)
1328  return diff > LOGRANGE? p : p + logtable[(int)(diff * LOGACCURACY)];
1329  return -diff > LOGRANGE? q : q + logtable[(int)(-diff * LOGACCURACY)];
1330  }
1331 
1333  static void init_log_table();
1334 
1336  static int32_t determine_logrange();
1337 
1339  static int32_t determine_logaccuracy(int32_t range);
1340 #else
1342  float64_t p, float64_t q)
1343  {
1344  float64_t diff;
1345 
1346  if (!CMath::is_finite(p))
1347  return q;
1348  if (!CMath::is_finite(q))
1349  return p;
1350  diff = p - q;
1351  if (diff > 0)
1352  return diff > LOGRANGE? p : p + log(1 + exp(-diff));
1353  return -diff > LOGRANGE? q : q + log(1 + exp(diff));
1354  }
1355 #endif
1356 #ifdef USE_LOGSUMARRAY
1357 
1362  static inline float64_t logarithmic_sum_array(
1363  float64_t *p, int32_t len)
1364  {
1365  if (len<=2)
1366  {
1367  if (len==2)
1368  return logarithmic_sum(p[0],p[1]) ;
1369  if (len==1)
1370  return p[0];
1371  return -INFTY ;
1372  }
1373  else
1374  {
1375  register float64_t *pp=p ;
1376  if (len%2==1) pp++ ;
1377  for (register int32_t j=0; j < len>>1; j++)
1378  pp[j]=logarithmic_sum(pp[j<<1], pp[1+(j<<1)]) ;
1379  }
1380  return logarithmic_sum_array(p,len%2 + (len>>1)) ;
1381  }
1382 #endif //USE_LOGSUMARRAY
1383 
1384 
1386  virtual const char* get_name() const { return "Math"; }
1387  public:
1390  static const float64_t NOT_A_NUMBER;
1393  static const float64_t INFTY;
1394  static const float64_t ALMOST_INFTY;
1395 
1398 
1400  static const float64_t PI;
1401 
1404 
1405  /* largest and smallest possible float64_t */
1408 
1409  /* Floating point Limits, Normalized */
1410  static const float32_t F_MAX_VAL32;
1412  static const float64_t F_MAX_VAL64;
1414 
1415  /* Floating point limits, Denormalized */
1416  static const float32_t F_MIN_VAL32;
1417  static const float64_t F_MIN_VAL64;
1418 
1419  protected:
1421  static int32_t LOGRANGE;
1422 
1424  static uint32_t seed;
1425 
1426 #ifdef USE_LOGCACHE
1427 
1429  static int32_t LOGACCURACY;
1431  static float64_t* logtable;
1433 #endif
1434 };
1435 
1436 //implementations of template functions
1437 template <class T1,class T2>
1439  {
1440  struct thread_qsort<T1,T2>* ps=(thread_qsort<T1,T2>*) p;
1441  T1* output=ps->output;
1442  T2* index=ps->index;
1443  int32_t size=ps->size;
1444  int32_t* qsort_threads=ps->qsort_threads;
1445  int32_t sort_limit=ps->sort_limit;
1446  int32_t num_threads=ps->num_threads;
1447 
1448  if (size<2)
1449  {
1450  return NULL;
1451  }
1452 
1453  if (size==2)
1454  {
1455  if (output[0] > output [1])
1456  {
1457  swap(output[0], output[1]);
1458  swap(index[0], index[1]);
1459  }
1460  return NULL;
1461  }
1462  //T1 split=output[random(0,size-1)];
1463  T1 split=output[size/2];
1464 
1465  int32_t left=0;
1466  int32_t right=size-1;
1467 
1468  while (left<=right)
1469  {
1470  while (output[left] < split)
1471  left++;
1472  while (output[right] > split)
1473  right--;
1474 
1475  if (left<=right)
1476  {
1477  swap(output[left], output[right]);
1478  swap(index[left], index[right]);
1479  left++;
1480  right--;
1481  }
1482  }
1483  bool lthread_start=false;
1484  bool rthread_start=false;
1485  pthread_t lthread;
1486  pthread_t rthread;
1487  struct thread_qsort<T1,T2> t1;
1488  struct thread_qsort<T1,T2> t2;
1489 
1490  if (right+1> 1 && (right+1< sort_limit || *qsort_threads >= num_threads-1))
1491  qsort_index(output,index,right+1);
1492  else if (right+1> 1)
1493  {
1494  (*qsort_threads)++;
1495  lthread_start=true;
1496  t1.output=output;
1497  t1.index=index;
1498  t1.size=right+1;
1499  t1.qsort_threads=qsort_threads;
1500  t1.sort_limit=sort_limit;
1501  t1.num_threads=num_threads;
1502  if (pthread_create(&lthread, NULL, parallel_qsort_index<T1,T2>, &t1) != 0)
1503  {
1504  lthread_start=false;
1505  (*qsort_threads)--;
1506  qsort_index(output,index,right+1);
1507  }
1508  }
1509 
1510 
1511  if (size-left> 1 && (size-left< sort_limit || *qsort_threads >= num_threads-1))
1512  qsort_index(&output[left],&index[left], size-left);
1513  else if (size-left> 1)
1514  {
1515  (*qsort_threads)++;
1516  rthread_start=true;
1517  t2.output=&output[left];
1518  t2.index=&index[left];
1519  t2.size=size-left;
1520  t2.qsort_threads=qsort_threads;
1521  t2.sort_limit=sort_limit;
1522  t2.num_threads=num_threads;
1523  if (pthread_create(&rthread, NULL, parallel_qsort_index<T1,T2>, &t2) != 0)
1524  {
1525  rthread_start=false;
1526  (*qsort_threads)--;
1527  qsort_index(&output[left],&index[left], size-left);
1528  }
1529  }
1530 
1531  if (lthread_start)
1532  {
1533  pthread_join(lthread, NULL);
1534  (*qsort_threads)--;
1535  }
1536 
1537  if (rthread_start)
1538  {
1539  pthread_join(rthread, NULL);
1540  (*qsort_threads)--;
1541  }
1542 
1543  return NULL;
1544  }
1545 
1546  template <class T1,class T2>
1547 void CMath::qsort_index(T1* output, T2* index, uint32_t size)
1548 {
1549  if (size<=1)
1550  return;
1551 
1552  if (size==2)
1553  {
1554  if (output[0] > output [1])
1555  {
1556  swap(output[0],output[1]);
1557  swap(index[0],index[1]);
1558  }
1559  return;
1560  }
1561  //T1 split=output[random(0,size-1)];
1562  T1 split=output[size/2];
1563 
1564  int32_t left=0;
1565  int32_t right=size-1;
1566 
1567  while (left<=right)
1568  {
1569  while (output[left] < split)
1570  left++;
1571  while (output[right] > split)
1572  right--;
1573 
1574  if (left<=right)
1575  {
1576  swap(output[left],output[right]);
1577  swap(index[left],index[right]);
1578  left++;
1579  right--;
1580  }
1581  }
1582 
1583  if (right+1> 1)
1584  qsort_index(output,index,right+1);
1585 
1586  if (size-left> 1)
1587  qsort_index(&output[left],&index[left], size-left);
1588 }
1589 
1590  template <class T1,class T2>
1591 void CMath::qsort_backward_index(T1* output, T2* index, int32_t size)
1592 {
1593  if (size<=1)
1594  return;
1595 
1596  if (size==2)
1597  {
1598  if (output[0] < output [1])
1599  {
1600  swap(output[0],output[1]);
1601  swap(index[0],index[1]);
1602  }
1603  return;
1604  }
1605 
1606  //T1 split=output[random(0,size-1)];
1607  T1 split=output[size/2];
1608 
1609  int32_t left=0;
1610  int32_t right=size-1;
1611 
1612  while (left<=right)
1613  {
1614  while (output[left] > split)
1615  left++;
1616  while (output[right] < split)
1617  right--;
1618 
1619  if (left<=right)
1620  {
1621  swap(output[left],output[right]);
1622  swap(index[left],index[right]);
1623  left++;
1624  right--;
1625  }
1626  }
1627 
1628  if (right+1> 1)
1629  qsort_backward_index(output,index,right+1);
1630 
1631  if (size-left> 1)
1632  qsort_backward_index(&output[left],&index[left], size-left);
1633 }
1634 
1635  template <class T>
1636 void CMath::nmin(float64_t* output, T* index, int32_t size, int32_t n)
1637 {
1638  if (6*n*size<13*size*CMath::log(size))
1639  for (int32_t i=0; i<n; i++)
1640  min(&output[i], &index[i], size-i) ;
1641  else
1642  qsort_index(output, index, size) ;
1643 }
1644 
1645 /* move the smallest entry in the array to the beginning */
1646  template <class T>
1647 void CMath::min(float64_t* output, T* index, int32_t size)
1648 {
1649  if (size<=1)
1650  return;
1651  float64_t min_elem=output[0];
1652  int32_t min_index=0;
1653  for (int32_t i=1; i<size; i++)
1654  {
1655  if (output[i]<min_elem)
1656  {
1657  min_index=i;
1658  min_elem=output[i];
1659  }
1660  }
1661  swap(output[0], output[min_index]);
1662  swap(index[0], index[min_index]);
1663 }
1664 
1665 #define COMPLEX128_ERROR_ONEARG_T(function) \
1666 template <> \
1667 inline complex128_t CMath::function<complex128_t>(complex128_t a) \
1668 { \
1669  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1670  #function);\
1671  return complex128_t(0.0, 0.0); \
1672 }
1673 
1674 #define COMPLEX128_ERROR_TWOARGS_T(function) \
1675 template <> \
1676 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b) \
1677 { \
1678  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1679  #function);\
1680  return complex128_t(0.0, 0.0); \
1681 }
1682 
1683 #define COMPLEX128_ERROR_THREEARGS_T(function) \
1684 template <> \
1685 inline complex128_t CMath::function<complex128_t>(complex128_t a, complex128_t b, complex128_t c) \
1686 { \
1687  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1688  #function);\
1689  return complex128_t(0.0, 0.0); \
1690 }
1691 
1692 #define COMPLEX128_ERROR_SORT_T(function) \
1693 template <> \
1694 inline void CMath::function<complex128_t>(complex128_t* output, int32_t b) \
1695 { \
1696  SG_SERROR("CMath::%s():: Not supported for complex128_t\n",\
1697  #function);\
1698 }
1699 
1702 
1705 
1708 
1710 // COMPLEX128_ERROR_ONEARG_T(sign)
1711 
1714 
1717 
1720 
1721 }
1722 #undef COMPLEX128_ERROR_ONEARG
1723 #undef COMPLEX128_ERROR_ONEARG_T
1724 #undef COMPLEX128_ERROR_TWOARGS_T
1725 #undef COMPLEX128_ERROR_THREEARGS_T
1726 #undef COMPLEX128_STDMATH
1727 #undef COMPLEX128_ERROR_SORT_T
1728 #endif
static float64_t sin(float64_t x)
Definition: Math.h:500
static float64_t normal_random(float64_t mean, float64_t std_dev)
Definition: Math.h:656
static const float32_t F_MAX_VAL32
Definition: Math.h:1410
static const float64_t MACHINE_EPSILON
Definition: Math.h:1403
static bool strtof(const char *str, float32_t *float_result)
Definition: Math.cpp:264
static int32_t binary_search(complex128_t *output, int32_t size, complex128_t elem)
binary_search not implemented for complex128_t
Definition: Math.h:1180
static void parallel_qsort_index(T1 *output, T2 *index, uint32_t size, int32_t n_threads, int32_t limit=262144)
Definition: Math.h:1083
static uint32_t seed
random generator seed
Definition: Math.h:1424
std::complex< float64_t > complex128_t
Definition: common.h:65
uint64_t random(uint64_t min_value, uint64_t max_value)
Definition: Random.h:100
float64_t std_normal_distrib() const
Definition: Random.cpp:238
static int is_finite(double f)
checks whether a float is finite
Definition: Math.cpp:249
static float64_t Align(char *seq1, char *seq2, int32_t l1, int32_t l2, float64_t gapCost)
Definition: Math.cpp:161
static float64_t sqrt(float64_t x)
x^0.5
Definition: Math.h:314
static floatmax_t random(floatmax_t min_value, floatmax_t max_value)
Definition: Math.h:625
static floatmax_t sqrt(floatmax_t x)
x^0.5
Definition: Math.h:320
static void linspace(float64_t *output, float64_t start, float64_t end, int32_t n=100)
Definition: Math.cpp:204
int32_t index_t
Definition: common.h:60
static float64_t ceil(float64_t d)
Definition: Math.h:276
static bool strtod(const char *str, float64_t *double_result)
Definition: Math.cpp:295
static complex128_t pow(complex128_t x, int32_t n)
x^n, x or n being a complex128_t
Definition: Math.h:402
virtual ~CMath()
Destructor - frees logtable.
Definition: Math.cpp:83
static const float64_t INFTY
infinity
Definition: Math.h:1393
static int32_t binary_search_helper(T *output, int32_t size, T elem)
Definition: Math.h:1138
#define COMPLEX128_ERROR_THREEARGS_T(function)
Definition: Math.h:1683
static void nmin(float64_t *output, T *index, int32_t size, int32_t n)
Definition: Math.h:1636
static void qsort_index(T1 *output, T2 *index, uint32_t size)
Definition: Math.h:1547
static void qsort_index(complex128_t *output, T *index, uint32_t size)
qsort_index not implemented for complex128_t
Definition: Math.h:1052
static float64_t log10(float64_t v)
tanh(x), x being a complex128_t
Definition: Math.h:466
#define COMPLEX128_ERROR_TWOARGS_T(function)
Definition: Math.h:1674
#define SG_SWARNING(...)
Definition: SGIO.h:180
static float32_t normal_random(float32_t mean, float32_t std_dev)
Definition: Math.h:633
uint64_t random_64() const
Definition: Random.cpp:135
static float64_t random(float64_t min_value, float64_t max_value)
Definition: Math.h:620
static int32_t binary_search(T *output, int32_t size, T elem)
Definition: Math.h:1171
static T sq(T x)
x^2
Definition: Math.h:302
static uint32_t get_seed()
returns number generator seed
Definition: Math.h:1275
static float32_t randn_float()
Definition: Math.h:663
static const float64_t MIN_REAL_NUMBER
Definition: Math.h:1407
static float64_t randn_double()
Definition: Math.h:670
static int32_t binary_search_max_lower_equal(T *output, int32_t size, T elem)
Definition: Math.h:1233
static float64_t atan(float64_t x)
tan(x), x being a complex128_t
Definition: Math.h:440
#define REQUIRE(x,...)
Definition: SGIO.h:208
static const float64_t F_MAX_VAL64
Definition: Math.h:1412
static float64_t cosh(float64_t x)
acos(x), x being a complex128_t not implemented
Definition: Math.h:540
static uint32_t get_log_range()
returns range of logtable
Definition: Math.h:1281
void split(v_array< ds_node< P > > &point_set, v_array< ds_node< P > > &far_set, int max_scale)
Definition: JLCoverTree.h:147
static float32_t random(float32_t min_value, float32_t max_value)
Definition: Math.h:615
static const float32_t F_MIN_VAL32
Definition: Math.h:1416
CRandom * sg_rand
Definition: init.cpp:31
static bool fequals_abs(const T &a, const T &b, const float64_t eps)
Definition: Math.h:210
#define RADIX_STACK_SIZE
Definition: Math.h:73
static float64_t imag(complex128_t c)
returns imag part of a complex128_t number
Definition: Math.h:1269
static float64_t floor(float64_t d)
Definition: Math.h:271
static uint64_t random()
Definition: Math.h:590
static int32_t get_num_nonzero(complex128_t *vec, int32_t len)
Definition: Math.h:685
static const float32_t F_MIN_NORM_VAL32
Definition: Math.h:1411
static float64_t real(complex128_t c)
returns real part of a complex128_t number
Definition: Math.h:1263
static uint8_t byte(complex128_t word, uint16_t p)
byte not implemented for complex128_t
Definition: Math.h:851
static void qsort(T *output, int32_t size)
Definition: Math.h:775
static int32_t LOGRANGE
range for logtable: log(1+exp(x)) -LOGRANGE <= x <= 0
Definition: Math.h:1421
static const float64_t ALMOST_NEG_INFTY
almost neg (log) infinity
Definition: Math.h:1397
static bool fequals(const T &a, const T &b, const float64_t eps, bool tolerant=false)
Definition: Math.h:227
static complex128_t pow(complex128_t x, complex128_t n)
Definition: Math.h:407
static float32_t invsqrt(float32_t x)
x^0.5, x being a complex128_t
Definition: Math.h:335
static uint8_t byte(T word, uint16_t p)
Definition: Math.h:845
float64_t normal_distrib(float64_t mu, float64_t sigma) const
Definition: Random.cpp:233
static complex128_t pow(float64_t x, complex128_t n)
Definition: Math.h:417
#define SG_SPRINT(...)
Definition: SGIO.h:182
CMath()
Constructor - initializes log-table.
Definition: Math.cpp:66
static float64_t pow(float64_t x, int32_t n)
Definition: Math.h:382
#define ASSERT(x)
Definition: SGIO.h:203
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:102
static int32_t random(int32_t min_value, int32_t max_value)
Definition: Math.h:610
static float64_t pow(float64_t x, float64_t n)
Definition: Math.h:396
static void qsort(T **vector, index_t length)
Definition: Math.h:972
static void init_random(uint32_t initseed=0)
Definition: Math.h:580
static int32_t pow(int32_t x, int32_t n)
Definition: Math.h:372
#define radix_pop(a, n, i)
Definition: Math.h:77
double float64_t
Definition: common.h:48
static void parallel_qsort_index(complex128_t *output, T *index, uint32_t size, int32_t n_threads, int32_t limit=0)
parallel_qsort_index not implemented for complex128_t
Definition: Math.h:1098
static void min(float64_t *output, complex128_t *index, int32_t size)
complex128_t cannot be used as index
Definition: Math.h:1115
long double floatmax_t
Definition: common.h:49
#define COMPLEX128_ERROR_ONEARG(function)
Definition: Math.h:112
static float64_t cos(float64_t x)
sinh(x), x being a complex128_t
Definition: Math.h:524
static void qsort_backword_index(complex128_t *output, T *index, uint32_t size)
qsort_backword_index not implemented for complex128_t
Definition: Math.h:1068
static int32_t get_num_nonzero(T *vec, int32_t len)
Definition: Math.h:676
static T log_sum_exp(SGVector< T > values)
Definition: Math.h:719
static T max(T a, T b)
return the maximum of two integers
Definition: Math.h:160
static float64_t area_under_curve(float64_t *xy, int32_t len, bool reversed)
cosh(x), x being a complex128_t
Definition: Math.h:548
static void display_bits(T word, int32_t width=8 *sizeof(T))
display bits (useful for debugging)
Definition: Math.h:1017
static uint64_t random(uint64_t min_value, uint64_t max_value)
Definition: Math.h:595
static int64_t factorial(int32_t n)
Definition: Math.h:572
static void qsort(complex128_t **vector, index_t length)
qsort not implemented for complex128_t
Definition: Math.h:1011
static float64_t atan2(float64_t x, float64_t y)
atan(x), x being a complex128_t not implemented
Definition: Math.h:449
static float64_t tan(float64_t x)
exp(x), x being a complex128_t
Definition: Math.h:431
static int32_t binary_search_max_lower_equal(complex128_t *output, int32_t size, complex128_t elem)
binary_search_max_lower_equal not implemented for complex128_t
Definition: Math.h:1246
float float32_t
Definition: common.h:47
static float64_t sinh(float64_t x)
asin(x), x being a complex128_t not implemented
Definition: Math.h:516
static int64_t nchoosek(int32_t n, int32_t k)
Definition: Math.h:693
static int32_t binary_search(T **vector, index_t length, T *elem)
Definition: Math.h:1195
static int is_infinity(double f)
checks whether a float is infinity
Definition: Math.cpp:230
void set_seed(uint32_t seed)
Definition: Random.cpp:57
static float64_t acos(float64_t x)
cos(x), x being a complex128_t
Definition: Math.h:532
static float64_t abs(complex128_t a)
return the absolute value of a complex number
Definition: Math.h:192
static float64_t log2(float64_t v)
log10(x), x being a complex128_t
Definition: Math.h:474
static void display_bits(complex128_t word, int32_t width=8 *sizeof(complex128_t))
disply_bits not implemented for complex128_t
Definition: Math.h:1036
#define COMPLEX128_STDMATH(function)
Definition: Math.h:120
static void radix_sort_helper(T *array, int32_t size, uint16_t i)
Definition: Math.h:858
static T sign(T a)
signum of type T variable a
Definition: Math.h:283
static float64_t tanh(float64_t x)
atan2(x), x being a complex128_t not implemented
Definition: Math.h:458
static int is_nan(double f)
checks whether a float is nan
Definition: Math.cpp:217
virtual const char * get_name() const
Definition: Math.h:1386
static float64_t asin(float64_t x)
sin(x), x being a complex128_t
Definition: Math.h:508
static T min(T a, T b)
return the minimum of two integers
Definition: Math.h:153
#define SG_SERROR(...)
Definition: SGIO.h:181
static float64_t exp(float64_t x)
Definition: Math.h:422
static const float64_t F_MIN_VAL64
Definition: Math.h:1417
static const float64_t F_MIN_NORM_VAL64
Definition: Math.h:1413
static float64_t log(float64_t v)
Definition: Math.h:483
static void radix_sort_helper(complex128_t *array, int32_t size, uint16_t i)
radix_sort_helper not implemented for complex128_t
Definition: Math.h:957
static const float64_t ALMOST_INFTY
Definition: Math.h:1394
Class which collects generic mathematical functions.
Definition: Math.h:133
static T log_mean_exp(SGVector< T > values)
Definition: Math.h:761
static void insertion_sort(T *output, int32_t size)
Definition: Math.h:817
static void swap(T &a, T &b)
swap e.g. floats a and b
Definition: Math.h:292
static complex128_t pow(complex128_t x, float64_t n)
Definition: Math.h:412
static void sort(int32_t *a, int32_t cols, int32_t sort_col=0)
Definition: Math.cpp:122
static int32_t binary_search_helper(complex128_t *output, int32_t size, complex128_t elem)
binary_search_helper not implemented for complex128_t
Definition: Math.h:1162
static float64_t round(float64_t d)
Definition: Math.h:266
static float32_t sqrt(float32_t x)
x^0.5
Definition: Math.h:308
static uint32_t generate_seed()
Definition: Random.cpp:315
static index_t floor_log(index_t n)
log(x), x being a complex128_t
Definition: Math.h:491
static void radix_sort(T *array, int32_t size)
Definition: Math.h:834
static uint32_t random(uint32_t min_value, uint32_t max_value)
Definition: Math.h:605
static floatmax_t powl(floatmax_t x, floatmax_t n)
x^n
Definition: Math.h:356
static float64_t logarithmic_sum(float64_t p, float64_t q)
Definition: Math.h:1341
#define radix_push(a, n, i)
Definition: Math.h:76
static T clamp(T value, T lb, T ub)
return the value clamped to interval [lb,ub]
Definition: Math.h:167
static int32_t binary_search(complex128_t **vector, index_t length, complex128_t *elem)
binary_search not implemented for complex128_t
Definition: Math.h:1226
#define COMPLEX128_ERROR_SORT_T(function)
Definition: Math.h:1692
static bool strtold(const char *str, floatmax_t *long_double_result)
Definition: Math.cpp:326
static int32_t pow(bool x, int32_t n)
Definition: Math.h:367
static const float64_t NOT_A_NUMBER
not a number
Definition: Math.h:1391
static void qsort_backward_index(T1 *output, T2 *index, int32_t size)
Definition: Math.h:1591
static const float64_t MAX_REAL_NUMBER
Definition: Math.h:1406
index_t vlen
Definition: SGVector.h:706
static T abs(T a)
return the absolute value of a number
Definition: Math.h:179
static void nmin(float64_t *output, complex128_t *index, int32_t size, int32_t n)
complex128_t cannot be used as index
Definition: Math.h:1127
static int64_t random(int64_t min_value, int64_t max_value)
Definition: Math.h:600
static const float64_t PI
Definition: Math.h:1400

SHOGUN 机器学习工具包 - 项目文档