Hi Csounders!

I haven't been doing anything with Csound in recent years, but I lurk on
this list.  Nearly 10 years ago I contributed a 31 bit pseudo-random
number generator to Csound for noise etc.  As far as I know, it is still
there.

However, I never understood how it worked.  After a week of research, I
know exactly how it works and have written some C, C++ and dsPIC (a
snappy microcontroller chip) code to implement the same algorithm.


Most of heavy research on PRNGs (pseudo-random number generators) is
for cryptographic or simulation applications.  These generators are
generally too slow or use too much memory if all that is required is
a long sequence of random sounding numbers, such as for noise, dither
etc. in audio DSP.  Trying to find a reliable, fast, algorithm for
this purpose is tricky, since many linear congruential generators have
failings.

I have a new page:

  http://www.firstpr.com.au/dsp/rand31/

which provides C, C++ and dsPIC assembly code for a fast, 32 bit
integer, no-division, implementation of the well-known Park Miller
(1988) "minimal standard" pseudo-random number generator.  This
follows the suggestion by David G. Carta in 1990 of how to implement
this algorithm without using division.  This generator produces 31
bit numbers, between 1 and (2^31-2) = 0x7FFFFFFE.

Park and Miller rejected Carta's work in 1993 - saying it didn't
work.  But it works fine!

I trace the history of this PRNG, look at some explanations, and give
an explanation of why Carta's approach works which is different to
and much simpler than Carta's explanation.

The dsPIC subroutine uses 18 CPU cycles, including the call and
return - so it can generate a million results a second at 20 to 30
MIPS (80 to 120 MHz).  With unoptimised C code, an 800 MHz Pentium
III can generate 13 million results a second.

 - Robin    http://www.firstpr.com.au    Melbourne Australia
-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Robin Whittle wrote:

> I haven't been doing anything with Csound in recent years, but I lurk on
> this list.  Nearly 10 years ago I contributed a 31 bit pseudo-random
> number generator to Csound for noise etc.  As far as I know, it is still
> there.

Do you mean the x = (x * 16807) % 0x7FFFFFFF algorithm ? It is
actually used in some newer opcodes as well. In fact, I think using
rand() is a bad idea for a few reasons: it may be 16 bit on some systems,
it may be 32 bit or more but slow on others, different implementations
generate different pseudo-random sequences with the same seed, and
it is not thread safe (even though Csound 5 would require that).
-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Hi,

I'm using Csound for some audio signal processing work.

I'd like to use the "display" opcode to show some variable values. The 
problem I'm having is that although the display window is being drawn, 
it is being overwritten by the ftables, so that all I see at the end of 
a render is the last ftable. How can make the earlier display values 
remain visible?

Also, is it possible to draw two sets of variables into one window, say, 
in different colours?

Cheers,

Chris
-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Hi Istvan,

This is the same algorithm:

> Do you mean the x = (x * 16807) % 0x7FFFFFFF algorithm ? 

This is probably the best known linear congruential PRNG - Park and
Miller's "minimal standard" from 1988.

What was not so clear was how to implement it with integer maths and
without divisions.  David G. Carta wrote a paper on this and Park and
Miller said his approach was wrong.  But Carta was right.

The full story and some C code is at:

  http://www.firstpr.com.au/dsp/rand31/

 - Robin


-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Michael Rempel wrote:

> No division is needed. Just clear the high order bit if it is set using an
> and. Also you need to ignore the constant setting of the overflow bit. If
> you dont care about the high order bit, which is the negative sign for
> signed integers, then you can just do the multiply. I dont know what that
> does to the 16807 constant though.
> 
> given mynum as a seed, the following pseudo-assembler will give you the next
> pseudo-random number stored in the 64 bit memory location mynum. x is a 64
> bit accumulator/register.
> 
> 	load	x, mynum
> 	mult	x, 16807
> 	and	x, 0x7FFFFFFF
> 	stor	x, mynum

This is wrong. Note that the algorithm is

x = (x * 16807) % 0x7FFFFFFF

