diff --git a/libsrc/common/rand.s b/libsrc/common/rand.s
index a272af80a..798d6343a 100644
--- a/libsrc/common/rand.s
+++ b/libsrc/common/rand.s
@@ -15,14 +15,14 @@
 ;  Multiplier must be 1 (mod 4)
 ;  Added value must be 1 (mod 2)
 ;  This guarantees max. period (2**32)
-;  The lowest bits have poor entropy and
-;  exhibit easily detectable patterns, so
-;  only the upper bits 16-22 and 24-31 of the
-;  4-byte state are returned.
+;  The quality of entropy in the bits of the seed are poorest in the lowest
+;  bits, and best in the highest bits.
 ;
-;  The best 8 bits, 24-31 are returned in the
-;  low byte A to provide the best entropy in the
-;  most commonly used part of the return value.
+;  The high 8 bits are used for the low byte A to provide the best entropy in
+;  the most commonly used part of the return value.
+;
+;  Finally XOR with the lower 2 bytes is used on the output, which breaks up
+;  some minor deficient sequential patterns. (#951)
 ;
 ;  Uses the following LCG values for ax + c (mod m)
 ;  a = $01010101
@@ -42,10 +42,15 @@
 
 ; The seed. When srand() is not called, the C standard says that that rand()
 ; should behave as if srand() was called with an argument of 1 before.
-rand:   .dword   1
+rand:   .dword $B5B5B4B4
 
 .code
 
+_srand: sta     rand+0          ; Store the seed
+        stx     rand+1
+        sta     rand+2          ; argument << 16 is convenient fill for MSW
+        stx     rand+3
+        ; fall through to rand() to sufficiently "shuffle" first rand() result
 _rand:  clc
         lda     rand+0
         adc     #$B3
@@ -54,18 +59,11 @@ _rand:  clc
         sta     rand+1
         adc     rand+2
         sta     rand+2
+        eor     rand+0
         and     #$7f            ; Suppress sign bit (make it positive)
         tax
         lda     rand+2
         adc     rand+3
         sta     rand+3
+        eor     rand+1
         rts                     ; return bit (16-22,24-31) in (X,A)
-
-_srand: sta     rand+0          ; Store the seed
-        stx     rand+1
-        lda     #0
-        sta     rand+2          ; Set MSW to zero
-        sta     rand+3
-        rts
-
-
diff --git a/test/val/rand.c b/test/val/rand.c
new file mode 100644
index 000000000..f2f604449
--- /dev/null
+++ b/test/val/rand.c
@@ -0,0 +1,110 @@
+/* This test verifies that the assembly implementation of rand() matches its
+ * theoretical high level equivalent.
+ *
+ * This does about 3000 tests from various starting srand() seeds.
+ * A more thorough test might visit the entire sequence with 2^32 tests, but
+ * that takes hours to simulate, and this should be a sufficient sampling.
+ *
+ * This will also fail if rand() is ever altered, which might be a warning to
+ * tread carefully. Some past discussion of RNG here:
+ * https://github.com/cc65/cc65/pull/951
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdint.h>
+
+/* for faster execution */
+#pragma static-locals (on)
+
+/* values tested per seed */
+#define SUBTESTS 50
+
+/* increments used between tested seeds */
+/* 653 is prime and divides 32768 by ~50 */
+#define TESTINC 653
+
+static uint32_t seed;
+
+int ref_rand()
+{
+	uint16_t output;
+	/* seed follows the LCG sequence * 0x01010101 + 0xB3B3B3B3 */
+	seed = seed * 0x01010101UL + 0xB3B3B3B3UL;
+	/* output uses the top two bytes (reversed) XOR with bottom two bytes */
+	{
+		uint16_t s0 = (seed >>  0) & 0xFF;
+		uint16_t s1 = (seed >>  8) & 0xFF;
+		uint16_t s2 = (seed >> 16) & 0xFF;
+		uint16_t s3 = (seed >> 24) & 0xFF;
+		uint16_t o0 = s3 ^ s1;
+		uint16_t o1 = s2 ^ s0;
+		output = o0 | (o1 << 8);
+	}
+	return (int)(output & 0x7FFF);
+}
+
+void ref_srand(int ax)
+{
+	uint32_t s = (unsigned int)ax;
+	seed = s | (s << 16); /* low 16 bits is convenient filler for high 16 bits */
+	ref_rand(); /* one pre-call "shuffles" the first rand() result so it isn't too predictable */
+}
+
+int main(void)
+{
+	unsigned int i,j;
+	int a,b;
+
+	/* test that startup state is equivalent to srand(1) */
+	{
+		//srand(1); // implied
+		ref_srand(1);
+		for (j=0; j<SUBTESTS; ++j)
+		{
+			a = rand();
+			b = ref_rand();
+			if (a != b)
+			{
+				printf("failed startup seed at test %d. rand()=%d reference=%d\n",j,a,b);
+				return EXIT_FAILURE;
+			}
+		}
+	}
+
+	/* test every power of 2 seed */
+	for (i = 0; i < 16; ++i)
+	{
+		srand(1<<i);
+		ref_srand(1<<i);
+		for (j=0; j<SUBTESTS; ++j)
+		{
+			a = rand();
+			b = ref_rand();
+			if (a != b)
+			{
+				printf("failed seed %d at test %d. rand()=%d reference=%d\n",(1<<i),j,a,b);
+				return EXIT_FAILURE;
+			}
+		}
+	}
+
+	/* test a sampling of seeds*/
+	for (i = 0; i < 32768UL; i += TESTINC)
+	{
+		srand(i);
+		ref_srand(i);
+		for (j=0; j<SUBTESTS; ++j)
+		{
+			a = rand();
+			b = ref_rand();
+			if (a != b)
+			{
+				printf("failed seed %d at test %d. rand()=%d reference=%d\n",(1<<i),j,a,b);
+				return EXIT_FAILURE;
+			}
+		}
+	}
+
+	return EXIT_SUCCESS;
+}