// Macros for doing double precision multiply
#define BPU ( 8*sizeof(unsigned) )       // Number of bits in an unsigned
#define lo(x) ( (x) & ((1<<(BPU/2))-1) ) // lower half of unsigned
#define hi(x) ( (x) >> (BPU/2) )         // upper half
#define lh(x) ( (x) << (BPU/2) )         // make upper half

void flex_unit::fast_mul( flex_unit &x, flex_unit &y, unsigned keep )
{
  // *this = (x*y) % (2**keep)
  unsigned i,limit = (keep+BPU-1)/BPU; // size of result in words
  reserve(limit); for (i=0; i<limit; i+=1) a[i] = 0;
  unsigned min = x.n; if (min>limit) min = limit;
  for (i=0; i<min; i+=1)
  {
    unsigned m = x.a[i];
    unsigned min = i+y.n; if (min>limit) min = limit;
#ifdef topspeed_asm
    unsigned c = dmul( m, min-i, y.a, a+i );
    unsigned j = min;
#else
    unsigned c = 0; // carry
    for ( unsigned j=i; j<min; j+=1 )
    {
      // This is the critical loop
      // Machine dependent code could help here
      // c:a[j] = a[j] + c + m*y.a[j-i];
      unsigned w, v = a[j], p = y.a[j-i];
      v += c; c = ( v < c );
      w = lo(p)*lo(m); v += w; c += ( v < w );
      w = lo(p)*hi(m); c += hi(w); w = lh(w); v += w; c += ( v < w );
      w = hi(p)*lo(m); c += hi(w); w = lh(w); v += w; c += ( v < w );
      c += hi(p) * hi(m);
      a[j] = v;
    }
#endif
    while ( c && j<limit )
    {
      a[j] += c;
      c = a[j] < c;
      j += 1;
    }
  }

  // eliminate unwanted bits
  keep %= BPU; if (keep) a[limit-1] &= (1<<keep)-1;

   // calculate n
  while (limit && a[limit-1]==0) limit-=1;
  n = limit;
};