448,589 Members | 1,202 Online Need help? Post your question and get tips & solutions from a community of 448,589 IT Pros & Developers. It's quick & easy.

# Random floats function

 P: n/a I tried making a function that would fill half of an array with random doubles that would go from -pi to +pi. Unfortunatly, all it ever returns is -pi. I haven't used an already existing random number generating function because I wasn't sure to know how to use them so I created my own based on the random number generator from the INMOS Transputer Development system which is x[n+1] = (1664525 x[n]) mod 2^32 double modfloat(double a, double b) { double c, d, x; c=a/b; d=(int32_t) c; x=c-d; return x; } double dblrand(double x) { return modfloat((1664525 * x * pow(2, 32)), pow(2, 32)) / pow(2, 32); } void phasenoise(double *s, int32_t M) { int32_t i; time_t seed=time(NULL); double pi=atan(1.0)*4.0; s[M/2+1]=dblrand((double) seed/pow(2, 32)) * (2*pi) - pi; for (i=M/2+2; i
15 Replies

 P: n/a Michel Rouzic wrote: I tried making a function that would fill half of an array with random doubles that would go from -pi to +pi. Unfortunatly, all it ever returns is -pi. I haven't used an already existing random number generating function because I wasn't sure to know how to use them so I created my own based on the random number generator from the INMOS Transputer Development system which is x[n+1] = (1664525 x[n]) mod 2^32 double modfloat(double a, double b) { double c, d, x; c=a/b; d=(int32_t) c; uint32_t, if anything. However, this is not guaranteed to work either. x=c-d; return x; } BTW: You do need to invent modfloat(), use fmod() instead. double dblrand(double x) { return modfloat((1664525 * x * pow(2, 32)), pow(2, 32)) / pow(2, 32); ^^^^^^^^^^ Lose that one. a*b%b == 0 So, you always get back something near 0.0, which leads with your scaling to -pi. } void phasenoise(double *s, int32_t M) { int32_t i; time_t seed=time(NULL); double pi=atan(1.0)*4.0; s[M/2+1]=dblrand((double) seed/pow(2, 32)) * (2*pi) - pi; for (i=M/2+2; i

 P: n/a Michael Mair wrote: Michel Rouzic wrote: I tried making a function that would fill half of an array with random doubles that would go from -pi to +pi. Unfortunatly, all it ever returns is -pi. I haven't used an already existing random number generating function because I wasn't sure to know how to use them so I created my own based on the random number generator from the INMOS Transputer Development system which is x[n+1] = (1664525 x[n]) mod 2^32 double modfloat(double a, double b) { double c, d, x; c=a/b; d=(int32_t) c; uint32_t, if anything. However, this is not guaranteed to work either. x=c-d; return x; } BTW: You do need to invent modfloat(), use fmod() instead. oh thanks i had never heard about it before. double dblrand(double x) { return modfloat((1664525 * x * pow(2, 32)), pow(2, 32)) / pow(2, 32); ^^^^^^^^^^ Lose that one. a*b%b == 0 So, you always get back something near 0.0, which leads with your scaling to -pi. yeah but no, x is a double between 0 and 1, if i do that it's to get back to the previous result before it's divided. } void phasenoise(double *s, int32_t M) { int32_t i; time_t seed=time(NULL); double pi=atan(1.0)*4.0; s[M/2+1]=dblrand((double) seed/pow(2, 32)) * (2*pi) - pi; for (i=M/2+2; i

 P: n/a Michael Mair wrote: Michel Rouzic wrote: I tried making a function that would fill half of an array with random doubles that would go from -pi to +pi. Unfortunatly, all it ever returns is -pi. I haven't used an already existing random number generating function because I wasn't sure to know how to use them so I created my own based on the random number generator from the INMOS Transputer Development system which is x[n+1] = (1664525 x[n]) mod 2^32 double modfloat(double a, double b) { double c, d, x; c=a/b; d=(int32_t) c; uint32_t, if anything. However, this is not guaranteed to work either. x=c-d; return x; } BTW: You do need to invent modfloat(), use fmod() instead. double dblrand(double x) { return modfloat((1664525 * x * pow(2, 32)), pow(2, 32)) / pow(2, 32); ^^^^^^^^^^ Lose that one. a*b%b == 0 So, you always get back something near 0.0, which leads with your scaling to -pi. } void phasenoise(double *s, int32_t M) { int32_t i; time_t seed=time(NULL); double pi=atan(1.0)*4.0; s[M/2+1]=dblrand((double) seed/pow(2, 32)) * (2*pi) - pi; for (i=M/2+2; i

 P: n/a "Michel Rouzic" writes: Michael Mair wrote: Michel Rouzic wrote: [...] > What's wrong with the way I do my pi? 1) The C standard makes no guarantee whatsoever about the precision of the math functions from the library. So, "your" pi could be way off. 2) Your pi has to be recomputed at every function call. oh, ok, well, I guess it's gonna be better if I stop doing that and put #define pi 3.1415926535897932 instead That's 17 significant digits; long double often has more precision than that. Using 40 decimal digits will cover anything with a mantissa up to 128 bits, which is enough for any hardware I've ever used. Use an 'L' suffix to make sure you get the full precision. It's very likely that you can get away with far less precision, depending on the application, but using more digits than you'll ever need isn't a burden, and it means one less thing to worry about if your program misbehaves. -- Keith Thompson (The_Other_Keith) ks***@mib.org San Diego Supercomputer Center <*> We must do something. This is something. Therefore, we must do this. Nov 15 '05 #7

 P: n/a Keith Thompson wrote: "Michel Rouzic" writes: Michael Mair wrote: Michel Rouzic wrote: [...] > What's wrong with the way I do my pi? 1) The C standard makes no guarantee whatsoever about the precision of the math functions from the library. So, "your" pi could be way off. 2) Your pi has to be recomputed at every function call. oh, ok, well, I guess it's gonna be better if I stop doing that and put #define pi 3.1415926535897932 instead That's 17 significant digits; long double often has more precision than that. Using 40 decimal digits will cover anything with a mantissa up to 128 bits, which is enough for any hardware I've ever used. Use an 'L' suffix to make sure you get the full precision. It's very likely that you can get away with far less precision, depending on the application, but using more digits than you'll ever need isn't a burden, and it means one less thing to worry about if your program misbehaves. I got all the precision of a double right here: /* BEGIN pi.c */ #include #include #include double pi(void); double pi2(void); int main(void) { static double Pi, Pi2; Pi = pi(); Pi2 = pi2(); printf("Pi is %f\n", Pi); printf("Pi2 is %f\n", Pi2); printf("Pi - 4 * atan(1) is %e\n", Pi - 4 * atan(1)); printf("Pi2 - 4 * atan(1) is %e\n", Pi2 - 4 * atan(1)); printf("Pi - Pi2 is %e\n", Pi - Pi2); return 0; } double pi(void) { double pi, b, first, second, numerator; unsigned denominator; pi = 0; numerator = 1 / 5.0; denominator = 1; do { first = numerator / denominator; denominator += 2; numerator /= 25.0; second = numerator / denominator; denominator += 2; numerator /= 25.0; b = first - second; pi += b; } while (b > DBL_EPSILON); pi *= 4; numerator = 1 / 239.0; denominator = 1; do { first = numerator / denominator; denominator += 2; numerator /= 57121.0; second = numerator / denominator; denominator += 2; numerator /= 57121.0; b = first - second; pi -= b; } while (b > DBL_EPSILON); return 4 * pi; } double pi2(void) { double pi, b, first, second, numerator; unsigned denominator; pi = 0; numerator = 1 / 2.0; denominator = 1; do { first = numerator / denominator; denominator += 2; numerator /= 4.0; second = numerator / denominator; denominator += 2; numerator /= 4.0; b = first - second; pi += b; } while (b > DBL_EPSILON); numerator = 1 / 3.0; denominator = 1; do { first = numerator / denominator; denominator += 2; numerator /= 9.0; second = numerator / denominator; denominator += 2; numerator /= 9.0; b = first - second; pi += b; } while (b > DBL_EPSILON); return 4 * pi; } /* END pi.c */ -- pete Nov 15 '05 #8

 P: n/a pete wrote: Based on the magnitude of the value of the pi variable when the loops complete, I think better epsilon numbers would be: double pi(void) { double pi, b, first, second, numerator; unsigned denominator; pi = 0; numerator = 1 / 5.0; denominator = 1; do { first = numerator / denominator; denominator += 2; numerator /= 25.0; second = numerator / denominator; denominator += 2; numerator /= 25.0; b = first - second; pi += b; } while (b > DBL_EPSILON); } while (b > 0.125 * DBL_EPSILON); pi *= 4; numerator = 1 / 239.0; denominator = 1; do { first = numerator / denominator; denominator += 2; numerator /= 57121.0; second = numerator / denominator; denominator += 2; numerator /= 57121.0; b = first - second; pi -= b; } while (b > DBL_EPSILON); } while (b > 0.5 * DBL_EPSILON); return 4 * pi; } double pi2(void) { double pi, b, first, second, numerator; unsigned denominator; pi = 0; numerator = 1 / 2.0; denominator = 1; do { first = numerator / denominator; denominator += 2; numerator /= 4.0; second = numerator / denominator; denominator += 2; numerator /= 4.0; b = first - second; pi += b; } while (b > DBL_EPSILON); } while (b > 0.25 * DBL_EPSILON); numerator = 1 / 3.0; denominator = 1; do { first = numerator / denominator; denominator += 2; numerator /= 9.0; second = numerator / denominator; denominator += 2; numerator /= 9.0; b = first - second; pi += b; } while (b > DBL_EPSILON); } while (b > 0.5 * DBL_EPSILON); return 4 * pi; } /* END pi.c */ -- pete Nov 15 '05 #9

 P: n/a pete wrote: Keith Thompson wrote: "Michel Rouzic" writes: Michael Mair wrote:> Michel Rouzic wrote: [...]> > What's wrong with the way I do my pi?>> 1) The C standard makes no guarantee whatsoever about the> precision of the math functions from the library. So, "your"> pi could be way off.> 2) Your pi has to be recomputed at every function call. oh, ok, well, I guess it's gonna be better if I stop doing that and put #define pi 3.1415926535897932 instead That's 17 significant digits; long double often has more precision than that. Using 40 decimal digits will cover anything with a mantissa up to 128 bits, which is enough for any hardware I've ever used. Use an 'L' suffix to make sure you get the full precision. It's very likely that you can get away with far less precision, depending on the application, but using more digits than you'll ever need isn't a burden, and it means one less thing to worry about if your program misbehaves. I got all the precision of a double right here: /* BEGIN pi.c */ Subsequent testing indicates maybe I don't got all the precision of a double. -- pete Nov 15 '05 #10

 P: n/a Keith Thompson wrote: "Michel Rouzic" writes: Michael Mair wrote: Michel Rouzic wrote: [...] > What's wrong with the way I do my pi? 1) The C standard makes no guarantee whatsoever about the precision of the math functions from the library. So, "your" pi could be way off. 2) Your pi has to be recomputed at every function call. oh, ok, well, I guess it's gonna be better if I stop doing that and put #define pi 3.1415926535897932 instead That's 17 significant digits; long double often has more precision than that. I'm only using double precision floats, and since I heard it could take between 15 and 17 signifiant digits... Using 40 decimal digits will cover anything with a mantissa up to 128 bits, which is enough for any hardware I've ever used. Use an 'L' suffix to make sure you get the full precision. It's very likely that you can get away with far less precision, depending on the application, but using more digits than you'll ever need isn't a burden, and it means one less thing to worry about if your program misbehaves. -- Keith Thompson (The_Other_Keith) ks***@mib.org San Diego Supercomputer Center <*> We must do something. This is something. Therefore, we must do this. Nov 15 '05 #11

 P: n/a Michel Rouzic wrote: Keith Thompson wrote:"Michel Rouzic" writes:Michael Mair wrote:Michel Rouzic wrote:[...]>What's wrong with the way I do my pi?1) The C standard makes no guarantee whatsoever about theprecision of the math functions from the library. So, "your"pi could be way off.2) Your pi has to be recomputed at every function call.oh, ok, well, I guess it's gonna be better if I stop doing that and put#define pi 3.1415926535897932 insteadThat's 17 significant digits; long double often has more precisionthan that. I'm only using double precision floats, and since I heard it could take between 15 and 17 signifiant digits... This may be the case for your implementation (platform/OS/compiler/standard library combination) but the C standard does not guarantee it. Whether or not double and long double are of different precision also is unspecified. Probably you have IEEE-like doubles, so this is moot but more digits do not hurt in _any_ respect. Whether or not it is a good idea to multiply with a long double constant may depend on what you want to do and your implementation. Cheers Michael -- E-Mail: Mine is an /at/ gmx /dot/ de address. Nov 15 '05 #12

 P: n/a pete wrote: pete wrote: Keith Thompson wrote: "Michel Rouzic" writes: > Michael Mair wrote: >> Michel Rouzic wrote: [...] >> > What's wrong with the way I do my pi? >> >> 1) The C standard makes no guarantee whatsoever about the >> precision of the math functions from the library. So, "your" >> pi could be way off. >> 2) Your pi has to be recomputed at every function call. > > oh, ok, well, I guess it's gonna be better > if I stop doing that and put > #define pi 3.1415926535897932 instead That's 17 significant digits; long double often has more precision than that. Using 40 decimal digits will cover anything with a mantissa up to 128 bits, which is enough for any hardware I've ever used. Use an 'L' suffix to make sure you get the full precision. It's very likely that you can get away with far less precision, depending on the application, but using more digits than you'll ever need isn't a burden, and it means one less thing to worry about if your program misbehaves. I got all the precision of a double right here: /* BEGIN pi.c */ Subsequent testing indicates maybe I don't got all the precision of a double. I'm pretty close though. My feeling is that since (double)3.1415926535897932384626433832795028841971 693993751 compares equal to 4 * atan(1) on my system, then, 4 * atan(1) is probably correct. The return value of pi(), seems to be low by 2 * DBL_EPSILON. I tried adding the positive half of an extra term in each loop in pi3, but it didn't help. /* BEGIN pi.c output */ PI = 3.141592653589793238462643383279502884197169399375 1; Pi = pi(); Pi2 = pi2(); Pi3 = pi3(); PI is 3.141593 Pi is 3.141593 Pi2 is 3.141593 Pi3 is 3.141593 PI - 4 * atan(1) is 0.000000e+000 Pi - 4 * atan(1) is -4.440892e-016 Pi2 - 4 * atan(1) is -4.440892e-016 Pi3 - 4 * atan(1) is -4.440892e-016 DBL_EPSILON * 2 is 4.440892e-016 Pi - Pi2 is 0.000000e+000 Pi - PI + DBL_EPSILON * 2 is 0.000000e+000 /* END pi.c output */ /* BEGIN pi.c */ #include #include #include double pi(void); double pi2(void); double pi3(void); int main(void) { double Pi, Pi2, Pi3, PI; PI = 3.141592653589793238462643383279502884197169399375 1; Pi = pi(); Pi2 = pi2(); Pi3 = pi3(); puts("/* BEGIN pi.c output */\n"); puts("PI = 3.1415926535897932384626" "433832795028841971693993751;"); puts("Pi = pi();\nPi2 = pi2();\nPi3 = pi3();\n"); printf("PI is %f\n", PI); printf("Pi is %f\n", Pi); printf("Pi2 is %f\n", Pi2); printf("Pi3 is %f\n", Pi3); printf("PI - 4 * atan(1) is %e\n", PI - 4 * atan(1)); printf("Pi - 4 * atan(1) is %e\n", Pi - 4 * atan(1)); printf("Pi2 - 4 * atan(1) is %e\n", Pi2 - 4 * atan(1)); printf("Pi3 - 4 * atan(1) is %e\n", Pi3 - 4 * atan(1)); printf("DBL_EPSILON * 2 is %e\n", DBL_EPSILON * 2); printf("Pi - Pi2 is %e\n", Pi - Pi2); printf("Pi - PI + DBL_EPSILON * 2 is %e\n", Pi - PI + DBL_EPSILON * 2); puts("\n/* END pi.c output */"); return 0; } double pi(void) { double pi, b, numerator; unsigned denominator; pi = 0; numerator = 1 / 5.0; denominator = 1; do { b = numerator / denominator; numerator /= 25; denominator += 2; b -= numerator / denominator; numerator /= 25; denominator += 2; pi += b; } while (b > DBL_EPSILON / 8); pi *= 4; numerator = 1 / 239.0; denominator = 1; do { b = numerator / denominator; numerator /= 57121; denominator += 2; b -= numerator / denominator; numerator /= 57121; denominator += 2; pi -= b; } while (b > DBL_EPSILON / 2); return 4 * pi; } double pi2(void) { double pi, b, numerator; unsigned denominator; pi = 0; numerator = 1 / 2.0; denominator = 1; do { b = numerator / denominator; numerator /= 4; denominator += 2; b -= numerator / denominator; numerator /= 4; denominator += 2; pi += b; } while (b > DBL_EPSILON / 4); numerator = 1 / 3.0; denominator = 1; do { b = numerator / denominator; numerator /= 9; denominator += 2; b -= numerator / denominator; numerator /= 9; denominator += 2; pi += b; } while (b > DBL_EPSILON / 2); return 4 * pi; } double pi3(void) { double pi, b, numerator; unsigned denominator; pi = 0; numerator = 1 / 2.0; denominator = 1; do { b = numerator / denominator; numerator /= 4; denominator += 2; b -= numerator / denominator; numerator /= 4; denominator += 2; pi += b; } while (b > DBL_EPSILON / 4); b = numerator / denominator; pi += b; numerator = 1 / 3.0; denominator = 1; do { b = numerator / denominator; numerator /= 9; denominator += 2; b -= numerator / denominator; numerator /= 9; denominator += 2; pi += b; } while (b > DBL_EPSILON / 2); b = numerator / denominator; pi += b; return 4 * pi; } /* END pi.c */ -- pete Nov 15 '05 #13

 P: n/a pete wrote: pete wrote: pete wrote: Keith Thompson wrote: > > "Michel Rouzic" writes: > > Michael Mair wrote: > >> Michel Rouzic wrote: > [...] > >> > What's wrong with the way I do my pi? > >> > >> 1) The C standard makes no guarantee whatsoever about the > >> precision of the math functions from the library. So, "your" > >> pi could be way off. > >> 2) Your pi has to be recomputed at every function call. > > > > oh, ok, well, I guess it's gonna be better > > if I stop doing that and put > > #define pi 3.1415926535897932 instead > > That's 17 significant digits; long double often has more precision > than that. > > Using 40 decimal digits will cover > anything with a mantissa up to 128 > bits, which is enough for any hardware I've ever used. Use an 'L' > suffix to make sure you get the full precision. > > It's very likely that you can get away with far less precision, > depending on the application, but using more digits than you'll ever > need isn't a burden, and it means one less thing to worry about if > your program misbehaves. I got all the precision of a double right here: /* BEGIN pi.c */ Subsequent testing indicates maybe I don't got all the precision of a double. Now, as far as my system is concerned, testing indicates that I *do* got all the precision of a double. /* BEGIN pi.c output */ PI = 3.141592653589793238462643383279502884197169399375 1; Pi = pi(); PI is 3.141593 Pi is 3.141593 PI - 4 * atan(1) is 0.000000e+000 Pi - 4 * atan(1) is 0.000000e+000 Pi - PI is 0.000000e+000 /* END pi.c output */ /* BEGIN pi.c */ #include #include #include double pi(void); int main(void) { double Pi, PI; PI = 3.141592653589793238462643383279502884197169399375 1; Pi = pi(); puts("/* BEGIN pi.c output */\n"); puts("PI = 3.1415926535897932384626" "433832795028841971693993751;"); puts("Pi = pi();\n"); printf("PI is %f\n", PI); printf("Pi is %f\n", Pi); printf("PI - 4 * atan(1) is %e\n", PI - 4 * atan(1)); printf("Pi - 4 * atan(1) is %e\n", Pi - 4 * atan(1)); printf("Pi - PI is %e\n", Pi - PI); puts("\n/* END pi.c output */"); return 0; } double pi(void) { double pi, b, numerator; unsigned denominator; pi = 0; numerator = 3; denominator = 1; do { numerator /= 9; b = numerator / denominator; numerator /= 9; denominator += 2; b -= numerator / denominator; denominator += 2; pi += b; } while (b > DBL_EPSILON / 2); numerator = 2; denominator = 1; do { numerator /= 4; b = numerator / denominator; numerator /= 4; denominator += 2; b -= numerator / denominator; denominator += 2; pi += b; } while (b > DBL_EPSILON / 4); return 4 * pi; } /* END pi.c */ -- pete Nov 15 '05 #14

 P: n/a pete wrote: pete wrote: pete wrote: pete wrote: > > Keith Thompson wrote: > > > > "Michel Rouzic" writes: > > > Michael Mair wrote: > > >> Michel Rouzic wrote: > > [...] > > >> > What's wrong with the way I do my pi? > > >> > > >> 1) The C standard makes no guarantee whatsoever about the > > >> precision of the math functions from the library. So, "your" > > >> pi could be way off. > > >> 2) Your pi has to be recomputed at every function call. > > > > > > oh, ok, well, I guess it's gonna be better > > > if I stop doing that and put > > > #define pi 3.1415926535897932 instead > > > > That's 17 significant digits; > > long double often has more precision > > than that. > > > > Using 40 decimal digits will cover > > anything with a mantissa up to 128 > > bits, which is enough for any hardware I've ever used. > > Use an 'L' > > suffix to make sure you get the full precision. > > > > It's very likely that you can > > get away with far less precision, > > depending on the application, > > but using more digits than you'll ever > > need isn't a burden, > > and it means one less thing to worry about if > > your program misbehaves. > > I got all the precision of a double right here: > > /* BEGIN pi.c */ Subsequent testing indicates maybe I don't got all the precision of a double. Now, as far as my system is concerned, testing indicates that I *do* got all the precision of a double. And even though it seems to work on my system, these next two lines should be switched because of the value of the pi variable when the loops terminate. } while (b > DBL_EPSILON / 2); } while (b > DBL_EPSILON / 4); -- pete Nov 15 '05 #15

 P: n/a Michael Mair wrote: Michel Rouzic wrote: Keith Thompson wrote:"Michel Rouzic" writes:Michael Mair wrote:>Michel Rouzic wrote:[...]>>What's wrong with the way I do my pi?>>1) The C standard makes no guarantee whatsoever about the>precision of the math functions from the library. So, "your">pi could be way off.>2) Your pi has to be recomputed at every function call.oh, ok, well, I guess it's gonna be better if I stop doing that and put#define pi 3.1415926535897932 insteadThat's 17 significant digits; long double often has more precisionthan that. I'm only using double precision floats, and since I heard it could take between 15 and 17 signifiant digits... This may be the case for your implementation (platform/OS/compiler/standard library combination) but the C standard does not guarantee it. Whether or not double and long double are of different precision also is unspecified. Probably you have IEEE-like doubles, so this is moot but more digits do not hurt in _any_ respect. Whether or not it is a good idea to multiply with a long double constant may depend on what you want to do and your implementation. I guess I could make it longer. Have an idea on how long is the longest mantissa in the different existing doubles? Anyways, if I'm using pi, it's for DSP, and since an imprecision of 10^-16 represents -160 dB, it's way good enough. But I guess I could make it long enough for the largest mantissa existing, doesn't cost anything but a few bytes in my code... Nov 15 '05 #16

### This discussion thread is closed

Replies have been disabled for this discussion. 