• 3

 点此查看本文代码

 相亲数(Amicable Pair)，又称亲和数、友爱数，指两个正整数中，彼此的全部约数之和（本身除外）与另一方相等。

寻找指定范围内的相亲数的一个方法如下：对于范围内的每个数，用试除法获得约数，之后求和：

    int GetAmicableNumber(int nNumber)
    ...{
            int nSumDivisors = SumDivisors(nNumber);
            if (nSumDivisors == nNumber) return 0;//perfect number
            if (nNumber == SumDivisors(nSumDivisors)) return nSumDivisors;
            return 0;
    }
    int SumDivisors(int n)
    ...{
            int sum = 0;                     //Declaring variable for sum of divisors
            for (int i = 1; i <= n / 2; i++)
                if (n % i == 0)
                    sum += i;
            return sum;                      //Return value of sum of divisors
    }

但是这个算法效率很低。一个优化的算法是在计算约数之前进行质因数分解，然后用质因数分解的结果来计算约数之和。

这里p_i为M的质因子，而k_i为质因子的幂，那么M的约数之和可以表示为

而上式中的求和可以用一个简单的除法来替代

这样做的好处是

1. 减少了运算的次数。
2. 质数表可以在分解不同的数时重用。由于质数表只生成一次，这里没做什么优化，就用了最古老的筛法。
   

   System::Boolean CAmicableNumberPairGenerator::StartWork() 
   ...{
       m_nRangeMax=Range.Max;
       m_nNumberInQuestion=m_nRangeMin=Range.Min;
       m_rAmicableNumberPairs=gcnew List<CAmicableNumberPair^>();
       m_rListPrimes=gcnew List<int>();;
       GeneratePrimes(m_nRangeMax+1);
       return true;
   }
   void CAmicableNumberPairGenerator::GeneratePrimes(int nMax) 
   ...{
       //should be able to use vector, but since this is not the bottleneck, use manageed array 
       array<Primality>^ arPrimality=gcnew array<Primality>(nMax);
       for each (Primality p in arPrimality)
       ...{
           p=Primality ::Unknown;
       }
       //not really, but it helps
       arPrimality[0]=Primality::Prime;
       arPrimality[1]=Primality::Prime;
       int i,j;
       //Loop through each unknown value
       for (i= 2; i<nMax; i++)
       ...{
           if(arPrimality[i]!=Primality::Unknown) continue;
           //The next FirstUnknown value must be prime
           arPrimality[i] = Primality::Prime;
           m_rListPrimes->Add(i);
           //mark out all multiples of this prime as being Composite
           for (j = i * 2; j < nMax; j += i)
           ...{
               arPrimality[j] = Primality::Composite;
           }
       }
   }

3. 质数表可以用来判断一个数是否是质数。质数不是相亲数。在这里我用的是遍历，但是可以用二叉树来进一步优化。
   

   //Primality check. Primes are not amicable
   for each (int i in m_rListPrimes)
   ...{
       if(i==nNumberInQuestion) return true;//next number
       if(i>nNumberInQuestion) break;
   }
   

   另外一些优化包括
4. 仅在计算约数和时使用64位整数
5. 跳过已经找出的相亲数对中的数。由于相亲数对并不多，这里不必用二叉树。
   

   //amicable pair check. skip if it is in an existing amicable pair
   for each (CAmicableNumberPair^ anp in m_rAmicableNumberPairs)
   ...{
       if(anp->n==nNumberInQuestion) return true;//next number
   }
   Int64 nRestrictedDivisorFunction =RestrictedDivisorFunction(nNumberInQuestion);
   
   // if nNumberInQuestion > RestrictedDivisorFunction (nNumberInQuestion) and nNumberInQuestion is amicable
   // then (nRestrictedDivisorFunction,nNumberInQuestion) can not pass the amicable pair check
   if(nRestrictedDivisorFunction< nNumberInQuestion) return true;//next number

跳过约数和小于284的数。284以内只有220一个是相亲数，其约数之和是284。

