Go Down

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

odometer

Jul 12, 2015, 12:33 pm Last 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

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


el_supremo

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

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

odometer

Just out of curiosity, what is "nuts" about the C++ modulo function?
The way it handles negative numbers.

el_supremo

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