自学教程：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 algorithmint 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 pPAlgebra::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 iterationsinline 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;//.........这里部分代码省略.........

// Reduce plaintext space to a divisor of the original plaintext spacevoid Ctxt::reducePtxtSpace(long newPtxtSpace){  long g = GCD(ptxtSpace, newPtxtSpace);  assert (g>1);  ptxtSpace = g;}

bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod){	assert(!equiv.IsNegative() && equiv < mod);	Integer gcd = GCD(equiv, mod);	if (gcd != Integer::One())	{		// the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv)		if (p <= gcd && gcd <= max && IsPrime(gcd))		{			p = gcd;			return true;		}		else			return false;	}	BuildPrimeTable();	if (p <= primeTable[primeTableSize-1])	{		word *pItr;		--p;		if (p.IsPositive())			pItr = std::upper_bound(primeTable, primeTable+primeTableSize, p.ConvertToLong());		else			pItr = primeTable;		while (pItr < primeTable+primeTableSize && *pItr%mod != equiv)			++pItr;		if (pItr < primeTable+primeTableSize)		{			p = *pItr;			return p <= max;		}		p = primeTable[primeTableSize-1]+1;	}	assert(p > primeTable[primeTableSize-1]);	if (mod.IsOdd())		return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1);	p += (equiv-p)%mod;	if (p>max)		return false;	PrimeSieve sieve(p, max, mod);	while (sieve.NextCandidate(p))	{		if (FastProbablePrimeTest(p) && IsPrime(p))			return true;	}	return false;}

// Encrypts plaintext, result returned in the ciphertext argument. The// returned value is the plaintext-space for that ciphertext. When called// with highNoise=true, returns a ciphertext with noise level~q/8.long FHEPubKey::Encrypt(Ctxt &ctxt, const ZZX& ptxt, long ptxtSpace,			bool highNoise) const{  FHE_TIMER_START;  assert(this == &ctxt.pubKey);  if (ptxtSpace != pubEncrKey.ptxtSpace) { // plaintext-space mistamtch    ptxtSpace = GCD(ptxtSpace, pubEncrKey.ptxtSpace);    if (ptxtSpace <= 1) Error("Plaintext-space mismatch on encryption");  }  // generate a random encryption of zero from the public encryption key  ctxt = pubEncrKey;  // already an encryption of zero, just not a random one  // choose a random small scalar r and a small random error vector e,  // then set ctxt = r*pubEncrKey + ptstSpace*e + (ptxt,0)  DoubleCRT e(context, context.ctxtPrimes);  DoubleCRT r(context, context.ctxtPrimes);  r.sampleSmall();  for (size_t i=0; i<ctxt.parts.size(); i++) {  // add noise to all the parts    ctxt.parts[i] *= r;    if (highNoise && i == 0) {      // we sample e so that coefficients are uniform over       // [-Q/(8*ptxtSpace)..Q/(8*ptxtSpace)]      ZZ B;      B = context.productOfPrimes(context.ctxtPrimes);      B /= ptxtSpace;      B /= 8;      e.sampleUniform(B);    }    else {       e.sampleGaussian();    }    e *= ptxtSpace;    ctxt.parts[i] += e;  }  // add in the plaintext  // FIXME: This relies on the first part, ctxt[0], to have handle to 1  if (ptxtSpace==2) ctxt.parts[0] += ptxt;  else { // The general case of ptxtSpace>2: for a ciphertext         // relative to modulus Q, we add ptxt * Q mod ptxtSpace.    long QmodP = rem(context.productOfPrimes(ctxt.primeSet), ptxtSpace);    ctxt.parts[0] += MulMod(ptxt,QmodP,ptxtSpace); // MulMod from module NumbTh  }  // fill in the other ciphertext data members  ctxt.ptxtSpace = ptxtSpace;  if (highNoise) {    // hack: we set noiseVar to Q^2/8, which is just below threshold     // that will signal an error    ctxt.noiseVar = xexp(2*context.logOfProduct(context.ctxtPrimes) - log(8.0));  }  else {    // We have <skey,ctxt>= r*<skey,pkey> +p*(e0+e1*s) +m, where VAR(<skey,pkey>)    // is recorded in pubEncrKey.noiseVar, VAR(ei)=sigma^2*phi(m), and VAR(s) is    // determined by the secret-key Hamming weight (skHwt).     // VAR(r)=phi(m)/2, hence the expected size squared is bounded by:    // E(X^2) <= pubEncrKey.noiseVar *phi(m) *stdev^2    //                               + p^2*sigma^2 *phi(m) *(skHwt+1) + p^2      long hwt = skHwts[0];    xdouble phim = to_xdouble(context.zMStar.getPhiM());    xdouble sigma2 = context.stdev * context.stdev;    xdouble p2 = to_xdouble(ptxtSpace) * to_xdouble(ptxtSpace);    ctxt.noiseVar = pubEncrKey.noiseVar*phim*0.5                     + p2*sigma2*phim*(hwt+1)*context.zMStar.get_cM() + p2;  }  return ptxtSpace;}

void TranslateBetweenGrids( const DistMatrix<T,MC,MR>& A, DistMatrix<T,MC,MR>& B ) {    DEBUG_ONLY(CSE cse("copy::TranslateBetweenGrids [MC,MR]"))    B.Resize( A.Height(), A.Width() );    // Just need to ensure that each viewing comm contains the other team's    // owning comm. Congruence is too strong.    // Compute the number of process rows and columns that each process    // needs to send to.    const Int colStride = B.ColStride();    const Int rowStride = B.RowStride();    const Int colRank = B.ColRank();    const Int rowRank = B.RowRank();    const Int colStrideA = A.ColStride();    const Int rowStrideA = A.RowStride();    const Int colGCD = GCD( colStride, colStrideA );    const Int rowGCD = GCD( rowStride, rowStrideA );    const Int colLCM = colStride*colStrideA / colGCD;    const Int rowLCM = rowStride*rowStrideA / rowGCD;    const Int numColSends = colStride / colGCD;    const Int numRowSends = rowStride / rowGCD;    const Int colAlign = B.ColAlign();    const Int rowAlign = B.RowAlign();    const Int colAlignA = A.ColAlign();    const Int rowAlignA = A.RowAlign();    const bool inBGrid = B.Participating();    const bool inAGrid = A.Participating();    if( !inBGrid && !inAGrid )        return;    const Int maxSendSize =        (A.Height()/(colStrideA*numColSends)+1) *        (A.Width()/(rowStrideA*numRowSends)+1);    // Translate the ranks from A's VC communicator to B's viewing so that    // we can match send/recv communicators. Since A's VC communicator is not    // necessarily defined on every process, we instead work with A's owning    // group and account for row-major ordering if necessary.    const int sizeA = A.Grid().Size();    vector<int> rankMap(sizeA), ranks(sizeA);    if( A.Grid().Order() == COLUMN_MAJOR )    {        for( int j=0; j<sizeA; ++j )            ranks[j] = j;    }    else    {        // The (i,j) = i + j*colStrideA rank in the column-major ordering is        // equal to the j + i*rowStrideA rank in a row-major ordering.        // Since we desire rankMap[i+j*colStrideA] to correspond to process        // (i,j) in A's grid's rank in this viewing group, ranks[i+j*colStrideA]        // should correspond to process (i,j) in A's owning group. Since the        // owning group is ordered row-major in this case, its rank is        // j+i*rowStrideA. Note that setting        // ranks[j+i*rowStrideA] = i+j*colStrideA is *NOT* valid.        for( int i=0; i<colStrideA; ++i )            for( int j=0; j<rowStrideA; ++j )                ranks[i+j*colStrideA] = j+i*rowStrideA;    }    mpi::Translate    ( A.Grid().OwningGroup(), sizeA, &ranks[0],      B.Grid().ViewingComm(), &rankMap[0] );    // Have each member of A's grid individually send to all numRow x numCol    // processes in order, while the members of this grid receive from all    // necessary processes at each step.    Int requiredMemory = 0;    if( inAGrid )        requiredMemory += maxSendSize;    if( inBGrid )        requiredMemory += maxSendSize;    vector<T> auxBuf( requiredMemory );    Int offset = 0;    T* sendBuf = &auxBuf[offset];    if( inAGrid )        offset += maxSendSize;    T* recvBuf = &auxBuf[offset];    Int recvRow = 0; // avoid compiler warnings...    if( inAGrid )        recvRow = Mod(Mod(A.ColRank()-colAlignA,colStrideA)+colAlign,colStride);    for( Int colSend=0; colSend<numColSends; ++colSend )    {        Int recvCol = 0; // avoid compiler warnings...        if( inAGrid )            recvCol=Mod(Mod(A.RowRank()-rowAlignA,rowStrideA)+rowAlign,                        rowStride);        for( Int rowSend=0; rowSend<numRowSends; ++rowSend )        {            mpi::Request sendRequest;            // Fire off this round of non-blocking sends            if( inAGrid )            {                // Pack the data                Int sendHeight = Length(A.LocalHeight(),colSend,numColSends);                Int sendWidth = Length(A.LocalWidth(),rowSend,numRowSends);//.........这里部分代码省略.........

int GCD(int num1, int num2) {    if ( num1 % num2 == 0)        return num2;    else        return GCD( num2, num1 % num2);}

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

void CRational::Normalize(){	const int gcd = GCD(abs(m_numerator), m_denominator);	m_numerator /= gcd;	m_denominator /= gcd;}

PAlgebra::PAlgebra(unsigned long mm, unsigned long pp,                     const vector<long>& _gens, const vector<long>& _ords ){  assert( ProbPrime(pp) );  assert( (mm % pp) != 0 );  assert( mm < NTL_SP_BOUND );  assert( mm > 1 );  cM  = 1.0; // default value for the ring constant  m = mm;  p = pp;  long k = NextPowerOfTwo(m);  if (mm == (1UL << k))    pow2 = k;  else    pow2 = 0;     // For dry-run, use a tiny m value for the PAlgebra tables  if (isDryRun()) mm = (p==3)? 4 : 3;  // Compute the generators for (Z/mZ)^* (defined in NumbTh.cpp)  if (_gens.size() == 0 || isDryRun())       ordP = findGenerators(this->gens, this->ords, mm, pp);  else {    assert(_gens.size() == _ords.size());    gens = _gens;    ords = _ords;    ordP = multOrd(pp, mm);  }  nSlots = qGrpOrd();  phiM = ordP * nSlots;  // Allocate space for the various arrays  T.resize(nSlots);  dLogT.resize(nSlots*gens.size());  Tidx.assign(mm,-1);    // allocate m slots, initialize them to -1  zmsIdx.assign(mm,-1);  // allocate m slots, initialize them to -1  long i, idx;  for (i=idx=0; i<(long)mm; i++) if (GCD(i,mm)==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.  // FIXME: is the comment above about reverse order true? It doesn't   // seem like it to me.  VJS.  // 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;  long ctr = 0;  do {    ctr++;    unsigned long t = exponentiate(buffer);    for (unsigned long j=0; j<buffer.size(); j++) dLogT[idx++] = buffer[j];    assert(GCD(t,mm) == 1); // sanity check for user-supplied gens    assert(Tidx[t] == -1);    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  assert(ctr == long(nSlots)); // sanity check for user-supplied gens  PhimX = Cyclotomic(mm); // 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];  }  //  pp_factorize(mFactors,mm); // prime-power factorization from NumbTh.cpp}

int GCD(int a,int b){    if(b==0)return a;    else return GCD(b,a%b);}

static CYTHON_INLINE struct ZZX* ZZX_gcd(struct ZZX* x, struct ZZX* y){    struct ZZX* g = new ZZX();    GCD(*g, *x, *y);    return g;}

ll GCD(ll a, ll b) { return b == 0 ? a : GCD(b, a%b); }

int main() {    int a, b;    printf("Input two numbers:/n");    scanf("%d %d", &a, &b);    printf("Greatest common divisor: %d/n", GCD(a, b));}

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

static unsigned long LCM(unsigned long x, unsigned long y){	unsigned long gcd = GCD(x, y);	return (gcd == 0) ? 0 : ((x / gcd) * y);}

int main(void) {	int num1, num2, denom1, denom2, num3, denom3, isLoop = 1;	char op, input[3], charContinue;	printf("Welcome to the fraction calculator program!/n");	while (isLoop == 1) {		printf("Please insert the first numerator:/n");		scanf("%s", input);		num1 = atoi(input);		fflush(NULL);		printf("Please insert the first denominator:/n");		scanf("%s", input);		denom1 = atoi(input);		fflush(NULL);		printf("This gives us a fraction of %i/%i./n", num1,denom1);                printf("Please insert the second numerator:/n");                scanf("%s", input);		num2 = atoi(input);                fflush(NULL);                printf("Please insert the second denominator:/n");                scanf("%s", input);		denom2 = atoi(input);                fflush(NULL);		printf("This gives us a second fraction of %i/%i./n", num2,denom2);		printf("Please insert the operand for the computation: [+ or - or * or /]/n");		scanf("%s", &op);		switch(op) {			case '+':				printf("You have chosen addition./n");				if (denom1 != denom2) {					denom3 = denom1 * denom2;					num1 = num1 * denom2;					num2 = num2 * denom1;					num3 = num1 + num2;					denom1 = denom3;					denom2 = denom3;				}				else {					denom3 = denom1;					num3 = num1 + num2;				}				break;			case '-':				printf("You have chosen substraction./n");				if (denom1 != denom2) {					denom3 = denom1 * denom2;					num1 = num1 * denom2;					num2 = num2 * denom1;					num3 = num1 - num2;					denom1 = denom3;					denom2 = denom3;				}				else {					denom3 = denom1;					num3 = num1 - num2;				}				break;			case '*':				printf("You have chosen multiplication./n");				num3 = num1 * num2;				denom3 = denom1 * denom2;						break;			case '/':				printf("You have chosen division./n");				num3 = num1 * denom2;				denom3 = denom1 * num2;				break;		}		printf("Thus, we have %i/%i %c %i/%i = %i/%i./n",num1,denom1,op,num2,denom2,num3,denom3);		printf("Would you like a simplified answer? [y/n]/n");		scanf("%s", &charContinue);		if (charContinue == 'y') {				int gcd;  			gcd = GCD(num3,denom3);			if (gcd != 0) {				if (num3 == (num3 / gcd)) {					printf("Your fraction could not be simplified any further./n");				}				else {					num3 = num3 / gcd;					denom3 = denom3 / gcd;					printf("Your simplified fraction is %i/%i./n", num3, denom3);				}			}			else {				printf("Your simplified fraction is 0./n");				}		}		printf("Would you like to continue? [y/n]/n");		scanf("%s", &charContinue);		if (charContinue != 'y') {			printf("Goodbye!/n");			isLoop = 0;		}	}	return 0;}