and not

x = (x * 16807) & 0x7FFFFFFF

Here is a fast (well, at least on x86) version from oscbnk.c:

static CS_PURE long oscbnk_rand31(long seed)
{
     uint64_t tmp1;
     uint32_t tmp2;

     /* x = (16807 * x) % 0x7FFFFFFF */
     tmp1 = (uint64_t) ((int32_t) seed * (int64_t) 16807);
     tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;
     tmp2 += (uint32_t) (tmp1 >> 31);
     if ((int32_t) tmp2 < (int32_t) 0)
       tmp2 -= (uint32_t) 0x7FFFFFFF;
     return (long) tmp2;
}

'seed' must be in the range 1 to 0x7FFFFFFE.
-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Still wrong. Note that the mod is not by 0x80000000 (which could really
be replaced with a simple AND), but one less and that is not so simple.
I know the tricks you describe, but you should read more carefully what
is to be implemented, as 2^31 and 2^31 - 1 are not the same.

Michael Rempel wrote:

> mod gives the same result as and. The division is perfect binary. tt is the
> reason 0x7FFFFFFF is choosen.
> 
> You can do the same thing for any power of 2 division. If you want the
> remainder of division by 16 for example, you can and your var with 15 to get
> the result.
> 
> in general for integer math y mod x for any x that is a power of 2 can be
> expressed as
> 
> y and (x-1)
> 
> Up to the overflow point anyway. This trick is why powers of 2 can be used
> to make cool tools.
> 
> Just for kicks, another one I like is the sum of integers from 1 to n. You
> can do a loop, or you can just do (n+1)*n/2 which is way faster. Proof for
> the equation, take any n you like. Write out the list from 1 to n. Then
> directly under that write out the list from n to 1. Now add each column.
> They are all n+1. n(n+1) is the total sum, divide by 2 to eliminate the
> duplicate list.
-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Michael Rempel wrote:

> In the line
>      tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;
> 
> what does the & operator do? In C it is an address operator, which makes no
> sense.

In C it is bitwise AND as well, but note the other lines that correct
the result so that it will be the same as if % was used, assuming
that the input value is in the valid range.
-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Nope. It is the bitwise AND in both languages (C++ supports all of C as a 
subset). && is logical AND.

& is also the address-of operator in both languages; as with the star, the 
context makes it clear what is intended.

Richard Dobson





Michael Rempel wrote:
> In the line
>      tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;
> 
> what does the & operator do? In C it is an address operator, which makes no
> sense.
> 
> In C++ it is logical and.
> 
> Michael

-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

No division is needed. Just clear the high order bit if it is set using an
and. Also you need to ignore the constant setting of the overflow bit. If
you dont care about the high order bit, which is the negative sign for
signed integers, then you can just do the multiply. I dont know what that
does to the 16807 constant though.

given mynum as a seed, the following pseudo-assembler will give you the next
pseudo-random number stored in the 64 bit memory location mynum. x is a 64
bit accumulator/register.

	load	x, mynum
	mult	x, 16807
	and	x, 0x7FFFFFFF
	stor	x, mynum

Michael Rempel

-----Original Message-----
From: Robin Whittle [mailto:rw@firstpr.com.au]
Sent: Thursday, September 22, 2005 4:29 AM
To: csound@lists.bath.ac.uk
Subject: Re: [Csnd] 31 bit pseudo-random number gen in C, C++ & dsPIC
assembly code


Hi Istvan,

This is the same algorithm:

> Do you mean the x = (x * 16807) % 0x7FFFFFFF algorithm ?

This is probably the best known linear congruential PRNG - Park and
Miller's "minimal standard" from 1988.

What was not so clear was how to implement it with integer maths and
without divisions.  David G. Carta wrote a paper on this and Park and
Miller said his approach was wrong.  But Carta was right.

The full story and some C code is at:

  http://www.firstpr.com.au/dsp/rand31/

 - Robin


--
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

