sin函数创建CORDIC表中的舍入误差

Question

我试图以rust为例在这c实现中实现CORDIC方法，但是在生成表时存在舍入错误问题。这是我的代码和结果。

fn generate_table () {
    let pi: f32 = 3.1415926536897932384626;
    let k1: f32 = 0.6072529350088812561694; // 1/k
    let num_bits: uint = 32;
    let num_elms: uint = num_bits;
    let mul: f32 = 1073741824.000000; // 1<<(bits-2)

    println!("Cordic sin in rust");
    println!("num bits {}", num_bits);
    println!("mul {}", mul);
    println!("pi is {}", pi);
    println!("k1 is {}", k1);

    let shift: f32 = 2.0;
    for ii in range(0, num_bits) {
        let ipow: f32 = 1f32/shift.powi(ii as i32);
        let cur: f32 = ipow.atan() * mul;
        //println!("table values {:.10f}", cur);
        println!("table values 0x{}", std::f32::to_str_hex(cur));
    }
}

fn main() {
    generate_table();
}

给出下表，请注意第一个和最后一个值，以查看最大的错误。

table values 0x3243f6c0
table values 0x1dac6700
table values 0xfadbb00
table values 0x7f56ea8
table values 0x3feab78
table values 0x1ffd55c
table values 0xfffaab
table values 0x7fff55.8
table values 0x3fffea.c
table values 0x1ffffd.6
table values 0xfffff.b
table values 0x7ffff.f8
table values 0x40000
table values 0x20000
table values 0x10000
table values 0x8000
table values 0x4000
table values 0x2000
table values 0x1000
table values 0x800
table values 0x400
table values 0x200
table values 0x100
table values 0x80
table values 0x40
table values 0x20
table values 0x10
table values 0x8
table values 0x4
table values 0x2
table values 0x1
table values 0x0.8

我为什么要得到这些(四舍五入)？错误以及我如何修复它们？

Answer 1

快速的答案是:你有多确定你在两种计算中输入相同的输入？特别是，您引用的C实现是这样的：int mul = (1<<(bits-2));

当你有硬编码时：let mul: f32 = 1073741824.000000; // 1<<(bits-2)

注1:您已经将mul的类型从int更改为f32。

注2:在我运行程序时得到的输出中，我看到了以下内容：mul 1073741844

值得注意的是，这与您上面所写的硬编码常量不同；它是关闭20。

在我注意到上面的问题之前，我通常是这样调试一个问题的，在这个例子中，我实际上是这样做的，就是用打印每个中间表达式的值来测试代码的C和Rust版本，以确定事情开始不同的地方，从而缩小了引入“错误”的操作。

在这种情况下，它涉及并行修改C代码和Rust代码，以打印一个表，表中不仅包括i (或Rust版本中的ii )和输出c，还包括每个中间结果。

下面是其中每一个的代码，以及我最后得到的输出表。(但只是分析了这些表，我才意识到这两个mul值是不同的！)

C代码：

#include <stdio.h>

#include <math.h>
#define PI 3.1415926536897932384626
#define K1 0.6072529350088812561694

int main(int argc, char **argv)
{
    int i;
    int bits = 32; // number of bits
    int mul = (1<<(bits-2));


    int n = bits; // number of elements. 
    int c;

    printf("Cordic sin in C\n");
    printf("num bits %d\n", bits);
    printf("mul %d\n", mul);
    printf("pi is %g\n", PI);
    printf("k1 is %g\n", K1);

    float shift = 2.0;
    printf("computing c = atan(pow(2, -i)) * mul\n");
    printf("       i \t       c           p           a       c2\n");
    for(i=0;i<n;i++)
    {
        c = (atan(pow(2, -i)) * mul);
        int neg_i = -i;
        double p = pow(2, neg_i);
        double a = atan(p);
        int c2 = a * mul;;
        printf("% 8d \t 0x%08X % 12g % 12g 0x%08X\n", i, c, p, a, c2);
    }
}

锈码：

fn generate_table () {
    let pi: f32 = 3.1415926536897932384626;
    let k1: f32 = 0.6072529350088812561694; // 1/k
    let num_bits: uint = 32;
    let num_elms: uint = num_bits;
    let mul: f32 = 1073741824.000000; // 1<<(bits-2)

    println!("Cordic sin in rust");
    println!("num bits {}", num_bits);
    println!("mul {}", mul);
    println!("1 << (bits - 2): {}", (1i << (num_bits-2)) as f32);
    println!("pi is {}", pi);
    println!("k1 is {}", k1);

    let shift: f32 = 2.0;
    println!("computing c = (1f32/shift.powi(ii as i32)).atan() * mul");
    println!("       i \t       c           p            a        c2\n");
    for ii in range(0, num_bits) {
        let ipow: f32 = 1f32/shift.powi(ii as i32);
        let cur: f32 = ipow.atan() * mul;
        let a = ipow.atan();
        let c2 = a * mul;
        //println!("table values {:.10f}", cur);
        // println!("table values 0x{}", std::f32::to_str_hex(cur));
        println!("{:8u} \t 0x{:8s} {:12f} {:12f} 0x{:8s}",
                 ii,
                 std::f32::to_str_hex(cur),
                 ipow,
                 a,
                 std::f32::to_str_hex(c2),
                 );
    }
}

