CORDIC trig functions

This commit is contained in:
bumbread 2022-06-18 13:53:07 +11:00
parent a16e2d5a6e
commit e79c192472
8 changed files with 866 additions and 30 deletions

View File

@ -0,0 +1,202 @@
double tK = 0x1.b7b2b62cef828p-1;
double hK = 0x1.a2166ada47ca4p-1;
double ttable[] = {
0x1.dac670561bb4fp-2,
0x1.f5b75f92c80ddp-3,
0x1.fd5ba9aac2f6ep-4,
0x1.ff55bb72cfdeap-5,
0x1.ffd55bba97625p-6,
0x1.fff555bbb729bp-7,
0x1.fffd555bbba97p-8,
0x1.ffff5555bbbb7p-9,
0x1.ffffd5555bbbcp-10,
0x1.fffff55555bbcp-11,
0x1.fffffd55555bcp-12,
0x1.ffffff555555cp-13,
0x1.ffffffd555556p-14,
0x1.fffffff555555p-15,
0x1.fffffffd55555p-16,
0x1.ffffffff55555p-17,
0x1.ffffffffd5555p-18,
0x1.fffffffff5555p-19,
0x1.fffffffffd555p-20,
0x1.ffffffffff555p-21,
0x1.ffffffffffd55p-22,
0x1.fffffffffff55p-23,
0x1.fffffffffffd5p-24,
0x1.ffffffffffff5p-25,
0x1.ffffffffffffdp-26,
0x1.fffffffffffffp-27,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
0x1.0000000000000p-32,
0x1.0000000000000p-33,
0x1.0000000000000p-34,
0x1.0000000000000p-35,
0x1.0000000000000p-36,
0x1.0000000000000p-37,
0x1.0000000000000p-38,
0x1.0000000000000p-39,
0x1.0000000000000p-40,
0x1.0000000000000p-41,
0x1.0000000000000p-42,
0x1.0000000000000p-43,
0x1.0000000000000p-44,
0x1.0000000000000p-45,
0x1.0000000000000p-46,
0x1.0000000000000p-47,
0x1.0000000000000p-48,
0x1.0000000000000p-49,
0x1.0000000000000p-50,
0x1.0000000000000p-51,
0x1.0000000000000p-52,
0x1.0000000000000p-53,
0x1.0000000000000p-54,
0x1.0000000000000p-55,
0x1.0000000000000p-56,
0x1.0000000000000p-57,
0x1.0000000000000p-58,
0x1.0000000000000p-59,
0x1.0000000000000p-60,
0x1.0000000000000p-61,
0x1.0000000000000p-62,
0x1.0000000000000p-63,
};
double htable[] = {
0x1.193ea7aad030bp-1,
0x1.058aefa811452p-2,
0x1.015891c9eaef8p-3,
0x1.005588ad375adp-4,
0x1.001558891aee2p-5,
0x1.000555888ad1dp-6,
0x1.000155588891bp-7,
0x1.000055558888bp-8,
0x1.0000155558889p-9,
0x1.0000055555889p-10,
0x1.0000015555589p-11,
0x1.0000005555559p-12,
0x1.0000001555556p-13,
0x1.0000000555555p-14,
0x1.0000000155555p-15,
0x1.0000000055555p-16,
0x1.0000000015555p-17,
0x1.0000000005555p-18,
0x1.0000000001555p-19,
0x1.0000000000555p-20,
0x1.0000000000155p-21,
0x1.0000000000055p-22,
0x1.0000000000015p-23,
0x1.0000000000005p-24,
0x1.0000000000001p-25,
0x1.0000000000000p-26,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
0x1.0000000000000p-32,
0x1.0000000000000p-33,
0x1.0000000000000p-34,
0x1.0000000000000p-35,
0x1.0000000000000p-36,
0x1.0000000000000p-37,
0x1.0000000000000p-38,
0x1.0000000000000p-39,
0x1.0000000000000p-40,
0x1.0000000000000p-41,
0x1.0000000000000p-42,
0x1.0000000000000p-43,
0x1.0000000000000p-44,
0x1.0000000000000p-45,
0x1.0000000000000p-46,
0x1.0000000000000p-47,
0x1.0000000000000p-48,
0x1.0000000000000p-49,
0x1.0000000000000p-50,
0x1.0000000000000p-51,
0x1.0000000000000p-52,
0x1.0000000000000p-53,
0x1.0000000000000p-54,
0x1.0000000000000p-55,
0x1.0000000000000p-56,
0x1.0000000000000p-57,
0x1.0000000000000p-58,
0x1.0000000000000p-59,
0x1.0000000000000p-60,
0x1.0000000000000p-61,
0x1.0000000000000p-62,
0x1.0000000000000p-63,
};
double ptable[] = {
0x1.0000000000000p-1,
0x1.0000000000000p-2,
0x1.0000000000000p-3,
0x1.0000000000000p-4,
0x1.0000000000000p-5,
0x1.0000000000000p-6,
0x1.0000000000000p-7,
0x1.0000000000000p-8,
0x1.0000000000000p-9,
0x1.0000000000000p-10,
0x1.0000000000000p-11,
0x1.0000000000000p-12,
0x1.0000000000000p-13,
0x1.0000000000000p-14,
0x1.0000000000000p-15,
0x1.0000000000000p-16,
0x1.0000000000000p-17,
0x1.0000000000000p-18,
0x1.0000000000000p-19,
0x1.0000000000000p-20,
0x1.0000000000000p-21,
0x1.0000000000000p-22,
0x1.0000000000000p-23,
0x1.0000000000000p-24,
0x1.0000000000000p-25,
0x1.0000000000000p-26,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
0x1.0000000000000p-32,
0x1.0000000000000p-33,
0x1.0000000000000p-34,
0x1.0000000000000p-35,
0x1.0000000000000p-36,
0x1.0000000000000p-37,
0x1.0000000000000p-38,
0x1.0000000000000p-39,
0x1.0000000000000p-40,
0x1.0000000000000p-41,
0x1.0000000000000p-42,
0x1.0000000000000p-43,
0x1.0000000000000p-44,
0x1.0000000000000p-45,
0x1.0000000000000p-46,
0x1.0000000000000p-47,
0x1.0000000000000p-48,
0x1.0000000000000p-49,
0x1.0000000000000p-50,
0x1.0000000000000p-51,
0x1.0000000000000p-52,
0x1.0000000000000p-53,
0x1.0000000000000p-54,
0x1.0000000000000p-55,
0x1.0000000000000p-56,
0x1.0000000000000p-57,
0x1.0000000000000p-58,
0x1.0000000000000p-59,
0x1.0000000000000p-60,
0x1.0000000000000p-61,
0x1.0000000000000p-62,
0x1.0000000000000p-63,
};

