# Crazy trig (8 or 9 good digits)

I can hardly believe I stayed up all night writing this wretched thing.

A lot of people say that Arduino-style floats aren’t good enough: they only give 7 good digits (8 if you’re lucky).

Here is this absurd thing I wrote in C++. It uses no floats or doubles, only integers (admittedly including some long long, i.e. 64 bit, integers).

The way I have it now, it only calculates sine.
The angle is expressed in units such that 1000 units = 1 arcsecond. Thus, 3600000 units make 1 degree, and 1296000000 units make 1 full turn.
It returns a (32-bit) long integer, in the range -1000000000 to 1000000000 (negative one billion to positive one billion). Sometimes the result is off by 1 or 2 units, I ought to investigate. But not now; I’m too tired.

Anyway, I got it to give me 8 good digits more or less consistently, and often 9 good digits.

Here you go:
EDIT: I found a stupid bug that was killing the accuracy on angles near 90 degrees. I squashed that bug, and now I get 9 good digits rather consistently.

``````long long one_billionth(long long x) {
long long y = (x >> 32) * 2305843009LL;
y >>= 29;
x -= 1000000000LL * y;
x += (x >> 4);
return (y + (((x >> 29) + 1) >> 1));
}

long long crazy_degree_table[] = {
//  angle        x         y         angle (dms)
0, 1000000000,        0, //  0: 0: 0
3599269, 999847757,  17448863, //  0:59:59.269
7198351, 999391106,  34891507, //  1:59:58.351
10799869, 998629568,  52335322, //  2:59:59.869
14392361, 997566633,  69719529, //  3:59:52.361
18001484, 996194071,  87162910, //  5: 0: 1.484
21603140, 994520304, 104543603, //  6: 0: 3.14
25198449, 992547068, 121861880, //  6:59:58.449
28800050, 990268035, 139173341, //  8: 0: 0.05
32397717, 987690072, 156423533, //  8:59:57.717
35999728, 984807982, 173646879, //  9:59:59.728
39601035, 981626226, 190813921, // 11: 0:1.035
43199592, 978148012, 207909756, // 11:59:59.592
46800693, 974369309, 224954328, // 13: 0: 0.693
50397751, 970298364, 241911316, // 13:59:57.751
53997701, 965928711, 258808279, // 14:59:57.701
57600865, 961260540, 275641387, // 16: 0: 0.865
61199837, 956304987, 292370949, // 16:59:59.837
64799302, 951057562, 309013776, // 17:59:59.302
68401066, 945516893, 325573041, // 19: 0: 1.066
72001124, 939690757, 342025264, // 20: 0: 1.124
75601453, 933577902, 358374526, // 21: 0: 1.453
79199676, 927184443, 374605137, // 21:59:59.676
82800308, 920504270, 390732503, // 23: 0: 0.308
86400832, 913543817, 406740328, // 24: 0: 0.832
89998876, 906310090, 422613323, // 24:59:58.876
93600202, 898793617, 438372027, // 26: 0: 0.202
97197971, 891010990, 453981735, // 26:59:57.971
100801394, 882944420, 469477530, // 28: 0: 1.394
104401002, 874617352, 484813869, // 29: 0: 1.002
108000141, 866025062, 500000592, // 30: 0: 0.141
111598780, 857170347, 515033005, // 30:59:58.78
115200617, 848046511, 529921801, // 32: 0: 0.617
118802334, 838664405, 544648525, // 33: 0: 2.334
122398991, 829040308, 559188848, // 33:59:58.991
125998685, 819155701, 573571214, // 34:59:58.685
129599798, 809017570, 587784460, // 35:59:59.798
133200048, 798635370, 601815209, // 37: 0: 0.048
136800553, 788009103, 615663588, // 38: 0: 0.553
140401537, 777141272, 629326182, // 39: 0: 1.537
143999915, 766044708, 642787294, // 39:59:59.915
147600516, 754707939, 656060917, // 41: 0: 0.516
151199457, 743146587, 669128650, // 41:59:59.457
154799036, 731356889, 681994942, // 42:59:59.036
158401431, 719334981, 694663361, // 44: 0: 1.431
161995261, 707123027, 707090535  // 44:59:55.261
};

long long crazy_minute_table [] = {
// angle        x         y     angle (dms)
0, 1000000000,      0, // 0: 0: 0
59784, 999999958,  289841, // 0: 0:59.784
119917, 999999831,  581374, // 0: 1:59.917
180056, 999999619,  872936, // 0: 3: 0.056
240017, 999999323, 1163635, // 0: 4: 0.017
300043, 999998942, 1454649, // 0: 5: 0.043
360108, 999998476, 1745852, // 0: 6: 0.108
421304, 999997914, 2042538, // 0: 7: 1.304
479938, 999997293, 2326803, // 0: 7:59.938
539611, 999996578, 2616105, // 0: 8:59.611
600013, 999995769, 2908941, // 0:10: 0.013
659920, 999994882, 3199377, // 0:10:59.92
720569, 999993898, 3493410, // 0:12: 0.569
779998, 999992850, 3781528, // 0:12:59.998
839932, 999991709, 4072094, // 0:13:59.932
900081, 999990479, 4363702, // 0:15: 0.081
959964, 999989170, 4654020, // 0:15:59.964
1020042, 999987772, 4945283, // 0:17: 0.042
1078867, 999986321, 5230471, // 0:17:58.867
1139920, 999984729, 5526460, // 0:18:59.92
1200101, 999983074, 5818221, // 0:20: 0.101
1260042, 999981341, 6108818, // 0:21: 0.042
1320065, 999979521, 6399812, // 0:22: 0.065
1380098, 999977616, 6690854, // 0:23: 0.098
1439844, 999975636, 6980504, // 0:23:59.844
1499876, 999973562, 7271540, // 0:24:59.876
1560081, 999971397, 7563414, // 0:26: 0.081
1619177, 999969189, 7849911, // 0:26:59.177
1679867, 999966836, 8144135, // 0:27:59.867
1740079, 999964416, 8436041, // 0:29: 0.079
1799931, 999961926, 8726201, // 0:29:59.931
1859870, 999959348, 9016782  // 0:30:59.87
};