fn main() {
    generate_table();
}

生成的表格：

% gcc gentable2.c && ./a.out
Cordic sin in C
num bits 32
mul 1073741824
pi is 3.14159
k1 is 0.607253
computing c = atan(pow(2, -i)) * mul
       i           c           p           a       c2
       0     0x3243F6A8            1     0.785398 0x3243F6A8
       1     0x1DAC6705          0.5     0.463648 0x1DAC6705
       2     0x0FADBAFC         0.25     0.244979 0x0FADBAFC
       3     0x07F56EA6        0.125     0.124355 0x07F56EA6
       4     0x03FEAB76       0.0625    0.0624188 0x03FEAB76
       5     0x01FFD55B      0.03125    0.0312398 0x01FFD55B
       6     0x00FFFAAA     0.015625    0.0156237 0x00FFFAAA
       7     0x007FFF55    0.0078125   0.00781234 0x007FFF55
       8     0x003FFFEA   0.00390625   0.00390623 0x003FFFEA
       9     0x001FFFFD   0.00195312   0.00195312 0x001FFFFD
      10     0x000FFFFF  0.000976562  0.000976562 0x000FFFFF
      11     0x0007FFFF  0.000488281  0.000488281 0x0007FFFF
      12     0x0003FFFF  0.000244141  0.000244141 0x0003FFFF
      13     0x0001FFFF   0.00012207   0.00012207 0x0001FFFF
      14     0x0000FFFF  6.10352e-05  6.10352e-05 0x0000FFFF
      15     0x00007FFF  3.05176e-05  3.05176e-05 0x00007FFF
      16     0x00003FFF  1.52588e-05  1.52588e-05 0x00003FFF
      17     0x00001FFF  7.62939e-06  7.62939e-06 0x00001FFF
      18     0x00000FFF   3.8147e-06   3.8147e-06 0x00000FFF
      19     0x000007FF  1.90735e-06  1.90735e-06 0x000007FF
      20     0x000003FF  9.53674e-07  9.53674e-07 0x000003FF
      21     0x000001FF  4.76837e-07  4.76837e-07 0x000001FF
      22     0x000000FF  2.38419e-07  2.38419e-07 0x000000FF
      23     0x0000007F  1.19209e-07  1.19209e-07 0x0000007F
      24     0x0000003F  5.96046e-08  5.96046e-08 0x0000003F
      25     0x0000001F  2.98023e-08  2.98023e-08 0x0000001F
      26     0x0000000F  1.49012e-08  1.49012e-08 0x0000000F
      27     0x00000008  7.45058e-09  7.45058e-09 0x00000008
      28     0x00000004  3.72529e-09  3.72529e-09 0x00000004
      29     0x00000002  1.86265e-09  1.86265e-09 0x00000002
      30     0x00000001  9.31323e-10  9.31323e-10 0x00000001
      31     0x00000000  4.65661e-10  4.65661e-10 0x00000000
% rustc gentable.rs && ./gentable 
gentable.rs:5:9: 5:17 warning: unused variable: `num_elms`, #[warn(unused_variables)] on by default
gentable.rs:5     let num_elms: uint = num_bits;
                      ^~~~~~~~
Cordic sin in rust
num bits 32
mul 1073741844
1 << (bits - 2): 1073741844
pi is 3.141593
k1 is 0.607253
computing c = (1f32/shift.powi(ii as i32)).atan() * mul
       i           c           p            a        c2

       0     0x3243f6c0            1     0.785398 0x3243f6c0
       1     0x1dac6700          0.5     0.463648 0x1dac6700
       2     0xfadbb00          0.25     0.244979 0xfadbb00 
       3     0x7f56ea8         0.125     0.124355 0x7f56ea8 
       4     0x3feab78        0.0625     0.062419 0x3feab78 
       5     0x1ffd55c       0.03125      0.03124 0x1ffd55c 
       6     0xfffaab       0.015625     0.015624 0xfffaab  
       7     0x7fff55.8     0.007813     0.007812 0x7fff55.8
       8     0x3fffea.c     0.003906     0.003906 0x3fffea.c
       9     0x1ffffd.6     0.001953     0.001953 0x1ffffd.6
      10     0xfffff.b      0.000977     0.000977 0xfffff.b 
      11     0x7ffff.f8     0.000488     0.000488 0x7ffff.f8
      12     0x40000        0.000244     0.000244 0x40000   
      13     0x20000        0.000122     0.000122 0x20000   
      14     0x10000        0.000061     0.000061 0x10000   
      15     0x8000         0.000031     0.000031 0x8000    
      16     0x4000         0.000015     0.000015 0x4000    
      17     0x2000         0.000008     0.000008 0x2000    
      18     0x1000         0.000004     0.000004 0x1000    
      19     0x800          0.000002     0.000002 0x800     
      20     0x400          0.000001     0.000001 0x400     
      21     0x200                 0            0 0x200     
      22     0x100                 0            0 0x100     
      23     0x80                  0            0 0x80      
      24     0x40                  0            0 0x40      
      25     0x20                  0            0 0x20      
      26     0x10                  0            0 0x10      
      27     0x8                   0            0 0x8       
      28     0x4                   0            0 0x4       
      29     0x2                   0            0 0x2       
      30     0x1                   0            0 0x1       
      31     0x0.8                 0            0 0x0.8     
%