EEPW论坛

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)
// 反復遞送