| 1 | /*
|
|---|
| 2 | * Copyright (c) 2012 Adam Hraska
|
|---|
| 3 | * All rights reserved.
|
|---|
| 4 | *
|
|---|
| 5 | * Redistribution and use in source and binary forms, with or without
|
|---|
| 6 | * modification, are permitted provided that the following conditions
|
|---|
| 7 | * are met:
|
|---|
| 8 | *
|
|---|
| 9 | * - Redistributions of source code must retain the above copyright
|
|---|
| 10 | * notice, this list of conditions and the following disclaimer.
|
|---|
| 11 | * - Redistributions in binary form must reproduce the above copyright
|
|---|
| 12 | * notice, this list of conditions and the following disclaimer in the
|
|---|
| 13 | * documentation and/or other materials provided with the distribution.
|
|---|
| 14 | * - The name of the author may not be used to endorse or promote products
|
|---|
| 15 | * derived from this software without specific prior written permission.
|
|---|
| 16 | *
|
|---|
| 17 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
|
|---|
| 18 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
|
|---|
| 19 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
|
|---|
| 20 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
|
|---|
| 21 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
|
|---|
| 22 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
|---|
| 23 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
|---|
| 24 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
|---|
| 25 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
|
|---|
| 26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|---|
| 27 | */
|
|---|
| 28 |
|
|---|
| 29 | #include "private/power_of_ten.h"
|
|---|
| 30 |
|
|---|
| 31 | #include <ieee_double.h>
|
|---|
| 32 | #include <assert.h>
|
|---|
| 33 |
|
|---|
| 34 | /** Precomputed normalized rounded-up powers of 10^k.
|
|---|
| 35 | *
|
|---|
| 36 | * The powers were computed using arbitrary precision arithmetic
|
|---|
| 37 | * and rounded up to the top 64 significant bits. Therefore:
|
|---|
| 38 | * 10^dec_exp == significand * 2^bin_exp +/- 0.5 ulp error
|
|---|
| 39 | *
|
|---|
| 40 | * The smallest interval of binary exponents computed by hand
|
|---|
| 41 | * is [-1083, 987]. Add 200 (exponent change > 3 * 64 bits)
|
|---|
| 42 | * to both bounds just to be on the safe side; ie [-1283, 1187].
|
|---|
| 43 | */
|
|---|
| 44 | static struct {
|
|---|
| 45 | uint64_t significand;
|
|---|
| 46 | int16_t bin_exp;
|
|---|
| 47 | int16_t dec_exp;
|
|---|
| 48 | } fp_powers_of_10[] = {
|
|---|
| 49 | { 0xb8e1cbc28bef0b69ULL, -1286, -368 },
|
|---|
| 50 | { 0x89bf722840327f82ULL, -1259, -360 },
|
|---|
| 51 | { 0xcd42a11346f34f7dULL, -1233, -352 },
|
|---|
| 52 | { 0x98ee4a22ecf3188cULL, -1206, -344 },
|
|---|
| 53 | { 0xe3e27a444d8d98b8ULL, -1180, -336 },
|
|---|
| 54 | { 0xa9c98d8ccb009506ULL, -1153, -328 },
|
|---|
| 55 | { 0xfd00b897478238d1ULL, -1127, -320 },
|
|---|
| 56 | { 0xbc807527ed3e12bdULL, -1100, -312 },
|
|---|
| 57 | { 0x8c71dcd9ba0b4926ULL, -1073, -304 },
|
|---|
| 58 | { 0xd1476e2c07286faaULL, -1047, -296 },
|
|---|
| 59 | { 0x9becce62836ac577ULL, -1020, -288 },
|
|---|
| 60 | { 0xe858ad248f5c22caULL, -994, -280 },
|
|---|
| 61 | { 0xad1c8eab5ee43b67ULL, -967, -272 },
|
|---|
| 62 | { 0x80fa687f881c7f8eULL, -940, -264 },
|
|---|
| 63 | { 0xc0314325637a193aULL, -914, -256 },
|
|---|
| 64 | { 0x8f31cc0937ae58d3ULL, -887, -248 },
|
|---|
| 65 | { 0xd5605fcdcf32e1d7ULL, -861, -240 },
|
|---|
| 66 | { 0x9efa548d26e5a6e2ULL, -834, -232 },
|
|---|
| 67 | { 0xece53cec4a314ebeULL, -808, -224 },
|
|---|
| 68 | { 0xb080392cc4349dedULL, -781, -216 },
|
|---|
| 69 | { 0x8380dea93da4bc60ULL, -754, -208 },
|
|---|
| 70 | { 0xc3f490aa77bd60fdULL, -728, -200 },
|
|---|
| 71 | { 0x91ff83775423cc06ULL, -701, -192 },
|
|---|
| 72 | { 0xd98ddaee19068c76ULL, -675, -184 },
|
|---|
| 73 | { 0xa21727db38cb0030ULL, -648, -176 },
|
|---|
| 74 | { 0xf18899b1bc3f8ca2ULL, -622, -168 },
|
|---|
| 75 | { 0xb3f4e093db73a093ULL, -595, -160 },
|
|---|
| 76 | { 0x8613fd0145877586ULL, -568, -152 },
|
|---|
| 77 | { 0xc7caba6e7c5382c9ULL, -542, -144 },
|
|---|
| 78 | { 0x94db483840b717f0ULL, -515, -136 },
|
|---|
| 79 | { 0xddd0467c64bce4a1ULL, -489, -128 },
|
|---|
| 80 | { 0xa54394fe1eedb8ffULL, -462, -120 },
|
|---|
| 81 | { 0xf64335bcf065d37dULL, -436, -112 },
|
|---|
| 82 | { 0xb77ada0617e3bbcbULL, -409, -104 },
|
|---|
| 83 | { 0x88b402f7fd75539bULL, -382, -96 },
|
|---|
| 84 | { 0xcbb41ef979346bcaULL, -356, -88 },
|
|---|
| 85 | { 0x97c560ba6b0919a6ULL, -329, -80 },
|
|---|
| 86 | { 0xe2280b6c20dd5232ULL, -303, -72 },
|
|---|
| 87 | { 0xa87fea27a539e9a5ULL, -276, -64 },
|
|---|
| 88 | { 0xfb158592be068d2fULL, -250, -56 },
|
|---|
| 89 | { 0xbb127c53b17ec159ULL, -223, -48 },
|
|---|
| 90 | { 0x8b61313bbabce2c6ULL, -196, -40 },
|
|---|
| 91 | { 0xcfb11ead453994baULL, -170, -32 },
|
|---|
| 92 | { 0x9abe14cd44753b53ULL, -143, -24 },
|
|---|
| 93 | { 0xe69594bec44de15bULL, -117, -16 },
|
|---|
| 94 | { 0xabcc77118461cefdULL, -90, -8 },
|
|---|
| 95 | { 0x8000000000000000ULL, -63, 0 },
|
|---|
| 96 | { 0xbebc200000000000ULL, -37, 8 },
|
|---|
| 97 | { 0x8e1bc9bf04000000ULL, -10, 16 },
|
|---|
| 98 | { 0xd3c21bcecceda100ULL, 16, 24 },
|
|---|
| 99 | { 0x9dc5ada82b70b59eULL, 43, 32 },
|
|---|
| 100 | { 0xeb194f8e1ae525fdULL, 69, 40 },
|
|---|
| 101 | { 0xaf298d050e4395d7ULL, 96, 48 },
|
|---|
| 102 | { 0x82818f1281ed44a0ULL, 123, 56 },
|
|---|
| 103 | { 0xc2781f49ffcfa6d5ULL, 149, 64 },
|
|---|
| 104 | { 0x90e40fbeea1d3a4bULL, 176, 72 },
|
|---|
| 105 | { 0xd7e77a8f87daf7fcULL, 202, 80 },
|
|---|
| 106 | { 0xa0dc75f1778e39d6ULL, 229, 88 },
|
|---|
| 107 | { 0xefb3ab16c59b14a3ULL, 255, 96 },
|
|---|
| 108 | { 0xb2977ee300c50fe7ULL, 282, 104 },
|
|---|
| 109 | { 0x850fadc09923329eULL, 309, 112 },
|
|---|
| 110 | { 0xc646d63501a1511eULL, 335, 120 },
|
|---|
| 111 | { 0x93ba47c980e98ce0ULL, 362, 128 },
|
|---|
| 112 | { 0xdc21a1171d42645dULL, 388, 136 },
|
|---|
| 113 | { 0xa402b9c5a8d3a6e7ULL, 415, 144 },
|
|---|
| 114 | { 0xf46518c2ef5b8cd1ULL, 441, 152 },
|
|---|
| 115 | { 0xb616a12b7fe617aaULL, 468, 160 },
|
|---|
| 116 | { 0x87aa9aff79042287ULL, 495, 168 },
|
|---|
| 117 | { 0xca28a291859bbf93ULL, 521, 176 },
|
|---|
| 118 | { 0x969eb7c47859e744ULL, 548, 184 },
|
|---|
| 119 | { 0xe070f78d3927556bULL, 574, 192 },
|
|---|
| 120 | { 0xa738c6bebb12d16dULL, 601, 200 },
|
|---|
| 121 | { 0xf92e0c3537826146ULL, 627, 208 },
|
|---|
| 122 | { 0xb9a74a0637ce2ee1ULL, 654, 216 },
|
|---|
| 123 | { 0x8a5296ffe33cc930ULL, 681, 224 },
|
|---|
| 124 | { 0xce1de40642e3f4b9ULL, 707, 232 },
|
|---|
| 125 | { 0x9991a6f3d6bf1766ULL, 734, 240 },
|
|---|
| 126 | { 0xe4d5e82392a40515ULL, 760, 248 },
|
|---|
| 127 | { 0xaa7eebfb9df9de8eULL, 787, 256 },
|
|---|
| 128 | { 0xfe0efb53d30dd4d8ULL, 813, 264 },
|
|---|
| 129 | { 0xbd49d14aa79dbc82ULL, 840, 272 },
|
|---|
| 130 | { 0x8d07e33455637eb3ULL, 867, 280 },
|
|---|
| 131 | { 0xd226fc195c6a2f8cULL, 893, 288 },
|
|---|
| 132 | { 0x9c935e00d4b9d8d2ULL, 920, 296 },
|
|---|
| 133 | { 0xe950df20247c83fdULL, 946, 304 },
|
|---|
| 134 | { 0xadd57a27d29339f6ULL, 973, 312 },
|
|---|
| 135 | { 0x81842f29f2cce376ULL, 1000, 320 },
|
|---|
| 136 | { 0xc0fe908895cf3b44ULL, 1026, 328 },
|
|---|
| 137 | { 0x8fcac257558ee4e6ULL, 1053, 336 },
|
|---|
| 138 | { 0xd6444e39c3db9b0aULL, 1079, 344 },
|
|---|
| 139 | { 0x9fa42700db900ad2ULL, 1106, 352 },
|
|---|
| 140 | { 0xede24ae798ec8284ULL, 1132, 360 },
|
|---|
| 141 | { 0xb13cc3832ef0c9acULL, 1159, 368 },
|
|---|
| 142 | { 0x840d57e2899d945fULL, 1186, 376 }
|
|---|
| 143 | };
|
|---|
| 144 |
|
|---|
| 145 | /**
|
|---|
| 146 | * Returns the smallest precomputed power of 10 such that
|
|---|
| 147 | * binary_exp <= power_of_10.bin_exp
|
|---|
| 148 | * where
|
|---|
| 149 | * 10^decimal_exp = power_of_10.significand * 2^bin_exp
|
|---|
| 150 | * with an error of 0.5 ulp in the significand.
|
|---|
| 151 | */
|
|---|
| 152 | void get_power_of_ten(int binary_exp, fp_num_t *power_of_10, int *decimal_exp)
|
|---|
| 153 | {
|
|---|
| 154 | const int powers_count = sizeof(fp_powers_of_10) / sizeof(fp_powers_of_10[0]);
|
|---|
| 155 | const int min_bin_exp = fp_powers_of_10[0].bin_exp;
|
|---|
| 156 | const int max_bin_exp = fp_powers_of_10[powers_count - 1].bin_exp;
|
|---|
| 157 | const int max_bin_exp_diff = 27;
|
|---|
| 158 |
|
|---|
| 159 | assert(min_bin_exp <= binary_exp && binary_exp <= max_bin_exp);
|
|---|
| 160 |
|
|---|
| 161 | /*
|
|---|
| 162 | * Binary exponent difference between adjacent powers of 10
|
|---|
| 163 | * is lg(10^8) = 26.575. The starting search index seed_idx
|
|---|
| 164 | * undershoots the actual position by less than 1.6%, ie it
|
|---|
| 165 | * skips 26.575/27 = 98.4% of all the smaller powers. This
|
|---|
| 166 | * translates to at most three extra tests.
|
|---|
| 167 | */
|
|---|
| 168 | int seed_idx = (binary_exp - min_bin_exp) / max_bin_exp_diff;
|
|---|
| 169 |
|
|---|
| 170 | assert(fp_powers_of_10[seed_idx].bin_exp < binary_exp);
|
|---|
| 171 |
|
|---|
| 172 | for (int i = seed_idx; i < powers_count; ++i) {
|
|---|
| 173 | /* Found the smallest power of 10 with bin_exp >= binary_exp. */
|
|---|
| 174 | if (binary_exp <= fp_powers_of_10[i].bin_exp) {
|
|---|
| 175 | assert(fp_powers_of_10[i].bin_exp <= binary_exp + max_bin_exp_diff);
|
|---|
| 176 |
|
|---|
| 177 | power_of_10->significand = fp_powers_of_10[i].significand;
|
|---|
| 178 | power_of_10->exponent = fp_powers_of_10[i].bin_exp;
|
|---|
| 179 | *decimal_exp = fp_powers_of_10[i].dec_exp;
|
|---|
| 180 | return;
|
|---|
| 181 | }
|
|---|
| 182 | }
|
|---|
| 183 |
|
|---|
| 184 | assert(false);
|
|---|
| 185 | }
|
|---|
