fork download
copy 
  1. // Cordic.cpp : Defines the entry point for the console application.
  2. //
  3.  
  4.  
  5. #include <stdio.h>
  6. #include <math.h>
  7.  
  8. /* working with IEEE doubles, means there are 53 bits
  9.  * of mantissa
  10.  */
  11. #define MAXBITS         53
  12. /* define this to perform all (non 2power) multiplications
  13.  * and divisions with the cordic linear method.
  14.  */
  15. #define CORDIC_LINEARx
  16. /* these are the constants needed */
  17. static double invGain1;
  18. static double invGain2;
  19. static double atanTable[MAXBITS];
  20. static double atanhTable[MAXBITS];
  21. static double gain1Cordic();
  22. static double gain2Cordic();
  23. void initCordic()
  24. {
  25.     /* must call this first to initialise the constants.
  26.      * of course, here i use the maths library, but the
  27.      * values would be precomputed.
  28.      */
  29.     double t = 1;
  30.     int i;
  31.     for (i = 0; i < MAXBITS; ++i) {
  32.         atanTable[i] = atan(t);
  33.         t /= 2;
  34.         atanhTable[i] = 0.5*log((1+t)/(1-t));
  35.     }
  36.     /* set constants */
  37.     invGain1 = 1/gain1Cordic();
  38.     invGain2 = 1/gain2Cordic();
  39. }
  40. void cordit1(double* x0, double* y0, double* z0, double vecmode)
  41. {
  42.     /* this is the circular method. 
  43.      * one slight change from the other methods is the y < vecmode 
  44.      * test. this is to implement arcsin, otherwise it can be
  45.      * y < 0 and you can compute arcsin from arctan using
  46.      * trig identities, so its not essential.
  47.      */
  48.     double t;
  49.     double x, y, z;
  50.     int i;
  51.     t = 1;
  52.     x = *x0; y = *y0; z = *z0;
  53.     for (i = 0; i < MAXBITS; ++i) {
  54.         double x1;
  55.         if (vecmode >= 0 && y < vecmode || vecmode<0  && z >= 0) {
  56.             x1 = x - y*t;
  57.             y = y + x*t;
  58.             z = z - atanTable[i];
  59.         }
  60.         else {
  61.             x1 = x + y*t;
  62.             y = y - x*t;
  63.             z = z + atanTable[i];
  64.         }
  65.         x = x1;
  66.         t /= 2;
  67.     }
  68.     *x0 = x;
  69.     *y0 = y;
  70.     *z0 = z;
  71. }
  72. void cordit2(double* x0, double* y0, double* z0, double vecmode)
  73. {
  74.     /* here's the hyperbolic methods. its very similar to
  75.      * the circular methods, but with some small differences.
  76.      *
  77.      * the `x' iteration have the opposite sign.
  78.      * iterations 4, 7, .. 3k+1 are repeated.
  79.      * iteration 0 is not performed.
  80.      */
  81.     double t;
  82.     double x, y, z;
  83.     int i;
  84.     t = 0.5;
  85.     x = *x0; y = *y0; z = *z0;
  86.     int k = 3;
  87.     for (i = 0; i < MAXBITS; ++i) {
  88.         double x1;
  89.         int j;
  90.         for (j = 0; j < 2; ++j) {
  91.             if (vecmode >= 0 && y < 0 || vecmode<0  && z >= 0) {
  92.                 x1 = x + y*t;
  93.                 y = y + x*t;
  94.                 z = z - atanhTable[i];
  95.             }
  96.             else {
  97.                 x1 = x - y*t;
  98.                 y = y - x*t;
  99.                 z = z + atanhTable[i];
  100.             }
  101.             x = x1;
  102.             if (k) {
  103.                 --k;
  104.                 break;
  105.             }
  106.             else k = 3;
  107.         }
  108.         t /= 2;
  109.     }
  110.     *x0 = x;
  111.     *y0 = y;
  112.     *z0 = z;
  113. }
  114. void cordit0(double* x0, double* y0, double* z0, double vecmode)
  115. {
  116.     /* the linear methods is the same as the circular but
  117.      * ive simplified out the x iteration as it doesnt change.
  118.      */
  119.     double t;
  120.     double x, y, z;
  121.     int i;
  122.     t = 1;
  123.     x = *x0; y = *y0; z = *z0;
  124.     for (i = 0; i < MAXBITS; ++i) {
  125.         if (vecmode >= 0 && y < 0 || vecmode<0  && z >= 0) {
  126.             y = y + x*t;
  127.             z = z - t;
  128.         }
  129.         else {
  130.             y = y - x*t;
  131.             z = z + t;
  132.         }
  133.         t /= 2;
  134.     }
  135.     *x0 = x;
  136.     *y0 = y;
  137.     *z0 = z;
  138. }
  139. /** Linear features ***********************************************/
  140. double mulCordic(double a, double b)
  141. {
  142.     double x, y, z;
  143.     x = a;
  144.     y = 0;
  145.     z = b;
  146.     cordit0(&x, &y, &z, -1);
  147.     return y;
  148. }
  149. double divCordic(double a, double b)
  150. {
  151.     double x, y, z;
  152.     x = b;
  153.     y = a;
  154.     z = 0;
  155.     cordit0(&x, &y, &z, 0);
  156.     return z;
  157. }
  158. #ifdef CORDIC_LINEAR
  159. #define MULT(_a, _b)    mulCordic(_a, _b)
  160. #define DIVD(_a, _b)    divCordic(_a, _b)
  161. #else
  162. #define MULT(_a, _b)    (_a)*(_b)
  163. #define DIVD(_a, _b)    (_a)/(_b)
  164. #endif 
  165. /** circular features ***********************************************/
  166. static double gain1Cordic()
  167. {
  168.     /* compute gain by evaluating cos(0) without inv gain */
  169.     double x, y, z;
  170.     x = 1;
  171.     y = 0;
  172.     z = 0;
  173.     cordit1(&x, &y, &z, -1);
  174.     return x;
  175. }
  176. double atanCordic(double a)
  177. {
  178.     /* domain: all a */
  179.     double x = 1;
  180.     double z = 0;
  181.     cordit1(&x, &a, &z, 0);
  182.     return z;
  183. }
  184. double sincosCordic(double a, double* cosp)
  185. {
  186.     /* |a| < 1.74 */
  187.     double sinp = 0;
  188.     *cosp = invGain1;
  189.     cordit1(cosp, &sinp, &a, -1);
  190.     return sinp;
  191. }
  192. double tanCordic(double a)
  193. {
  194.     /* |a| < 1.74 */
  195.     double sinp = 0;
  196.     double cosp = invGain1;
  197.     cordit1(&cosp, &sinp, &a, -1);
  198.     return DIVD(sinp, cosp);
  199. }
  200. double asinCordic(double a)
  201. {
  202.     /* |a| < 0.98 */
  203.     double x, y, z;
  204.  
  205.     x = invGain1;
  206.     y = 0;
  207.     z = 0;
  208.     int neg = 1;
  209.     if (a < 0) {
  210.         a = -a;
  211.         neg = 0;
  212.     }
  213.  
  214.     cordit1(&x, &y, &z, a);
  215.     if (neg) z = -z;
  216.     return z;
  217. }
  218. /** hyperbolic features ********************************************/
  219. double gain2Cordic()
  220. {
  221.     /* calculate hyperbolic gain */
  222.     double x, y, z;
  223.     x = 1;
  224.     y = 0;
  225.     z = 0;
  226.     cordit2(&x, &y, &z, -1);
  227.     return x;
  228. }
  229. double sinhcoshCordic(double a, double* coshp)
  230. {
  231.     /* |a| < 1.13 */
  232.     double y;
  233.     *coshp = invGain2;
  234.     y = 0;
  235.     cordit2(coshp, &y, &a, -1);
  236.     return y;
  237. }
  238. double tanhCordic(double a)
  239. {
  240.     /* |a| < 1.13 */
  241.    double sinhp, coshp;
  242.    coshp = invGain2;
  243.    sinhp = 0;
  244.    cordit2(&coshp, &sinhp, &a, -1);
  245.    return DIVD(sinhp,coshp);
  246. }
  247. double atanhCordic(double a)
  248. {
  249.     /* |a| < 1.13 */
  250.     double x, z;
  251.     x = 1;
  252.     z = 0;
  253.     cordit2(&x, &a, &z, 0);
  254.     return z;
  255. }
  256. double logCordic(double a)
  257. {
  258.     /* 0.1 < a < 9.58 */
  259.     double x, y, z;
  260.     x = a + 1;
  261.     y = a - 1;
  262.     z = 0;
  263.     cordit2(&x, &y, &z, 0);
  264.     return 2*z;
  265. }
  266. double sqrtCordic(double a)
  267. {
  268.     /* 0.03 < a < 2 */
  269.     double x, y, z;
  270.     x = a + 0.25;
  271.     y = a - 0.25;
  272.     z = 0;
  273.     cordit2(&x, &y, &z, 0);
  274.     return MULT(x, invGain2);
  275. }
  276. double expCordic(double a)
  277. {
  278.     double sinhp, coshp;
  279.     coshp = invGain2;
  280.     sinhp = 0;
  281.     cordit2(&coshp, &sinhp, &a, -1);
  282.     return sinhp + coshp;
  283. }
  284.  
  285. double asinCordic1 (double a)
  286. {
  287. 	return 2.0 * atanCordic(a) * a / (1.0 + sqrt(1.0 - a*a));
  288. }
  289.  
  290. int main()
  291. {
  292.     /* just run a few tests */
  293.     double v;
  294.     double x;
  295.     double c;
  296.     initCordic();
  297.     x = 1;
  298.     v = atanCordic(x);
  299.     printf("atan %f = %.18e\n", x, v);
  300.     x = 1;
  301.     v = sincosCordic(x, &c);
  302.     printf("sin %f = %.18e\n", x, v);
  303.     printf("cos %f = %.18e\n", x, c);
  304.     x = 1;
  305.     v = tanCordic(x);
  306.     printf("tan %f = %.18e\n", x, v);
  307.     x = 0.5;
  308.     v = asinCordic(x);
  309.     printf("asin %f = %.18e\n", x, v);
  310. 	x = 0.999999;
  311.     v = asinCordic(x);
  312.     printf("asin %f = %.18e\n", x, v);
  313. 	v = asinCordic1(x);
  314.     printf("asin %f = %.18e\n", x, v);
  315. 	x = 0.61;
  316.     v = asinCordic(x);
  317.     printf("asin %f = %.18e\n", x, v);
  318. 	v = asinCordic1(x);
  319.     printf("asin %f = %.18e\n", x, v);
  320. 	x = 0.65;
  321.     v = asinCordic(x);
  322.     printf("asin %f = %.18e\n", x, v);
  323. 	v = asinCordic1(x);
  324.     printf("asin %f = %.18e\n", x, v);
  325.     x = 1;
  326.     v = sinhcoshCordic(x, &c);
  327.     printf("sinh %f = %.18e\n", x, v);
  328.     printf("cosh %f = %.18e\n", x, c);
  329.     x = 1;
  330.     v = tanhCordic(x);
  331.     printf("tanhh %f = %.18e\n", x, v);
  332.     x = 0.5;
  333.     v = atanhCordic(x);
  334.     printf("atanh %f = %.18e\n", x, v);
  335.     x = 0.8;
  336.     v = logCordic(x);
  337.     printf("log %f = %.18e\n", x, v);
  338.     x = 2;
  339.     v = sqrtCordic(x);
  340.     printf("sqrt %f = %.18e\n", x, v);
  341.     x = 1;
  342.     v = expCordic(x);
  343.     printf("exp %f = %.18e\n", x, v);
  344.  
  345.     return 0;
  346. }
  347.  
Success #stdin #stdout 0s 3300KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
atan 1.000000 = 7.853981633974483900e-01
sin 1.000000 = 8.414709848078963939e-01
cos 1.000000 = 5.403023058681397650e-01
tan 1.000000 = 1.557407724654901848e+00
asin 0.500000 = 5.235987755982988157e-01
asin 0.999999 = 1.743286620472339621e+00
asin 0.999999 = 1.568575455870315771e+00
asin 0.610000 = 7.548049243241690132e-01
asin 0.610000 = 3.728198447779662583e-01
asin 0.650000 = 7.548049243241690132e-01
asin 0.650000 = 4.257476123251776046e-01
sinh 1.000000 = 1.175201193643801378e+00
cosh 1.000000 = 1.543080634815243712e+00
tanhh 1.000000 = 7.615941559557648510e-01
atanh 0.500000 = 5.493061443340548911e-01
log 0.800000 = -2.231435513142098759e-01
sqrt 2.000000 = 1.414213562373094923e+00
exp 1.000000 = 2.718281828459045091e+00
https://ideone.com/XkUmf9
language:
C++14 (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!