summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--packages/base/hmatrix.cabal3
-rw-r--r--packages/base/src/C/vector-aux.c57
2 files changed, 59 insertions, 1 deletions
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index c007509..b3559a9 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -1,5 +1,5 @@
1Name: hmatrix 1Name: hmatrix
2Version: 0.16.1.1 2Version: 0.16.1.2
3License: BSD3 3License: BSD3
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: Alberto Ruiz
@@ -118,6 +118,7 @@ library
118 if arch(i386) 118 if arch(i386)
119 cc-options: -arch i386 119 cc-options: -arch i386
120 frameworks: Accelerate 120 frameworks: Accelerate
121 cpp-options: -DOSX
121 122
122 if os(freebsd) 123 if os(freebsd)
123 extra-lib-dirs: /usr/local/lib 124 extra-lib-dirs: /usr/local/lib
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index 54cc413..946d069 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -700,6 +700,61 @@ int saveMatrix(char * file, char * format, KDMAT(a)){
700 700
701//////////////////////////////////////////////////////////////////////////////// 701////////////////////////////////////////////////////////////////////////////////
702 702
703#ifdef OSX
704
705#pragma message "randomVector is not thread-safe in OSX"
706
707double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
708{
709 double V1=*pV1, V2=*pV2, S=*pS;
710 double X;
711
712 if(*phase == 0) {
713 do {
714 double U1 = (double)random() / (double)RAND_MAX;
715 double U2 = (double)random() / (double)RAND_MAX;
716
717 V1 = 2 * U1 - 1;
718 V2 = 2 * U2 - 1;
719 S = V1 * V1 + V2 * V2;
720 } while(S >= 1 || S == 0);
721
722 X = V1 * sqrt(-2 * log(S) / S);
723 } else
724 X = V2 * sqrt(-2 * log(S) / S);
725
726 *phase = 1 - *phase;
727 *pV1=V1; *pV2=V2; *pS=S;
728
729 return X;
730
731}
732
733int random_vector(unsigned int seed, int code, DVEC(r)) {
734 int phase = 0;
735 double V1,V2,S;
736
737 int k;
738 switch (code) {
739 case 0: { // uniform
740 for (k=0; k<rn; k++) {
741 rp[k] = (double)random() / (double)RAND_MAX;
742 }
743 OK
744 }
745 case 1: { // gaussian
746 for (k=0; k<rn; k++) {
747 rp[k] = gaussrand(&phase,&V1,&V2,&S);
748 }
749 OK
750 }
751
752 default: ERROR(BAD_CODE);
753 }
754}
755
756#else
757
703inline double urandom(struct random_data * buffer) { 758inline double urandom(struct random_data * buffer) {
704 int32_t res; 759 int32_t res;
705 random_r(buffer,&res); 760 random_r(buffer,&res);
@@ -766,6 +821,8 @@ int random_vector(unsigned int seed, int code, DVEC(r)) {
766 } 821 }
767} 822}
768 823
824#endif
825
769//////////////////////////////////////////////////////////////////////////////// 826////////////////////////////////////////////////////////////////////////////////
770 827
771int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) { 828int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {