Ошибка округления в Rust при попытке создать CORDIC таблицу для функции sin
Я пытаюсь реализовать метод CORDIC в rust используя эту реализацию c в качестве примера, однако у меня возникают проблемы с ошибками округления при создании таблицы. Вот мой код и результаты.
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
Почему я получаю эти (округления?) Ошибки и как их исправить?
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 (или же ii в версии Rust) и вывод 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
%