mod gives the same result as and. The division is perfect binary. tt is the
reason 0x7FFFFFFF is choosen.

You can do the same thing for any power of 2 division. If you want the
remainder of division by 16 for example, you can and your var with 15 to get
the result.

in general for integer math y mod x for any x that is a power of 2 can be
expressed as

y and (x-1)

Up to the overflow point anyway. This trick is why powers of 2 can be used
to make cool tools.

Just for kicks, another one I like is the sum of integers from 1 to n. You
can do a loop, or you can just do (n+1)*n/2 which is way faster. Proof for
the equation, take any n you like. Write out the list from 1 to n. Then
directly under that write out the list from n to 1. Now add each column.
They are all n+1. n(n+1) is the total sum, divide by 2 to eliminate the
duplicate list.

Cheers,

Michael


-----Original Message-----
From: Istvan Varga [mailto:istvan_v@fibermail.hu]
Sent: Thursday, September 22, 2005 7:49 AM
To: csound@lists.bath.ac.uk
Subject: Re: [Csnd] 31 bit pseudo-random number gen in C, C++ & dsPIC
assembly code


Michael Rempel wrote:

> No division is needed. Just clear the high order bit if it is set using an
> and. Also you need to ignore the constant setting of the overflow bit. If
> you dont care about the high order bit, which is the negative sign for
> signed integers, then you can just do the multiply. I dont know what that
> does to the 16807 constant though.
>
> given mynum as a seed, the following pseudo-assembler will give you the
next
> pseudo-random number stored in the 64 bit memory location mynum. x is a 64
> bit accumulator/register.
>
> 	load	x, mynum
> 	mult	x, 16807
> 	and	x, 0x7FFFFFFF
> 	stor	x, mynum

This is wrong. Note that the algorithm is

x = (x * 16807) % 0x7FFFFFFF

and not

x = (x * 16807) & 0x7FFFFFFF

Here is a fast (well, at least on x86) version from oscbnk.c:

static CS_PURE long oscbnk_rand31(long seed)
{
     uint64_t tmp1;
     uint32_t tmp2;

     /* x = (16807 * x) % 0x7FFFFFFF */
     tmp1 = (uint64_t) ((int32_t) seed * (int64_t) 16807);
     tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;
     tmp2 += (uint32_t) (tmp1 >> 31);
     if ((int32_t) tmp2 < (int32_t) 0)
       tmp2 -= (uint32_t) 0x7FFFFFFF;
     return (long) tmp2;
}

'seed' must be in the range 1 to 0x7FFFFFFE.
--
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

In the line
     tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;

what does the & operator do? In C it is an address operator, which makes no
sense.

In C++ it is logical and.

Michael

-----Original Message-----
From: Michael Rempel [mailto:michael_rempel@shaw.ca]
Sent: Thursday, September 22, 2005 10:11 AM
To: csound@lists.bath.ac.uk
Subject: Re: [Csnd] 31 bit pseudo-random number gen in C, C++ & dsPIC
assembly code


mod gives the same result as and. The division is perfect binary. tt is the
reason 0x7FFFFFFF is choosen.

You can do the same thing for any power of 2 division. If you want the
remainder of division by 16 for example, you can and your var with 15 to get
the result.

in general for integer math y mod x for any x that is a power of 2 can be
expressed as

y and (x-1)

Up to the overflow point anyway. This trick is why powers of 2 can be used
to make cool tools.

Just for kicks, another one I like is the sum of integers from 1 to n. You
can do a loop, or you can just do (n+1)*n/2 which is way faster. Proof for
the equation, take any n you like. Write out the list from 1 to n. Then
directly under that write out the list from n to 1. Now add each column.
They are all n+1. n(n+1) is the total sum, divide by 2 to eliminate the
duplicate list.

Cheers,

Michael


-----Original Message-----
From: Istvan Varga [mailto:istvan_v@fibermail.hu]
Sent: Thursday, September 22, 2005 7:49 AM
To: csound@lists.bath.ac.uk
Subject: Re: [Csnd] 31 bit pseudo-random number gen in C, C++ & dsPIC
assembly code


