From c5fcf7179a83ef65c86c6a4a390029149e518649 Mon Sep 17 00:00:00 2001 From: Franciszek Malinka Date: Tue, 5 Oct 2021 21:49:54 +0200 Subject: Duzy commit ze smieciami --- semestr-3/anm/cordic/CORDICtable.c | 21 +++++ semestr-3/anm/cordic/cordic.cpp | 189 +++++++++++++++++++++++++++++++++++++ semestr-3/anm/cordic/cordic.jl | 81 ++++++++++++++++ semestr-3/anm/cordic/csin.cpp | 49 ++++++++++ semestr-3/anm/cordic/taylor.jl | 103 ++++++++++++++++++++ semestr-3/anm/cordic/trig.jl | 43 +++++++++ 6 files changed, 486 insertions(+) create mode 100644 semestr-3/anm/cordic/CORDICtable.c create mode 100644 semestr-3/anm/cordic/cordic.cpp create mode 100644 semestr-3/anm/cordic/cordic.jl create mode 100644 semestr-3/anm/cordic/csin.cpp create mode 100644 semestr-3/anm/cordic/taylor.jl create mode 100644 semestr-3/anm/cordic/trig.jl (limited to 'semestr-3/anm/cordic') diff --git a/semestr-3/anm/cordic/CORDICtable.c b/semestr-3/anm/cordic/CORDICtable.c new file mode 100644 index 0000000..f222d08 --- /dev/null +++ b/semestr-3/anm/cordic/CORDICtable.c @@ -0,0 +1,21 @@ +//CORDIC, 21 bits, 18 iterations +// 1.0 = 262144.000000 multiplication factor +// A = 1.743279 convergence angle (limit is 1.7432866 = 99.9deg) +// F = 1.646760 gain (limit is 1.64676025812107) +// 1/F = 0.607253 inverse gain (limit is 0.607252935008881) +// pi = 3.141593 (3.1415926536897932384626) + +#define CORDIC_A 1.743279 // CORDIC convergence angle A +#define CORDIC_F 0x00069648 // CORDIC gain F +#define CORDIC_1F 0x00026DD4 // CORDIC inverse gain 1/F +#define CORDIC_HALFPI 0x0006487F +#define CORDIC_PI 0x000C90FE +#define CORDIC_TWOPI 0x001921FB +#define CORDIC_MUL 262144.000000 // CORDIC multiplication factor M = 2^18 +#define CORDIC_MAXITER 18 + +int CORDIC_ZTBL[] = { + 0x0003243F, 0x0001DAC6, 0x0000FADC, 0x00007F57, 0x00003FEB, 0x00001FFD, 0x00001000, 0x00000800, + 0x00000400, 0x00000200, 0x00000100, 0x00000080, 0x00000040, 0x00000020, 0x00000010, 0x00000008, + 0x00000004, 0x00000002 }; + diff --git a/semestr-3/anm/cordic/cordic.cpp b/semestr-3/anm/cordic/cordic.cpp new file mode 100644 index 0000000..b0dfd39 --- /dev/null +++ b/semestr-3/anm/cordic/cordic.cpp @@ -0,0 +1,189 @@ +#include +#include +//#define M_PI 3.1415926536897932384626 +int main(int argc, char **argv) +{ + FILE *f; + char tname[50], cname[10]; + int n, n2, mp2, niter, bits, t; + double F, A, mul, tmul; // CORDIC gain, convergence angle, multiplication factor + printf("0)circular, 1)linear, 2)hyperbolic? "); + scanf("%d", &t); + switch (t) + { + case 0: + sprintf(cname, "%s", ""); + break; + case 1: + sprintf(cname, "%s", "_LIN"); + break; + case 2: + sprintf(cname, "%s", "_HYPER"); + break; + } + sprintf(tname, "CORDICtable%s.c", cname); + if (NULL == (f = fopen(tname, "wt"))) + { + printf("cannot write to %s\n", tname); + return 0; + } + + printf("number of bits for mantissa (e.g. 30)? "); + scanf("%d", &bits); + printf("0) mul factor is 2^n (easier output scaling), or\n" + "1) 2pi is 2^n (easier implementation)\n ? "); + scanf("%d", &mp2); + printf("suggested multiplication factor "); + if (mp2 == 0) + { + tmul = (double)(1 << (bits - 3)); + printf("2^%d = %f\n", bits - 3, tmul); + } + else + { + tmul = (double)(1 << (bits - 2)) / M_PI; + printf("2^%d/pi = %f\n", bits - 2, tmul); + } + printf("multiplication factor (0 for suggested)? "); + scanf("%lf", &mul); + if (mul < 0.1) + { + mul = tmul; + printf("%f\n", mul); + } + else + mp2 = -1; // custom mul factor + switch (t) + { + case 0: + for (n = 0; n < bits; n++) + if ((int)round(atan(pow(2.0, (double)(-n))) * mul) == 0) + break; + break; + case 1: + for (n = 0; n < bits; n++) + if ((int)round((pow(2.0, (double)(-n))) * mul) == 0) + break; + break; + case 2: + for (n = 1, n2 = 4; n < bits;) + { + if ((int)round(atanh(pow(2.0, (double)(-n))) * mul) == 0) + break; + if (n == n2) + n2 = 3 * n + 1; + else + n++; + } + break; + } + printf("iterations (up to %d)? ", n); + scanf("%d", &niter); + + F = 1.0; + A = 0.0; + switch (t) + { + case 0: + for (n = 0; n < niter; n++) + { + F = F * sqrt(1 + pow(2.0, -2.0 * n)); + A += atan(pow(2.0, (double)(-n))); + } + break; + case 1: + for (n = 0; n < niter; n++) + { + F = F * sqrt(1); + A += (pow(2.0, (double)(-n))); + } + break; + case 2: + for (n = 1, n2 = 4; n < niter;) + { + F = F * sqrt(1 - pow(2.0, -2.0 * n)); + A += atanh(pow(2.0, (double)(-n))); + if (n == n2) + n2 = 3 * n + 1; + else + n++; + } + break; + } + + fprintf(f, "//CORDIC%s, %d bits, %d iterations\n", cname, bits, niter); + fprintf(f, "// 1.0 = %f multiplication factor\n", mul); + switch (t) + { + case 0: + fprintf(f, "// A = %lf convergence angle " + "(limit is 1.7432866 = 99.9deg)\n", + A); + fprintf(f, "// F = %lf gain (limit is 1.64676025812107)\n", F); + fprintf(f, "// 1/F = %lf inverse gain (limit is 0.607252935008881)\n", 1.0 / F); + break; + case 1: + fprintf(f, "// A = %lf convergence angle (limit is 2)\n", A); + fprintf(f, "// F = %lf gain (limit is 1.0)\n", F); + fprintf(f, "// 1/F = %lf inverse gain (limit is 1.0)\n", 1.0 / F); + break; + case 2: + fprintf(f, "// A = %lf convergence angle " + "(limit is 1.1181730 = 64.0deg)\n", + A); + fprintf(f, "// F = %lf gain (limit is 0.82978162013890)\n", F); + fprintf(f, "// 1/F = %lf inverse gain (limit is 1.20513635844646)\n", 1.0 / F); + break; + } + fprintf(f, "// pi = %lf (3.1415926536897932384626)\n", M_PI); + fprintf(f, "\n"); + fprintf(f, "#define CORDIC%s_A %f // CORDIC convergence angle A\n", cname, A); + fprintf(f, "#define CORDIC%s_F 0x%08X // CORDIC gain F\n", + cname, (int)round(mul * F)); + fprintf(f, "#define CORDIC%s_1F 0x%08X // CORDIC inverse gain 1/F\n", + cname, (int)round(mul / F)); + fprintf(f, "#define CORDIC%s_HALFPI 0x%08X\n", cname, (int)round(mul * (M_PI / 2.0))); + fprintf(f, "#define CORDIC%s_PI 0x%08X\n", cname, (int)round(mul * (M_PI))); + fprintf(f, "#define CORDIC%s_TWOPI 0x%08X\n", cname, (int)round(mul * (2.0 * M_PI))); + fprintf(f, "#define CORDIC%s_MUL %f // CORDIC multiplication factor M", cname, mul); + switch (mp2) + { + case 0: + fprintf(f, " = 2^%d\n", bits - 3); + break; + case 1: + fprintf(f, " = 2^%d/pi\n", bits - 2); + break; + default: + fprintf(f, "\n"); + break; + } + fprintf(f, "#define CORDIC%s_MAXITER %d\n\n", cname, niter); + fprintf(f, "int CORDIC%s_ZTBL[] = {", cname); + for (n = 0; n < niter; n++) + { + if ((n % 8) == 0) + fprintf(f, "\n "); + switch (t) + { + case 0: + fprintf(f, "0x%08X", (int)round(atan(pow(2.0, (double)(-n))) * mul)); + break; + case 1: + fprintf(f, "0x%08X", (int)round((pow(2.0, (double)(-n))) * mul)); + break; + case 2: + n = n == 0 ? 1 : n; + fprintf(f, "0x%08X", (int)round(atanh(pow(2.0, (double)(-n))) * mul)); + break; + } + if (n < (niter - 1)) + fprintf(f, ", "); + else + fprintf(f, " "); + } + fprintf(f, "};\n\n"); + fclose(f); + printf("table written to %s\n", tname); + return 0; +} diff --git a/semestr-3/anm/cordic/cordic.jl b/semestr-3/anm/cordic/cordic.jl new file mode 100644 index 0000000..9509fb2 --- /dev/null +++ b/semestr-3/anm/cordic/cordic.jl @@ -0,0 +1,81 @@ +using Printf + +global ITERATIONS = 30 +global CORDIC_MUL_POW = 30 +global CORDIC_MUL = 2.0^CORDIC_MUL_POW +global CORDIC_ATANS +global CORDIC_F +global CORDIC_F_INV + +CORDIC_ATANS = [843314857, 497837829, 263043837, 133525159, 67021687, 33543516, 16775851, 8388437, 4194283, 2097149, 1048576, 524288, 262144, 131072, 65536, 32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2] +CORDIC_F = 1768195363 +CORDIC_F_INV = 652032874 + +function preprocess_atan(iterations) + global CORDIC_MUL + atan2pow = Array{Float64}(undef, iterations) + @printf("CORDIC_ATANS = [") + for i in 1:iterations + atan2pow[i] = round(atan(1.0 / Float64(BigInt(2)^(i - 1))) * CORDIC_MUL) + @printf("%d", atan2pow[i]) + if i < iterations + @printf(", ") + end + end + @printf("]\n") +end + +function preprocess_scaling_factor(iterations) + # @printf("Preprocessing scaling factor, %d iterations\n", iterations) + CORDIC_F = 1.0 + for i in 0:iterations + CORDIC_F *= sqrt(1. + 1. / Float64(BigInt(2)^(2 * i))) + end + # @printf("Scaling factor: %.16f\n", F) + @printf("CORDIC_F = %d\nCORDIC_F_INV = %d\n", round(CORDIC_F * CORDIC_MUL), round(CORDIC_MUL / CORDIC_F)) +end + +function approx_trig(x, iterations) + global CORDIC_ATANS + global CORDIC_F_INV + X = CORDIC_F_INV + Y = 0 + Z = round(x * CORDIC_MUL) + s = 1 + for i in 0:(iterations - 1) + println(X, " ", Y) + + tempX = X + if Z == 0 + break + end + if Z >= 0 + X -= s * (Y >> i) + Y += s * (tempX >> i) + Z -= s * CORDIC_ATANS[i + 1] + else + X += s * (Y >> i) + Y -= s * (tempX >> i) + Z += s * CORDIC_ATANS[i + 1] + end + end + + println(X, " ", Y) + return (Float64(X) / CORDIC_MUL, Float64(Y) / CORDIC_MUL) +end + +# works only for real numbers +function approx_sin(x, y) + return (approx_trig(x, ITERATIONS)[2], 0) +end + +# works only for real numbers +function approx_cos(x, y) + return (approx_trig(x, ITERATIONS)[1], 0) +end + +function main() + println("Preprocessing CORDIC constants.") + preprocess_atan(CORDIC_MUL_POW) + preprocess_scaling_factor(CORDIC_MUL_POW) +end diff --git a/semestr-3/anm/cordic/csin.cpp b/semestr-3/anm/cordic/csin.cpp new file mode 100644 index 0000000..bcc9af4 --- /dev/null +++ b/semestr-3/anm/cordic/csin.cpp @@ -0,0 +1,49 @@ +#include "CORDICtable.c" +#include +// angle is radians multiplied by CORDIC multiplication factor M +// modulus can be set to CORDIC inverse gain 1/F to avoid post-division +void CORDICsincos(int a, int m, int *s, int *c) +{ + int k, tx, x = m, y = 0, z = a, fl = 0; + if (z > +CORDIC_HALFPI) + { + fl = +1; + z = (+CORDIC_PI) - z; + } + else if (z < -CORDIC_HALFPI) + { + fl = +1; + z = (-CORDIC_PI) - z; + } + for (k = 0; k < CORDIC_MAXITER; k++) + { + std::cout << x << " " << y << " " << z << "\n"; + tx = x; + if (z >= 0) + { + x -= (y >> k); + y += (tx >> k); + z -= CORDIC_ZTBL[k]; + } + else + { + x += (y >> k); + y -= (tx >> k); + z += CORDIC_ZTBL[k]; + } + } + if (fl) + x = -x; + *c = x; // m*cos(a) multiplied by gain F and factor M + *s = y; // m*sin(a) multiplied by gain F and factor M +} + +int main() +{ + double x; + std::cin >> x; + int sinus, cosinus; + CORDICsincos(x * CORDIC_MUL, CORDIC_1F, &sinus, &cosinus); + std::cout << sinus << " " << cosinus << "\n"; + std::cout << (double)sinus / (double)CORDIC_MUL << " " << (double)cosinus / (double)CORDIC_MUL << "\n"; +} \ No newline at end of file diff --git a/semestr-3/anm/cordic/taylor.jl b/semestr-3/anm/cordic/taylor.jl new file mode 100644 index 0000000..5d2d26d --- /dev/null +++ b/semestr-3/anm/cordic/taylor.jl @@ -0,0 +1,103 @@ +module taylor +using Base +using Printf + +ITER = 12 +HYPERBOLIC_MAX = 0.2 + +function series(x, parity, change_sign, iterations) + elements = ones(Float64, 2 * iterations) + res = 0.0 + i = 2 + while i <= 2 * iterations + parity - 2 + elements[i + 1] = elements[i] / i + if change_sign && (i % 2 == parity) + elements[i + 1] = -elements[i + 1] + end - + i += 1 + end + i = 2 * iterations + parity - 2 + while i >= parity + res *= x * x + res += elements[i + 1] + i -= 2 + end + if parity == 1 + res *= x + end + return res +end + +function real_sin(r, iterations) + r = r - floor(r / (2 * pi)) * 2 * pi + if r > pi + return -real_sin(r - pi, iterations) + end + if r > pi / 2 + return real_cos(r - pi / 2, iterations) + end + if r > pi / 4 + return real_cos(pi / 2 - r, iterations) + end + + return series(r, 1, true, iterations) +end + +function real_cos(r, iterations) + r = r - floor(r / (2 * pi)) * 2 * pi + if r > pi + return -real_cos(r - pi, iterations) + end + if r > pi / 2 + return -real_sin(r - pi / 2, iterations) + end + if r > pi / 4 + return real_sin(pi / 2 - r, iterations) + end + + return series(r, 0, true, iterations) +end + +function real_sinh(r, iterations) + if abs(r) > HYPERBOLIC_MAX + return 2 * real_sinh(r / 2, iterations) * real_cosh(r / 2, iterations) + end + return series(r, 1, false, iterations) +end + +function real_cosh(r, iterations) + if abs(r) > HYPERBOLIC_MAX + s = real_sinh(r / 2, iterations) + c = real_cosh(r / 2, iterations) + return s * s + c * c + end + return series(r, 0, false, iterations) +end + +function complex_sin(a, b, iterations) + return (real_sin(a, iterations) * real_cosh(b, iterations), + real_cos(a, iterations) * real_sinh(b, iterations)) +end + +function complex_cos(a, b, iterations) + return (real_cos(a, iterations) * real_cosh(b, iterations), + -real_sin(a, iterations) * real_sinh(b, iterations)) +end + +# c = a + bi +function csin(a, b) + return complex_sin(a, b, ITER) +end + +function ccos(a, b) + return complex_cos(a, b, ITER) +end + +function rsinh(r) + return real_sinh(r, ITER) +end + +function rcosh(r) + return real_cosh(r, ITER) +end +end \ No newline at end of file diff --git a/semestr-3/anm/cordic/trig.jl b/semestr-3/anm/cordic/trig.jl new file mode 100644 index 0000000..bf1dc2c --- /dev/null +++ b/semestr-3/anm/cordic/trig.jl @@ -0,0 +1,43 @@ +include("cordic.jl") +include("taylor.jl") +using Printf + +# preprocessing +cordic.main() + +function rel_error(x, y) + return abs((x - y) / y) +end + +function test_real_sin(arg, trig_func) + @printf("Testing relative error for function %s and argument %f.\n", String(Symbol(trig_func)), arg) + res = trig_func(arg, 0) + real_res = sin(arg) + @printf("Result: %.50f\n", res[1]) + @printf("LibRes: %.50F\n", real_res) + @printf("Relative error: %e\n", rel_error(res[1], real_res)) +end + +function test_complex_sin(arg_real, arg_imag, trig_func) + @printf("Testing relative error for function %s and argument %f + %f i.\n", String(Symbol(trig_func)), arg_real, arg_imag) + res = trig_func(arg_real, arg_imag) + real_res = sin(arg_real + arg_imag * im) + @printf("Result: %.50f + %.50f i\n", res[1], res[2]) + @printf("LibRes: %.50f + %.50f i\n", real(real_res), imag(real_res)) + @printf("Relative error: %e, %e\n", rel_error(res[1], real(real_res)), rel_error(res[2], imag(real_res))) +end + +test_real_sin(0.5, taylor.csin) +test_real_sin(0.5, cordic.approx_sin) + + +test_real_sin(0.001, taylor.csin) +test_real_sin(0.001, cordic.approx_sin) + + +test_real_sin(0.1, taylor.csin) +test_real_sin(0.1, cordic.approx_sin) + + +# test_complex_sin(0.5 + 2pi, 0.5, taylor.csin) +# test_complex_sin(100, 0.5, taylor.csin) \ No newline at end of file -- cgit v1.2.3