summaryrefslogtreecommitdiff
path: root/packages/base
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base')
-rw-r--r--packages/base/src/C/vector-aux.c21
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
703inline 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
711double gaussrand(struct drand48_data *buffer, 704double 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
737int random_vector(unsigned int seed, int code, DVEC(r)) { 730int 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 }