Thread: Issue with math

Code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

int xPosCalc(int sigmaMatter, int muMatter, int sigmaAntimatter, int muAntimatter, int xRand)
{
    int seed = time(NULL);
    srand(seed);
    xRand = rand() % 36;
    int xPosMatter = 1/(sigmaMatter*sqrt(2*3.14159))*exp(-((xRand-muMatter)*(xRand-muMatter))/(2*(sigmaMatter*sigmaMatter)));
    int xPosAntimatter = 1/(sigmaAntimatter*sqrt(2*3.14159))*exp(-pow((xRand-muAntimatter), 2)/(2*(pow(sigmaAntimatter, 2))));
    int xProbability = xPosMatter * xPosAntimatter;
    //printf("\nMatter: %d\nAntimatter: %d", xPosMatter, xPosAntimatter);
    //printf("\nProbability is: %d", xProbability);
}

int main(void)
{
    int sigmaMatter = 4.2;
    int muMatter = 27.0;
    int sigmaAntimatter = 6.0;
    int muAntimatter = 17.0;
    int xRand = 0.0;
    xPosCalc(sigmaMatter, muMatter, sigmaAntimatter, muAntimatter, xRand);
    int i, p;
    int count = 0; 
    int pTotal = 0;
    for(i =0; i <= 4999; i=i+1)
    {
        p = 1/(4.2*sqrt(2*3.14159))*exp(-((27-xRand)*(27-xRand))/((2*4.2)*(2*4.2)))*1/(6*sqrt(2*3.14159))*exp(-((17-xRand)*(17-xRand))/((2*6)*(2*6)));
        count++;
        pTotal = pTotal + p;
    }
    int pFinal = pTotal/count;
    printf("\npTotal: %d", pTotal);
    printf("\nTotal iterations: %d", count);
    printf("\nFinal probability: %d\n", pFinal);
}

Alright, p in the main routine is supposed to be the first attachment,

and xPosmatter/xPosantimatter are supposed to be the second attachment.

I'm getting the wrong output, so I just want to make sure I translated the math right before I go any further.

I'll separate out the steps with some temporary variables and see if I can get anywhere.

Also, if you look at the attachments I posted, the square is on the outside of the parentheses. (2*4)^2 does equal (2*4)*(2*4).

Originally Posted by Paul Omans

Also, if you look at the attachments I posted, the square is on the outside of the parentheses. (2*4)^2 does equal (2*4)*(2*4).

Well, I don't see (2*4)^2 in your attachments. I only see 2 * 4.2^2 which is 2 * 4.2 * 4.2.

Bye, Andreas

Since the summation does not contain the summation variable i in (using Maple syntax)

Code:

P := 1/5000 + sum( (1/(4.2*sqrt(2*Pi))) * exp(-(27-x)^2/(2*4.2^2)) / (6*sqrt(2*Pi)) * exp(-(17-x)^2/(2*6^2)), i=0..4999 );

the sum is equivalent to

Code:

P = 0.0000000000001211344564478819687921177341468858685481 * exp(x*(2.002834467120181405895691609977324263038 - 0.04223356009070294784580498866213151927438 * x);

I suspect both x in the summation should have been x_i? In that case, P can be written in pseudocode as

Code:

P = 0.0
for i = [0, 4999]:
    P = P + exp( -4.013888888888888888888888888888888888889
               + x[i] * ( 2.002834467120181405895691609977324263038
                 - x[i] * ( 0.04223356009070294784580498866213151927438 )))
end for
P = P * 0.00000000000000134124305433139340656790345171524598923

---

The second one is a simple Gaussian probability density function. Using Maple syntax,

Code:

p := 1/(sigma * sqrt(2*Pi)) * exp( - (x - mu)^2 / (2*sigma^2) );

it boils down to

Code:

p = 0.3989422804014326779399460599343818684758 / sigma * exp(-0.5*((x-mu)/sigma)*((x-mu)/sigma));

---

I'm a very firm believer in first simplifying the formulas with a computer algebra system (I like Maple, but Maxima is free and open source), then writing them in pseudocode (often awk for quick tests!), and only then implementing them in C or other low-level language. As long as you check your results on each step, you always have a reliable basis you can rely on.

The OP is in a difficult position: it is very hard to say where the problem lies. Could be the math (notation error in original equations), could be the formula, could be the C implementation.

Originally Posted by AndiPersti

Well, I don't see (2*4)^2 in your attachments. I only see 2 * 4.2^2 which is 2 * 4.2 * 4.2.

Bye, Andreas

Oh shoot, you're right. Now I feel pretty stupid. I was looking at the wrong section.

Thanks for all the helpful replies people!! I'm still working on the code, I'll post my progress later tonight (it's difficult to sit indoors and code when it's 85 F outside and breezy!)

Originally Posted by Paul Omans

Thanks for all the helpful replies people!! I'm still working on the code, I'll post my progress later tonight (it's difficult to sit indoors and code when it's 85 F outside and breezy!)

So code outside.

Originally Posted by anduril462

So code outside.

I've heard of people coding while longboarding, but personally I would find that a little dangerous.

Here's what I got:

Code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

int xPosCalc(float sigmaMatter, float muMatter, float sigmaAntimatter, float muAntimatter, float xRand)
{
    float seed = time(NULL);
    srand(seed);
    xRand = rand() % 36;
    float xPosMatter = 1/(sigmaMatter*2.506627216)*exp(-((xRand-muMatter)*(xRand-muMatter))/(2*(sigmaMatter*sigmaMatter)));
    float xPosAntimatter = 1/(sigmaAntimatter*2.506627216)*exp(-pow((xRand-muAntimatter), 2)/(2*(pow(sigmaAntimatter, 2))));
    float xProbability = xPosMatter * xPosAntimatter;
    //printf("\nMatter: %d\nAntimatter: %d", xPosMatter, xPosAntimatter);
    //printf("\nProbability is: %d", xProbability);
}

int main(void)
{
    float sigmaMatter = 4.2;
    float muMatter = 27.0;
    float sigmaAntimatter = 6.0;
    float muAntimatter = 17.0;
    float xRand = 0.0;
    xPosCalc(sigmaMatter, muMatter, sigmaAntimatter, muAntimatter, xRand);
    int i;
    int count = 0; 
    float p = 0;
    float pTotal = 0;
    for(i =0; i <= 4999; i=i+1) //in step 2, addition between should be changed to multiplication
    {
        p = 0.09498629735*exp(-((27-xRand)*(27-xRand))/25.28)*0.06649040814*exp(-((17-xRand)*(17-xRand))/72);
        count++;
        pTotal = pTotal + p;
    }
    float pFinal = pTotal/count;
    printf("\npTotal: %d", pTotal);
    printf("\nTotal iterations: %d", count);
    printf("\nFinal probability: %d\n", pFinal);
}

I condensed the math, and I now get numbers for pTotal and pFinal instead of just zeros. I'll continue working on it tomorrow, I still need to split the functions into separate .c/.h files and create a makefile.