From: M DOT A DOT Bukin AT inp DOT nsk DOT su
To: djgpp AT delorie DOT com
Subject: Re: An inline assembler RNG for C
References: <36CB0BF1 DOT 23FF9E2 AT stat DOT fsu DOT edu>
Date: 18 Feb 1999 19:18:05 +0600
In-Reply-To: George Marsaglia's message of "Wed, 17 Feb 1999 13:35:29 -0500"
Message-ID: <20yalvsvea.fsf@Sky.inp.nsk.su>
Lines: 105
X-Mailer: Gnus v5.5/Emacs 19.34
Reply-To: djgpp AT delorie DOT com

George Marsaglia <geo AT stat DOT fsu DOT edu> writes:

> The method is Multiply-with-Carry for 32-bit integers.
> This is the idea.  Suppose you have
>           a current 32-bit integer x
>           a current 32-bit carry   c
>           a fixed multiplier, say, a=2083801278
> Now form  a*x+c in 64 bits.
> Return the bottom 32 bits as the new x,
> and keep the top 32 bits as the new carry, c.

Can it be proved that there is no situation when x(i+1)=x(i) and
c(i+1)=c(i), and generated sequence will not degrade to x(i)?

> If you are familiar with how inline assembler works
> for C compilers you use, and can provide a version of
> the above for those compilers, a posting might be of
> interest to many potential users.

It is of interest to me, and I tried to make one

<prng.c>
static unsigned long x = 1;
static unsigned long c = 0;

#if !defined(__GNUC__)
static void
addto (unsigned long *rh, unsigned long *rl,
       const unsigned long ah, const unsigned long al)
{
  *rl = (*rl + al) & 0xFFFFFFFFUL;
  *rh = (*rh + ah + ((*rl < al) ? 1 : 0)) & 0xFFFFFFFFUL;
}
#endif

unsigned long
prng (void)
{
  const unsigned long a = 2083801278UL;

#if defined(__GNUC__)
#if defined(__i386__)
  __asm__ __volatile__
    ("mull %%edx\n\t"
     "addl %4, %%eax\n\t"
     "adcl $0, %%edx\n\t"
     : "=a" (x), "=d" (c)
     : "0" (x), "1" (a), "r" (c));

  return x;
#else
  unsigned long long tmp;

  tmp = (unsigned long long) x * a + c;
  c = (tmp >> 32) & 0xFFFFFFFFUL;
  x = (tmp & 0xFFFFFFFFUL);

  return x;
#endif
#else
  const unsigned long xl = (x & 0xFFFFUL);
  const unsigned long xh = (x >> 16) & 0xFFFFUL;
  const unsigned long al = (a & 0xFFFFUL);
  const unsigned long ah = (a >> 16) & 0xFFFFUL;

  unsigned long tmp;

  x = c;
  c = 0;
  tmp = xl * al;
  addto (&c, &x, 0, tmp);
  tmp = xl * ah;
  addto (&c, &x, (tmp >> 16), (tmp & 0xFFFFUL) << 16);
  tmp = xh * ah;
  addto (&c, &x, tmp, 0);
  tmp = xh * al;
  addto (&c, &x, (tmp >> 16), (tmp & 0xFFFFUL) << 16);

  return x;
#endif
}

void
seed_prng (unsigned long seed)
{
  x = 1;
  c = seed;
}
</prng.c>

With full optimizations (Pentium-gcc with `-O6') code generated for
long long version and inline assembly is the same.

Maybe there are better ways to make portable implementation of this
RNG, using only ANSI-C.  If const are gcc-isms, then they should be
removed.

The choice for `seed_prng' -- seed is assigned to carry and initial x
is 1, which makes it possible to obtain multiplier constant by seeding
RNG with 0.

If I made any mistakes, please, point them out.

-- 
Michael Bukin