Re: Individually, each of the 31 bits are equally random.

You don't have to use the 16-bit rand() routine.
Below is some VC++ 8.0 code to generate 31-bit random numbers.

I received the original code came from a Computer Science professor
at the University of Washington... make it 63-bit, if you want.

#pragma warning( disable: \
 4007 4101 4127 4189 4201 4244 4305 4430 4508 4706 4709 4996 )

Deal(); // float Rand() { return Deal() / float( Rand_Mask ); }

WinMain( int, int, int, int ) {
 int Random = Deal();
} // Here, " Random " is a random number from 0 to " 2^31 - 1 ".
  // Individually, each of the 31 bits are equally random.

#define Tics ( QueryPerformanceCounter( & _Tics ), _Tics )

#define LOOP while ( 1 )
#define Loop( N ) int J = -1, LLL = N ; while ( ++ J < LLL )
#define LoopD( N ) int J = N ; while( -- J )

__int64 _Tics ;

typedef unsigned int uint ; typedef unsigned char uchar ;
typedef int * pInt ; typedef void * pVoid ; typedef uchar * ucharP ;
const int SzInt = sizeof( int );

extern "C" pVoid memmove( pVoid, pVoid, uint );
extern "C" pVoid memset( pVoid, int, uint );
extern "C" int __stdcall QueryPerformanceCounter( __int64 * );

const int HighResQual = 0x400, DecM = HighResQual * 4- 1
, LongLeg = 63, Short_Leg = 37, BothLegs = LongLeg + Short_Leg
, Seed_Buf_Sz = 2 * BothLegs - 1, StopBit = 1<<30
, EvenBits = StopBit - 2, Rand_Mask = StopBit - 1 ;

int LegsBuf[ BothLegs ], HighResBuf[ HighResQual ];
pInt RanE = HighResBuf + BothLegs, RanP ;
ucharP Dec = ucharP( HighResBuf );

Shuffle () {
 int StreamSperation = 70 - 1, Seed_Buf[ Seed_Buf_Sz ]
 , Seed = Tics, Seed_2 = Seed + 2 & EvenBits; RanP = RanE ;
 { Loop( BothLegs ) {
     Seed_Buf[ J ] = Seed_2, Seed_2 <<= 1 ;
     if ( Seed_2 >= StopBit ) Seed_2 -= EvenBits ; } }
 memset( Seed_Buf + BothLegs, 0, ( BothLegs - 1 ) * SzInt );
 Seed_Buf[ 1 ] ++, Seed_2 = Seed & Rand_Mask ;
 LOOP { { LoopD( BothLegs ) Seed_Buf[ 2 * J ] = Seed_Buf[ J ]; }
   { pInt P = Seed_Buf + 1, R = Seed_Buf + Seed_Buf_Sz - 1;
     Loop( BothLegs - LongLeg / 2 - 1 )
       P[ 2 * J ] = R[ -2 * J ] & EvenBits; }
   { pInt B = Seed_Buf, P = B + Seed_Buf_Sz, R = B + Short_Leg - 1 ;
     LoopD( BothLegs ) { if ( ! ( * -- P & 1 ) ) continue;
     R[ J ] = R[ J ] - * P & Rand_Mask ;
     B[ J - 1 ] = B[ J - 1 ] - * P & Rand_Mask ; } }
   if ( Seed_2 & 1 ) {
     memmove( Seed_Buf + 1, Seed_Buf, BothLegs * SzInt );
     * Seed_Buf = Seed_Buf[ BothLegs ];
     if ( Seed_Buf[ BothLegs ] & 1 )
       Seed_Buf [ Short_Leg ]
       = Seed_Buf [ Short_Leg ]
         - Seed_Buf[ BothLegs ] & Rand_Mask ; }
   if ( ! Seed_2 && ! -- StreamSperation ) break;
   Seed_2 >>= 1 ; }
 memmove( LegsBuf, Seed_Buf + Short_Leg, LongLeg * SzInt );
 memmove( LegsBuf + LongLeg, Seed_Buf, Short_Leg * SzInt );
 Deal(); }

int Deal() { if ( ! RanP ) Shuffle();
 if ( RanP && RanP < RanE ) return * RanP ++ ;
 RanP = HighResBuf, memmove( HighResBuf, LegsBuf, sizeof LegsBuf );
 { Loop( HighResQual - BothLegs )
     HighResBuf [ J + BothLegs ]
     = HighResBuf [ J ] - HighResBuf [ J + LongLeg ] & Rand_Mask ; }
 pInt P = HighResBuf + HighResQual - Short_Leg ;
 { Loop( Short_Leg )
     LegsBuf [ J ] = P [ J - LongLeg ] - P [ J ] & Rand_Mask ; }
 P = HighResBuf + HighResQual - LongLeg ;
 Loop( LongLeg )
   LegsBuf [ J + Short_Leg ] = P [ J ] - LegsBuf [ J ] & Rand_Mask ;
 return * RanP ++ ; }