View File

@ -0,0 +1,106 @@
float tKf = 0x1.b7b2b62cef828p-1;
float hKf = 0x1.a2166ada47ca4p-1;
float ttablef[] = {
0x1.dac670561bb4fp-2,
0x1.f5b75f92c80ddp-3,
0x1.fd5ba9aac2f6ep-4,
0x1.ff55bb72cfdeap-5,
0x1.ffd55bba97625p-6,
0x1.fff555bbb729bp-7,
0x1.fffd555bbba97p-8,
0x1.ffff5555bbbb7p-9,
0x1.ffffd5555bbbcp-10,
0x1.fffff55555bbcp-11,
0x1.fffffd55555bcp-12,
0x1.ffffff555555cp-13,
0x1.ffffffd555556p-14,
0x1.fffffff555555p-15,
0x1.fffffffd55555p-16,
0x1.ffffffff55555p-17,
0x1.ffffffffd5555p-18,
0x1.fffffffff5555p-19,
0x1.fffffffffd555p-20,
0x1.ffffffffff555p-21,
0x1.ffffffffffd55p-22,
0x1.fffffffffff55p-23,
0x1.fffffffffffd5p-24,
0x1.ffffffffffff5p-25,
0x1.ffffffffffffdp-26,
0x1.fffffffffffffp-27,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
};
float htablef[] = {
0x1.193ea7aad030bp-1,
0x1.058aefa811452p-2,
0x1.015891c9eaef8p-3,
0x1.005588ad375adp-4,
0x1.001558891aee2p-5,
0x1.000555888ad1dp-6,
0x1.000155588891bp-7,
0x1.000055558888bp-8,
0x1.0000155558889p-9,
0x1.0000055555889p-10,
0x1.0000015555589p-11,
0x1.0000005555559p-12,
0x1.0000001555556p-13,
0x1.0000000555555p-14,
0x1.0000000155555p-15,
0x1.0000000055555p-16,
0x1.0000000015555p-17,
0x1.0000000005555p-18,
0x1.0000000001555p-19,
0x1.0000000000555p-20,
0x1.0000000000155p-21,
0x1.0000000000055p-22,
0x1.0000000000015p-23,
0x1.0000000000005p-24,
0x1.0000000000001p-25,
0x1.0000000000000p-26,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
};
float ptablef[] = {
0x1.0000000000000p-1,
0x1.0000000000000p-2,
0x1.0000000000000p-3,
0x1.0000000000000p-4,
0x1.0000000000000p-5,
0x1.0000000000000p-6,
0x1.0000000000000p-7,
0x1.0000000000000p-8,
0x1.0000000000000p-9,
0x1.0000000000000p-10,
0x1.0000000000000p-11,
0x1.0000000000000p-12,
0x1.0000000000000p-13,
0x1.0000000000000p-14,
0x1.0000000000000p-15,
0x1.0000000000000p-16,
0x1.0000000000000p-17,
0x1.0000000000000p-18,
0x1.0000000000000p-19,
0x1.0000000000000p-20,
0x1.0000000000000p-21,
0x1.0000000000000p-22,
0x1.0000000000000p-23,
0x1.0000000000000p-24,
0x1.0000000000000p-25,
0x1.0000000000000p-26,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
};

