非恢复浮点数平方根算法

Question

我正在尝试使用非恢复算法计算浮点数的平方根。

例如，假设x = 1001，则平方根为31.6386。

我想使用非恢复方法计算此平方根。

我尝试按照论文中的方法进行： Implementation of Single Precision Floating Point Square Root on FPGAs 但似乎我的结果略有偏差。我无法弄清原因。

例如，我编写的程序将生成以下结果：

correct_result =
  41FD1BD2

myresult =  
  41FD1BD1

error =    
    1.192093e-007

C++代码的版本 :

#include <iostream>
#include <cmath>

using namespace std;

  union newfloat{
    float f;
    int i;
  };

int main () {
// Input number
newfloat x;
cout << "Enter Number: ";
cin >> x.f;

// Pull out exponent and mantissa
int exponent = (x.i >> 23) & 0xFF;
int mantissa = (x.i & 0x7FFFFF) | ((exponent && exponent) << 23);

// Calculate new exponent
int new_exponent = (exponent >> 1) + 63 + (exponent & 1);


// Shift right (paper says shift left but shift left doesn't work?)
if (exponent & 1) {
    mantissa = mantissa  >> 1;
    cout << " Shifted right " << endl;
}

// Create an array with the bits of the mantissa
unsigned int D [48];
for (int i = 47; i >= 0; i--) {
  if (i >= 24) {
    D[i] = (mantissa >> (i-24)) & 1;
  } else {
    D[i] = 0;
  }
}


// == Perform square root ==
// Set q24 = 0, r24 = 0 and then iterate from k = 23 to 0
int q[25] = {0}; // 25 element array, indexing ends at 24
int r[25] = {0};

for (int k = 23; k >= 0; k--) {
    if (r[k+1] >= 0) {
        r[k] = ((r[k+1] << 2) | (D[2*k+1] << 1) | D[2*k] ) - (q[k+1] << 2 | 1 );
        } else {
        r[k] = ((r[k+1] << 2) | (D[2*k+1] << 1) | D[2*k] ) + (q[k+1] << 2 | 0x3 );
        } 

    if (r[k] >= 0) {
        q[k] = (q[k+1] << 1) | 1;
        } else {
        q[k] = q[k+1] << 1;
    }

    if (k == 0) {
        if (r[0] < 0) {
            r[0] = r[0] + (q[0] << 1) | 1;
        }
    }
}

// Create quotient from LSBs of q[]
int Q = 0;
for (int i = 0; i <= 23; i++) {
    Q = Q | ((q[i] & 1) << i);
}

// Option 1 Rounding
//if (r[0] > 0) // Works for 10, 1001, 1021, but not 1012
// Q = Q + 1;

// Option 2 Rounding (No rounding)
// Works for 1012, Doesn't work for 10, 1001, 1021

// Option 3 Rounding (Calculate the next 3 Quotient bits to get a guard round and sticky bit)

// Calculate correct result:
newfloat correct_result;
correct_result.f = sqrt(x.f);

// Form my result into a single number
newfloat myresult;
myresult.i = (new_exponent << 23) | (Q & 0x7FFFFF);

// Print results
cout << hex << "My result: " << myresult.i << endl;
cout << hex << "Correct:   " <<  correct_result.i << endl;
return 0;
}

- Veridian

---

我不知道答案，但强烈怀疑检查您的代码中的错误并帮助您处理该部分会比编写关于技术的长篇教程容易得多（而您已经熟悉论文中的基础知识）。请展示您代码的相关部分。 - Neil Slater

你的工作基础是哪篇论文还不清楚。你看过这篇论文吗：Anuja Nanhe, Gaurav Gawali, Shashank Ahire和K. Sivasankaran：Implementation of Fixed and Floating Point Square Root Using Nonrestoring Algorithm on FPGA。《国际计算机与电气工程杂志》2013年10月第5卷第5期，533-537页。 - njuffa

@njuffa，我已经更新了论文链接，它在问题中有说明。我看过那篇论文，但由于糟糕的英语，算法有些令人困惑。此外，我不想开始关注另一篇论文。 - Veridian

“convertnum2fpparts”、“combinefpparts”等函数的代码在哪里？你的代码无法运行。你检查过该函数的输出是否正确吗？另一件需要注意的事情是，如果使用double进行计算而你想要单精度，这看起来像是你尝试过了，但这可能会很棘手。在其中添加一些对“whos”或“class”的调用以进行检查。另一个需要注意的问题是“0”与“1”索引错误。你的代码将“a=1001”，但你的问题涉及到一个值为“1000”，那么用哪个值生成所示的十六进制值？ - horchler

@horchler，我添加了缺失的代码。我已经检查了你提到的那些输出，它们是正确的。我也想让代码设置a = 1001。一些值提供了正确的结果，一些错误，所以代码将a = 1001，因为它是显示错误结果的值。代码确实可以运行，请告诉我是否遇到问题。我正在运行MATLAB R2010a。更新版本可能会出现问题，因为它在新版本中使用不同的参数进行bitor命令。 - Veridian

Answer 1

首先让我从论文中突出相关部分：

你需要重新审视加减法的处理方式。你的代码在使用普通的双精度数字进行计算，但我认为该算法是针对整数模运算设计的。

因此，如果你查看论文后面列出的示例，0011 - 0101 的计算会回绕到 1110。

这可能可以解释为什么你得到了错误的结果，我想 :)

- Amro

---

我添加了一个使用整数算术的 C++ 版本的代码，但仍然没有得到正确的结果。也许是舍入问题？有什么想法吗？ - Veridian

@starbox：我认为你误解了论文中使用的符号（我不怪你，它不容易阅读）... 首先，r(k)和R(k)之间存在差异（q和Q同理）。大写字母表示单个位，而小写字母表示一系列位... - Amro

@starbox：我尝试按照我理解的方式实现了该算法（整数情况，而非浮点数情况）。它对许多输入给出了正确的结果，但也有错误的结果... 我也不确定我做错了什么！如果您愿意，我可以发布我的版本。 - Amro

Answer 2

我今天查看了你的程序的c++版本并阅读了相关文档。据我所知，该算法旨在提供商和余数两个结果。就提供的示例而言，他使用自己的算法来获得127的平方根，提供了11+R6的结果。11²+6=127。

虽然这是对整数的操作，但每种数据类型都有其精度限制。这让我相信你的程序正在按预期执行，只是你已经失去了精度，至少对于计算平方根的方式和所使用的数据类型而言。我认为你会发现你失去的精度在r[0]中。

我从代码的注释中看到，你打算或尝试计算额外的精度。这似乎是一个合理的尝试路径。请注意，除了需要进行的其他更改之外，您还必须删除（或移动）检查k == 0；因为它修改了余数，这将破坏循环。

我认为真正的问题是你能接受多大的精度误差。例如，C++的sqrt函数（以及你的函数）在sqrt(2)上的误差为0.00000002。没有人似乎在意。考虑到你编写的程序与C++的sqrt函数相比，在不匹配的情况下的误差要小于这个值。我花了大部分时间拆分它，测试各个部分，并审查了主题，但没有发现明显的错误。对我来说，它似乎足够接近政府工作的标准。

- Taekahn