From e79c192472e8968c505e91b3a9b2cfb2c7488afb Mon Sep 17 00:00:00 2001 From: bumbread Date: Sat, 18 Jun 2022 13:53:07 +1100 Subject: [PATCH] CORDIC trig functions --- src/code/math/cordic/cordic_data.c | 202 +++++++++++++++++++++++ src/code/math/cordic/cordic_dataf.c | 106 ++++++++++++ src/code/math/cordic/cordic_datal.c | 202 +++++++++++++++++++++++ src/code/math/cordic/maketab.py | 58 +++++++ src/code/math/generic.c | 8 + src/code/math/generic.h | 247 ++++++++++++++++++++++++---- test/plot.py | 47 ++++++ test/test_trig.c | 26 +++ 8 files changed, 866 insertions(+), 30 deletions(-) create mode 100644 src/code/math/cordic/cordic_data.c create mode 100644 src/code/math/cordic/cordic_dataf.c create mode 100644 src/code/math/cordic/cordic_datal.c create mode 100644 src/code/math/cordic/maketab.py create mode 100644 test/plot.py create mode 100644 test/test_trig.c diff --git a/src/code/math/cordic/cordic_data.c b/src/code/math/cordic/cordic_data.c new file mode 100644 index 0000000..a208623 --- /dev/null +++ b/src/code/math/cordic/cordic_data.c @@ -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, +}; + diff --git a/src/code/math/cordic/cordic_dataf.c b/src/code/math/cordic/cordic_dataf.c new file mode 100644 index 0000000..ccfa432 --- /dev/null +++ b/src/code/math/cordic/cordic_dataf.c @@ -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, +}; + diff --git a/src/code/math/cordic/cordic_datal.c b/src/code/math/cordic/cordic_datal.c new file mode 100644 index 0000000..7e23664 --- /dev/null +++ b/src/code/math/cordic/cordic_datal.c @@ -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, +}; + diff --git a/src/code/math/cordic/maketab.py b/src/code/math/cordic/maketab.py new file mode 100644 index 0000000..d46700e --- /dev/null +++ b/src/code/math/cordic/maketab.py @@ -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(); diff --git a/src/code/math/generic.c b/src/code/math/generic.c index b61b39b..6bf666e 100644 --- a/src/code/math/generic.c +++ b/src/code/math/generic.c @@ -5,21 +5,29 @@ #include #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 suffix(name) name ## f +#include "cordic/cordic_dataf.c" #include "generic.h" #undef ftype #undef suffix #define ftype double #define suffix(name) name +#include "cordic/cordic_data.c" #include "generic.h" #undef ftype #undef suffix #define ftype long double #define suffix(name) name ## l +#include "cordic/cordic_datal.c" #include "generic.h" #undef ftype #undef suffix diff --git a/src/code/math/generic.h b/src/code/math/generic.h index 872f5c0..e94c78f 100644 --- a/src/code/math/generic.h +++ b/src/code/math/generic.h @@ -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) { + if(x == 0) return x; + if(isinf(x)) return NAN; 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 t = x*LOG2E; + int64_t k = (int64_t)t; + x = (t - (ftype)k) / LOG2E; + ftype xx = x*x; + ftype xxx = xx*x; + ftype xxxx = xx*xx; + ftype xxxxx = xxx*xx; + ftype expx = 1+x + xx/2.0 + xxx/6.0+ xxxx/24.0 + xxxxx/720.0; + if(k>0) expx *= 2.0; + if(k>0) while(k-- > 0) expx *= 2.0; + if(k<0) while(k++ < 0) expx /= 2.0; + return expx; } +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) { return suffix(exp)(x * ln2); } diff --git a/test/plot.py b/test/plot.py new file mode 100644 index 0000000..2623566 --- /dev/null +++ b/test/plot.py @@ -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() + + diff --git a/test/test_trig.c b/test/test_trig.c new file mode 100644 index 0000000..b57fb9b --- /dev/null +++ b/test/test_trig.c @@ -0,0 +1,26 @@ + +#include +#include +#include +#include +#include +#include + +// 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; +}