View File

@ -0,0 +1,202 @@
long double tKl = 0x1.b7b2b62cef828p-1;
long double hKl = 0x1.a2166ada47ca4p-1;
long double ttablel[] = {
0x1.dac670561bb4fp-2,
0x1.f5b75f92c80ddp-3,
0x1.fd5ba9aac2f6ep-4,
0x1.ff55bb72cfdeap-5,
0x1.ffd55bba97625p-6,
0x1.fff555bbb729bp-7,
0x1.fffd555bbba97p-8,
0x1.ffff5555bbbb7p-9,
0x1.ffffd5555bbbcp-10,
0x1.fffff55555bbcp-11,
0x1.fffffd55555bcp-12,
0x1.ffffff555555cp-13,
0x1.ffffffd555556p-14,
0x1.fffffff555555p-15,
0x1.fffffffd55555p-16,
0x1.ffffffff55555p-17,
0x1.ffffffffd5555p-18,
0x1.fffffffff5555p-19,
0x1.fffffffffd555p-20,
0x1.ffffffffff555p-21,
0x1.ffffffffffd55p-22,
0x1.fffffffffff55p-23,
0x1.fffffffffffd5p-24,
0x1.ffffffffffff5p-25,
0x1.ffffffffffffdp-26,
0x1.fffffffffffffp-27,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
0x1.0000000000000p-32,
0x1.0000000000000p-33,
0x1.0000000000000p-34,
0x1.0000000000000p-35,
0x1.0000000000000p-36,
0x1.0000000000000p-37,
0x1.0000000000000p-38,
0x1.0000000000000p-39,
0x1.0000000000000p-40,
0x1.0000000000000p-41,
0x1.0000000000000p-42,
0x1.0000000000000p-43,
0x1.0000000000000p-44,
0x1.0000000000000p-45,
0x1.0000000000000p-46,
0x1.0000000000000p-47,
0x1.0000000000000p-48,
0x1.0000000000000p-49,
0x1.0000000000000p-50,
0x1.0000000000000p-51,
0x1.0000000000000p-52,
0x1.0000000000000p-53,
0x1.0000000000000p-54,
0x1.0000000000000p-55,
0x1.0000000000000p-56,
0x1.0000000000000p-57,
0x1.0000000000000p-58,
0x1.0000000000000p-59,
0x1.0000000000000p-60,
0x1.0000000000000p-61,
0x1.0000000000000p-62,
0x1.0000000000000p-63,
};
long double htablel[] = {
0x1.193ea7aad030bp-1,
0x1.058aefa811452p-2,
0x1.015891c9eaef8p-3,
0x1.005588ad375adp-4,
0x1.001558891aee2p-5,
0x1.000555888ad1dp-6,
0x1.000155588891bp-7,
0x1.000055558888bp-8,
0x1.0000155558889p-9,
0x1.0000055555889p-10,
0x1.0000015555589p-11,
0x1.0000005555559p-12,
0x1.0000001555556p-13,
0x1.0000000555555p-14,
0x1.0000000155555p-15,
0x1.0000000055555p-16,
0x1.0000000015555p-17,
0x1.0000000005555p-18,
0x1.0000000001555p-19,
0x1.0000000000555p-20,
0x1.0000000000155p-21,
0x1.0000000000055p-22,
0x1.0000000000015p-23,
0x1.0000000000005p-24,
0x1.0000000000001p-25,
0x1.0000000000000p-26,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
0x1.0000000000000p-32,
0x1.0000000000000p-33,
0x1.0000000000000p-34,
0x1.0000000000000p-35,
0x1.0000000000000p-36,
0x1.0000000000000p-37,
0x1.0000000000000p-38,
0x1.0000000000000p-39,
0x1.0000000000000p-40,
0x1.0000000000000p-41,
0x1.0000000000000p-42,
0x1.0000000000000p-43,
0x1.0000000000000p-44,
0x1.0000000000000p-45,
0x1.0000000000000p-46,
0x1.0000000000000p-47,
0x1.0000000000000p-48,
0x1.0000000000000p-49,
0x1.0000000000000p-50,
0x1.0000000000000p-51,
0x1.0000000000000p-52,
0x1.0000000000000p-53,
0x1.0000000000000p-54,
0x1.0000000000000p-55,
0x1.0000000000000p-56,
0x1.0000000000000p-57,
0x1.0000000000000p-58,
0x1.0000000000000p-59,
0x1.0000000000000p-60,
0x1.0000000000000p-61,
0x1.0000000000000p-62,
0x1.0000000000000p-63,
};
long double ptablel[] = {
0x1.0000000000000p-1,
0x1.0000000000000p-2,
0x1.0000000000000p-3,
0x1.0000000000000p-4,
0x1.0000000000000p-5,
0x1.0000000000000p-6,
0x1.0000000000000p-7,
0x1.0000000000000p-8,
0x1.0000000000000p-9,
0x1.0000000000000p-10,
0x1.0000000000000p-11,
0x1.0000000000000p-12,
0x1.0000000000000p-13,
0x1.0000000000000p-14,
0x1.0000000000000p-15,
0x1.0000000000000p-16,
0x1.0000000000000p-17,
0x1.0000000000000p-18,
0x1.0000000000000p-19,
0x1.0000000000000p-20,
0x1.0000000000000p-21,
0x1.0000000000000p-22,
0x1.0000000000000p-23,
0x1.0000000000000p-24,
0x1.0000000000000p-25,
0x1.0000000000000p-26,
0x1.0000000000000p-27,
0x1.0000000000000p-28,
0x1.0000000000000p-29,
0x1.0000000000000p-30,
0x1.0000000000000p-31,
0x1.0000000000000p-32,
0x1.0000000000000p-33,
0x1.0000000000000p-34,
0x1.0000000000000p-35,
0x1.0000000000000p-36,
0x1.0000000000000p-37,
0x1.0000000000000p-38,
0x1.0000000000000p-39,
0x1.0000000000000p-40,
0x1.0000000000000p-41,
0x1.0000000000000p-42,
0x1.0000000000000p-43,
0x1.0000000000000p-44,
0x1.0000000000000p-45,
0x1.0000000000000p-46,
0x1.0000000000000p-47,
0x1.0000000000000p-48,
0x1.0000000000000p-49,
0x1.0000000000000p-50,
0x1.0000000000000p-51,
0x1.0000000000000p-52,
0x1.0000000000000p-53,
0x1.0000000000000p-54,
0x1.0000000000000p-55,
0x1.0000000000000p-56,
0x1.0000000000000p-57,
0x1.0000000000000p-58,
0x1.0000000000000p-59,
0x1.0000000000000p-60,
0x1.0000000000000p-61,
0x1.0000000000000p-62,
0x1.0000000000000p-63,
};

