divtab.c 6.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204
  1. /* Copyright (C) 2003-2015 Free Software Foundation, Inc.
  2. This file is free software; you can redistribute it and/or modify it
  3. under the terms of the GNU General Public License as published by the
  4. Free Software Foundation; either version 3, or (at your option) any
  5. later version.
  6. This file is distributed in the hope that it will be useful, but
  7. WITHOUT ANY WARRANTY; without even the implied warranty of
  8. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  9. General Public License for more details.
  10. Under Section 7 of GPL version 3, you are granted additional
  11. permissions described in the GCC Runtime Library Exception, version
  12. 3.1, as published by the Free Software Foundation.
  13. You should have received a copy of the GNU General Public License and
  14. a copy of the GCC Runtime Library Exception along with this program;
  15. see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
  16. <http://www.gnu.org/licenses/>. */
  17. /* Calculate division table for SH5Media integer division
  18. Contributed by Joern Rennecke
  19. joern.rennecke@superh.com */
  20. #include <stdio.h>
  21. #include <math.h>
  22. #define BITS 5
  23. #define N_ENTRIES (1 << BITS)
  24. #define CUTOFF_BITS 20
  25. #define BIAS (-330)
  26. double max_defect = 0.;
  27. double max_defect_x;
  28. double min_defect = 1e9;
  29. double min_defect_x;
  30. double max_defect2 = 0.;
  31. double max_defect2_x;
  32. double min_defect2 = 0.;
  33. double min_defect2_x;
  34. double min_defect3 = 01e9;
  35. double min_defect3_x;
  36. int min_defect3_val;
  37. double max_defect3 = 0.;
  38. double max_defect3_x;
  39. int max_defect3_val;
  40. static double
  41. note_defect3 (int val, double d2, double y2d, double x)
  42. {
  43. int cutoff_val = val >> CUTOFF_BITS;
  44. double cutoff;
  45. double defect;
  46. if (val < 0)
  47. cutoff_val++;
  48. cutoff = (cutoff_val * (1<<CUTOFF_BITS) - val) * y2d;
  49. defect = cutoff + val * d2;
  50. if (val < 0)
  51. defect = - defect;
  52. if (defect > max_defect3)
  53. {
  54. max_defect3 = defect;
  55. max_defect3_x = x;
  56. max_defect3_val = val;
  57. }
  58. if (defect < min_defect3)
  59. {
  60. min_defect3 = defect;
  61. min_defect3_x = x;
  62. min_defect3_val = val;
  63. }
  64. }
  65. /* This function assumes 32-bit integers. */
  66. static double
  67. calc_defect (double x, int constant, int factor)
  68. {
  69. double y0 = (constant - (int) floor ((x * factor * 64.))) / 16384.;
  70. double y1 = 2 * y0 -y0 * y0 * (x + BIAS / (1.*(1LL<<30)));
  71. double y2d0, y2d;
  72. int y2d1;
  73. double d, d2;
  74. y1 = floor (y1 * (1024 * 1024 * 1024)) / (1024 * 1024 * 1024);
  75. d = y1 - 1 / x;
  76. if (d > max_defect)
  77. {
  78. max_defect = d;
  79. max_defect_x = x;
  80. }
  81. if (d < min_defect)
  82. {
  83. min_defect = d;
  84. min_defect_x = x;
  85. }
  86. y2d0 = floor (y1 * x * (1LL << 60-16));
  87. y2d1 = (int) (long long) y2d0;
  88. y2d = - floor ((y1 - y0 / (1<<30-14)) * y2d1) / (1LL<<44);
  89. d2 = y1 + y2d - 1/x;
  90. if (d2 > max_defect2)
  91. {
  92. max_defect2 = d2;
  93. max_defect2_x = x;
  94. }
  95. if (d2 < min_defect2)
  96. {
  97. min_defect2 = d2;
  98. min_defect2_x = x;
  99. }
  100. /* zero times anything is trivially zero. */
  101. note_defect3 ((1 << CUTOFF_BITS) - 1, d2, y2d, x);
  102. note_defect3 (1 << CUTOFF_BITS, d2, y2d, x);
  103. note_defect3 ((1U << 31) - (1 << CUTOFF_BITS), d2, y2d, x);
  104. note_defect3 ((1U << 31) - 1, d2, y2d, x);
  105. note_defect3 (-1, d2, y2d, x);
  106. note_defect3 (-(1 << CUTOFF_BITS), d2, y2d, x);
  107. note_defect3 ((1U << 31) - (1 << CUTOFF_BITS) + 1, d2, y2d, x);
  108. note_defect3 (-(1U << 31), d2, y2d, x);
  109. return d;
  110. }
  111. int
  112. main ()
  113. {
  114. int i;
  115. unsigned char factors[N_ENTRIES];
  116. short constants[N_ENTRIES];
  117. int steps = N_ENTRIES / 2;
  118. double step = 1. / steps;
  119. double eps30 = 1. / (1024 * 1024 * 1024);
  120. for (i = 0; i < N_ENTRIES; i++)
  121. {
  122. double x_low = (i < steps ? 1. : -3.) + i * step;
  123. double x_high = x_low + step - eps30;
  124. double x_med;
  125. int factor, constant;
  126. double low_defect, med_defect, high_defect, max_defect;
  127. factor = (1./x_low- 1./x_high) / step * 256. + 0.5;
  128. if (factor == 256)
  129. factor = 255;
  130. factors[i] = factor;
  131. /* Use minimum of error function for x_med. */
  132. x_med = sqrt (256./factor);
  133. if (x_low < 0)
  134. x_med = - x_med;
  135. low_defect = 1. / x_low + x_low * factor / 256.;
  136. high_defect = 1. / x_high + x_high * factor / 256.;
  137. med_defect = 1. / x_med + x_med * factor / 256.;
  138. max_defect
  139. = ((low_defect > high_defect) ^ (x_med < 0)) ? low_defect : high_defect;
  140. constant = (med_defect + max_defect) * 0.5 * 16384. + 0.5;
  141. if (constant < -32768 || constant > 32767)
  142. abort ();
  143. constants[i] = constant;
  144. calc_defect (x_low, constant, factor);
  145. calc_defect (x_med, constant, factor);
  146. calc_defect (x_high, constant, factor);
  147. }
  148. printf ("/* This table has been generated by divtab.c .\n");
  149. printf ("Defects for bias %d:\n", BIAS);
  150. printf (" Max defect: %e at %e\n", max_defect, max_defect_x);
  151. printf (" Min defect: %e at %e\n", min_defect, min_defect_x);
  152. printf (" Max 2nd step defect: %e at %e\n", max_defect2, max_defect2_x);
  153. printf (" Min 2nd step defect: %e at %e\n", min_defect2, min_defect2_x);
  154. printf (" Max div defect: %e at %d:%e\n", max_defect3, max_defect3_val,
  155. max_defect3_x);
  156. printf (" Min div defect: %e at %d:%e\n", min_defect3, min_defect3_val,
  157. min_defect3_x);
  158. printf (" Defect at 1: %e\n",
  159. calc_defect (1., constants[0], factors[0]));
  160. printf (" Defect at -2: %e */\n",
  161. calc_defect (-2., constants[steps], factors[steps]));
  162. printf ("\t.section\t.rodata\n");
  163. printf ("\t.balign 2\n");
  164. printf ("/* negative division constants */\n");
  165. for (i = steps; i < 2 * steps; i++)
  166. printf ("\t.word\t%d\n", constants[i]);
  167. printf ("/* negative division factors */\n");
  168. for (i = steps; i < 2*steps; i++)
  169. printf ("\t.byte\t%d\n", factors[i]);
  170. printf ("\t.skip %d\n", steps);
  171. printf ("\t.global GLOBAL(div_table):\n");
  172. printf ("GLOBAL(div_table):\n");
  173. printf ("\t.skip %d\n", steps);
  174. printf ("/* positive division factors */\n");
  175. for (i = 0; i < steps; i++)
  176. printf ("\t.byte\t%d\n", factors[i]);
  177. printf ("/* positive division constants */\n");
  178. for (i = 0; i < steps; i++)
  179. printf ("\t.word\t%d\n", constants[i]);
  180. exit (0);
  181. }