Thread: Vector buildings and Box Muller transform

  1. #1
    Registered User DeliriumCordia's Avatar
    Join Date
    Mar 2012
    Posts
    59

    Vector buildings and Box Muller transform

    Hi guys, I got 2 problems:

    • I need to build an array that has a function similar to a triangle. My function goes from -40 to 40 with steps of 0.002. So it has to be a 40000 elements array, made like that:

    var filler goes from -40 to 40 with step of 0.002


    function is, for each element:
    (40.0-fabs(filler))/40.0

    so I expect my triangle starting from o (40-|-40|)/40=0/40=0, with an apex of (40-0)/40=1, and getting back to 0 at the end.

    Here is how I implemented it:

    Code:
    filler=-40.0;
        for(i=0;i<40000;i++){
            VettoreIniziale[i].img=0.0;
            VettoreIniziale[i].real=(40.0-fabs(filler))/40.0;
            filler+=0.002;
              }


    It starts from 0.0000 but it has an apex of about 0.998 and end with a negative number, like -0.004.. Why it happens? What did I do wrong?


    • Second problem. I need to fill an array with random numbers in the range [-1,1] with mean = 0 and variance =1. I tried to use Box Muller transform, but it gives me strange results... Here the main and the method.


    Box Muller
    Code:
    #include <stdlib.h>double rand_normal(double mean, double stddev) {
        static double n2 = 0.0;
        static int n2_cached = 0;
        if (!n2_cached) {
            double x, y, r;
    do {
       x = 2.0*rand()/RAND_MAX - 1;
    y = 2.0*rand()/RAND_MAX - 1;
        r = x*x + y*y;
    } while (r == 0.0 || r > 1.0);
            {
            double d = sqrt(-2.0*log(r)/r);
    double n1 = x*d;
    n2 = y*d;
            double result = n1*stddev + mean;
            n2_cached = 1;
            return result;
            }
        } else {
            n2_cached = 0;
            return n2*stddev + mean;
        }
    }

    This is the main code regarding it:
    Code:
    #include <stdio.h>
    #include <stdlib.h>
    #include "MISC.h"
    #include <math.h>
    #include <time.h>
    
    ...
    srand(time(NULL));
    ...
    
    for(i=0;i<LENGTH;i++){        GaussianWhiteNoise[i].real=(float)rand_normal(0.0,1.0);
            GaussianWhiteNoise[i].img=0.0;
            printf("%d %3.4f \n",i,GaussianWhiteNoise[i].real);
    I tried not to use (float) to convert, but it happens the same.

    This is the output...

    0 4.0000
    1 2665544.0000
    2 4.0000
    3 2665544.0000
    4 0.0000
    5 2665544.0000
    6 0.0000
    7 2665544.0000
    8 0.0000
    9 2665544.0000
    10 0.0000
    11 2665544.0000
    12 0.0000
    13 2665544.0000
    14 0.0000
    15 2665544.0000
    16 0.0000


    etc etc

  2. #2
    Registered User DeliriumCordia's Avatar
    Join Date
    Mar 2012
    Posts
    59
    Don't care about 2nd problem. I forgot to put function prototype on the main -.-

  3. #3
    - - - - - - - - oogabooga's Avatar
    Join Date
    Jan 2008
    Posts
    2,808
    It happens because fractional values are not stored exactly in floating-point numbers.
    So every time you add in 0.002 you accumulate some error.

    Try this:
    Code:
        for (i = 0; i < 40000; i++){
            filler = -40 + 0.002 * i;
            VettoreIniziale[i].img = 0.0;
            VettoreIniziale[i].real = (40.0 - fabs(filler)) / 40.0;
        }
    You'll still get some error, but it won't accumulate.
    The cost of software maintenance increases with the square of the programmer's creativity. - Robert D. Bliss

  4. #4
    Registered User DeliriumCordia's Avatar
    Join Date
    Mar 2012
    Posts
    59
    Thank you very much, you've been very kind

    It can be closed if necessary.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Transform help
    By fadlan12 in forum C++ Programming
    Replies: 7
    Last Post: 04-27-2010, 07:32 PM
  2. how to transform this loop..
    By transgalactic2 in forum C Programming
    Replies: 17
    Last Post: 01-17-2009, 02:10 AM
  3. Ever heard of a Grapple? (No, not for scaling buildings.)
    By LuckY in forum A Brief History of Cprogramming.com
    Replies: 7
    Last Post: 01-24-2006, 03:53 PM
  4. transform help!!!
    By what3v3r in forum C++ Programming
    Replies: 7
    Last Post: 01-16-2006, 10:27 PM
  5. Fourier Transform
    By srikanthreddy in forum C++ Programming
    Replies: 2
    Last Post: 04-09-2005, 08:36 AM