6

I am trying to implement a random-number generator with Mersenne prime (231-1) as the modulus. The following working code was based on several related posts:

  1. How do I extract specific 'n' bits of a 32-bit unsigned integer in C?
  2. Fast multiplication and subtraction modulo a prime
  3. Fast multiplication modulo 2^16 + 1

However,

It does not work with uint32_t hi, lo;, which means I do not understand signed vs. unsigned aspect of the problem.

Based on #2 above, I was expecting the answer to be (hi+lo). Which means, I do not understand why the following statement is needed.

 if (x1 > r) x1 += r + 2;

  • Can someone please clarify the source of my confusion?

  • Can the code itself be improved?

  • Should the generator avoid 0 or 231-1 as a seed?

  • How would the code change for a prime (2p-k)?

Original code

#include <inttypes.h> // x1 = a*x0 (mod 2^31-1) int32_t lgc_m(int32_t a, int32_t x) { printf("x %"PRId32"\n", x); if (x == 2147483647){ printf("x1 %"PRId64"\n", 0); return (0); } uint64_t c, r = 1; c = (uint64_t)a * (uint64_t)x; if (c < 2147483647){ printf("x1 %"PRId64"\n", c); return (c); } int32_t hi=0, lo=0; int i, p = 31;//2^31-1 for (i = 1; i < p; ++i){ r |= 1 << i; } lo = (c & r) ; hi = (c & ~r) >> p; uint64_t x1 = (uint64_t ) (hi + lo); // NOT SURE ABOUT THE NEXT STATEMENT if (x1 > r) x1 += r + 2; printf("c %"PRId64"\n", c); printf("r %"PRId64"\n", r); printf("\tlo %"PRId32"\n", lo); printf("\thi %"PRId32"\n", hi); printf("x1 %"PRId64"\n", x1); printf("\n" ); return((int32_t) x1); } int main(void) { int32_t r; r = lgc_m(1583458089, 1); r = lgc_m(1583458089, 2000000000); r = lgc_m(1583458089, 2147483646); r = lgc_m(1583458089, 2147483647); return(0); }
Improve this question
7
  • Look at the assembler output of the sgined operations. You will generally see Logical Shifts to the Left, but Arithmetic Shifts to the Right with signed values. Other bitwise operators should behave the same. The ASL operator will preserve the sign bits when shifting downwards. Commented Jun 24, 2015 at 22:51
  • Thanks. The first attempt was signed and therefore added to the confusion. The second solution used unsigned and revealed the real problem. Commented Jun 24, 2015 at 23:50
  • uint64_t x1 = (uint64_t ) ((hi + lo) ); Hmmm maybe uint64_t x1 = (uint64_t ) hi + lo;? Commented Jun 25, 2015 at 3:40
  • Minor: r |= 1 << i; --> r |= (uint32_t)1 << i; (Cope with 16-bit int. and signed overflow.) Commented Jun 25, 2015 at 3:43
  • 1
    It would be nicer to remove the solution from the question, post it as an answer (what it is !) and accept it. Currently, the question appears as a not answered one. Commented Jun 25, 2015 at 5:40

2 Answers 2

Reset to default
2

The following if statement

if (x1 > r) x1 += r + 2;

should be written as

if (x1 > r) x1 -= r;

Both results are the same modulo 2^31:

x1 + r + 2 = x1 + 2^31 - 1 + 2 = x1 + 2^31 + 1 x1 - r = x1 - (2^31 - 1) = x1 - 2^31 + 1

The first solution overflows an int32_t and assumes that conversion from uint64_t to int32_t is modulo 2^31. While many C compilers handle the conversion this way, this is not mandated by the C standard. The actual result is implementation-defined.

The second solution avoids the overflow and works with both int32_t and uint32_t.

You can also use an integer constant for r:

uint64_t r = 0x7FFFFFFF; // 2^31 - 1

Or simply

uint64_t r = INT32_MAX;

EDIT: For primes of the form 2^p-k, you have to use masks with p bits and calculate the result with

uint32_t x1 = (k * hi + lo) % ((1 << p) - k)

If k * hi + lo can overflow a uint32_t (that is (k + 1) * (2^p - 1) >= 2^32), you have to use 64-bit arithmetic:

uint32_t x1 = ((uint64_t)a * x) % ((1 << p) - k)

Depending on the platform, the latter might be faster anyway.

Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

Thanks. Even though, the code 'worked' I did not really understand what was going on. Your answer addresses all points of confusion.
1

Sue provided this as a solution:

With some experimentation (new code at the bottom), I was able to use uint32_t, which further suggests that I do not understand how the signed integers work with bit operations.

The following code uses uint32_t for input as well as hi and lo.

 #include <inttypes.h> // x1 = a*x0 (mod 2^31-1) uint32_t lgc_m(uint32_t a, uint32_t x) { printf("x %"PRId32"\n", x); if (x == 2147483647){ printf("x1 %"PRId64"\n", 0); return (0); } uint64_t c, r = 1; c = (uint64_t)a * (uint64_t)x; if (c < 2147483647){ printf("x1 %"PRId64"\n", c); return (c); } uint32_t hi=0, lo=0; int i, p = 31;//2^31-1 for (i = 1; i < p; ++i){ r |= 1 << i; } hi = c >> p; lo = (c & r) ; uint64_t x1 = (uint64_t ) ((hi + lo) ); // NOT SURE ABOUT THE NEXT STATEMENT if (x1 > r){ printf("x1 - r = %"PRId64"\n", x1- r); x1 -= r; } printf("c %"PRId64"\n", c); printf("r %"PRId64"\n", r); printf("\tlo %"PRId32"\n", lo); printf("\thi %"PRId32"\n", hi); printf("x1 %"PRId64"\n", x1); printf("\n" ); return((uint32_t) x1); } int main(void) { uint32_t r; r = lgc_m(1583458089, 1583458089); r = lgc_m(1583458089, 2147483645); return(0); }

The issue was that my assumption that the reduction will be complete after one pass. If (x > 231-1), then by definition the reduction has not occurred and a second pass is necessary. Subtracting 231-1, in that case does the trick. In the second attempt above, r = 2^31-1 and is therefore the modulus. x -= r achieves the final reduction.

Perhaps someone with expertise in random numbers or modular reduction could explain it better.

Cleaned function without printf()s.

uint32_t lgc_m(uint32_t a, uint32_t x){ uint64_t c, x1, m = 2147483647; //modulus: m = 2^31-1 if (x == m) return (0); c = (uint64_t)a * (uint64_t)x; if (c < m)//no reduction necessary return (c); uint32_t hi, lo, p = 31;//2^p-1, p = 31 hi = c >> p; lo = c & m; x1 = (uint64_t)(hi + lo); if (x1 > m){//one more pass needed //this block can be replaced by x1 -= m; hi = x1 >> p; lo = (x1 & m); x1 = (uint64_t)(hi + lo); } return((uint32_t) x1); }
Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.