Go Down

### Topic: Crazy trig (8 or 9 good digits) (Read 1 time)previous topic - next topic

#### odometer

##### Jul 12, 2015, 12:33 pmLast Edit: Jul 12, 2015, 02:34 pm by odometer Reason: squashing a bug
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.
Code: [Select]
`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: https://ideone.com/qNzyYE

#### robtillaart

#1
##### Jul 12, 2015, 02:21 pm
check - http://forum.arduino.cc/index.php?topic=69723.0 -

might be inspiring to extend your code
Rob Tillaart

Nederlandse sectie - http://arduino.cc/forum/index.php/board,77.0.html -
(Please do not PM for private consultancy)

#### odometer

#2
##### Jul 12, 2015, 02:44 pm
I was inspired by CORDIC and Gal's accurate tables.

#### el_supremo

#3
##### Jul 12, 2015, 05:57 pm
Just out of curiosity, what is "nuts" about the C++ modulo function?

Pete
Don't send me technical questions via Private Message.

#### odometer

#4
##### Jul 13, 2015, 02:18 am
Just out of curiosity, what is "nuts" about the C++ modulo function?
The way it handles negative numbers.

#### el_supremo

#5
##### Jul 13, 2015, 02:48 am
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
Don't send me technical questions via Private Message.

Go Up