天天看點

多項式相乘的高效分治算法一

        多項式相乘最容易想到的方法就是采用2個for循環,依次周遊相乘累加即可。但時間複雜度為O(n2)。那麼是否可用其他的方法來降低周遊比較的次數呢?可用分治法來進行對應切分,然後相乘,再相加。但是這好像并不能減少周遊的次數,因為與之前的合并排序不同的是,合并排序中,如果左邊的數若比右邊的數小,則左邊的左邊就不需要與右邊的所有數進行比較,因為數是sorted,是以預設比右邊的小,無需再比較;而兩個多項式相乘,不存在排序,而且每個數都必須與對面的每個數相乘,因為數都是不同的,是以必須每個數與另一個多項式的每一項相乘(注:一個加号連接配接兩個項),那麼根據以上的分析,無論如何都要進行N2次乘法,這樣即使采用分治法,最後的時間複雜度也為O(n2)。這個問題讓我不得其解,因為不管怎麼運用何種算法都不可能減少相乘的次數。然後高斯采用了一個出其不意的思想方法來減少了相乘的次數,變4次大相乘為3次大相乘,而每個數又都相乘了。

         該方法的理論依據是:

  1. 在計算上,乘法比加法更耗時,是以可用加法來代替乘法;
  2. 該方法來源于反思考宏觀的推理。

        是以,多項式相乘高效方法來源于高斯的推理。那麼首先從多項式的4次分治相乘的程式設計開始。如何采用分治法進行分割相乘呢?基本的步驟如下:

      例如有多項式A=[2,3,5,4]和B=[3,1,3,5]  ,注:多項式隻存系數,系數的位序代表了該項的幂次。位序==幂次是程式設計實作的關鍵。如何關鍵呢?

     Step1. 将A,B分别從中間分割成A0和A1,B0和B1;

     Step2.  A*B=A0*B0+A0*B1+A1*B0+A1*B1,對這四項進行遞歸分割相乘,然後疊加;

      首先,根據分治法,編寫一個主體的遞歸程式如下:

vector<double> DivideMerge::polynomialsMultiply(vector<double>&numA, int lowA, int highA, vector<double>&numB, int lowB, int highB)
{
	vector<double>numXY,numXYZ;
	vector<double>numX, numY, numZ;                                                            

	if (lowA < highA && lowB<highB)
	{
		int midA = (highA + lowA) / 2;
		int midB = (highB + lowB) / 2;

		numX = polynomialsMultiply(numA, lowA, midA, numB, lowB, midB);

		vector<double> numa, numb;
		int dy;                 //交叉相乘遞歸
		numa = polynomialsMultiply(numA, lowA, midA, numB, midB + 1, highB);
		numb = polynomialsMultiply(numA, midA + 1, highA, numB, lowB, midB);
		numY = crossAdd(numa, midA, highB, numb, highA, midB,dy);  //A0*B1+A1*B0

		int dxy;
		numXY = crossAdd(numX, midA, midB, numY, 0, dy, dxy);     //A0B0+(A0B1+A1B0)
		numZ = polynomialsMultiply(numA, midA + 1, highA, numB, midB + 1, highB);              

		int dxyz;
		numXYZ = crossAdd(numXY, 0, dxy, numZ, highA, highB,dxyz);
	}
else
	{
	 numXYZ.push_back(numA[lowA] * numB[lowB]);                                                 //調用的這個函數傳回的是一個void類型
	}
}
           

      這個遞歸中隻考慮if (lowA < highA && lowB<highB),當遞歸到終止條件lowA==highA或lowB==highB時,若用此一個if語句,當時出現了兩個bug,一是位址越界:即numa = polynomialsMultiply(numA, lowA, midA, numB, midB + 1, highB); 其中midB+1==1,highB==0,則遞歸不會出錯,但是在終止條件傳回時則出現位址越界,因為這時的lowB>hgihB>numB.size()-1,這個錯誤是由于采用一個統籌的if語句引起的,前面numA、numB都是low<high,而mid肯定是比high小,是以在分治遞歸上沒有問題,但是當一個滿足if條件而另一個不滿足時,這時仍然采用之前的mid計算方法,則就會使得mid==high,而分治遞歸時,mid+1>high導緻終止條件後面的位址越界。當然可用條件判斷來根據實際的low和high控制是否調用交叉相乘遞歸。但是本文采用了另外一種邏輯更清晰的結構,即歸并排序的while循環結構,整個遞歸程式設計如下:

