Function DivStep62

IntMod.cpp:110–239

// ------------------------------------------------

inline void DivStep62(int64_t u0,int64_t v0,
  int64_t* eta,
  int64_t* uu,int64_t* uv,
  int64_t* vu,int64_t* vv) {

  // u' = (uu*u + uv*v) >> bitCount
  // v' = (vu*u + vv*v) >> bitCount
  // Do not maintain a matrix for r and s, the number of
  // 'added P' can be easily calculated

  int bitCount;

#if 0

  #define SWAP_ADD(x,y) x+=y;y-=x;
  #define SWAP_SUB(x,y) x-=y;y+=x;

  // Former divstep62 (using __builtin_ctzll)
  // Do not use eta, u and v have an exponential decay in worst case
  // but with low probability to reach this worst case complexity
  // Avg: 581 Kinv/s

  bitCount = 62;
  int64_t nb0;

  while(true) {

    int zeros = __builtin_ctzll(v0 | (UINT64_MAX << bitCount));
    v0 >>= zeros;
    *uu <<= zeros;
    *uv <<= zeros;
    bitCount -= zeros;

    if(bitCount <= 0)
      break;

    nb0 = (v0 + u0) & 0x3;
    if(nb0 == 0) {
      SWAP_ADD(*vv,*uv);
      SWAP_ADD(*vu,*uu);
      SWAP_ADD(v0,u0);
    } else {
      SWAP_SUB(*vv,*uv);
      SWAP_SUB(*vu,*uu);
      SWAP_SUB(v0,u0);
    }

  }


#endif

#if 1

  #define SWAP_NEG(tmp,x,y) tmp = x; x = y; y = -tmp;

  int64_t m,w,x,y,z;
  bitCount = 62;