1. 中值定理
设abc为顺序排列间距为P的3个整数,ABC是他们的平方,则有:
b^2 = ( a^2 + c^2 ) / 2 - R,
即: B = (A+C)/2 - R
其中, 修正值 R = P^2
证明:
a = b + P, c = b- P
=> a^2 = (b+P)^2 , c^2 = (b-P)^2
=> a^2 + b^2 = 2b^2 + 2P^2
=> b^2 = ( a^2 + c^2)/2 - P^2
特别的,
当 P=2^n = 2*2^(n-1) = 2Pn-1,
则 R = P^2 = 2^(2n) = (2Pn-1)^2=4(Pn-1)^2=4Rn-1
2. 笔算开均方
被开放数 A , A = (10a + b)^2
(10a + b)^2 = 100a^2 + 20ab + b^2 = 100a^2 + b(20a+b)
a 代表已经开方出来的结果的高位,b代表当前需要计算的位上的数据
计算过程中,100a^2被减掉,剩下b(20a+b), 这里需要做的就是找到最大的整数b', 使得b'(20a+b')=<b(20a+b)
<<九章算术>>中讲了开方,开立方的方法, 计算步骤和现在的基本一样.
例如:
2 3 4 _____________________ ^_/ 5, 47, 56 2 4 ------------------------ 1 4 7 20*2+3 1 2 9 ------------------------- 1 8 56 20*23+4 1 8 56 ------------------------ 0
二进制开5^2 = 25
1 0 1 _____________________ ^_/ 1, 10, 01 1*1 1 ------------------------ 10 (2*1+0)*0 00 ------------------------- 1 0 01 (2*10+1)*1 1 0 01 ------------------------ 0
3. 2进制数的开均方根
一个32bit的整数开方后的数为16bit, 这16bit由1和0组成, 因此可以从高位到低位逐个逐个求出每一bit,即从16bit开始试根1,然后比较大小,小于或等于被开放数,则试根1正确,否则试根为0
相关程序参考:
程序1
// 0x 04 00 00 00 0l -> 0x40 00 00 00 #define step(shift) \ if((0x40000000 >> shift) + sqrtVal <= val) \ { \ val -= (0x40000000 >> shift) + sqrtVal; \ sqrtVal = (sqrtVal >> 1) | (0x40000000 >> shift); \ } \ else \ { \ sqrtVal = sqrtVal >> 1; \ }
unsigned long xxgluSqrtFx(unsigned long val) { // Note: This fast square root // only works with an even Q_FACTOR unsigned long sqrtVal = 0; step(0); step(2); step(4); step(6); step(8); step(10); step(12); step(14); step(16); step(18); step(20); step(22); step(24); step(26); step(28); step(30); if(sqrtVal < val) { ++sqrtVal; } // sqrtVal <<= (Q_FACTOR)/2; return(sqrtVal); }
程序2
// (root << 1 + b)<< bshft-- static unsigned julery_isqrt(unsigned long val) { unsigned long temp, g=0, b = 0x8000, bshft = 15; // g 均方根 do { // if ( val >= ( temp = ( ((g << 1) + b)<<(bshft--) ) ) ) // bshft 先參與運算,然後自減1 // 等效與下式 temp = g << 1; // divisor = ROOT*2 temp += b; // divisor += B (1000 0000) temp <<= bshft; // 左邊移動 bshft位 衝15開始 bshft--; // if( val >= temp) // VAL ? divisor { // val -= temp; // g += b; // ROOT += B; } // } while (b >>= 1); // B = B/2 return g; }
程序3
// q = (x^2 - 4*p^2)/(4*p+q) (4) // root = p, rem=yushu, divsior = 4*p+1 // a >>30 get high 2bit // a <<=2 cut off high 2bit // root 2 timers <<1, = 4*P unsigned short sqrt_1(unsigned long a) { unsigned long rem = 0; // 餘數 unsigned long root = 0; // 均方根 unsigned long divisor = 0; unsigned char i; for(i=0; i<16; i++) { root <<= 1; // 均方根*2 rem = ((rem << 2) + (a >> 30)); // 餘數 = 上次結果 + 本次取被開方數高2位 a <<= 2; // 移除高2位 divisor = (root<<1) + 1; // 均方根*2+1 if(divisor <= rem) // 均方根試商1 與 被開方數比較 { // >=, 試商1正確 rem -= divisor; // 新的餘數 root++; // 根 + 1 } } return (unsigned short)(root); }
程序4
// 如果對一個32位元無符號數開方,那麼結果一定是一個16位元無符號數。 // 現設被開方的數為a,開方結果: // b = b[15] * 2^15 + b[14] * 2^14 + ... + b[0] * 2^0 // 上式中的b[15]表示16位元無符號數b的第15位(最高位),2^15表示2的15次方;以此類推。 // 開方的基本原理如下: // 首先假設b[15] = 1,此時如果(b[15] * 2^15)^2 <= a成立,則b[15]=1,否則b[15]=0; // 然後再假設b[14]=1,此時如果(b[15] * 2^15 + b[14] * 2^14)^2 <= a成立,則b[14]=1,否則b[14]=0。 // 以此類推,經過16次比較之後可以得出最終的開方結果。 unsigned int MUXfixed_sqrt(unsigned long a) { unsigned char i; unsigned int ROOT=0; unsigned long c; for(i=0; i<16; i++) { c = (ROOT | (1 << (15 - i))); if((c * c) <= a) { ROOT = (unsigned int)c; } } return ROOT; }
程序5
unsigned int SFsqrt_fixed(unsigned long a) { unsigned long i,c; unsigned long ROOT = 0; for(i = 0x40000000; i != 0; i >>= 2) { c = i + ROOT; ROOT >>= 1; if(c <= a) { a -= c; ROOT += i; } } return (unsigned int)ROOT; }
程序6
// 遞送發+尋找最高位 // A = X ^ 2 // Xn+1 = ( Xn + A/Xn)/2 unsigned int Sqrt_RC(unsigned long a) { unsigned int ROOT = 0; unsigned long b = 65536; unsigned char i; for(i=0; i<32; i++) { if( a & (1 << (32 - i))) { ROOT = (b + a/b)/2; ROOT = (ROOT + a/ROOT)/2; ROOT = (ROOT + a/ROOT)/2; ROOT = (ROOT + a/ROOT)/2; ROOT = (ROOT + a/ROOT)/2; break; } else { b >= 1; } } return (unsigned int)ROOT; }
// 若使用 二分法,设置合适的初始值,可以缩短递送次数 // f(x) = x^2 - A, a 就是 f(x) = 0的根 // f(n) < 0, f(m) >0, a = (n , m) // 若 f[(m+m)/2] <0 a = ( (m+m)/2, m) // 否則 f[(m+m)/2] >0 a= (n, (m+m)/2) // 反復遞送
|