const long ONE_CRAZY_SECOND = 1000L;
const long ONE_CRAZY_MINUTE = 60L * ONE_CRAZY_SECOND;
const long ONE_CRAZY_DEGREE = 60L * ONE_CRAZY_MINUTE;
const long ONE_CRAZY_OCTANT = 45L * ONE_CRAZY_DEGREE;
const long ONE_CRAZY_TURN   =  8L * ONE_CRAZY_OCTANT;
//
// note that 2 * ONE_CRAZY_TURN WILL NOT FIT in a long!!

long crazy_sin (long theta) {
// normalize theta (range: -180  to +180 degrees)
while (theta < (-4 * ONE_CRAZY_OCTANT)) theta += ONE_CRAZY_TURN;
while (theta > (4 * ONE_CRAZY_OCTANT)) theta -= ONE_CRAZY_TURN;
// use an identity to reduce range to: -90 to +90 degrees
if (theta < (-2 * ONE_CRAZY_OCTANT)) {
theta = (-4 * ONE_CRAZY_OCTANT) - theta;
}
if (theta > (2 * ONE_CRAZY_OCTANT)) {
theta = (4 * ONE_CRAZY_OCTANT) - theta;
}
// where do we look in the table?
// this is complicated because C++ modulo function is nuts
int where = theta / ONE_CRAZY_DEGREE;
if ((theta-(where*ONE_CRAZY_DEGREE))>(30*ONE_CRAZY_MINUTE)) where++;
else if ((theta-(where*ONE_CRAZY_DEGREE))<(-30*ONE_CRAZY_MINUTE)) where--;
// now do the table lookup
long long x; long long y;
if (where > 90) return -1111111111L;
else if (where > 45) {
where = 90 - where;
theta -= ((90 * ONE_CRAZY_DEGREE) - crazy_degree_table[where*3]);
y = crazy_degree_table[(where*3)+1];
x = crazy_degree_table[(where*3)+2];
}
else if (where >= 0) {
theta -= crazy_degree_table[where*3];
x = crazy_degree_table[(where*3)+1];
y = crazy_degree_table[(where*3)+2];
}
else if (where >= -45) {
where = -where;
theta += crazy_degree_table[where*3];
x = crazy_degree_table[(where*3)+1];
y = -crazy_degree_table[(where*3)+2];
}
else if (where >= -90) {
where += 90;
theta += ((90 * ONE_CRAZY_DEGREE) - crazy_degree_table[where*3]);
y = -crazy_degree_table[(where*3)+1];
x = crazy_degree_table[(where*3)+2];
}
else return -1111111111L;
// lather, rinse, repeat
where = theta / ONE_CRAZY_MINUTE;
if ((theta-(where*ONE_CRAZY_MINUTE))>(30*ONE_CRAZY_SECOND)) where++;
else if ((theta-(where*ONE_CRAZY_MINUTE))<(-30*ONE_CRAZY_SECOND)) where--;
// okay, now it gets complicated
// let's try to squeeze in a few "guard bits"
x <<= 2;
y <<= 2;
// and we need some new temporary storage
long long x_2;
long long y_2;
long long t;
// so let's take care of the minutes ...
if (where > 31) return -1111111111L;
else if (where >= 0) {
theta -= crazy_minute_table[where*3];
x_2 = crazy_minute_table[(where*3)+1];
y_2 = crazy_minute_table[(where*3)+2];
}
else if (where >= -31) {
where = -where;
theta += crazy_minute_table[where*3];
x_2 = crazy_minute_table[(where*3)+1];
y_2 = -crazy_minute_table[(where*3)+2];
}
else return -1111111111L;
// I told you it would be complicated ...
t = one_billionth((x * x_2) - (y * y_2));
y = one_billionth((x * y_2) + (y * x_2));
x = t;
// next we take care of, that's right, the seconds
// our angle unit is 1/1000 arcsecond,
// which is 4.848136811e-9 of a radian
t = (((635455LL * theta) >> 15)+1)>>1;
x_2 = 2000000000LL - one_billionth((t * t) >> 2); // bug squashed
y_2 = t;
// wiggle it... just a little bit...
// t = one_billionth((x * x_2) - (y * y_2));
y = one_billionth((x * y_2) + (y * x_2));
// x = t;
// x += 4;
y += 4;
// x >>= 3;
y >>= 3;
return y;
}
``````

Link to my silly test here: qNzyYE - Online C++0x Compiler & Debugging Tool - Ideone.com

might be inspiring to extend your code

I was inspired by CORDIC and Gal's accurate tables.

Just out of curiosity, what is "nuts" about the C++ modulo function?

Pete

el_supremo: Just out of curiosity, what is "nuts" about the C++ modulo function?

The way it handles negative numbers.

The problem IIRC is that if either or both of the operands of % are negative then the operation is undefined and will give a machine-dependent result. On the Arduino it appears that the result of a%b has the correct magnitude and the same sign as a. If my quick tests can be extrapolated to all possible a and b, you can just calculate abs(a%b) and the result will be correct.

Pete