cordic methods describe essentially the same algorithm that with suitably chosen inputs can be used to calculate a whole range of scientific functions including; sin, cos, tan, arctan, arcsin, arccos, sinh, cosh, tanh, arctanh, log, exp, square root and even multiply and divide.
the method dates back to volder [1959], and due to its versatility and compactness, it made possible the microcoding of the hp35 pocket scientific calculator in 1972.
here is some code to illustrate the techniques. ive split the methods into three parts; linear, circular and hyperbolic. in the hp35 microcode these would be unified into one function (for space reasons). because the linear mode can perform multiply and divide, you only need add/subtract and shift to complete the implementation.
you can select in the code whether to do the multiples and divides also by cordic means. other multiplies and divides are all powers of 2 (these dont count). to eliminate these too, would involve ieee hackery.
#include <stdio.h> #include <math.h>
/* working with IEEE doubles, means there are 53 bits * of mantissa */ #define MAXBITS 53
/* define this to perform all (non 2power) multiplications * and divisions with the cordic linear method. */ #define CORDIC_LINEARx
/* these are the constants needed */ static double invGain1; static double invGain2; static double atanTable[MAXBITS]; static double atanhTable[MAXBITS];
static double gain1Cordic(); static double gain2Cordic();
void initCordic()
{
/* must call this first to initialise the constants.
* of course, here i use the maths library, but the
* values would be precomputed.
*/
double t = 1;
int i;
for (i = 0; i < MAXBITS; ++i) {
atanTable[i] = atan(t);
t /= 2;
atanhTable[i] = 0.5*log((1+t)/(1-t));
} /* set constants */
invGain1 = 1/gain1Cordic();
invGain2 = 1/gain2Cordic();
}void cordit1(double* x0, double* y0, double* z0, double vecmode)
{
/* this is the circular method.
* one slight change from the other methods is the y < vecmode
* test. this is to implement arcsin, otherwise it can be
* y < 0 and you can compute arcsin from arctan using
* trig identities, so its not essential.
*/ double t;
double x, y, z;
int i; t = 1;
x = *x0; y = *y0; z = *z0; for (i = 0; i < MAXBITS; ++i) {
double x1; if (vecmode >= 0 && y < vecmode || vecmode<0 && z >= 0) {
x1 = x - y*t;
y = y + x*t;
z = z - atanTable[i];
}
else {
x1 = x + y*t;
y = y - x*t;
z = z + atanTable[i];
} x = x1;
t /= 2;
} *x0 = x;
*y0 = y;
*z0 = z;
}void cordit2(double* x0, double* y0, double* z0, double vecmode)
{
/* here's the hyperbolic methods. its very similar to
* the circular methods, but with some small differences.
*
* the `x' iteration have the opposite sign.
* iterations 4, 7, .. 3k+1 are repeated.
* iteration 0 is not performed.
*/ double t;
double x, y, z;
int i; t = 0.5;
x = *x0; y = *y0; z = *z0; int k = 3;
for (i = 0; i < MAXBITS; ++i) {
double x1;
int j; for (j = 0; j < 2; ++j) {
if (vecmode >= 0 && y < 0 || vecmode<0 && z >= 0) {
x1 = x + y*t;
y = y + x*t;
z = z - atanhTable[i];
}
else {
x1 = x - y*t;
y = y - x*t;
z = z + atanhTable[i];
}
x = x1; if (k) {
--k;
break;
}
else k = 3;
} t /= 2;
} *x0 = x;
*y0 = y;
*z0 = z;
}void cordit0(double* x0, double* y0, double* z0, double vecmode)
{ /* the linear methods is the same as the circular but
* ive simplified out the x iteration as it doesnt change.
*/
double t;
double x, y, z;
int i; t = 1;
x = *x0; y = *y0; z = *z0; for (i = 0; i < MAXBITS; ++i) { if (vecmode >= 0 && y < 0 || vecmode<0 && z >= 0) {
y = y + x*t;
z = z - t;
}
else {
y = y - x*t;
z = z + t;
}
t /= 2;
} *x0 = x;
*y0 = y;
*z0 = z;
}/** Linear features ***********************************************/
double mulCordic(double a, double b)
{
double x, y, z;
x = a;
y = 0;
z = b;
cordit0(&x, &y, &z, -1);
return y;
}double divCordic(double a, double b)
{
double x, y, z;
x = b;
y = a;
z = 0;
cordit0(&x, &y, &z, 0);
return z;
}#ifdef CORDIC_LINEAR #define MULT(_a, _b) mulCordic(_a, _b) #define DIVD(_a, _b) divCordic(_a, _b) #else #define MULT(_a, _b) (_a)*(_b) #define DIVD(_a, _b) (_a)/(_b) #endif
/** circular features ***********************************************/
static double gain1Cordic()
{
/* compute gain by evaluating cos(0) without inv gain */
double x, y, z; x = 1;
y = 0;
z = 0;
cordit1(&x, &y, &z, -1);
return x;
}double atanCordic(double a)
{
/* domain: all a */
double x = 1;
double z = 0;
cordit1(&x, &a, &z, 0);
return z;
}double sincosCordic(double a, double* cosp)
{
/* |a| < 1.74 */
double sinp = 0;
*cosp = invGain1;
cordit1(cosp, &sinp, &a, -1);
return sinp;
}double tanCordic(double a)
{
/* |a| < 1.74 */
double sinp = 0;
double cosp = invGain1;
cordit1(&cosp, &sinp, &a, -1);
return DIVD(sinp, cosp);
}double asinCordic(double a)
{
/* |a| < 0.98 */
double x, y, z;
x = invGain1;
y = 0;
z = 0; int neg = 1;
if (a < 0) {
a = -a;
neg = 0;
}
cordit1(&x, &y, &z, a);
if (neg) z = -z;
return z;
}/** hyperbolic features ********************************************/
double gain2Cordic()
{
/* calculate hyperbolic gain */
double x, y, z;
x = 1;
y = 0;
z = 0;
cordit2(&x, &y, &z, -1);
return x;
}double sinhcoshCordic(double a, double* coshp)
{
/* |a| < 1.13 */
double y;
*coshp = invGain2;
y = 0;
cordit2(coshp, &y, &a, -1);
return y;
}double tanhCordic(double a)
{
/* |a| < 1.13 */
double sinhp, coshp;
coshp = invGain2;
sinhp = 0;
cordit2(&coshp, &sinhp, &a, -1);
return DIVD(sinhp,coshp);
}double atanhCordic(double a)
{
/* |a| < 1.13 */
double x, z;
x = 1;
z = 0;
cordit2(&x, &a, &z, 0);
return z;
}double logCordic(double a)
{
/* 0.1 < a < 9.58 */
double x, y, z; x = a + 1;
y = a - 1;
z = 0;
cordit2(&x, &y, &z, 0);
return 2*z;
}double sqrtCordic(double a)
{
/* 0.03 < a < 2 */
double x, y, z; x = a + 0.25;
y = a - 0.25;
z = 0;
cordit2(&x, &y, &z, 0);
return MULT(x, invGain2);
}double expCordic(double a)
{
double sinhp, coshp;
coshp = invGain2;
sinhp = 0;
cordit2(&coshp, &sinhp, &a, -1);
return sinhp + coshp;
}main()
{
/* just run a few tests */
double v;
double x;
double c;initCordic();
x = 1;
v = atanCordic(x);
printf("atan %f = %.18e\n", x, v); x = 1;
v = sincosCordic(x, &c);
printf("sin %f = %.18e\n", x, v);
printf("cos %f = %.18e\n", x, c); x = 1;
v = tanCordic(x);
printf("tan %f = %.18e\n", x, v); x = 0.5;
v = asinCordic(x);
printf("asin %f = %.18e\n", x, v); x = 1;
v = sinhcoshCordic(x, &c);
printf("sinh %f = %.18e\n", x, v);
printf("cosh %f = %.18e\n", x, c); x = 1;
v = tanhCordic(x);
printf("tanhh %f = %.18e\n", x, v); x = 0.5;
v = atanhCordic(x);
printf("atanh %f = %.18e\n", x, v); x = 0.8;
v = logCordic(x);
printf("log %f = %.18e\n", x, v); x = 2;
v = sqrtCordic(x);
printf("sqrt %f = %.18e\n", x, v); x = 1;
v = expCordic(x);
printf("exp %f = %.18e\n", x, v);
return 0;
}
the more i look at it, the more im impressed with the hp35. dont forget that it performed all this and a lot more in only 768 instructions and that included code to store the constants as well as the basic arithmetic.
you'll notice that the example functions above have somewhat limited domains. to finish the implementation, trig identities need to be used to bring the input value into the quoted range. so, for example, sin(x) converges for |x| < pi/2, and to make it correct for any x, you must write x = d(pi/2) + d, where |d| < pi/2, and deal separately with the `d' quadrant. similar treatment is needed for logs and square root.
what is interesting is that the above methods require sinh and cosh to compute log and exp. if classic hp models used the same cordic approach, they must have had hyperbolic functionally internally but never put it on keyboard. my theory is they either didnt have the keyboard space, didnt think them wanted or, possibly, their methods did not use the convergence 3k+1 acceleration for the hyperbolic case but used a different trick to ensure the final accuracy of log and exp (eg a single final newton iteration) this being the case would render sinh and cosh too inaccurate to offer to users directly.