Talk:Multiply-with-carry pseudorandom number generator

OR

Somebody has tagged this page as "original research" and "conflict of interests".

The Wikipedia policy on citing oneself does not prevent a scientist from citing his own work. In fact, this may be the only way to get articles on highly specialized fields.

I do agree, however, that there was a bias in the original article in the sense that published criticism of the MWC generator was omitted. I have added a reference to the article by Couture and l'Ecuyer to counter this bias.

We need a general discussion of the pros and cons of different RNGs. I am not qualified to write such an article. Therefore, I have replaced the tags with the "need expert attention" tag.

I think that the MWC generator has certain advantages in terms of high bifurcation, but I have no reliable reference for this claim.

Arnold90 13:39, 20 June 2007 (UTC)

It's a bad title

A piano uses strings and information about constructing and playing pianos is valuable, but unwanted on a page titled string. Likewise Multiply-with-carry in isolation is too general a title for this specific application of a common mathematical function that is widely found in FFTs, digital filters and general-purpose microprocessors. Cuddlyable3 (talk) 14:17, 27 August 2008 (UTC) perhaps "Multiply-with-carry(PRNG)" would be a better title? 118.90.129.46 (talk) 08:46, 18 September 2010 (UTC)

• Strongly agree - there needs to be "PRNG" in the title somewhere, as that's what this is, it's a PRNG — Preceding unsigned comment added by Fatphil (talk • contribs) 11:02, 19 December 2011 (UTC)
• Agreed, "Multiply-with-carry" absolutely doesn't look right for this article. I actually landed on the page searching for MWC8222; there certainly should be a redirect for that abbreviation (if not the primary title). --77.91.170.13 (talk) 11:50, 17 December 2015 (UTC)

mod 2^32-1 operation

In article one mentions that modulo for b=2^32-1, is slow operation. It is true. But article mentions something about speeding it be converting i somehow to b+1 (2^32). How this is achived. I don't see any simple way for this. —Preceding unsigned comment added by 149.156.82.207 (talk) 20:34, 13 January 2010 (UTC)

Not sure if this is the right place for such a discussion, but I believe you'll find clarification at http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm#Shift_optimizations (unsigned)

That optimisation does mod 2^n+1.

For mod(2^n-1) use addition instead of subtraction. This is analogous to mod operation performed by Casting out nines but done in arithmetic base 2^n instead of base 10 thus yielding a modulus base 2^n-1 instead of base 9. Jasen betts (talk) 09:06, 18 September 2010 (UTC)

It also says using mod ((2^ n) -1) is superior in "avoiding nagging difficulties", but makes no claim as to what they might be or why this is so. — Preceding unsigned comment added by 68.203.10.112 (talk) 23:13, 15 August 2013 (UTC)

Linear?

Is it linear? That would be useful for the article to state. Cmcqueen1975 (talk) 11:19, 2 March 2011 (UTC)

incorrect formula

it says...

 If we choose an a of 28 to 31 bits such that ab−1 is a safe prime, that is both ab − 1 and ab/2 − 1 are prime, 

but the definition of a safe prime disagrees with this statement.

99.155.186.249 (talk) 15:49, 28 March 2011 (UTC)

• Set p = ab/2-1. Then 2*p+1 = 2*(ab/2-1)+1 = (ab-2)+1 = ab-1. The formula is correct. p=ab/2-1 is the Sophie Germain prime, 2p+1=ab-1 is the safe prime. Fatphil (talk) 10:57, 19 December 2011 (UTC)

Implementation example

I think the used example is not a good one. The article explain how to use base b=2^32, but this example uses b=(2^32)-1. Marsaglia's article "On the Randomness of Pi and Other Decimal Expansions" shows some examples with b=2^32 and larger periods. With b=2^32-1, the algorithm can't return 32 unbiased bits, so it's necessary take 16 bits from two sucessive output to return an integer from 0 to 2^32-1 (maybe wikipedia article should mention it). In my oppinion, a better example, from Marsaglia's article, is b=2^32, a=2169967, r=25000. It has period > 10^240833 > 2^800029. I can put the code here, but I do not know how to write this on wikipedia, and my English is not very good. — Preceding unsigned comment added by Edua1978 (talk • contribs) 02:42, 8 January 2012 (UTC)

Updated version of code (is it, the original, or neither correct?)

The example code seems to be a from George's 2003 comp lang c posting, which provided the following code

static unsigned long Q[4096], c = 362436;

unsigned long CMWC4096(void) {
        unsigned long long t, a = 18782LL;
        static unsigned long i = 4095;
        unsigned long x, r = 0xfffffffe;

        i = (i + 1) & 4095;
        t = a * Q[i] + c;
        c = (t >> 32);
        x = t + c;
        if (x < c) {
                x++;
                c++;
        }
        return (Q[i] = r - x);
}

He also posted another version in the sci crypt group in 2005

static unsigned long Q[4096], c = 362436;

