Python code from C code does not give the same results

Question

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?

Answer 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