fork(1) download
copy 
  1. #include <iostream>
  2. #include <iomanip>
  3. #include <cmath>
  4. using namespace std;
  5.  
  6. long long one_billionth(long long x) {
  7. 	long long y = (x >> 32) * 2305843009LL;
  8. 	y >>= 29;
  9. 	x -= 1000000000LL * y;
  10. 	x += (x >> 4);
  11. 	return (y + (((x >> 29) + 1) >> 1));
  12. }
  13.  
  14. long long crazy_degree_table[] = {
  15.   //  angle        x         y         angle (dms)
  16.           0, 1000000000,        0, //  0: 0: 0
  17.     3599269, 999847757,  17448863, //  0:59:59.269
  18.     7198351, 999391106,  34891507, //  1:59:58.351
  19.    10799869, 998629568,  52335322, //  2:59:59.869
  20.    14392361, 997566633,  69719529, //  3:59:52.361
  21.    18001484, 996194071,  87162910, //  5: 0: 1.484
  22.    21603140, 994520304, 104543603, //  6: 0: 3.14
  23.    25198449, 992547068, 121861880, //  6:59:58.449
  24.    28800050, 990268035, 139173341, //  8: 0: 0.05
  25.    32397717, 987690072, 156423533, //  8:59:57.717
  26.    35999728, 984807982, 173646879, //  9:59:59.728
  27.    39601035, 981626226, 190813921, // 11: 0:1.035
  28.    43199592, 978148012, 207909756, // 11:59:59.592
  29.    46800693, 974369309, 224954328, // 13: 0: 0.693
  30.    50397751, 970298364, 241911316, // 13:59:57.751
  31.    53997701, 965928711, 258808279, // 14:59:57.701
  32.    57600865, 961260540, 275641387, // 16: 0: 0.865
  33.    61199837, 956304987, 292370949, // 16:59:59.837
  34.    64799302, 951057562, 309013776, // 17:59:59.302
  35.    68401066, 945516893, 325573041, // 19: 0: 1.066
  36.    72001124, 939690757, 342025264, // 20: 0: 1.124
  37.    75601453, 933577902, 358374526, // 21: 0: 1.453
  38.    79199676, 927184443, 374605137, // 21:59:59.676
  39.    82800308, 920504270, 390732503, // 23: 0: 0.308
  40.    86400832, 913543817, 406740328, // 24: 0: 0.832
  41.    89998876, 906310090, 422613323, // 24:59:58.876
  42.    93600202, 898793617, 438372027, // 26: 0: 0.202
  43.    97197971, 891010990, 453981735, // 26:59:57.971
  44.   100801394, 882944420, 469477530, // 28: 0: 1.394
  45.   104401002, 874617352, 484813869, // 29: 0: 1.002
  46.   108000141, 866025062, 500000592, // 30: 0: 0.141
  47.   111598780, 857170347, 515033005, // 30:59:58.78
  48.   115200617, 848046511, 529921801, // 32: 0: 0.617
  49.   118802334, 838664405, 544648525, // 33: 0: 2.334
  50.   122398991, 829040308, 559188848, // 33:59:58.991
  51.   125998685, 819155701, 573571214, // 34:59:58.685
  52.   129599798, 809017570, 587784460, // 35:59:59.798
  53.   133200048, 798635370, 601815209, // 37: 0: 0.048
  54.   136800553, 788009103, 615663588, // 38: 0: 0.553
  55.   140401537, 777141272, 629326182, // 39: 0: 1.537
  56.   143999915, 766044708, 642787294, // 39:59:59.915
  57.   147600516, 754707939, 656060917, // 41: 0: 0.516
  58.   151199457, 743146587, 669128650, // 41:59:59.457
  59.   154799036, 731356889, 681994942, // 42:59:59.036
  60.   158401431, 719334981, 694663361, // 44: 0: 1.431
  61.   161995261, 707123027, 707090535  // 44:59:55.261
  62. };
  63.  
  64. long long crazy_minute_table [] = {
  65.   // angle        x         y     angle (dms)
  66.         0, 1000000000,      0, // 0: 0: 0
  67.     59784, 999999958,  289841, // 0: 0:59.784
  68.    119917, 999999831,  581374, // 0: 1:59.917
  69.    180056, 999999619,  872936, // 0: 3: 0.056
  70.    240017, 999999323, 1163635, // 0: 4: 0.017
  71.    300043, 999998942, 1454649, // 0: 5: 0.043
  72.    360108, 999998476, 1745852, // 0: 6: 0.108
  73.    421304, 999997914, 2042538, // 0: 7: 1.304
  74.    479938, 999997293, 2326803, // 0: 7:59.938
  75.    539611, 999996578, 2616105, // 0: 8:59.611
  76.    600013, 999995769, 2908941, // 0:10: 0.013
  77.    659920, 999994882, 3199377, // 0:10:59.92
  78.    720569, 999993898, 3493410, // 0:12: 0.569
  79.    779998, 999992850, 3781528, // 0:12:59.998
  80.    839932, 999991709, 4072094, // 0:13:59.932
  81.    900081, 999990479, 4363702, // 0:15: 0.081
  82.    959964, 999989170, 4654020, // 0:15:59.964
  83.   1020042, 999987772, 4945283, // 0:17: 0.042
  84.   1078867, 999986321, 5230471, // 0:17:58.867
  85.   1139920, 999984729, 5526460, // 0:18:59.92
  86.   1200101, 999983074, 5818221, // 0:20: 0.101
  87.   1260042, 999981341, 6108818, // 0:21: 0.042
  88.   1320065, 999979521, 6399812, // 0:22: 0.065
  89.   1380098, 999977616, 6690854, // 0:23: 0.098
  90.   1439844, 999975636, 6980504, // 0:23:59.844
  91.   1499876, 999973562, 7271540, // 0:24:59.876
  92.   1560081, 999971397, 7563414, // 0:26: 0.081
  93.   1619177, 999969189, 7849911, // 0:26:59.177
  94.   1679867, 999966836, 8144135, // 0:27:59.867
  95.   1740079, 999964416, 8436041, // 0:29: 0.079
  96.   1799931, 999961926, 8726201, // 0:29:59.931
  97.   1859870, 999959348, 9016782  // 0:30:59.87 
  98. };
  99.  
  100. const long ONE_CRAZY_SECOND = 1000L;
  101. const long ONE_CRAZY_MINUTE = 60L * ONE_CRAZY_SECOND;
  102. const long ONE_CRAZY_DEGREE = 60L * ONE_CRAZY_MINUTE;
  103. const long ONE_CRAZY_OCTANT = 45L * ONE_CRAZY_DEGREE;
  104. const long ONE_CRAZY_TURN   =  8L * ONE_CRAZY_OCTANT;
  105. //
  106. // note that 2 * ONE_CRAZY_TURN WILL NOT FIT in a long!!
  107.  
  108. long crazy_sin (long theta) {
  109.   // normalize theta (range: -180  to +180 degrees)
  110.   while (theta < (-4 * ONE_CRAZY_OCTANT)) theta += ONE_CRAZY_TURN;
  111.   while (theta > (4 * ONE_CRAZY_OCTANT)) theta -= ONE_CRAZY_TURN;
  112.   // use an identity to reduce range to: -90 to +90 degrees
  113.   if (theta < (-2 * ONE_CRAZY_OCTANT)) {
  114.     theta = (-4 * ONE_CRAZY_OCTANT) - theta;
  115.   }
  116.   if (theta > (2 * ONE_CRAZY_OCTANT)) {
  117.     theta = (4 * ONE_CRAZY_OCTANT) - theta;
  118.   }
  119.   // where do we look in the table?
  120.   // this is complicated because C++ modulo function is nuts
  121.   int where = theta / ONE_CRAZY_DEGREE;
  122.   if ((theta-(where*ONE_CRAZY_DEGREE))>(30*ONE_CRAZY_MINUTE)) where++;
  123.   else if ((theta-(where*ONE_CRAZY_DEGREE))<(-30*ONE_CRAZY_MINUTE)) where--;
  124.   // now do the table lookup
  125.   long long x; long long y;
  126.   if (where > 90) return -1111111111L;
  127.   else if (where > 45) {
  128.     where = 90 - where;
  129.     theta -= ((90 * ONE_CRAZY_DEGREE) - crazy_degree_table[where*3]);
  130.     y = crazy_degree_table[(where*3)+1];
  131.     x = crazy_degree_table[(where*3)+2];
  132.   } 
  133.   else if (where >= 0) {
  134.     theta -= crazy_degree_table[where*3];
  135.     x = crazy_degree_table[(where*3)+1];
  136.     y = crazy_degree_table[(where*3)+2]; 
  137.   }
  138.   else if (where >= -45) {
  139.     where = -where;
  140.     theta += crazy_degree_table[where*3];
  141.     x = crazy_degree_table[(where*3)+1];
  142.     y = -crazy_degree_table[(where*3)+2];     
  143.   }
  144.   else if (where >= -90) {
  145.     where += 90;
  146.     theta += ((90 * ONE_CRAZY_DEGREE) - crazy_degree_table[where*3]);
  147.     y = -crazy_degree_table[(where*3)+1];
  148.     x = crazy_degree_table[(where*3)+2]; 
  149.   }
  150.   else return -1111111111L;
  151.   // lather, rinse, repeat
  152.   where = theta / ONE_CRAZY_MINUTE;
  153.   if ((theta-(where*ONE_CRAZY_MINUTE))>(30*ONE_CRAZY_SECOND)) where++;
  154.   else if ((theta-(where*ONE_CRAZY_MINUTE))<(-30*ONE_CRAZY_SECOND)) where--;
  155.   // okay, now it gets complicated
  156.   // let's try to squeeze in a few "guard bits"
  157.   x <<= 2;
  158.   y <<= 2;
  159.   // and we need some new temporary storage
  160.   long long x_2;
  161.   long long y_2;
  162.   long long t; 
  163.   // so let's take care of the minutes ...
  164.   if (where > 31) return -1111111111L;
  165.   else if (where >= 0) {
  166.     theta -= crazy_minute_table[where*3];
  167.     x_2 = crazy_minute_table[(where*3)+1];
  168.     y_2 = crazy_minute_table[(where*3)+2]; 
  169.   }
  170.   else if (where >= -31) {
  171.     where = -where;
  172.     theta += crazy_minute_table[where*3];
  173.     x_2 = crazy_minute_table[(where*3)+1];
  174.     y_2 = -crazy_minute_table[(where*3)+2];     
  175.   }
  176.   else return -1111111111L;
  177.   // I told you it would be complicated ...
  178.   t = one_billionth((x * x_2) - (y * y_2));
  179.   y = one_billionth((x * y_2) + (y * x_2));
  180.   x = t;
  181.   // next we take care of, that's right, the seconds
  182.   // our angle unit is 1/1000 arcsecond,
  183.   // which is 4.848136811e-9 of a radian
  184.   t = (((635455LL * theta) >> 15)+1)>>1;
  185.   x_2 = 2000000000LL - one_billionth((t * t) >> 2); // bug squashed
  186.   y_2 = t;
  187.   // wiggle it... just a little bit...
  188.   // t = one_billionth((x * x_2) - (y * y_2));
  189.   y = one_billionth((x * y_2) + (y * x_2));
  190.   // x = t;
  191.   // x += 4;
  192.   y += 4;
  193.   // x >>= 3;
  194.   y >>= 3;
  195.   return y; 
  196. }
  197.  
  198. int main() {
  199. 	long one_deg = 3600000L;
  200. 	long one_step = one_deg/200;
  201. 	double maxerr = 0.0;
  202. 	for (long i = -360*one_deg; i <= 360*one_deg; i += one_step) {
  203. 	    double th = (i*1.0) / one_deg; 
  204. 	    double cs = (crazy_sin(i)*1.0e-9);
  205. 	    double rs = sin(th * 0.0174532925199433);
  206. 	    double err = abs(cs-rs);
  207. 	    if (err > maxerr) maxerr = err;
  208. 	    // cout << th << "  " ;
  209. 	    // cout << setprecision(9) << cs << "  " << rs << "  ";
  210. 	    // cout << err << "\n";
  211. 	}
  212. 	cout << maxerr << "\n";
  213. 	// your code goes here
  214. 	return 0;
  215. }
Success #stdin #stdout 0.03s 3100KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
9.89892e-10
https://ideone.com/qNzyYE
language:
C++14 (gcc 8.3)
created:
8 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!