 ## "DSP chips can produce random numbers using proven algorithm"

### Paul Mennen, Tektronix Inc EDN January 21, 1991, page 141

Psuedorandom number sequences prove useful for many measurement techniques. The basic idea is to excite the system under test with a signal that simultaneously has all frequencies of interest instead of at one frequency at a time as in the swept-sine method. You then compare the Fourier transforms of the system input and output to generate a frequency response.

Of the hundreds of methods in use for generating psuedorandom sequences, this listing for the WE DSP16 use the Linear Congruential Sequence" (LCS). One advantage of the LCS is that it is simple enough that mathematicians have been able to prove many of its properties.

The defining equation for the LCS is:

R' = (aR + c) mod m

where R is the random number, a and c are constants, and mod m is used to signify the operation of taking the remainder after dividing by m (modulus). First pick a starting value for R, called the seed. (i.e. assign the seed as the first value for R). Then multiply R by a constant, add a different constant to the product, divide the sum by m, and assign the remainder resulting from this division to be the next value of R. These steps are repeated every time a new random number is required.

A little thought will reveal that the numbers produced (R) are between 0 and m-1 and that the sequence can not be longer than m. (assuming that a, c, m, and the seed are positive). Now you are only left with the non-trivial task of picking appropriate values for a, c, and m. (The choice of the seed turns out to be unimportant). I chose these numbers to satisfy the following 6 criteria. (Additionally define n = sqrt(m)).

1.) c and m are relatively prime

2.) a-1 is a multiple of every prime factor of m

3.) a-1 is a multiple of 4 if m is a multiple of 4

4.) a and c are both greater than n

5.) trunc(a/n) is a power of 2

6.) a mod n is less than n/2

