diff options
author | Kiwamu Ishikura <ishikura.kiwamu@gmail.com> | 2014-12-29 14:09:46 +0900 |
---|---|---|
committer | Kiwamu Ishikura <ishikura.kiwamu@gmail.com> | 2014-12-29 14:09:46 +0900 |
commit | 2a04e4e69e6c7179a6bc92f7ae150baac4553297 (patch) | |
tree | c6aac7bfc6e5e12e35d857f5c2c93884fb62dc52 /packages/base/src | |
parent | cc7d4003ae9ccd7f8729de58f9fbfad481f94c60 (diff) | |
parent | cc38f9a8e6d51b2e96037375d338546cc606640d (diff) |
merge master
Diffstat (limited to 'packages/base/src')
-rw-r--r-- | packages/base/src/C/vector-aux.c | 80 |
1 files changed, 76 insertions, 4 deletions
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c index d01e788..0f1b3c7 100644 --- a/packages/base/src/C/vector-aux.c +++ b/packages/base/src/C/vector-aux.c | |||
@@ -700,7 +700,10 @@ int saveMatrix(char * file, char * format, KDMAT(a)){ | |||
700 | 700 | ||
701 | //////////////////////////////////////////////////////////////////////////////// | 701 | //////////////////////////////////////////////////////////////////////////////// |
702 | 702 | ||
703 | // http://c-faq.com/lib/gaussian.html | 703 | #ifdef OSX |
704 | |||
705 | #pragma message "randomVector is not thread-safe in OSX" | ||
706 | |||
704 | double gaussrand(int *phase, double *pV1, double *pV2, double *pS) | 707 | double gaussrand(int *phase, double *pV1, double *pV2, double *pS) |
705 | { | 708 | { |
706 | double V1=*pV1, V2=*pV2, S=*pS; | 709 | double V1=*pV1, V2=*pV2, S=*pS; |
@@ -730,10 +733,8 @@ double gaussrand(int *phase, double *pV1, double *pV2, double *pS) | |||
730 | int random_vector(unsigned int seed, int code, DVEC(r)) { | 733 | int random_vector(unsigned int seed, int code, DVEC(r)) { |
731 | int phase = 0; | 734 | int phase = 0; |
732 | double V1,V2,S; | 735 | double V1,V2,S; |
736 | |||
733 | int k; | 737 | int k; |
734 | |||
735 | srandom(seed); | ||
736 | |||
737 | switch (code) { | 738 | switch (code) { |
738 | case 0: { // uniform | 739 | case 0: { // uniform |
739 | for (k=0; k<rn; k++) { | 740 | for (k=0; k<rn; k++) { |
@@ -752,6 +753,77 @@ int random_vector(unsigned int seed, int code, DVEC(r)) { | |||
752 | } | 753 | } |
753 | } | 754 | } |
754 | 755 | ||
756 | #else | ||
757 | |||
758 | inline double urandom(struct random_data * buffer) { | ||
759 | int32_t res; | ||
760 | random_r(buffer,&res); | ||
761 | return (double)res/RAND_MAX; | ||
762 | } | ||
763 | |||
764 | |||
765 | // http://c-faq.com/lib/gaussian.html | ||
766 | double gaussrand(struct random_data *buffer, | ||
767 | int *phase, double *pV1, double *pV2, double *pS) | ||
768 | { | ||
769 | double V1=*pV1, V2=*pV2, S=*pS; | ||
770 | double X; | ||
771 | |||
772 | if(*phase == 0) { | ||
773 | do { | ||
774 | double U1 = (double)random() / (double)RAND_MAX; | ||
775 | double U2 = (double)random() / (double)RAND_MAX; | ||
776 | |||
777 | V1 = 2 * U1 - 1; | ||
778 | V2 = 2 * U2 - 1; | ||
779 | S = V1 * V1 + V2 * V2; | ||
780 | } while(S >= 1 || S == 0); | ||
781 | |||
782 | X = V1 * sqrt(-2 * log(S) / S); | ||
783 | } else | ||
784 | X = V2 * sqrt(-2 * log(S) / S); | ||
785 | |||
786 | *phase = 1 - *phase; | ||
787 | *pV1=V1; *pV2=V2; *pS=S; | ||
788 | |||
789 | return X; | ||
790 | |||
791 | } | ||
792 | |||
793 | int random_vector(unsigned int seed, int code, DVEC(r)) { | ||
794 | struct random_data buffer; | ||
795 | char random_state[128]; | ||
796 | memset(&buffer, 0, sizeof(struct random_data)); | ||
797 | memset(random_state, 0, sizeof(random_state)); | ||
798 | |||
799 | initstate_r(seed,random_state,sizeof(random_state),&buffer); | ||
800 | // setstate_r(random_state,&buffer); | ||
801 | // srandom_r(seed,&buffer); | ||
802 | |||
803 | int phase = 0; | ||
804 | double V1,V2,S; | ||
805 | |||
806 | int k; | ||
807 | switch (code) { | ||
808 | case 0: { // uniform | ||
809 | for (k=0; k<rn; k++) { | ||
810 | rp[k] = urandom(&buffer); | ||
811 | } | ||
812 | OK | ||
813 | } | ||
814 | case 1: { // gaussian | ||
815 | for (k=0; k<rn; k++) { | ||
816 | rp[k] = gaussrand(&buffer,&phase,&V1,&V2,&S); | ||
817 | } | ||
818 | OK | ||
819 | } | ||
820 | |||
821 | default: ERROR(BAD_CODE); | ||
822 | } | ||
823 | } | ||
824 | |||
825 | #endif | ||
826 | |||
755 | //////////////////////////////////////////////////////////////////////////////// | 827 | //////////////////////////////////////////////////////////////////////////////// |
756 | 828 | ||
757 | int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) { | 829 | int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) { |