Stare semestry, niepoukladane

8 files changed, 486 insertions, 0 deletions

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 b/Semestr 3/anm/cordic/cordic
new file mode 100644
index 0000000..cb593c0
--- /dev/null
+++ b/Semestr 3/anm/cordic/cordic
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 <stdio.h>

+#include <math.h>

+//#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 b/Semestr 3/anm/cordic/csin
new file mode 100644
index 0000000..49d3cfb
--- /dev/null
+++ b/Semestr 3/anm/cordic/csin
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 <iostream>

+// 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