The first three criteria insure that the sequence will be maximal length, ie the sequence will contain all the numbers between 0 and m-1 before repeating. (This result can be proved, but it's beyond the scope of the article.) Criterion 4 insures sufficient spectral flatness and distribution uniformity. Criteria 5 and 6 will turn out to simplify the computation without compromising the sequence quality.

Note that you can simplify the computation by choosing m to be a power of 2, since the division can be implemented by shifting. The DSP16 includes a 16-bit multiplier, so it would seem that m = 2**16 would be a logical choice. Assuming a maximal length sequence and a 512-kHz output rate, the sequence will repeat every .128 seconds (65536 / 512000). Too short for most practical applications.

The next logical choice then is m = 2**32, resulting in a sequence length of 2.33 hours (65536 * .128 seconds). This is long enough to meet even the most demanding situations.

For this choice of m, you need to do a 32-bit multiply. This operation requires multiple-precision arithmetic which complicates the code considerably. Define alo and ahi to be the lower and higher 16 bits of the constant a. Similarly denote the two halves of the random number as Rlo and Rhi. Now the LCS using multiple-precision arithmetic is

EQUATION 1:

R = c + alo*Rlo + ((alo*Rhi + ahi*Rlo + zz) << 16)

where

* signifies a product of two 16-bit numbers (32-bit result)
+ signifies a sum of two 32-bit numbers
<< signifies a logical left shift operation
zz = alo if Rlo is greater than 32768
zz = 0 otherwise

To see where this comes from first look at the multiple-precision product of the 32-bit constant "a" and the 32-bit random number "R". (Fig 5 in article).

Because the remainder after dividing by 2**32 (m) is simply the low order half of the 64 bit product, we don't even need to compute the upper half. Thus we need only save the lower half of the 2nd and 3rd partials and we don't need to compute the 4th partial at all. Because of criterion #5, the 3rd partial is a simple left shift, leaving only 2 required multiplies.

By now you should be able to verify equation 1 except for the mysterious zz term. This term stems from the fact that the partial products must be computed using unsigned multiplies for the 64-bit multiple-precision product to be correct. However, most multipliers (including the DSP16's) are 2's complement signed multipliers. Using a signed multiply for the 2nd, 3rd, and 4th partial products produces errors only in the unused high-order half of the product, so only the 1st partial remains a concern. Although alo will always be interpreted by the multiplier as a positive number (because of criterion #6), Rlo can be (and will be) any combination of 16 bits. Whenever the most significant bit of Rlo is 1, we need to correct the result of the signed multiply (with the zz term) to fool the multiplier into producing an unsigned result.

The actual values I chose for constants a and c are:

a = 107465 (hex) c = 234567 (hex)

You can verify that these constants satisfy the above 6 criteria. Note however that there is a large number of equally fine choices. Using these values for a and c equation 1 in psuedo code is

TEMP = 7465 * Rhi ; 2nd partial product (TEMP is 32 bits)
TEMP = TEMP + Rlo << 4
; add 3rd partial product (ahi = 10 hex) if (Rlo > 7fff)
; if Rlo > 32767 correct the 1st partial TEMP = TEMP + 7465
; product for using a signed multiply
TEMP = TEMP << 10 ; shift left 16 (discard 16 high order bits)
TEMP = TEMP + 7465 * Rlo ; add 1st partial product
R = TEMP + 234567 ; 32 bit add produces the new random number

All constants in the psuedo code are in hexidecimal. As before, Rlo and Rhi are the lower and upper 16 bits respectively of the random number R

For those familiar with the AT&T DSP16, the code segment which computes the LCS sequence follows. The computations have been carefully ordered to maximize the number of parallel operations. Note that the routine computes a 32-bit random number, but my application needed only 16 bits to send to a DAC. I chose to send the upper 16 bits the DAC since it turns out to be "more random" than the lower half.

Note that it takes only 15 cycles to output each random number. This means that if you were to use the fastest version of the DSP16A (25 ns) and the DSP chip was doing nothing else, you could output a new random number every 375 ns.

```__________________________________________________________________________

Code sequence for computing the LCS on the DSP16
__________________________________________________________________________
.ram
ALO:	int			/* storage for a(lo) - Multiplier	*/
CHI:	int			/* storage for c(hi) - Add constant	*/
CLO:	int			/* storage for c(lo)			*/
.endram
alo:	int 0x7465		/* a(lo)				*/
chi:	int 0x0023		/* c(hi)				*/
clo:	int 0x4567		/* c(lo)				*/
i = 0			/* used as table increment		*/
a0 = i			/* initialize a0 (seed)  		*/
r0 = ALO		/* points to alo			*/
rb = ALO		/* beginning of circular memory buffer 	*/
re = CLO		/* end of circular memory buffer	*/
pt = alo		/* point to table of constants		*/
do 3 {
y=a0 x=*pt++	/* save constants alo, chi, and clo to	*/
*r0++ = x	/* memory locations ALO, CHI, and CLO  */
}
pt = alo		/* point to alo				*/
y=a0 x=*pt++i		/* x will always contain constant alo	*/

/* NOW WE ARE SET UP, AND READY TO ROLL OUT THE RANDOM NUMBERS QUICKLY */

do n {	/* Execute the loop in the instruction cache (needed for speed) */
y = a0			/* y = Rhi				*/
a1l = *r0++		/* a1l = alo (zz for Rlo >= 32768)	*/
a0 = a0 << 16		/* a0h = Rlo				*/
if pl a1=a1<<16		/* clear zz if Rlo < 32768		*/
p=x*y y=a0 x=*pt++i	/* p = alo * Rhi; y = Rlo		*/
a1 = a1 + p		/* 2nd partial product + zz		*/
a1 = a1 << 16		/* save only 16 LS bits in a1h		*/
a0 = a0 << 4		/* a0 = ahi * Rlo (3rd partial product)	*/
p=x*y y=a0 x=*pt++i	/* p = alo*Rlo; yh = ahi*Rlo = 16*Rlo	*/
a0=a1+y	 y =*r0++	/* add 3rd partial product, yh = cH	*/
a0=a0+p	 yl=*r0++	/* add 1st partial product, yl = cL	*/
pdx0 = a1		/* output high part of R		*/
}
__________________________________________________________________________

For the PIC:
In the example below it computes: Rnew = (Rold * 221 + 53) mod 256

;221 = 256 - 32 - 4 + 1
;256 can be eliminated
;so we need to calculate Rnew = Rold * (1 - 32 - 4) + 53 using
;truncating arithmetic
;or Rnew = Rold * (-32 - 3) + 53

clrc
rlf RANDOM, f
swapf RANDOM, w
andlw 0xE0
rrf RANDOM, f
sublw 0x35
movwf RANDOM

```

 file: /Techref/microchip/Linear Congruential Sequence.htm, 8KB, , updated: 2000/6/22 15:58, local time: 2023/9/27 04:06, TOP NEW HELP FIND:  44.192.115.114:LOG IN

 ©2023 These pages are served without commercial sponsorship. (No popup ads, etc...).Bandwidth abuse increases hosting cost forcing sponsorship or shutdown. This server aggressively defends against automated copying for any reason including offline viewing, duplication, etc... Please respect this requirement and DO NOT RIP THIS SITE. Questions?Please DO link to this page! Digg it! / MAKE! microchip Linear Congruential Sequence

After you find an appropriate page, you are invited to your to this massmind site! (posts will be visible only to you before review) Just type a nice message (short messages are blocked as spam) in the box and press the Post button. (HTML welcomed, but not the <A tag: Instead, use the link box to link to another page. A tutorial is available Members can login to post directly, become page editors, and be credited for their posts.

Attn spammers: All posts are reviewed before being made visible to anyone other than the poster.
 Did you find what you needed? "No. I'm looking for: " "No. Take me to the search page." "No. Take me to the top so I can drill down by catagory" "No. I'm willing to pay for help, please refer me to a qualified consultant"

### Welcome to sxlist.com!

Site supported by
& kind contributors
just like you!

(here's why

Copies of the site on CD
are available at minimal cost.

.