diff options
Diffstat (limited to 'packages')
-rw-r--r-- | packages/base/src/C/vector-aux.c | 21 |
1 files changed, 6 insertions, 15 deletions
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c index f8feb37..51bff73 100644 --- a/packages/base/src/C/vector-aux.c +++ b/packages/base/src/C/vector-aux.c | |||
@@ -700,24 +700,16 @@ int saveMatrix(char * file, char * format, KDMAT(a)){ | |||
700 | 700 | ||
701 | //////////////////////////////////////////////////////////////////////////////// | 701 | //////////////////////////////////////////////////////////////////////////////// |
702 | 702 | ||
703 | inline double urandom(struct drand48_data * buffer) { | ||
704 | double res; | ||
705 | drand48_r(buffer,&res); | ||
706 | return res; | ||
707 | } | ||
708 | |||
709 | |||
710 | // http://c-faq.com/lib/gaussian.html | 703 | // http://c-faq.com/lib/gaussian.html |
711 | double gaussrand(struct drand48_data *buffer, | 704 | double gaussrand(int *phase, double *pV1, double *pV2, double *pS) |
712 | int *phase, double *pV1, double *pV2, double *pS) | ||
713 | { | 705 | { |
714 | double V1=*pV1, V2=*pV2, S=*pS; | 706 | double V1=*pV1, V2=*pV2, S=*pS; |
715 | double X; | 707 | double X; |
716 | 708 | ||
717 | if(*phase == 0) { | 709 | if(*phase == 0) { |
718 | do { | 710 | do { |
719 | double U1 = urandom(buffer); | 711 | double U1 = (double)random() / (double)RAND_MAX; |
720 | double U2 = urandom(buffer); | 712 | double U2 = (double)random() / (double)RAND_MAX; |
721 | 713 | ||
722 | V1 = 2 * U1 - 1; | 714 | V1 = 2 * U1 - 1; |
723 | V2 = 2 * U2 - 1; | 715 | V2 = 2 * U2 - 1; |
@@ -732,11 +724,10 @@ double gaussrand(struct drand48_data *buffer, | |||
732 | *pV1=V1; *pV2=V2; *pS=S; | 724 | *pV1=V1; *pV2=V2; *pS=S; |
733 | 725 | ||
734 | return X; | 726 | return X; |
727 | |||
735 | } | 728 | } |
736 | 729 | ||
737 | int random_vector(unsigned int seed, int code, DVEC(r)) { | 730 | int random_vector(unsigned int seed, int code, DVEC(r)) { |
738 | struct drand48_data buffer; | ||
739 | srand48_r(seed,&buffer); | ||
740 | int phase = 0; | 731 | int phase = 0; |
741 | double V1,V2,S; | 732 | double V1,V2,S; |
742 | 733 | ||
@@ -744,13 +735,13 @@ int random_vector(unsigned int seed, int code, DVEC(r)) { | |||
744 | switch (code) { | 735 | switch (code) { |
745 | case 0: { // uniform | 736 | case 0: { // uniform |
746 | for (k=0; k<rn; k++) { | 737 | for (k=0; k<rn; k++) { |
747 | rp[k] = urandom(&buffer); | 738 | rp[k] = (double)random() / (double)RAND_MAX; |
748 | } | 739 | } |
749 | OK | 740 | OK |
750 | } | 741 | } |
751 | case 1: { // gaussian | 742 | case 1: { // gaussian |
752 | for (k=0; k<rn; k++) { | 743 | for (k=0; k<rn; k++) { |
753 | rp[k] = gaussrand(&buffer,&phase,&V1,&V2,&S); | 744 | rp[k] = gaussrand(&phase,&V1,&V2,&S); |
754 | } | 745 | } |
755 | OK | 746 | OK |
756 | } | 747 | } |