View File

@ -0,0 +1,58 @@
import math;
import os;
import sys;
abspath = os.path.abspath(sys.argv[0])
dname = os.path.dirname(abspath)
os.chdir(dname)
data=[
('cordic_dataf.c', 'float', 'f', 32),
('cordic_data.c', 'double', '', 64),
('cordic_datal.c', 'long double', 'l', 64),
]
for f in data:
fname = f[0]
ftype = f[1]
s = f[2]
N = f[3]
f = open(fname, 'w')
tK = 1
for i in range(1, N):
tK *= math.cos(math.atan(2**(-i)));
f.write(ftype + ' tK' + s + ' = ' + float.hex(tK) + ';\n')
f.write('\n')
hK = 1
for i in range(1, N):
hK *= math.cos(math.atanh(2**(-i)));
f.write(ftype + ' hK' + s + ' = ' + float.hex(hK) + ';\n')
f.write('\n')
f.write(ftype + ' ttable' + s + '[] = {\n');
for i in range(1, N):
v = math.atan(2**(-i));
f.write(' ' + float.hex(v) + ',\n')
f.write('};\n')
f.write('\n')
f.write(ftype + ' htable' + s + '[] = {\n');
for i in range(1, N):
v = math.atanh(2**(-i));
f.write(' ' + float.hex(v) + ',\n')
f.write('};\n')
f.write('\n')
f.write(ftype + ' ptable' + s + '[] = {\n');
for i in range(1, N):
v = 2**(-i);
f.write(' ' + float.hex(v) + ',\n')
f.write('};\n')
f.write('\n')
f.close();

