1
2 // demo-is-self-contained
3
4 #include "aux0/print-fixed.h" // print_fixed()
5 #include "nextarg.h"
6
7 #include "fxttypes.h" // ulong
8 #include <cmath>
9
10 //% CORDIC computation of hyperbolic sine and cosine.
11
12 double *cordic_htab;
13 //const double cordic_Kp = 0.8281593609602156270761983277591751468694538376908425291;
14 const double cordic_1Kp = 1.207497067763072128877721011310915836812783221769813422;
15
16
17 void
make_cordic_htab(ulong na)18 make_cordic_htab(ulong na)
19 {
20 double s = 1.0;
21 cordic_htab[0] = 0;
22 for (ulong k=1; k<na; ++k)
23 {
24 s *= 0.5;
25 cordic_htab[k] = atanh(s); // cordic_htab[0] unused
26 }
27 }
28 // -------------------------
29
30 //#define TEX
31 #ifdef TEX
32 #define COLSEP " & "
33 #define ENDL " \\\\"
34 #else // TEX
35 #define COLSEP " "
36 #define ENDL ""
37 #endif // TEX
38 const long nd = 9; // Number of Digits displayed
39
40 void
print_xyza(long n,double x,double y,double z,double a)41 print_xyza(long n, double x, double y, double z, double a)
42 {
43 cout << setw(2) << n;
44 print_fixed(COLSEP, x, nd, 0);
45 print_fixed(COLSEP, y, nd, 0);
46 print_fixed(COLSEP, z, nd, 1);
47 print_fixed(COLSEP, a, nd, 1);
48 cout << ENDL << endl;
49 }
50 // -------------------------
51
52
53 void
cordic_hyp(double theta,double & s,double & c,ulong n)54 cordic_hyp(double theta, double &s, double &c, ulong n)
55 {
56 double x = cordic_1Kp;
57 double y = 0;
58 double z = theta;
59 double v = 1.0;
60 print_xyza(-1, x, y, z, 0.0);
61 ulong i = 4;
62 for (ulong k=1; k<n; ++k)
63 {
64 v *= 0.5;
65 again:
66 double d = ( z>=0 ? +1 : -1 );
67 double tx = x + d * v * y;
68 double ty = y + d * v * x;
69 double tz = z - d * cordic_htab[k];
70 x = tx; y = ty; z = tz;
71 print_xyza((long)k, x, y, z, -d*cordic_htab[k]);
72 if ( k==i ) { i=3*i+1; goto again; }
73 }
74 c = x;
75 s = y;
76 }
77 // -------------------------
78
79 int
main(int argc,char ** argv)80 main(int argc, char **argv)
81 {
82 double t = 1.0;
83 NXARGFL(t, "Argument for CORDIC sinh/cosh computation");
84
85 ulong n = 16;
86 NXARG(n, "Number of iterations");
87
88 ulong na = 64;
89 // cordic_ctab = new double[na];
90 // make_cordic_ctab(na);
91 cordic_htab = new double[na];
92 make_cordic_htab(na);
93
94 double s, c;
95
96 cordic_hyp(t, s, c, n);
97
98
99 double zs = sinh(t), zc = cosh(t);
100 print_xyza(99, zc, zs, 0.0, 0.0);
101
102 cout << endl;
103 cout << setprecision(19);
104 cout << "dc = " << setw(12) << c-zc << endl;
105 cout << "ds = " << setw(12) << s-zs << endl;
106 cout << "2**(-n) = " << pow(2.0, -(double)n) << endl;
107 cout << endl;
108
109 delete [] cordic_htab;
110
111 return 0;
112 }
113 // -------------------------
114
115 /// Emacs:
116 /// Local Variables:
117 /// MyRelDir: "demo/arith"
118 /// makefile-dir: "../../"
119 /// make-target: "1demo DSRC=demo/arith/cordic-hyp-demo.cc"
120 /// make-target2: "1demo DSRC=demo/arith/cordic-hyp-demo.cc DEMOFLAGS=-DTIMING"
121 /// End:
122
123