//the first amicable partner is 284, and perfect numbers are not our target
if( nRestrictedDivisorFunction < 284 || nRestrictedDivisorFunction ==nNumberInQuestion) return true;//next number

如果一个数可以被6整除，且其约数之和是奇数，那么这个数不是相亲数，可以跳过一步求约数和运算

//if nNumberInQuestion mod 6==0 and nRestrictedDivisorFunction is even
//then they can not be amicable pairs
//Lee, E. J. "On Divisibility of the Sums of Even Amicable Pairs." Math. Comput. 23, 545-548, 1969.
if(nNumberInQuestion % 6 == 0 && nRestrictedDivisorFunction % 2 == 1) return true;//next number

质因数分解中，不用遍历整个质数表，在商小于当前质数的平方时就可以终止分解

Int64 CAmicableNumberPairGenerator::RestrictedDivisorFunction(Int64 nNumber)
...{
    if(nNumber<4) return 1; //4 is the first composite number

    //decomposition results
    stack<CFactor> stackFactors;

    //temp vars
    int nPrime;
    CFactor f; 
    Int64 nNumberToFactor=nNumber;
    for each (nPrime in m_rListPrimes)
    ...{
        if(nNumberToFactor==1) break;//done

        //if nPrime> sqrt(nNumberToFactor), then nNumberToFactor is a prime
        if(nPrime > nNumberToFactor/nPrime) 
        ...{
            //nNumberToFactor is a prime
            f.exponent=1;f.prime=(int)nNumberToFactor;
            stackFactors.push(f);
            break;//done
        }
        //composite
        while(nNumberToFactor % nPrime==0)
        ...{
            //increase the exponent if it is already a factor
            if(stackFactors.size()>0 && nPrime==stackFactors.top().prime)
            ...{
                f = stackFactors.top();
                stackFactors.pop();
                f.exponent = f.exponent + 1;
                stackFactors.push(f);
            }
            else
            ...{
                //a new factor
                f.exponent=1;f.prime=nPrime;
                stackFactors.push(f);
            }
            nNumberToFactor /= nPrime;
        }
    }
    //calculate the sum of the aliquot divisors
    Int64 nSumDivisor=1;
    Int64 nPrime64;//cache for 64 bit operations
    //for each factor f, calculate the sum of the divisors of pow( f.nPrime, f.exponent)
    //and multiply with existing result, nSumDivisor
    while(stackFactors.size()>0)
    ...{
        f = stackFactors.top();
        stackFactors.pop();
        nPrime64=f.prime;
        //the sum of the divisors of pow( nPrime, exponent)
        //is 1 + nPrime + nPrime ^ 2 + nPrime ^ 3 +...+ nPrime ^ exponent
        //and can be calculated by the fomula 
        // ( pow ( nPrime, exponent + 1) - 1 ) / ( nPrime - 1 )
        nSumDivisor*= ( pow ( nPrime64, f.exponent + 1) - 1 );
        nSumDivisor/=(nPrime64- 1);
    }
    //calculate the restricted divisor function 
    return nSumDivisor-nNumber;
}

计算乘幂时，使用平方乘幂法

Int64 CAmicableNumberPairGenerator::pow(Int64 m, int n)
...{
    //uses
    //Local Data
    Int64 p = 1; //working storage for the power.
    if( n<0 )...{
        p=0;
    }
    else // n>=0
    ...{
        while(n>0) // here pow(m0,n0) is p * pow(m, n)
        // where m0 and n0 were the initial values of m,n.
        ...{ //n>0
            if( n % 2 == 0 ) // so that n is 2*k or equivalently
            // k is n/2.
            ...{
                m*= m;
                n = n/2;
            }
            else // n is odd and 2*k+1
            ...{
                p = p * m;
                n--;
            }
        } // n==0
        // So pow(m0,n0) is p*pow(m,0) is p.
    } 
    return p;
}

尽可能使用STL和基本数据类型来替换范型和托管数据类型。

尽管编译器可以对代码进行优化，但是真正高效的算法还是要手工编写的。