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 ++ ; }