summaryrefslogtreecommitdiff
path: root/packages/base/src/C
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/C')
-rw-r--r--packages/base/src/C/vector-aux.c49
1 files changed, 49 insertions, 0 deletions
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index f1bb371..5b9c171 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -695,3 +695,52 @@ int saveMatrix(char * file, char * format, KDMAT(a)){
695 OK 695 OK
696} 696}
697 697
698////////////////////////////////////////////////////////////////////////////////
699
700// http://c-faq.com/lib/gaussian.html
701double gaussrand()
702{
703 static double V1, V2, S;
704 static int phase = 0;
705 double X;
706
707 if(phase == 0) {
708 do {
709 double U1 = (double)rand() / RAND_MAX;
710 double U2 = (double)rand() / RAND_MAX;
711
712 V1 = 2 * U1 - 1;
713 V2 = 2 * U2 - 1;
714 S = V1 * V1 + V2 * V2;
715 } while(S >= 1 || S == 0);
716
717 X = V1 * sqrt(-2 * log(S) / S);
718 } else
719 X = V2 * sqrt(-2 * log(S) / S);
720
721 phase = 1 - phase;
722
723 return X;
724}
725
726int random_vector(int seed, int code, DVEC(r)) {
727 srand(seed);
728 int k;
729 switch (code) {
730 case 0: { // uniform
731 for (k=0; k<rn; k++) {
732 rp[k] = (double)rand()/RAND_MAX;
733 }
734 OK
735 }
736 case 1: { // gaussian
737 for (k=0; k<rn; k++) {
738 rp[k] = gaussrand();
739 }
740 OK
741 }
742
743 default: ERROR(BAD_CODE);
744 }
745}
746