unsigned long CMWC4096(void){
        unsigned long long t, a = 18782LL, b = 4294967295LL;
        static unsigned long i = 4095;
        unsigned long x, r = b - 1;

        i = (i + 1) & 4095;
        t = a * Q[i] + c;
        c = (t >> 32);
        t = (t & b) + c;
        if (t > r) {
                c++;
                t = t - b;
        }
        return (Q[i] = r - t);
}

Clean up the code a bit in the style of the first one, we get

static unsigned long Q[4096], c = 362436;

unsigned long CMWC4096(void){
        unsigned long long t, a = 18782LL;
        static unsigned long i = 4095;
        unsigned long x, r = 0xfffffffe;

        i = (i + 1) & 4095;
        t = a * Q[i] + c;
        c = (t >> 32);
        x = (t & 0xffffffff) + c;
        if (x > r) {
                x = x - 0xffffffff;
                c++;
        }
        return (Q[i] = r - x);
}

This begs the question of why has the last little bit changed for the earlier version? Was it incorrect? Is the if statement suppose to be making sure $r-x>=0$? Does this mean the routine returns uniform values on $[0,2^{32}-2]$? It all seems a bit odd. — Preceding unsigned comment added by 129.100.171.100 (talk) 20:27, 13 January 2012 (UTC)

I would add that based on the code and math on the page, I believe the intent was to implement a complementary multiply with carry generator with base $b=2^{32}-1$ (so it does indeed return uniform numbers on $[0,2^{32}-2]$). The following should then be correct

uint32_t Q[0x1000], c;

uint32_t CMWC4096(void){
        const uint32_t a = 18782;
        static uint i = 0x0fff;
        uint64_t t;

        i = (i + 1) & 0x0fff;
        t = a * Q[i] + c;
        c = t >> 32;
        x = t & 0xffffffff;
        if (x > 0xfffffffe) {
                x = x - 0xffffffff;
                c++;
        }
        return (Q[i] = 0xfffffffe - x);
}

The if statement is an adjustment for the base being $2^{32}-1$ and not $2^{32}$.

Initialization code for $Q$ should also ensure that each element is in the range $[0,2^{32}-2]$. I don't believe the sample code does this for $Q$. Initialization code for $c$ should ensure it is in the range $[0,a-1]$ (as per both figuring out what the maximum possible carry value is and George's own comments in his emails to the newsgroups group). Neither the sample code nor the two codes he provides do this.

References

Licensing and use in GPL'd program

Hello, I'm working on a program that requires a fast PRNG that is simple to use. I found this one but I haven't found any information regarding licensing. Who would I ask about using the implementation provided on the main page or on this talk page? 99.246.88.56 (talk) 17:12, 23 February 2012 (UTC)

lag-r for r>1

I doubt this is accurate:

(but apparently only for lag-1 MWC sequences)

Each round in the algorithm represents one step in a long multiplication. If you consider r rounds, then r inputs (say x₀ to x_r-1) have been read in as a single value in the range [0,b^r) and multiplied by a, with the result over b^r given in c and the remainder in x_r to x_2r-1. It's the same algorithm as if b had actually been b^r at the outset and r was 1.

If some analysis has suggested that r>1 doesn't work, it's probably through failing to regard sets of r outputs as components of a single integer. --217.140.96.21 (talk) 12:01, 27 June 2012 (UTC)

Python alternative

I spent a couple hours trying to get both the original and the Talk-updated versions running in TinyC. Both failed. For some reason, both t and c stayed zero in successive passes through the generator. This could be due to either TinyC or the fact that I was running on a (Intel 32-bit) CPU, which the main page seems to think is OK.

Here is a Python translation which seems to work. (I make no clams for the correctness of the random algorithm itself; I simply implemented what I think is the algorithm as presented.)

# Python implementation of rand_cmwc from MW: Multiply-with-carry
phi = 0x9e3779b9
q = []
for dummy in range(4096): q.append(0)
ones16 = 0x0ffff
ones32 = 0xffffffff
two32 = ones32 + 1
def init_rand(seed):
    q[0] = seed
    q[1] = (seed + phi) & ones32
    q[2] = (seed + phi + phi) & ones32
    for j in range(3,4096):
        q[j] = q[j-3] ^ q[j-2] ^ phi ^ i
        # The "^" operator is bitwise exclusive or.
global i, c
i = ones16        
c = 362436
def randcmwc():
    global i, c
    a = 18782
    r = ones32 - 1
    i = (i+1) & ones16
    t = a * q[i] + c
    c = t // ones32
    x = t & ones32
    if x > r:
        x -= ones32
        c += 1   
    q[i] = (r - x) & ones32
    return q[i]
if __name__ == "__main__":
    init_rand(12345)
    for j in range(2):
        random = randcmwc()
        print "%2d %10d" % (j+1, random)

I used & rather than % in the remaindering operations, because I think it's faster, though that's academic in an interpreted compiler. ;-) Docduke (talk) 02:59, 16 July 2013 (UTC)

Quality of generator, performance in test suites

dear Wikipedians,
This article has no clear statement of the quality of the MWC generator.

1. How does it perform in TestU01 and other test suites?
2. How well does it comply with US, UK, German, and other standards ?

Who is willing to provide some answers to these questions ? -- jw (talk) 20:53, 24 February 2020 (UTC)