Apple 的定点平方根算法背后的数学理论是什么

【问题标题】：What's the math theory behind Apple's square root of fixed point algorithmApple 的定点平方根算法背后的数学理论是什么
【发布时间】：2017-10-29 12:41:04
【问题描述】：

我正在一个项目中工作，我需要计算一个 16 位数字的平方根。我正在使用不支持浮点操作的 NXP 微控制器。在网上查找我发现了以下方法Fixed-Point Square Root by Ken Turkowsky from Apple。我还在 WEB 上找到了相同的算法，并进行了一些修改，包括：Looking for an efficient integer square root algorithm for ARM Thumb2 我在这里复制代码，以便在链接断开时可用：

    void put_space(unsigned int number)
{
    int i;
    char spaces = 12;
    char buffer[15];
    char number_of_chars = sprintf(buffer,"%u",number);
    for(i=0;i<spaces-number_of_chars;i++)
    {
        printf(" ");    
    }   
}
void print_binary(unsigned int number)
{
    int i;
    unsigned char bit;
    for(i=0;i<32;i++)
    {
        if(i%8 == 0 && i!=0)
        {
            printf(" ");
        }
        bit = (number & (2147483648 >> i)) == (2147483648 >> i);
        printf("%u",bit);
    }
}
void print_number(unsigned int number)
{
    print_binary(number);
    printf(" ");
    printf("(%u)",number);
    put_space(number);
}


typedef signed int TFract; /* 2 integer bits, 30 fractional bits */
    TFract FFracSqrt(TFract x)
    {
        register unsigned int root, remHi, remLo, testDiv, count;
        root = 0; /* Clear root */ printf("root: "); print_number(root);        
        remHi = 0; /* Clear high part of partial remainder */  printf("remHi: "); print_number(remHi);
        remLo = x; /* Get argument into low part of partial remainder */ printf("remLo: "); print_number(remLo);
        count = 30; /* Load loop counter */

        printf("\n\n");
    do 
    {
        remHi = (remHi<<2) | (remLo>>30); printf("remHi: "); print_number(remHi);
        remLo <<= 2; /* get 2 bits of arg */ printf("remLo: "); print_number(remLo);
        root <<= 1; /* Get ready for the next bit in the root */ printf("root: "); print_number(root);
        testDiv = (root << 1) + 1; /* Test radical */ printf("testDiv: "); print_number(testDiv);       
        if (remHi >= testDiv) 
        {
            remHi -= testDiv; printf("remHi: "); print_number(remHi);
            root++; printf("root: "); print_number(root);
        }
        printf("\n");
    } while (count-- != 0);

    printf("\n\nResult: %u\n\n",root);  
    return(root);
}

我可以理解算法中发生了什么，包括移位和 OR 操作。但是，我无法理解这段代码背后的数学原理。 FFracSqrt 接收一个“x”整数作为参数并返回 floor(sqrt(x)*32768)。这是怎么发生的？ 请帮助我和其他人理解它，因为到目前为止我找不到任何关于该理论的信息。 真诚的，

我在这里添加输入为 2 时的输出： 在 do-while 循环之前 根：00000000 00000000 00000000 00000000 (0)
remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000000 00000010 (2)

进入do-while循环 remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000000 00001000 (8)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000000 00100000 (32)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000000 10000000 (128)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000010 00000000 (512)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00001000 00000000 (2048)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00100000 00000000 (8192)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 10000000 00000000 (32768)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000010 00000000 00000000 (131072)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00001000 00000000 00000000 (524288)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00100000 00000000 00000000 (2097152)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 10000000 00000000 00000000 (8388608)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000010 00000000 00000000 00000000 (33554432)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0) remLo: 00001000 00000000 00000000 00000000 (134217728)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00100000 00000000 00000000 00000000 (536870912)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 10000000 00000000 00000000 00000000 (2147483648)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000010 (2)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1) 如果条件为真
remHi: 00000000 00000000 00000000 00000001 (1)
根：00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000100 (4)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000000 00000010 (2)
testDiv: 00000000 00000000 00000000 00000101 (5)

remHi: 00000000 00000000 00000000 00010000 (16)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000000 00000100 (4)
testDiv: 00000000 00000000 00000000 00001001 (9) 如果条件为 TRUE remHi: 00000000 00000000 00000000 00000111 (7)
根：00000000 00000000 00000000 00000101 (5)

remHi: 00000000 00000000 00000000 00011100 (28)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000000 00001010 (10)
testDiv: 00000000 00000000 00000000 00010101 (21) 如果条件为真
remHi: 00000000 00000000 00000000 00000111 (7)
根：00000000 00000000 00000000 00001011 (11)