View File

@ -5,21 +5,29 @@
#include <stdint.h> #include <stdint.h>
#define ln2 0.69314718055994530941723212145817 #define ln2 0.69314718055994530941723212145817
static double HALF_PI = 1.570796326794896619231321691639751;
static double PI = 3.141592653589793238462643383279502;
static double LOG2E = 1.442695040888963407359924681001892;
#define countof(arr) (sizeof arr/sizeof arr[0])
#define ftype float #define ftype float
#define suffix(name) name ## f #define suffix(name) name ## f
#include "cordic/cordic_dataf.c"
#include "generic.h" #include "generic.h"
#undef ftype #undef ftype
#undef suffix #undef suffix
#define ftype double #define ftype double
#define suffix(name) name #define suffix(name) name
#include "cordic/cordic_data.c"
#include "generic.h" #include "generic.h"
#undef ftype #undef ftype
#undef suffix #undef suffix
#define ftype long double #define ftype long double
#define suffix(name) name ## l #define suffix(name) name ## l
#include "cordic/cordic_datal.c"
#include "generic.h" #include "generic.h"
#undef ftype #undef ftype
#undef suffix #undef suffix

View File

@ -1,38 +1,225 @@
static ftype suffix(cordic)(
ftype x0, // x initial value
ftype y0, // y initial value
ftype z0, // initial 'angle'
int m, // system: 1 for trig, -1 for hyperbolic
int v, // mode: 1 for vectoring, 0 for rotation
ftype *c, // output x
ftype *s // output y
) {
int tab_count = 0;
ftype *tab = NULL;
if(m == 1) {
tab = suffix(ttable);
tab_count = countof(suffix(ttable));
}
else if(m == -1) {
tab = suffix(htable);
tab_count = countof(suffix(htable));
}
ftype x = x0;
ftype y = y0;
ftype z = z0;
for(int i = 0; i != tab_count; ++i) {
ftype sign;
if(v) sign = (y < 0)? 1 : -1;
else sign = (z < 0)? -1 : 1;
ftype x1 = suffix(fma)(suffix(ptable)[i], -m*sign*y, x);
ftype y1 = suffix(fma)(suffix(ptable)[i], x*sign, y);
x = x1;
y = y1;
z = suffix(fma)(tab[i], -sign, z);
}
if(c!=NULL) *c = x;
if(s!=NULL) *s = y;
return z;
}
ftype suffix(sin)(ftype x) {
if(isinf(x)) return NAN;
if(isnan(x)) return NAN;
int k;
ftype alpha = suffix(remquo)(x, HALF_PI, &k);
if(x == 0) return x;
k = (((k % 4) +4) %4);
ftype sinx;
ftype cosx;
suffix(cordic)(suffix(tK), 0, alpha, 1, 0, &cosx, &sinx);
switch(k) {
case 0: return sinx;
case 1: return cosx;
case 2: return -sinx;
case 3: return -cosx;
}
return 0;
}
ftype suffix(cos)(ftype x) {
if(isinf(x)) return NAN;
if(isnan(x)) return NAN;
int k;
ftype alpha = suffix(remquo)(x, HALF_PI, &k);
k = (((k % 4) +4) %4);
ftype sinx;
ftype cosx;
suffix(cordic)(suffix(tK), 0, alpha, 1, 0, &cosx, &sinx);
switch(k) {
case 0: return cosx;
case 1: return -sinx;
case 2: return -cosx;
case 3: return sinx;
}
return 0;
}
ftype suffix(tan)(ftype x) {
if(isinf(x)) return NAN;
if(isnan(x)) return NAN;
int k;
ftype alpha = suffix(remquo)(x, HALF_PI, &k);
if(x == 0) return x;
k = (((k % 2) +2) %2);
ftype sinx;
ftype cosx;
suffix(cordic)(suffix(tK), 0, alpha, 1, 0, &cosx, &sinx);
switch(k) {
case 0: return sinx/cosx;
case 1: return -cosx/sinx;
}
return 0;
}
ftype suffix(cot)(ftype x) {
if(x == 0) return x;
if(isinf(x)) return NAN;
if(isnan(x)) return NAN;
int k;
ftype alpha = suffix(remquo)(x, HALF_PI, &k);
k = (((k % 2) +2) %2);
ftype sinx;
ftype cosx;
suffix(cordic)(suffix(tK), 0, alpha, 1, 0, &cosx, &sinx);
switch(k) {
case 0: return cosx/sinx;
case 1: return -sinx/cosx;
}
return 0;
}
ftype suffix(exp)(ftype x) { ftype suffix(exp)(ftype x) {
if(x == 0) return x;
if(isinf(x)) return NAN;
if(isnan(x)) return NAN; if(isnan(x)) return NAN;
if(x > 0 && isinf(x)) return +INFINITY; ftype t = x*LOG2E;
if(x < 0 && isinf(x)) return +0.0; int64_t k = (int64_t)t;
if(x == 0) return 1.0; x = (t - (ftype)k) / LOG2E;
if(x > 709.8) { ftype xx = x*x;
#if math_errhandling & MATH_ERREXCEPT ftype xxx = xx*x;
feraiseexcept(FE_OVERFLOW); ftype xxxx = xx*xx;
#endif ftype xxxxx = xxx*xx;
#if math_errhandling & MATH_ERRNO ftype expx = 1+x + xx/2.0 + xxx/6.0+ xxxx/24.0 + xxxxx/720.0;
errno = ERANGE; if(k>0) expx *= 2.0;
#endif if(k>0) while(k-- > 0) expx *= 2.0;
return +INFINITY; if(k<0) while(k++ < 0) expx /= 2.0;
} return expx;
if(x < -708.4) {
#if math_errhandling & MATH_ERREXCEPT
feraiseexcept(FE_OVERFLOW);
#endif
#if math_errhandling & MATH_ERRNO
errno = ERANGE;
#endif
return 0;
}
ftype e = 1.0;
ftype xp = 1.0;
ftype f = 1;
for(uint64_t i = 1; i != 10; ++i) {
f *= i;
xp *= x;
e += xp / f;
}
return e;
} }
ftype suffix(atan)(ftype x) {
if(x == 0) return x;
if(isinf(x)) return INFINITY;
if(isnan(x)) return NAN;
ftype atan;
if(x > 1) {
atan = HALF_PI - suffix(cordic)(x, 1, 0, 1, 1, NULL, NULL);
}
else if(x < -1) {
atan = -HALF_PI + suffix(cordic)(-x, 1, 0, 1, 1, NULL, NULL);
}
else {
atan = suffix(cordic)(1, x, 0, 1, 1, NULL, NULL);
}
return atan;
}
ftype suffix(acos)(ftype x) {
if(x == 0) return HALF_PI;
if(x == -1) return PI;
if(x == 1) return 0;
if(isinf(x)) return INFINITY;
if(isnan(x)) return NAN;
if(fabs(x) > 1) return NAN;
ftype atan;
x = sqrt(fma(-x,x,1))/x;
if(x > 1) {
atan = HALF_PI - suffix(cordic)(x, 1, 0, 1, 1, NULL, NULL);
}
else if(x < -1) {
atan = -HALF_PI + suffix(cordic)(-x, 1, 0, 1, 1, NULL, NULL);
}
else {
atan = suffix(cordic)(1, x, 0, 1, 1, NULL, NULL);
}
if(x < 0) atan += PI;
return atan;
}
ftype suffix(asin)(ftype x) {
if(x == 0) return 0;
if(x == -1) return -HALF_PI;
if(x == 1) return HALF_PI;
if(isinf(x)) return INFINITY;
if(isnan(x)) return NAN;
if(fabs(x) > 1) return NAN;
ftype atan;
x /= sqrt(fma(-x,x,1));
if(x > 1) {
atan = HALF_PI - suffix(cordic)(x, 1, 0, 1, 1, NULL, NULL);
}
else if(x < -1) {
atan = -HALF_PI + suffix(cordic)(-x, 1, 0, 1, 1, NULL, NULL);
}
else {
atan = suffix(cordic)(1, x, 0, 1, 1, NULL, NULL);
}
return atan;
}
// ftype suffix(exp)(ftype x) {
// if(isnan(x)) return NAN;
// if(x > 0 && isinf(x)) return +INFINITY;
// if(x < 0 && isinf(x)) return +0.0;
// if(x == 0) return 1.0;
// if(x > 709.8) {
// #if math_errhandling & MATH_ERREXCEPT
// feraiseexcept(FE_OVERFLOW);
// #endif
// #if math_errhandling & MATH_ERRNO
// errno = ERANGE;
// #endif
// return +INFINITY;
// }
// if(x < -708.4) {
// #if math_errhandling & MATH_ERREXCEPT
// feraiseexcept(FE_OVERFLOW);
// #endif
// #if math_errhandling & MATH_ERRNO
// errno = ERANGE;
// #endif
// return 0;
// }
// ftype e = 1.0;
// ftype xp = 1.0;
// ftype f = 1;
// for(uint64_t i = 1; i != 10; ++i) {
// f *= i;
// xp *= x;
// e += xp / f;
// }
// return e;
// }
ftype suffix(exp2)(ftype x) { ftype suffix(exp2)(ftype x) {
return suffix(exp)(x * ln2); return suffix(exp)(x * ln2);
} }

47
test/plot.py Normal file
View File

@ -0,0 +1,47 @@
import matplotlib.pyplot as plt
import numpy as np
import subprocess;
import math;
import csv;
plt.style.use('_mpl-gallery')
pi=math.pi
with open('bin/data.csv', newline='') as csvfile:
reader = csv.reader(csvfile, delimiter=' ', quotechar='|')
rows = [row for row in reader];
# make data
x = [float.fromhex(xi[:-1]) for xi in rows[0][:-1]]
y = [float.fromhex(yi[:-1]) for yi in rows[1][:-1]]
yex = [math.asin(x0) for x0 in x]
err = np.subtract(y, yex);
maxerrv = max(err)
minerrv = min(err)
maxerr = [maxerrv for xi in rows[0][:-1]]
minerr = [minerrv for xi in rows[0][:-1]]
errscale = 10**-14
lo = x[0]
hi = x[-1]
fig, a = plt.subplots(2, 1, constrained_layout=True)
# a.step(x, y, linewidth=1)
a[0].set_xlim(lo, hi);
# a[0].set_ylim(-1.000001, 1.000001);
a[0].set_xlabel('x');
a[0].set_ylabel('f(x)');
a[0].plot(x, y, x, yex)
a[1].set_xlim(lo, hi);
a[1].ticklabel_format(useOffset=False)
a[1].set_ylim(-errscale, errscale);
a[1].set_xlabel('x');
a[1].set_ylabel('absolute error');
a[1].plot(x, err, x, maxerr, x, minerr);
plt.show()

26
test/test_trig.c Normal file
View File

@ -0,0 +1,26 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <math.h>
// TODO: won't work until %a is implemented
int main(int argc, char **argv) {
double lo = -1;
double hi = 1;
double step = 0.0001;
FILE *out = fopen("bin/data.csv", "w");
for(double i = lo; i < hi; i += step) {
fprintf(out, "%f", i);
fprintf(out, ", ");
}
fprintf(out, "\n");
for(double i = lo; i < hi; i += step) {
fprintf(out, "%a", asin(i));
fprintf(out, ", ");
}
fclose(out);
return 0;
}