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