2

I have C code that was used for a paper. I wanted to write the exact code in Python.

Here is everything that is needed:

The prime function that is used in C code is:

#include "mrand_seeds.h"
#define norm   2.328306549295728e-10  /* 1.0/(m1+1) */
#define norm2  2.328318825240738e-10  /* 1.0/(m2+1) */
#define m1     4294967087.0
#define m2     4294944443.0

double mrand(int stream)
{
    long k;
    double p,
           s10 = drng[stream][0], s11 = drng[stream][1], s12 = drng[stream][2],
           s20 = drng[stream][3], s21 = drng[stream][4], s22 = drng[stream][5];

    p = 1403580.0 * s11 - 810728.0 * s10;
    k = p / m1;  p -= k*m1;  if (p < 0.0) p += m1;
    s10 = s11;   s11 = s12;  s12 = p;

    p = 527612.0 * s22 - 1370589.0 * s20;
    k = p / m2;  p -= k*m2;  if (p < 0.0) p += m2;
    s20 = s21;   s21 = s22;  s22 = p; 

    drng[stream][0] = s10;  drng[stream][1] = s11;  drng[stream][2] = s12;
    drng[stream][3] = s20;  drng[stream][4] = s21;  drng[stream][5] = s22;

    if (s12 <= s22) return ((s12 - s22 + m1) * norm);
    else return ((s12 - s22) * norm);
}

And the drng is a list of 60000 integers in mrand_seeds.h: (below list is not the complete list from file)

static double drng[][6] =
{
           0,           0,           1,           0,           0,           1, 
  1772212344,  1374954571,  2377447708,   540628578,  1843308759,   549575061, 
  2602294560,  1764491502,  3872775590,  4089362440,  2683806282,   437563332, 
   376810349,  1545165407,  3443838735,  3650079346,  1898051052,  2606578666, 
  1847817841,  3038743716,  2014183350,  2883836363,  3242147124,  1955620878, 
  1075987441,  3468627582,  2694529948,   368150488,  2026479331,  2067041056, 
   134547324,  4246812979,  1700384422,  2358888058,    83616724,  3045736624, 
  2816844169,   885735878,  1824365395,  2629582008,  3405363962,  1835381773, 
   675808621,   434584068,  4021752986,  3831444678,  4193349505,  2833414845, 
  2876117643,  1466108979,   163986545,  1530526354,    68578399,  1111539974, 
   411040508,   544377427,  2887694751,   702892456,   758163486,  2462939166};

Now I wrote the mrand function in Python:

M1 = 4294967087
M2 = 4294944443
NORM1 = 2.328306549295728e-10
NORM2 = 2.328318825240738e-10

def mrand(stream):
    s10 = drng1[stream][0]
    s11 = drng1[stream][1]
    s12 = drng1[stream][2]
    s20 = drng1[stream][3]
    s21 = drng1[stream][4]
    s22 = drng1[stream][5]

    p = 1403580.0 * s11 - 810728.0 * s10
    k = p / M1
    p -= k * M1

    if p < 0:
        p += M1
    s10 = s11
    s11 = s12
    s12 = p

    p = 527612.0 * s22 - 1370589.0 * s20
    k = p / M2
    p -= k*M2

    if p < 0.0:
        p += M2
    s20 = s21
    s21 = s22
    s22 = p

    drng1[stream][0] = s10
    drng1[stream][1] = s11
    drng1[stream][2] = s12
    drng1[stream][3] = s20
    drng1[stream][4] = s21
    drng1[stream][5] = s22
    
    if s12 <= s22 :
        return ((s12 - s22 + M1) * NORM1)
    else:
        return ((s12 - s22) * NORM1)

And defined the list in another .py file and imported to Python code and converted it to a 2D array.

import myfunc
drng = myfunc.retlist()
drng1 = [drng[i:i+6] for i in range(0, len(drng), 6)]

The rtlist function simply defines the list and returns it.

Now my problem is when I'm executing C code I get different output with different parameter but in Python I always get 0.9999999997671695 is output even with different parameter.

What am I doing wrong?

Improve this question
2
  • 1
    Why don't you simply debug it step by step until you find the variable state in Python doesn't match that in C so you can investigate what causes that difference? Commented Jun 5, 2021 at 15:21
  • Actually, I did that, But since I haven't worked with C or C++ for about 5 or 6 years I didn't know how to debug the code through the library since all this stuff happens in the library. I was searching for how to do this but the answer was great and solved my problem. Commented Jun 5, 2021 at 16:41

1 Answer 1

Reset to default
1

In the C code, the k variable is declared as an integral type, so the following code sequence acts as a sort of fmod operation, leaving in p the remainder after the division:

    k = p / M1;  // Here, k will be TRUNCATED to the integral part of the division
    p -= k * M1; // So re-multiplying and then subtracting will leave the remainder

Thus, given initial values for p and M1 of 9.0 and 7.0, respectively, after those two lines of code, p will be 2.0.

However, in your Python code, the k variable will be a floating-point type, the division will be 'exact' and the value of p (given the same starting values for p and M1 as above) will be zero after the divide-multiply-subtract operation sequence:

p = 9.0
M1 = 7.0
k = p / M1   # k here is NOT an integer ...
p -= k * M1  # ... so this will reduce p to zero
print(p)

In fact, p will always be (very near) zero after those two lines of code, when k is of the same (real) type as p and M1 (and the problem is repeated with the k = p / M2; operation).

There are likely various ways to fix this problem (I'm no Python expert), but a simple solution is to convert the result of the division to a integer:

k = int(p / M1)
p -= k * M1

Alternatively, a perhaps more 'Pythonic' way to achieve the same result is to use the floor division operator (//):

k = p // M1 # Floor division - returns the integral part of the quotient 
p -= k * M1
Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.