sin 1° на калькуляторе
Придумал этот алгоритм Джек Волдер, который тогда работал в компании Конвэйр над навигационным вычислителем вышеупомянутого бомбардировщика.
Главное преимущество метода «цифра за цифрой» в том, что он использует только операции сложения и деления на два (которое легко реализовать сдвигом вправо).
Кроме того, алгоритм можно заставить работать прямо в двоично-десятичном коде, который используется в большинстве калькуляторов, но в приведённом ниже примере мы в эти дебри не полезем.
Алгоритм итеративный и использует таблицу арктангенсов, по одному на итерацию. Таблицу нужно посчитать заранее:
#include
#include
int main(int argc, char **argv)
{
int bits = 32;
int cordic_one = 1 << (bits - 2);
printf("// Число с фиксированной точкой, соответствующее единице с плавающей точкой\n");
printf("static const int cordic_one = 0x%08x;\n", cordic_one);
printf("static const int cordic_table[] = {\n");
double k = 1;
for (int i = 0; i < bits; i++) {
printf("0x%08x, // 0x%08x * atan(1/%.0f) \n", (int)(atan(pow(2, -i)) * cordic_one), cordic_one, pow(2, i));
k /= sqrt(1 + pow(2, -2 * i));
}
printf("};\n");
printf("static const int cordic_k = 0x%08x; // %.16f * 0x%08x\n", (int)(k * cordic_one), k, cordic_one);
}
Заодно считается масштабирующий коэффициент cordic_k.
После этого посчитать пресловутый sin 1° можно так:
#include
#include
// Число с фиксированной точкой, соответствующее единице с плавающей точкой
static const int cordic_one = 0x40000000;
static const int cordic_table[] = {
0x3243f6a8, // 0x40000000 * atan(1/1)
0x1dac6705, // 0x40000000 * atan(1/2)
0x0fadbafc, // 0x40000000 * atan(1/4)
0x07f56ea6, // 0x40000000 * atan(1/8)
0x03feab76, // 0x40000000 * atan(1/16)
0x01ffd55b, // 0x40000000 * atan(1/32)
0x00fffaaa, // 0x40000000 * atan(1/64)
0x007fff55, // 0x40000000 * atan(1/128)
0x003fffea, // 0x40000000 * atan(1/256)
0x001ffffd, // 0x40000000 * atan(1/512)
0x000fffff, // 0x40000000 * atan(1/1024)
0x0007ffff, // 0x40000000 * atan(1/2048)
0x0003ffff, // 0x40000000 * atan(1/4096)
0x0001ffff, // 0x40000000 * atan(1/8192)
0x0000ffff, // 0x40000000 * atan(1/16384)
0x00007fff, // 0x40000000 * atan(1/32768)
0x00003fff, // 0x40000000 * atan(1/65536)
0x00001fff, // 0x40000000 * atan(1/131072)
0x00000fff, // 0x40000000 * atan(1/262144)
0x000007ff, // 0x40000000 * atan(1/524288)
0x000003ff, // 0x40000000 * atan(1/1048576)
0x000001ff, // 0x40000000 * atan(1/2097152)
0x000000ff, // 0x40000000 * atan(1/4194304)
0x0000007f, // 0x40000000 * atan(1/8388608)
0x0000003f, // 0x40000000 * atan(1/16777216)
0x0000001f, // 0x40000000 * atan(1/33554432)
0x0000000f, // 0x40000000 * atan(1/67108864)
0x00000008, // 0x40000000 * atan(1/134217728)
0x00000004, // 0x40000000 * atan(1/268435456)
0x00000002, // 0x40000000 * atan(1/536870912)
0x00000001, // 0x40000000 * atan(1/1073741824)
0x00000000, // 0x40000000 * atan(1/2147483648)
};
static const int cordic_k = 0x26dd3b6a; // 0.6072529350088813 * 0x40000000
void cordic(int theta, int& s, int& c)
{
c = cordic_k;
s = 0;
for (int k = 0; k < 32; ++k) {
int d = (theta >= 0) ? 0 : -1;
int tx = c - (((s >> k) ^ d) - d);
int ty = s + (((c >> k) ^ d) - d);
c = tx; s = ty;
theta -= ((cordic_table[k] ^ d) - d);
}
}
int main(void)
{
double alpha = M_PI / 180;
int sine, cosine;
cordic(alpha * cordic_one, sine, cosine);
printf("CORDIC: %.8f\nExpected: %.8f\n", (double)sine / cordic_one, sin(alpha));
}
Результат:
CORDIC: 0.01745240
Expected: 0.01745241
Тут 32 итерации, поэтому осталась небольшая ошибка. Калькуляторы обычно используют 40 итераций.
