找回密码
 欢迎注册
楼主: 无心人

[讨论] 二进制32位整数快速平方根

[复制链接]
 楼主| 发表于 2009-2-12 17:54:17 | 显示全部楼层
FWAIT是
    Causes the processor to check for and handle pending, unmasked, floating-point exceptions before proceeding.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-2-12 19:54:37 | 显示全部楼层
在386时代
FWAIT是等待FST FSTP FISTP写入内存完毕的指令

我想,P4时代应该不必了
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-2-12 20:01:26 | 显示全部楼层
重新回到整数算法上来

假设bsr方法得到初始值

因为,初始值是个2^k类型的值
所以首次的牛顿迭代是无需乘法和除法的

现在的问题是,首次迭代的值和精确值的差距是多少
二次迭代是否能得到准确值??

如果不能得到
那不能得到精确值的输入是否有范围?
范围多大?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-2-12 20:19:23 | 显示全部楼层
网上找到的,似乎有用处
http://tech.ddvip.com/2008/08/121928805856259.html

为工作的需要,要在单片机上实现开根号的操作。目前开平方的方法大部分是用牛顿

  迭代法。我在查了一些资料以后找到了一个比牛顿迭代法更加快速的方法。

  1.原理

  因为排版的原因,用pow(X,Y)表示X的Y次幂,用B[0],B[1],...,B[m-1]表示一个序列,

  其中[x]为下标。

  假设:

  B[x],b[x]都是二进制序列,取值0或1。

  M = B[m-1]*pow(2,m-1) + B[m-2]*pow(2,m-2) + ... + B[1]*pow(2,1) + B[0]*pow

  (2,0)

  N = b[n-1]*pow(2,n-1) + b[n-2]*pow(2,n-2) + ... + b[1]*pow(2,1) + n[0]*pow

  (2,0)

  pow(N,2) = M

  (1) N的最高位b[n-1]可以根据M的最高位B[m-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=k-1,n=k=m/2

  所以b[n-1]完全由B[m-1]决定。

  余数 M[1] = M - b[n-1]*pow(2, 2*n-2)

  (2) N的次高位b[n-2]可以采用试探法来确定。

  因为b[n-1]=1,假设b[n-2]=1,则 pow(b[n-1]*pow(2,n-1) + b[n-1]*pow(2,n-2),

  2) = b[n-1]*pow(2,2*n-2) + (b[n-1]*pow(2,2*n-2) + b[n-2]*pow(2,2*n-4)),

  然后比较余数M[1]是否大于等于 (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4)。这种

  比较只须根据B[m-1]、B[m-2]、...、B[2*n-4]便可做出判断,其余低位不做比较。

  若 M[1] >= (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4), 则假设有效,b[n-2] =

  1;

  余数 M[2] = M[1] - pow(pow(2,n-1)*b[n-1] + pow(2,n-2)*b[n-2], 2) = M[1] -

  (pow(2,2)+1)*pow(2,2*n-4);

  若 M[1] < (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4), 则假设无效,b[n-2] =

  0;余数 M[2] = M[1]。

  (3) 同理,可以从高位到低位逐位求出M的平方根N的各位。

  使用这种算法计算32位数的平方根时最多只须比较16次,而且每次比较时不必把M的各位逐

  一比较,尤其是开始时比较的位数很少,所以消耗的时间远低于牛顿迭代法。

  2. 流程图

  (制作中,稍候再上)

  3. 实现代码

  这里给出实现32位无符号整数开方得到16位无符号整数的C语言代码。
  1. /****************************************/
  2. /*Function: 开根号处理         */
  3. /*入口参数:被开方数,长整型      */
  4. /*出口参数:开方结果,整型       */
  5. /****************************************/
  6. unsigned int sqrt_16(unsigned long M)
  7. {
  8.   unsigned int N, i;
  9.   unsigned long tmp, ttp; // 结果、循环计数
  10.   if (M == 0)       // 被开方数,开方结果也为0
  11.     return 0;
  12.   N = 0;
  13.   tmp = (M >> 30);     // 获取最高位:B[m-1]
  14.   M <<= 2;
  15.   if (tmp > 1)       // 最高位为1
  16.   {
  17.     N ++;        // 结果当前位为1,否则为默认的0
  18.     tmp -= N;
  19.   }
  20.   for (i=15; i>0; i--)   // 求剩余的15位
  21.   {
  22.     N <<= 1;       // 左移一位
  23.     tmp <<= 2;
  24.     tmp += (M >> 30);  // 假设
  25.     ttp = N;
  26.     ttp = (ttp<<1)+1;
  27.     M <<= 2;
  28.     if (tmp >= ttp)   // 假设成立
  29.     {
  30.       tmp -= ttp;
  31.       N ++;
  32.     }
  33.   }
  34.   return N;
  35. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-2-12 20:43:29 | 显示全部楼层
/****************************************/
/*Function: 开根号处理         */
/*入口参数:被开方数,长整型      */
/*出口参数:开方结果,整型       */
/****************************************/
#include <stdio.h>
#include <stdlib.h>

unsigned int sqrt_16(unsigned int M)
{
  unsigned int N, i;
  unsigned int tmp, ttp;  // 结果、循环计数
  if (M == 0)    // 被开方数,开方结果也为0
    return 0;
  N = 0;
  tmp = (M >> 30);  // 获取最高位:B[m-1]
  M <<= 2;
  if (tmp > 1)    //最高位为1
  {
    N ++;  // 结果当前位为1,否则为默认的0
    tmp -= N;
  }

  for (i=15; i>0; i--)   // 求剩余的15位
  {
    N <<= 1;   // 左移一位
    tmp <<= 2;
    tmp += (M >> 30);  // 假设
    ttp = N;
    ttp = (ttp<<1)+1;
    M <<= 2;
    if (tmp >= ttp)  // 假设成立
    {
      tmp -= ttp;
      N ++;
    }
  }
  return N;
}

int main(void)
{
  unsigned int M, N;
  printf("Input M:");
  scanf("%u", &M);
  N = sqrt_16(M);
  printf("\nSQRT(%u) = %u\n", M, N);
  return 0;
}

写了个测试程序
测试了几个例子
E:\MinGW\MSYS\math>sqrt
Input M:1023

SQRT(1023) = 31

E:\MinGW\MSYS\math>sqrt
Input M:4000000000

SQRT(4000000000) = 63245

E:\MinGW\MSYS\math>sqrt
Input M:0

SQRT(0) = 0

E:\MinGW\MSYS\math>sqrt
Input M:1

SQRT(1) = 1

E:\MinGW\MSYS\math>sqrt
Input M:23

SQRT(23) = 4

E:\MinGW\MSYS\math>sqrt
Input M:400000000

SQRT(400000000) = 20000
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2009-2-12 21:21:49 | 显示全部楼层
很奇怪的现象:
虽说我的43#代码够简洁,
但测试结果速度并不理想,基本上比无心人和liangbch写的函数要慢一倍。

我又测试了HugeCalc内部采用的纯C/C++版本(不含浮点计算),
虽然语句比较多,还有跳转语句(不见得编译后就有跳转指令),
但测试起来速度反而比用fsqrt的都快很多很多!

不信大家可以测试测试:
  1. UINT32 Sqrt( UINT32 n )
  2. {
  3.     UINT32 t;
  4.     UINT32 ret = 0;

  5.     #define SQRT_CORE(x)        \
  6.         t = ret + (1UL<<(x));   \
  7.         ret >>= 1;              \
  8.         if ( n >= t )           \
  9.         {                       \
  10.             n -= t;             \
  11.             ret += (1UL<<(x));  \
  12.         }

  13.     SQRT_CORE(30);
  14.     SQRT_CORE(28);
  15.     SQRT_CORE(26);
  16.     SQRT_CORE(24);
  17.     SQRT_CORE(22);
  18.     SQRT_CORE(20);
  19.     SQRT_CORE(18);
  20.     SQRT_CORE(16);
  21.     SQRT_CORE(14);
  22.     SQRT_CORE(12);
  23.     SQRT_CORE(10);
  24.     SQRT_CORE(8);
  25.     SQRT_CORE(6);
  26.     SQRT_CORE(4);
  27.     SQRT_CORE(2);
  28.     SQRT_CORE(0);

  29.     #undef SQRT_CORE

  30.     return ret;
  31. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-2-12 21:30:28 | 显示全部楼层
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-2-12 21:34:53 | 显示全部楼层
56#的似乎类似于54#的算法

呵呵

因为不存在乘法和除法
只有逻辑移位和比较
所以比较快
===================
fsqrt的大概58-60个时钟
郭老大的代码最少32个时钟吧

有时间依照54#的逻辑写个汇编
呵呵
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2009-2-12 22:28:53 | 显示全部楼层
他这个算法是固定15次循环,我想这样做。
  1.做一个256 byte 的表。
  2. 如果被开放数小于8bit 直接查表得到结果。
    否则如果被开方数是k(<8k<=32)bit,则首先通过一次得到最高结果的最高4bit,然后循环 (k-8)/2次循环得到结果的其他bit.
  3.程序仍然采用先前的结果,bsr 加跳表。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-2-13 09:29:53 | 显示全部楼层
测试了

FPU2, 54#, 56#代码
2100000000 -- 2200000000

FPU2 Version: 3587.139 ms
iSqrt_16 Version: 26195.432 ms
iSqrt_GxQ Version: 16059.504 ms
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-29 07:39 , Processed in 0.049778 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表