vector<double> DivideMerge::polynomialsMultiply(vector<double>&numA, int lowA, int highA, vector<double>&numB, int lowB, int highB)
{
	vector<double>numXY,numXYZ;
	vector<double>numX, numY, numZ;              

	if (lowA < highA && lowB<highB)
	{
		int midA = (highA + lowA) / 2;
		int midB = (highB + lowB) / 2;

		numX = polynomialsMultiply(numA, lowA, midA, numB, lowB, midB);

		vector<double> numa, numb;
		int dy;
		numa = polynomialsMultiply(numA, lowA, midA, numB, midB + 1, highB);
		numb = polynomialsMultiply(numA, midA + 1, highA, numB, lowB, midB);
		numY = crossAdd(numa, midA, highB, numb, highA, midB,dy);

		int dxy;
		numXY = crossAdd(numX, midA, midB, numY, 0, dy, dxy);
		numZ = polynomialsMultiply(numA, midA + 1, highA, numB, midB + 1, highB);               

		int dxyz;
		numXYZ = crossAdd(numXY, 0, dxy, numZ, highA, highB,dxyz);
	
	}
	else if (lowA<highA)
	{
		int midA = (highA + lowA) / 2;
		numX = polynomialsMultiply(numA, lowA, midA, numB, lowB, highB);                          //隻有一個就不用交叉
		numZ = polynomialsMultiply(numA, midA + 1, highA, numB, lowB, highB);

		int dxz;
		numXYZ = crossAdd(numX, midA, highB, numZ, highA, highB,dxz);
	}
	else if (lowB < highB)
	{
		int midB = (lowB+highB) / 2;
		numX = polynomialsMultiply(numA, lowA, highA, numB, lowB, midB);                          //隻有一個就不用交叉
		numZ = polynomialsMultiply(numA, lowA, highA, numB, midB+1, highB);

		int dxz;
		numXYZ = crossAdd(numX, highA, midB, numZ, highA, highB,dxz);
	}
	else
	{
		numXYZ.push_back(numA[lowA] * numB[lowB]);                                                 //調用的這個函數傳回的是一個void類型
	}

	return numXYZ;
}
           

        其次,4次相乘依次結束後,就需要對同幂次的項進行疊加,如何對所有遞歸回來的項進行正确的疊加是實作該算法的核心和關鍵,前前後後花了2個星期的思考才徹底解決了這個問題:(numX,numY和numZ 分别表示A0*B0,A0B1+B0A1和A1*B1遞歸相乘傳回的結果)

        (思考1)最先想到的是對numX的最後一項與numY的第一項相加,這樣僅僅試用對3×3的多項式,但當多項式較多時,numX的最後一項的幂次就與numY的幂次就不對應了。

        (思考2)根據4×4的例子,采用n/2的劃分來疊加,即元素項的總個數n/2對應numC = numX+numY,n/2+n/2對應于numC+numZ,測試用例得出,該方法也不能普遍試用所有的情況,當項數較多時,結果疊加不正确,也不适用于1×3的情況。

       那麼,到底應該如果疊加呢,當多項式的項數增多時,且兩個多項式的項數不一緻時,得出的numX,numY和numZ的情況僅僅通過n/2的方法進行疊加,你會發現會出現疊加錯誤,根本原因是這裡根據相乘結果的位置來疊加,而numX 和numY位置上的元素的幂次不能保證對應相等,比如numX最後一位為3次幂,而numY的第1項可能使2次幂,這樣就沒有對應幂次相加。到底該如何通用的疊加呢?這個問題一直思考了1個星期。

         最終回歸到問題的本身,依據什麼疊加,依據是同幂次才能相加,那麼隻要保證相加的數是同幂次,則就能保證疊加的正确性。那麼如何保證相加的數是同幂次呢?就要根據幂次在程式中的定義來程式設計。還記得本文開頭加黑的位序==幂次是程式設計實作的關鍵,多項式的存儲是按照幂次的高低依次存儲在vector中的。是以隻有根據總體的low和high的值就能準确計算出整個遞歸過程中的結果項和中間項的幂次,即通過vector的下标就能周遊該項,這樣再次印證了“資料善置,算法自成”,對資料能掌控的存和取,就能對資料進行心中想要的計算。

      幂次的資料疊加步驟如下:

      Step1. 對于給定的兩個vector資料項numa 和numb 的最高次幂,就能推到出numa中前面項的幂次,因為幂次都是按順序存儲在numa中的,若最高次幂為n,則前一項為n-1,而第一項則是n-numa.size()+1,這樣程式設計就能很容易實作:若numa的首項的幂次小于相對應的numb的首項。則當i的循環值還沒有到達numa和numb的同幂次項時,則先讓numa[i]進vector,這裡有一個同幂次基準下标就是numb的首項幂次,是以可以通過if判斷下标i(i定義為numa的下标)若小于低次幂的元素個數(低次幂的個數=numa的首項幂次-numa的首項幂次),就直接numa[i]進vector;而兩項相加則注意i的值小于numa.size(),避免下标越界。

