Random floats function

Discussion in 'C Programming' started by Michel Rouzic, Oct 12, 2005.

  1. 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<M; i++)
    s=dblrand((s[i-1]+pi) / (2*pi)) * (2*pi) - pi;
    }
     
    Michel Rouzic, Oct 12, 2005
    #1
    1. Advertising

  2. Michel Rouzic

    Michael Mair Guest

    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<M; i++)
    > s=dblrand((s[i-1]+pi) / (2*pi)) * (2*pi) - pi;
    > }
    >


    Notes:
    If the number of random digits does not play that much of a role,
    then use (double)rand()/RANDMAX instead (with an appropriate seed
    set by srand()).
    Otherwise, use another off-the-shelf generator -- it probably has
    better mathematical properties than your generator. This is not
    the right place to discuss this, though.

    I would rather work with symbolic constants instead of
    pow(2, 32) and pi.
    A convenient portable way for the former is
    #define POW2TO32_DBL ((1UL << 31)*2.0)
    For the latter, just take enough digits of pi.


    Cheers
    Michael
    --
    E-Mail: Mine is an /at/ gmx /dot/ de address.
     
    Michael Mair, Oct 12, 2005
    #2
    1. Advertising

  3. 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<M; i++)
    > > s=dblrand((s[i-1]+pi) / (2*pi)) * (2*pi) - pi;
    > > }
    > >

    >
    > Notes:
    > If the number of random digits does not play that much of a role,
    > then use (double)rand()/RANDMAX instead (with an appropriate seed
    > set by srand()).


    ok, so how do I do that (that's where I sound like a newbie), i do
    something like srand(time(NULL)); and then (double) rand()/RANDMAX??
    And I admit it kinda annoys me to use a random generator with only 2^15
    possibilities (right?) but I guess it should be alright.

    > Otherwise, use another off-the-shelf generator -- it probably has
    > better mathematical properties than your generator. This is not
    > the right place to discuss this, though.
    >
    > I would rather work with symbolic constants instead of
    > pow(2, 32) and pi.
    > A convenient portable way for the former is
    > #define POW2TO32_DBL ((1UL << 31)*2.0)


    I have no idea what ((1UL << 31)*2.0) can mean..

    > For the latter, just take enough digits of pi.


    What's wrong with the way I do my pi?
     
    Michel Rouzic, Oct 12, 2005
    #3
  4. 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<M; i++)
    > > s=dblrand((s[i-1]+pi) / (2*pi)) * (2*pi) - pi;
    > > }
    > >

    >
    > Notes:
    > If the number of random digits does not play that much of a role,
    > then use (double)rand()/RANDMAX instead (with an appropriate seed
    > set by srand()).
    > Otherwise, use another off-the-shelf generator -- it probably has
    > better mathematical properties than your generator. This is not
    > the right place to discuss this, though.
    >
    > I would rather work with symbolic constants instead of
    > pow(2, 32) and pi.
    > A convenient portable way for the former is
    > #define POW2TO32_DBL ((1UL << 31)*2.0)
    > For the latter, just take enough digits of pi.
    >
    >
    > Cheers
    > Michael
    > --
    > E-Mail: Mine is an /at/ gmx /dot/ de address.


    I used fmod() instead of my modfloat() and now it works. Well, now the
    values seem to be quite random and in the right range, so I guess I'll
    keep using my dblrand() function. Thanks alot for help!
     
    Michel Rouzic, Oct 12, 2005
    #4
  5. Michel Rouzic

    Michael Mair Guest

    Michel Rouzic wrote:
    > Michael Mair wrote:
    >
    >>Michel Rouzic wrote:

    <snip!>
    >>
    >>>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.


    You are right, I did not really think about that one...
    Question: Why are you not just using
    return fmod(1664525.0 * x, 1.0);
    instead?

    >>Notes:
    >>If the number of random digits does not play that much of a role,
    >>then use (double)rand()/RANDMAX instead (with an appropriate seed
    >>set by srand()).

    >
    > ok, so how do I do that (that's where I sound like a newbie), i do
    > something like srand(time(NULL)); and then (double) rand()/RANDMAX??
    > And I admit it kinda annoys me to use a random generator with only 2^15
    > possibilities (right?) but I guess it should be alright.


    RANDMAX is _at_least_ 2^15-1 -- your implementation's standard
    library is allowed to do better...

    >>I would rather work with symbolic constants instead of
    >>pow(2, 32) and pi.
    >>A convenient portable way for the former is
    >>#define POW2TO32_DBL ((1UL << 31)*2.0)

    >
    > I have no idea what ((1UL << 31)*2.0) can mean..


    Okay, a << N is a*2^N and has the same type as a. As int is
    only guaranteed to have at least 16 bits, I use unsigned long.
    1UL is a constant of type unsigned long with value 1.
    Unfortunately, unsigned long is only guaranteed to have 32 bits
    (and can hold 2^31 but not necessarily 2^32).
    As we need a constant of type double, I just multiply the best
    I can do by an appropriate double factor.
    2^31 as double consequently is (1UL<<31)*1.0, as float
    (1UL<<31)*1.0F and so on.
    I prefer
    ((1UL << 31)*2.0)
    to actually writing
    4294967296.0
    as it is easier to see that the former is a power of 2.

    So, why bother at all with symbolic constants?
    My solution gives you a compile time constant whereas
    your solution gives you something which may be computed each
    and every time. In addition, you can emphasize the _role_
    of the constant to tell one nameless 2^32 from another which
    comes in handy when you have to change it...


    >>For the latter, just take enough digits of pi.

    >
    > 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.


    Cheers
    Michael
    --
    E-Mail: Mine is an /at/ gmx /dot/ de address.
     
    Michael Mair, Oct 12, 2005
    #5
  6. Michael Mair wrote:
    > Michel Rouzic wrote:
    > > Michael Mair wrote:
    > >
    > >>Michel Rouzic wrote:

    > <snip!>
    > >>
    > >>>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.

    >
    > You are right, I did not really think about that one...
    > Question: Why are you not just using
    > return fmod(1664525.0 * x, 1.0);
    > instead?


    oh yeah of course, it's much simpler.

    > >>Notes:
    > >>If the number of random digits does not play that much of a role,
    > >>then use (double)rand()/RANDMAX instead (with an appropriate seed
    > >>set by srand()).

    > >
    > > ok, so how do I do that (that's where I sound like a newbie), i do
    > > something like srand(time(NULL)); and then (double) rand()/RANDMAX??
    > > And I admit it kinda annoys me to use a random generator with only 2^15
    > > possibilities (right?) but I guess it should be alright.

    >
    > RANDMAX is _at_least_ 2^15-1 -- your implementation's standard
    > library is allowed to do better...
    >
    > >>I would rather work with symbolic constants instead of
    > >>pow(2, 32) and pi.
    > >>A convenient portable way for the former is
    > >>#define POW2TO32_DBL ((1UL << 31)*2.0)

    > >
    > > I have no idea what ((1UL << 31)*2.0) can mean..

    >
    > Okay, a << N is a*2^N and has the same type as a. As int is
    > only guaranteed to have at least 16 bits, I use unsigned long.
    > 1UL is a constant of type unsigned long with value 1.
    > Unfortunately, unsigned long is only guaranteed to have 32 bits
    > (and can hold 2^31 but not necessarily 2^32).
    > As we need a constant of type double, I just multiply the best
    > I can do by an appropriate double factor.
    > 2^31 as double consequently is (1UL<<31)*1.0, as float
    > (1UL<<31)*1.0F and so on.
    > I prefer
    > ((1UL << 31)*2.0)
    > to actually writing
    > 4294967296.0
    > as it is easier to see that the former is a power of 2.
    >
    > So, why bother at all with symbolic constants?
    > My solution gives you a compile time constant whereas
    > your solution gives you something which may be computed each
    > and every time. In addition, you can emphasize the _role_
    > of the constant to tell one nameless 2^32 from another which
    > comes in handy when you have to change it...


    OK, so it's not going to compute it everytime?

    > >>For the latter, just take enough digits of pi.

    > >
    > > 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
     
    Michel Rouzic, Oct 13, 2005
    #6
  7. "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) <http://www.ghoti.net/~kst>
    San Diego Supercomputer Center <*> <http://users.sdsc.edu/~kst>
    We must do something. This is something. Therefore, we must do this.
     
    Keith Thompson, Oct 13, 2005
    #7
  8. Michel Rouzic

    pete Guest

    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 <stdio.h>
    #include <float.h>
    #include <math.h>

    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
     
    pete, Oct 13, 2005
    #8
  9. Michel Rouzic

    pete Guest

    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
     
    pete, Oct 13, 2005
    #9
  10. Michel Rouzic

    pete Guest

    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
     
    pete, Oct 13, 2005
    #10
  11. 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) <http://www.ghoti.net/~kst>
    > San Diego Supercomputer Center <*> <http://users.sdsc.edu/~kst>
    > We must do something. This is something. Therefore, we must do this.
     
    Michel Rouzic, Oct 14, 2005
    #11
  12. Michel Rouzic

    Michael Mair Guest

    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 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...


    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.
     
    Michael Mair, Oct 14, 2005
    #12
  13. Michel Rouzic

    pete Guest

    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.1415926535897932384626433832795028841971693993751
    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.1415926535897932384626433832795028841971693993751;
    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 <stdio.h>
    #include <float.h>
    #include <math.h>

    double pi(void);
    double pi2(void);
    double pi3(void);

    int main(void)
    {
    double Pi, Pi2, Pi3, PI;

    PI = 3.1415926535897932384626433832795028841971693993751;
    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
     
    pete, Oct 14, 2005
    #13
  14. Michel Rouzic

    pete Guest

    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.1415926535897932384626433832795028841971693993751;
    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 <stdio.h>
    #include <float.h>
    #include <math.h>

    double pi(void);

    int main(void)
    {
    double Pi, PI;

    PI = 3.1415926535897932384626433832795028841971693993751;
    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
     
    pete, Oct 14, 2005
    #14
  15. Michel Rouzic

    pete Guest

    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
     
    pete, Oct 14, 2005
    #15
  16. 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 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...

    >
    > 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...
     
    Michel Rouzic, Oct 14, 2005
    #16
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Kosio

    Floats to chars and chars to floats

    Kosio, Sep 16, 2005, in forum: C Programming
    Replies:
    44
    Views:
    1,311
    Tim Rentsch
    Sep 23, 2005
  2. globalrev
    Replies:
    4
    Views:
    784
    Gabriel Genellina
    Apr 20, 2008
  3. VK
    Replies:
    15
    Views:
    1,238
    Dr J R Stockton
    May 2, 2010
  4. Norah Jones
    Replies:
    7
    Views:
    137
    llanitedave
    Mar 13, 2013
  5. Boris FELD
    Replies:
    0
    Views:
    107
    Boris FELD
    Mar 12, 2013
Loading...

Share This Page