C++ GCD函数代码示例

void CBitPatternTreeMethod::convertToIntegers(CMatrix< C_FLOAT64 > & values)
{
  bool Problems = false;

  static const C_FLOAT64 limit = 10.0 / std::numeric_limits< C_INT32 >::max();

  size_t Size = values.size();
  size_t Columns = values.numCols();

  C_FLOAT64 * pColumn = values.array();
  C_FLOAT64 * pColumnEnd = pColumn + Columns;
  C_FLOAT64 * pValue = values.array();
  C_FLOAT64 * pValueEnd = pColumn + Size;

  for (; pColumn < pColumnEnd; ++pColumn)
    {
      unsigned C_INT32 Multiplier = 1;
      unsigned C_INT32 m00, m01, m10, m11;
      unsigned C_INT32 maxden = 10000000;
      C_INT32 GCD1, GCD2;
      unsigned C_INT32 ai;

      C_FLOAT64 x;

      for (pValue = pColumn; pValue < pValueEnd; pValue += Columns)
        {
          x = fabs(*pValue);

          if (x < 100 * std::numeric_limits< C_FLOAT64 >::epsilon())
            {
              continue;
            }

          /*
           * Find rational approximation to given real number
           * David Eppstein / UC Irvine / 8 Aug 1993
           *
           * With corrections from:
           *   Arno Formella, May 2008
           *   Stefan Hoops, Sept 2009
           *
           * Based on the theory of continued fractions
           * if x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...)))
           * then best approximation is found by truncating this series
           * (with some adjustments in the last term).
           *
           * Note the fraction can be recovered as the first column of the matrix
           *  (a1 1 ) (a2 1 ) (a3 1 ) ...
           *  (1  0 ) (1  0 ) (1  0)
           * Instead of keeping the sequence of continued fraction terms,
           * we just keep the last partial product of these matrices.
           */

          /* initialize matrix */
          m00 = m11 = 1;
          m01 = m10 = 0;

          /* loop finding terms until denom gets too big */
          while (m10 *(ai = (unsigned C_INT32) x) + m11 <= maxden)
            {
              C_INT32 t;
              t = m00 * ai + m01;
              m01 = m00;
              m00 = t;

              t = m10 * ai + m11;
              m11 = m10;
              m10 = t;

              if (fabs(x - (C_FLOAT64) ai) < limit)
                break;     // SH: We reached the numerical precision of the machine;

              x = 1 / (x - (C_FLOAT64) ai);
            }

          if ((C_FLOAT64) m10 *(C_FLOAT64) ai + (C_FLOAT64) m11 > (C_FLOAT64) maxden)
            {
              Problems = true;
            }

          if (fabs(fabs(*pValue) - ((C_FLOAT64) m00) / ((C_FLOAT64) m10)) > 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon())
            {
              ai = (maxden - m11) / m10;
              m00 = m00 * ai + m01;
              m10 = m10 * ai + m11;
            }

          // Find the greatest common divisor (GCD) of the multiplier and the current denominator.
          // Euclidean algorithm
          GCD1 = m10;
          GCD2 = Multiplier;

          GCD(GCD1, GCD2);

          // Calculate the least common multiplier: LCM = v1 * v2 / GCD(v1, v2)
          Multiplier *= m10 / GCD1;
        }

      for (pValue = pColumn; pValue < pValueEnd; pValue += Columns)
        {
//.........这里部分代码省略.........

// Euclid's algorithm
int GCD(int x, int y) {
    if (y == 0)
        return x;
    return GCD(y, x%y);
}

LL LCM(LL a, LL b)
{
    return a * b / GCD(a, b);
}

static int LCMDenoCalculator( int n, int m )
{
	int gcd = GCD( n, m );
	return m/gcd;
}

int main()
{
   long n;
   GF2X a, b, c, c1, ss, ss1, tt, tt1;
   double t;
   long iter, i;

   cout << WD(12,"n") << WD(12,"OldGCD") <<  WD(12,"GCD") << WD(12,"OldXGCD")
        << WD(12, "XGCD") << "\n";

   cout.precision(3);
   cout.setf(ios::scientific);


   for (n = 32; n <= (1L << 18); n = n << 3) {
      random(a, n);
      random(b, n);
      OldGCD(c, a, b);
      GCD(c1, a, b);
      OldXGCD(c, ss, tt, a, b);
      XGCD(c1, ss1, tt1, a, b);
      if (c1 != c || ss1 != ss || tt1 != tt) {
         cerr << "**** GF2XTest FAILED!\n";
         return 1;
      }

      cout << WD(12,n); 

      iter = 0;
      do {
         iter = iter ? (2*iter) : 1;
         t = GetTime();
         for (i = 0; i < iter; i++)
            OldGCD(c, a, b);
         t = GetTime()-t;
      } while (t < 0.5);

      cout << WD(12,t/iter);

      iter = 0;
      do {
         iter = iter ? (2*iter) : 1;
         t = GetTime();
         for (i = 0; i < iter; i++)
            GCD(c, a, b);
         t = GetTime()-t;
      } while (t < 0.5);

      cout << WD(12,t/iter);

      iter = 0;
      do {
         iter = iter ? (2*iter) : 1;
         t = GetTime();
         for (i = 0; i < iter; i++)
            OldXGCD(c, ss, tt, a, b);
         t = GetTime()-t;
      } while (t < 0.5);

      cout << WD(12,t/iter);

      iter = 0;
      do {
         iter = iter ? (2*iter) : 1;
         t = GetTime();
         for (i = 0; i < iter; i++)
            XGCD(c, ss, tt, a, b);
         t = GetTime()-t;
      } while (t < 0.5);

      cout << WD(12,t/iter);

      cout << "\n";
   }

   return 0;
}

int LCM(int x, int y){    
    return x * y / GCD(x, y);
}

long IterIrredTest(const ZZ_pEX& f)
{
   if (deg(f) <= 0) return 0;
   if (deg(f) == 1) return 1;

   ZZ_pEXModulus F;

   build(F, f);
   
   ZZ_pEX h;

   FrobeniusMap(h, F);

   long CompTableSize = 2*SqrRoot(deg(f));

   ZZ_pEXArgument H;

   build(H, h, F, CompTableSize);

   long i, d, limit, limit_sqr;
   ZZ_pEX g, X, t, prod;


   SetX(X);

   i = 0;
   g = h;
   d = 1;
   limit = 2;
   limit_sqr = limit*limit;

   set(prod);


   while (2*d <= deg(f)) {
      sub(t, g, X);
      MulMod(prod, prod, t, F);
      i++;
      if (i == limit_sqr) {
         GCD(t, f, prod);
         if (!IsOne(t)) return 0;

         set(prod);
         limit++;
         limit_sqr = limit*limit;
         i = 0;
      }

      d = d + 1;
      if (2*d <= deg(f)) {
         CompMod(g, g, H, F);
      }
   }

   if (i > 0) {
      GCD(t, f, prod);
      if (!IsOne(t)) return 0;
   }

   return 1;
}

unsigned int LCM( unsigned int a, unsigned int b ) {

  unsigned int tmp=a/GCD(a,b);
  return tmp*b;

}

//.........这里部分代码省略.........
    {14112, 18415, 28, 11278},
    {15004, 15709, 22,  3867},
    {15360, 20485, 24, 12767},
 // {16384, 21845, 16, 12798},
    {17208 ,21931, 24, 18387},
    {18000, 18631, 25,  4208},
    {18816, 24295, 28, 16360},
    {19200, 21607, 40, 35633},
    {21168, 27305, 28, 15407},
    {23040, 23377, 48,  5292},
    {24576, 24929, 48,  5612},
    {27000, 32767, 15, 20021},
    {31104, 31609, 71,  5149},
    {42336, 42799, 21,  5952},
    {46080, 53261, 24, 33409},
    {49140, 57337, 39,  2608},
    {51840, 59527, 72, 21128},
    {61680, 61681, 40,  1273},
    {65536, 65537, 32,  1273},
    {75264, 82603, 56, 36484},
    {84672, 92837, 56, 38520}
  };

#if 0

  for (long i = 0; i < 25; i++) {
    long m = ms[i][1];
    PAlgebra alg(m);
    alg.printout();
    cout << "\n";
    // compute phi(m) directly
    long phim = 0;
    for (long j = 0; j < m; j++)
      if (GCD(j, m) == 1) phim++;

    if (phim != alg.phiM()) cout << "ERROR\n";
  }

  exit(0);


#endif



  // find the first m satisfying phi(m)>=N and d | ord(2) in Z_m^*
  long m = 0;
  for (unsigned i=0; i<sizeof(ms)/sizeof(long[3]); i++) 
    if (ms[i][0]>=N && (ms[i][2] % d) == 0) {
      m = ms[i][1];
      c_m = 0.001 * (double) ms[i][3];
      break;
    }
  if (m==0) Error("Cannot support this L,d combination");
#endif
  //  m = 257;
  FHEcontext context(m);
#if 0
  context.stdev = to_xdouble(0.5); // very low error
#endif
  activeContext = &context; // Mark this as the "current" context

  context.zMstar.printout();
  cout << endl;

// Generate the representation of Z_m^* for a given odd integer m
// and plaintext base p
PAlgebra::PAlgebra(unsigned long mm, unsigned long pp)
{
    m = mm;
    p = pp;

    assert( (m&1) == 1 );
    assert( ProbPrime(p) );
    // replaced by Mahdi after a conversation with Shai
    // assert( m > p && (m % p) != 0 );	// original line
    assert( (m % p) != 0 );
    // end of replace by Mahdi
    assert( m < NTL_SP_BOUND );

    // Compute the generators for (Z/mZ)^*
    vector<unsigned long> classes(m);
    vector<long> orders(m);

    unsigned long i;
    for (i=0; i<m; i++) { // initially each element in its own class
        if (GCD(i,m)!=1)
            classes[i] = 0; // i is not in (Z/mZ)^*
        else
            classes[i] = i;
    }

    // Start building a representation of (Z/mZ)^*, first use the generator p
    conjClasses(classes,p,m);  // merge classes that have a factor of 2

    // The order of p is the size of the equivalence class of 1
    ordP = (unsigned long) count (classes.begin(), classes.end(), 1);

    // Compute orders in (Z/mZ)^*/<p> while comparing to (Z/mZ)^*
    long idx, largest;
    while (true) {
        compOrder(orders,classes,true,m);
        idx = argmax(orders);      // find the element with largest order
        largest = orders[idx];

        if (largest <= 0) break;   // stop comparing to order in (Z/mZ)^*

        // store generator with same order as in (Z/mZ)^*
        gens.push_back(idx);
        ords.push_back(largest);
        conjClasses(classes,idx,m); // merge classes that have a factor of idx
    }
    // Compute orders in (Z/mZ)^*/<p> without comparing to (Z/mZ)^*
    while (true) {
        compOrder(orders,classes,false,m);
        idx = argmax(orders);      // find the element with largest order
        largest = orders[idx];

        if (largest <= 0) break;   // we have the trivial group, we are done

        // store generator with different order than (Z/mZ)^*
        gens.push_back(idx);
        ords.push_back(-largest);  // store with negative sign
        conjClasses(classes,idx,m);  // merge classes that have a factor of idx
    }

    nSlots = qGrpOrd();
    phiM = ordP * nSlots;

    // Allocate space for the various arrays
    T.resize(nSlots);
    dLogT.resize(nSlots*gens.size());
    Tidx.assign(m,-1);    // allocate m slots, initialize them to -1
    zmsIdx.assign(m,-1);  // allocate m slots, initialize them to -1
    for (i=idx=0; i<m; i++) if (GCD(i,m)==1) zmsIdx[i] = idx++;

    // Now fill the Tidx and dLogT translation tables. We identify an element
    // t\in T with its representation t = \prod_{i=0}^n gi^{ei} mod m (where
    // the gi's are the generators in gens[]) , represent t by the vector of
    // exponents *in reverse order* (en,...,e1,e0), and order these vectors
    // in lexicographic order.

    // buffer is initialized to all-zero, which represents 1=\prod_i gi^0
    vector<unsigned long> buffer(gens.size()); // temporaty holds exponents
    i = idx = 0;
    do {
        unsigned long t = exponentiate(buffer);
        for (unsigned long j=0; j<buffer.size(); j++) dLogT[idx++] = buffer[j];
        T[i] = t;       // The i'th element in T it t
        Tidx[t] = i++;  // the index of t in T is i

        // increment buffer by one (in lexigoraphic order)
    } while (nextExpVector(buffer)); // until we cover all the group

    PhimX = Cyclotomic(m); // compute and store Phi_m(X)

    // initialize prods array
    long ndims = gens.size();
    prods.resize(ndims+1);
    prods[ndims] = 1;
    for (long j = ndims-1; j >= 0; j--) {
        prods[j] = OrderOf(j) * prods[j+1];
    }
}

/**************************
//Paillier HE system
//len: length of params p,q
**************************/
void Paillier(int len=512){
  ZZ n, n2, p, q, g, lamda;
  ZZ p1, q1, p1q1;
  ZZ miu;
  
  ZZ m1, m2;
  ZZ BSm, HEm; //baseline and HE result
  ZZ c, c1, c2, cm1, cm2, r;

  //key gen
  start = std::clock();

  GenPrime(p, len);
  GenPrime(q, len);
  mul(n, p, q);
  mul(n2, n, n);

  sub(p1,p,1);
  sub(q1,q,1);
  GCD(lamda,p1,q1);
  mul(p1q1,p1,q1);
  div(lamda, p1q1, lamda);

  RandomBnd(g, n2);

  PowerMod(miu,g,lamda,n2);
  sub(miu, miu, 1);
  div(miu,miu,n); //should add 1?
  InvMod(miu, miu, n);
  
  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
  cout<<"Pailler Setup:"<< duration <<'\n';

  //Enc
  start = std::clock();
  RandomBnd(m1,n);
  RandomBnd(m2,n);

  RandomBnd(r,n); //enc m1
  PowerMod(c1, g,m1,n2);
  PowerMod(c2, r,n,n2);
  MulMod(cm1, c1,c2, n2);

  RandomBnd(r,n); //enc m2
  PowerMod(c1, g,m2,n2);
  PowerMod(c2, r,n,n2);
  MulMod(cm2, c1,c2, n2);

  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
  cout<<"Pailler Enc:"<< duration/2 <<'\n';

  //Evaluation
  start = std::clock();
  c=cm1;
  for(int i=0; i<TIMES; i++)
  	MulMod(c,c,cm2,n2);
  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
  cout<<"Pailler Eval:"<< duration <<'\n';

  //c=cm2;
  //Dec  
  start = std::clock();
  PowerMod(c,c,lamda,n2);
  sub(c,c,1);
  div(c,c,n); //should add 1?
  MulMod(HEm, c, miu, n);  

  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
  cout<<"Pailler Dec:"<< duration <<'\n';

  //baseline
  BSm=m1;
  for(int i=0; i<TIMES; i++)
  	AddMod(BSm,BSm,m2,n);

  assert(BSm==HEm);
}

int GCD(int a, int b) {
	//printf(">> %d %d\n", a, b);
	if (b==0) return a;
	else return GCD(b, a%b);
};

void main(void){
  printf("%d\n",GCD(6,12));
}

//.........这里部分代码省略.........

          C_INT64 TMP = Identity[CurrentRowIndex];
          Identity[CurrentRowIndex]  = Identity[NonZeroIndex];
          Identity[NonZeroIndex] = TMP;
        }

      if (*(pActiveRowStart + CurrentColumnIndex) < 0)
        {
          for (pRow = pActiveRowStart; pRow < pActiveRowEnd; ++pRow)
            {
              *pRow *= -1;
            }

          Identity[CurrentRowIndex] *= -1;
        }

      // For each row
      pRow = pActiveRowStart + NumCols;
      pIdentity = Identity.array() + CurrentRowIndex + 1;

      C_INT64 ActiveRowValue = *(pActiveRowStart + CurrentColumnIndex);
      *(pActiveRowStart + CurrentColumnIndex) = Identity[CurrentRowIndex];

      for (; pRow < pRowEnd; pRow += NumCols, ++pIdentity)
        {
          C_INT64 RowValue = *(pRow + CurrentColumnIndex);

          if (RowValue == 0)
            continue;

          *(pRow + CurrentColumnIndex) = 0;

          // compute GCD(*pActiveRowStart, *pRow)
          C_INT64 GCD1 = abs64(ActiveRowValue);
          C_INT64 GCD2 = abs64(RowValue);

          GCD(GCD1, GCD2);

          C_INT64 alpha = ActiveRowValue / GCD1;
          C_INT64 beta = RowValue / GCD1;

          // update rest of row
          pActiveRow = pActiveRowStart;
          pCurrent = pRow;
          *pIdentity *= alpha;

          GCD1 = abs64(*pIdentity);

          for (; pActiveRow < pActiveRowEnd; ++pActiveRow, ++pCurrent)
            {
              // Assert that we do not have a numerical overflow.
              assert(fabs(((C_FLOAT64) alpha) *((C_FLOAT64) * pCurrent) - ((C_FLOAT64) beta) *((C_FLOAT64) * pActiveRow)) < std::numeric_limits< C_INT64 >::max());

              *pCurrent = alpha * *pCurrent - beta * *pActiveRow;

              // We check that the row values do not have any common divisor.
              if (GCD1 > 1 &&
                  (GCD2 = abs64(*pCurrent)) > 0)
                {
                  GCD(GCD1, GCD2);
                }
            }

          if (GCD1 > 1)
            {
              *pIdentity /= GCD1;

// TODO: Add the ability to set a maximum number of iterations
inline BigInt FindFactor
( const BigInt& n,
  Int a,
  const PollardRhoCtrl& ctrl )
{
    if( a == 0 || a == -2 )
        Output("WARNING: Problematic choice of Pollard rho shift");
    BigInt tmp, gcd;
    BigInt one(1);

    auto xAdvance =
      [&]( BigInt& x )
      {
        if( ctrl.numSteps == 1 )
        {
            // TODO: Determine if there is a penalty to x *= x
            /*
            tmp = x;
            tmp *= x;
            tmp += a;
            x = tmp;
            x %= n;
            */
            x *= x;
            x += a;
            x %= n;
        }
        else
        {
            PowMod( x, 2*ctrl.numSteps, n, x );
            x += a;
            x %= n;
        }
      };

    auto QAdvance =
      [&]( const BigInt& x, const BigInt& x2, BigInt& Q )
      {
        tmp = x2;
        tmp -= x;
        Q *= tmp;
        Q %= n;
      };

    Int gcdDelay = ctrl.gcdDelay;
    BigInt xi=ctrl.x0;
    BigInt x2i(xi);
    BigInt xiSave=xi, x2iSave=x2i;
    BigInt Qi(1);
    Int k=1, i=1; // it is okay for i to overflow since it is just for printing
    while( true )
    {
        // Advance xi once
        xAdvance( xi );

        // Advance x2i twice
        xAdvance( x2i );
        xAdvance( x2i );

        // Advance Qi
        QAdvance( xi, x2i, Qi );

        if( k >= gcdDelay )
        {
            GCD( Qi, n, gcd );
            if( gcd > one )
            {
                // NOTE: This was not suggested by Pollard's original paper
                if( gcd == n )
                {
                    if( gcdDelay == 1 )
                    {
                        RuntimeError("(x) converged before (x mod p) at i=",i);
                    }
                    else
                    {
                        if( ctrl.progress )
                            Output("Backtracking at i=",i);
                        i = Max( i-(gcdDelay+1), 0 );
                        gcdDelay = 1;
                        xi = xiSave;
                        x2i = x2iSave;
                    }
                }
                else
                {
                    if( ctrl.progress )
                        Output("Found factor ",gcd," at i=",i); 
                    return gcd;
                }
            }

            // NOTE: This was not suggested by Pollard's original paper
            k = 0;
            xiSave = xi;
            x2iSave = x2i;
            Qi = 1;
        }
        ++k;
//.........这里部分代码省略.........