Michael Rempel wrote:

> No division is needed. Just clear the high order bit if it is set using an
> and. Also you need to ignore the constant setting of the overflow bit. If
> you dont care about the high order bit, which is the negative sign for
> signed integers, then you can just do the multiply. I dont know what that
> does to the 16807 constant though.
>
> given mynum as a seed, the following pseudo-assembler will give you the
next
> pseudo-random number stored in the 64 bit memory location mynum. x is a 64
> bit accumulator/register.
>
> 	load	x, mynum
> 	mult	x, 16807
> 	and	x, 0x7FFFFFFF
> 	stor	x, mynum

This is wrong. Note that the algorithm is

x = (x * 16807) % 0x7FFFFFFF

and not

x = (x * 16807) & 0x7FFFFFFF

Here is a fast (well, at least on x86) version from oscbnk.c:

static CS_PURE long oscbnk_rand31(long seed)
{
     uint64_t tmp1;
     uint32_t tmp2;

     /* x = (16807 * x) % 0x7FFFFFFF */
     tmp1 = (uint64_t) ((int32_t) seed * (int64_t) 16807);
     tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;
     tmp2 += (uint32_t) (tmp1 >> 31);
     if ((int32_t) tmp2 < (int32_t) 0)
       tmp2 -= (uint32_t) 0x7FFFFFFF;
     return (long) tmp2;
}

'seed' must be in the range 1 to 0x7FFFFFFE.
--
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

--
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Yes, now we understand each other. I was simplifying.

We know that x mod y is approximately x mod (y-1) for very big y values,
the messing with bits after the & fixes the difference.

cheers,

Michael

-----Original Message-----
From: Istvan Varga [mailto:istvan_v@fibermail.hu]
Sent: Thursday, September 22, 2005 8:38 AM
To: csound@lists.bath.ac.uk
Subject: Re: [Csnd] 31 bit pseudo-random number gen in C, C++ & dsPIC
assembly code


Michael Rempel wrote:

> In the line
>      tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;
>
> what does the & operator do? In C it is an address operator, which makes
no
> sense.

In C it is bitwise AND as well, but note the other lines that correct
the result so that it will be the same as if % was used, assuming
that the input value is in the valid range.
--
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

Thanks. It has been a while since I used c/c++. I use Delphi most of the
time.

cheers,

Michael

-----Original Message-----
From: Richard Dobson [mailto:richarddobson@blueyonder.co.uk]
Sent: Thursday, September 22, 2005 8:53 AM
To: csound@lists.bath.ac.uk
Subject: Re: [Csnd] 31 bit pseudo-random number gen in C, C++ & dsPIC
assembly code


Nope. It is the bitwise AND in both languages (C++ supports all of C as a
subset). && is logical AND.

& is also the address-of operator in both languages; as with the star, the
context makes it clear what is intended.

Richard Dobson





Michael Rempel wrote:
> In the line
>      tmp2 = (uint32_t) tmp1 & (uint32_t) 0x7FFFFFFF;
>
> what does the & operator do? In C it is an address operator, which makes
no
> sense.
>
> In C++ it is logical and.
>
> Michael

--
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk

what version of csound are you using?

[ i use display opcode a lot in MacCsound ]

-m

On Thu, 22 Sep 2005, Chris Share wrote:

> Hi,
>
> I'm using Csound for some audio signal processing work.
>
> I'd like to use the "display" opcode to show some variable values. The 
> problem I'm having is that although the display window is being drawn, it is 
> being overwritten by the ftables, so that all I see at the end of a render is 
> the last ftable. How can make the earlier display values remain visible?
>
> Also, is it possible to draw two sets of variables into one window, say, in 
> different colours?
>
> Cheers,
>
> Chris
> -- 
> Send bugs reports to this list.
> To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk
>
-- 
Send bugs reports to this list.
To unsubscribe, send email to csound-unsubscribe@lists.bath.ac.uk