最後的同幂次疊加算法如下:(注:該算法不針對特定的某個項,是一個普适的隻要知道元素項和最高次幂,是以不管是累加還是交叉疊加都可以)

std::vector<double> DivideMerge::crossAdd(vector<double> &numa, int midA, int highB, vector<double>&numb, int highA, int midB,int &d)  //這裡的int類型的參數都是最高次幂
{
	vector<double> numc;
	int i = 0;
	int j = 0;
	int na = midA + highB;                             //numa的最高幂次 2
	int nb = highA + midB;                             //numb的最高幂次 3
	int diffA = na - numa.size() + 1;
	int diffB = nb - numb.size() + 1;
	int flag = 0;                                     //相等就為0

	d = na > nb ? na : nb;                            //原則:傳回最高次幂

	if (diffA > diffB)
	{
		flag = 1;
	}
	else if (diffA < diffB)
	{
		flag = 2;
	}

	while (i < numa.size() || j < numb.size())
	{
		switch (flag)
		{
			case 0:{
					   if (i < numa.size() && j < numb.size()) 
					   {
							numc.push_back(numa[i] + numb[j]);
							i++;
							j++;
					   }
					   else if (i<numa.size())
					   {
						   numc.push_back(numa[i]);
						   i++;
					   }
					   else
					   {
						   numc.push_back(numb[j]);
						   j++;
					   }
					  
			}continue;

			case 1:{
					   if (j < diffA - diffB)
					   {
						   numc.push_back(numb[j]);
						   j++;
					   }
					   else if (j < numb.size())
					   {
						   numc.push_back(numa[i] + numb[j]);
						   i++;
						   j++;
					   }
					   else
					   {
						   numc.push_back(numa[i]);
						   i++;
					   }

			}continue;

			default:{

						if (i < diffB - diffA)
						{
							numc.push_back(numa[i]);
							i++;
						}
						else if (i<numa.size())
						{
							numc.push_back(numa[i] + numb[j]);
							i++;
							j++;
						}
						else
						{
							numc.push_back(numb[j]);
							j++;
						}
			}continue;
		}//switch
	}//while

	return numc;
}
           

       總結:本文基于分治法的多項式相乘算法的時間複雜度仍為O(n2),而降低時間複雜度的方法就是采用高斯方法,即減少相乘的次數,這個将在下一篇文中寫到。

      在疊加算法中,第一次采用了switch... case ... continue的搭配,這是由于外層有循環條件,是以具體用break 還是continue 要根據實際的功能需求而定,break跳出整個循環,而continue則僅僅是終止接下來語句的執行,并不一定是循環,思維一定要開闊。在測試中,僅僅針對于項數2×3,4×4,1×3來測試 或8×8等,這樣的測試用例沒有全考慮涵蓋,是以最後根據1×2、1×3、1×4... 2×2、2×3、2×4 等等,根據周遊的方法既能都考慮到,又避免重複的全測試方法。

繼續閱讀