1、.單片機C語言求平方根函數2009年07月12日 星期日 11:01轉自：在單片機中要開平方.可以用到下面算法: 算法1: 本算法只采用移位、加減法、判斷和循環實現，因為它不需要浮點運算，也不需要乘除運算，因此可以很方便地運用到各種芯片上去。 我們先來看看10進制下是如何手工計算開方的。 先看下面兩個算式， x = 10*p + q (1) 公式(1)左右平方之后得： x2 = 100*p2 + 20pq + q2 (2) 現在假設我們知道x2和p，希望求出q來，求出了q也就求出了x2的開方x了。 我們把公式(2)改寫為如下格式： q = (x2 - 100*p2)/(20*p+q) (3)

2、這個算式左右都有q，因此無法直接計算出q來，因此手工的開方算法和手工除法算法一樣有一步需要猜值。 我們來一個手工計算的例子：計算1234567890的開方 首先我們把這個數兩位兩位一組分開，計算出最高位為3。也就是(3)中的p，最下面一行的334為余數，也就是公式(3)中的(x2 - 100*p2)近似值 3 - | 12 34 56 78 90 9 - | 3 34 下面我們要找到一個0-9的數q使它最接近滿足公式(3)。我們先把p乘以20寫在334左邊： 3 q - | 12 34 56 78 90 9 - 6q| 3 34 我們看到q為5時(60+q*q)的值最接近334，而且不超過33

3、4。于是我們得到： 3 5 - | 12 34 56 78 90 9 - 65| 3 34 | 3 25 - 9 56 接下來就是重復上面的步驟了，這里就不再啰嗦了。 這個手工算法其實和10進制關系不大，因此我們可以很容易的把它改為二進制，改為二進制之后，公式(3)就變成了： q = (x2 - 4*p2)/(4*p+q) (4) 我們來看一個例子，計算100(二進制1100100)的開方： 1 0 1 0 - | 1 10 01 00 1 - 100| 0 10 | 0 00 - | 10 011001| 10 01 - 0 00 這里每一步不再是把p乘以20了，而是把p乘以4，也就是把p右

4、移兩位，而由于q的值只能為0或者1，所以我們只需要判斷余數(x2 - 4*p2)和(4*p+1)的大小關系，如果余數大于等于(4*p+q)那么該上一個1，否則該上一個0。 下面給出完成的C語言程序，其中root表示p，rem表示每步計算之后的余數，divisor表示(4*p+1)，通過a>>30取a的最高 2位，通過a<<=2將計算后的最高2位剔除。其中root的兩次<<1相當于4*p。程序完全是按照手工計算改寫的，應該不難理解。 unsigned short sqrt(unsigned long a) unsigned long rem = 0; unsi

5、gned long root = 0; unsigned long divisor = 0; for(int i=0; i<16; i+) root <<= 1; rem = (rem << 2) + (a >> 30); a <<= 2; divisor = (root<<1) + 1; if(divisor <= rem) rem -= divisor; root+; return (unsigned short)(root); 算法2 :單片機開平方的快速算法 因為工作的需要，要在單片機上實現開根號的操作。目前開平方

6、的方法大部分是用牛頓 迭代法。我在查了一些資料以后找到了一個比牛頓迭代法更加快速的方法。不敢獨享，介 紹給大家，希望會有些幫助。 1.原理 因為排版的原因，用pow(X,Y)表示X的Y次冪，用B0，B1，.，Bm-1表示一個序列， 其中x為下標。 假設： Bx,bx都是二進制序列,取值0或1。 M = Bm-1*pow(2,m-1) + Bm-2*pow(2,m-2) + . + B1*pow(2,1) + B0*pow (2,0) N = bn-1*pow(2,n-1) + bn-2*pow(2,n-2) + . + b1*pow(2,1) + n0*pow (2,0) pow(N,2) =

7、 M (1) N的最高位bn-1可以根據M的最高位Bm-1直接求得。 設 m 已知,因為 pow(2, m-1) <= M <= pow(2, m)，所以 pow(2, (m-1)/2) <= N <= pow(2, m/2) 如果 m 是奇數，設m=2*k+1, 那么 pow(2,k) <= N < pow(2, 1/2+k) < pow(2, k+1), n-1=k, n=k+1=(m+1)/2 如果 m 是偶數，設m=2k, 那么 pow(2,k) > N >= pow(2, k-1/2) > pow(2, k-1), n-1=

8、k-1,n=k=m/2 所以bn-1完全由Bm-1決定。 余數 M1 = M - bn-1*pow(2, 2*n-2) (2) N的次高位bn-2可以采用試探法來確定。 因為bn-1=1，假設bn-2=1，則 pow(bn-1*pow(2,n-1) + bn-1*pow(2,n-2), 2) = bn-1*pow(2,2*n-2) + (bn-1*pow(2,2*n-2) + bn-2*pow(2,2*n-4), 然后比較余數M1是否大于等于 (pow(2,2)*bn-1 + bn-2) * pow(2,2*n-4)。這種 比較只須根據Bm-1、Bm-2、.、B2*n-4便可做出判斷，其余低位

9、不做比較。 若 M1 >= (pow(2,2)*bn-1 + bn-2) * pow(2,2*n-4), 則假設有效，bn-2 = 1； 余數 M2 = M1 - pow(pow(2,n-1)*bn-1 + pow(2,n-2)*bn-2, 2) = M1 - (pow(2,2)+1)*pow(2,2*n-4)； 若 M1 < (pow(2,2)*bn-1 + bn-2) * pow(2,2*n-4), 則假設無效，bn-2 = 0；余數 M2 = M1。 (3) 同理，可以從高位到低位逐位求出M的平方根N的各位。 使用這種算法計算32位數的平方根時最多只須比較16次，而且每次比較

10、時不必把M的各位逐 一比較，尤其是開始時比較的位數很少，所以消耗的時間遠低于牛頓迭代法。 2. 實現代碼 這里給出實現32位無符號整數開方得到16位無符號整數的C語言代碼。 - - /*/ /*Function: 開根號處理 */ /*入口參數：被開方數，長整型 */ /*出口參數：開方結果，整型 */ /*/ unsigned int sqrt_16(unsigned long M) unsigned int N, i; unsigned long tmp, ttp; / 結果、循環計數 if (M = 0) / 被開方數，開方結果也為0 return 0; N = 0; tmp = (M >> 30); / 獲取最高位：Bm-1 M <<= 2; if (tmp > 1) / 最高位為1 N +; / 結果當前位為1，否則為默認的0 tmp -= N; for (i=15; i>0; i-) / 求剩余的15位 N <