remHi: 00000000 00000000 00000000 00011100 (28)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000000 00010110 (22)
testDiv: 00000000 00000000 00000000 00101101 (45)

remHi: 00000000 00000000 00000000 01110000 (112)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000000 00101100 (44)
testDiv: 00000000 00000000 00000000 01011001 (89) 如果条件为 TRUE remHi: 00000000 00000000 00000000 00010111 (23)
根：00000000 00000000 00000000 00101101 (45)

remHi: 00000000 00000000 00000000 01011100 (92)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000000 01011010 (90)
testDiv: 00000000 00000000 00000000 10110101 (181)

remHi: 00000000 00000000 00000001 01110000 (368)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000000 10110100 (180)
testDiv: 00000000 00000000 00000001 01101001 (361) 如果条件为真
remHi: 00000000 00000000 00000000 00000111 (7)
根：00000000 00000000 00000000 10110101 (181)

remHi: 00000000 00000000 00000000 00011100 (28)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000001 01101010 (362)
testDiv: 00000000 00000000 00000010 11010101 (725)

remHi: 00000000 00000000 00000000 01110000 (112)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000010 11010100 (724)
testDiv: 00000000 00000000 00000101 10101001 (1449)

remHi: 00000000 00000000 00000001 11000000 (448)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00000101 10101000 (1448)
testDiv: 00000000 00000000 00001011 01010001 (2897)

remHi: 00000000 00000000 00000111 00000000 (1792)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00001011 01010000 (2896)
testDiv: 00000000 00000000 00010110 10100001 (5793)

remHi: 00000000 00000000 00011100 00000000 (7168)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00010110 10100000 (5792)
testDiv: 00000000 00000000 00101101 01000001 (11585)

remHi: 00000000 00000000 01110000 00000000 (28672)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 00101101 01000000 (11584)
testDiv: 00000000 00000000 01011010 10000001 (23169) 如果条件为真
remHi: 00000000 00000000 00010101 01111111 (5503)
根：00000000 00000000 00101101 01000001 (11585)

remHi: 00000000 00000000 01010101 11111100 (22012)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 01011010 10000010 (23170)
testDiv: 00000000 00000000 10110101 00000101 (46341)

remHi: 00000000 00000001 01010111 11110000 (88048)
remLo: 00000000 00000000 00000000 00000000 (0)
根：00000000 00000000 10110101 00000100 (46340)
testDiv: 00000000 00000001 01101010 00001001 (92681)

结果：46340

【问题讨论】：

  • 你用铅笔和纸做数学题了吗？
  • 我做得更好。我修改了算法以显示二进制数字。正如我所说，我知道发生了什么，但我无法理解数学。
  • 维基百科有a good collection of square root methods。以上看起来很像二进制逐位方法，详细here。
  • 正如 John McFarlane 所指出的，此代码实现了在手持计算器出现之前曾在学校教授的具有数百年历史的长平方根算法。此变体以 2 为基数运行，因此每次迭代生成一位结果。

标签： microcontroller fixed-point square-root


【解决方案1】：

这是我对 cme​​ts 的尝试

TFract FFracSqrt(TFract x)
{
   unsigned root = 0; /* Clear root */

   //remHi & remLo forming a 64bit unsigned integer
   //remHiLo containing the remainder from the radicant not yet transfered into root
   unsigned remHi = 0;
   unsigned remLo = x;

   unsigned bitcount = 30; /* Load loop counter */

   do 
   {
       //shift left the remHiLo package by two bit positions
       remHi = (remHi<<2) | (remLo>>30);
       remLo <<= 2;

       //As we multiplied remHiLo with 4, it is fair to multiply the root with sqrt(4)
       root <<= 1;

       //We just shifted in a zero bit in root. Should it be set?

       //We should check whether (root+1) is too large

       //(root+1) * (root+1)
       //root^2 + 2*root + 1

       //We always deduct root^2 from remHi, so we need to check only the delta
       //between root^2 and (root+1)^2

       unsigned testDiv = (root << 1) + 1;

       if (remHi < testDiv)
       {
          //root+1 would overshoot the radicant, so leave everything as it is
       } 
       else
       {
           root++; //set LSB
           remHi -= testDiv; //update remHi to reflect the new bit in root
       }
    } while (bitcount-- != 0);
    return root;

【讨论】：

    猜你喜欢
    • 2018-04-11
    • 2018-03-19
    • 2011-11-04
    • 2020-02-15
    • 1970-01-01
    • 2019-01-04
    • 1970-01-01
    • 1970-01-01
    • 2017-07-06
    相关资源
    最近更新 更多
    热门标签
    Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode