Skip to content

Instantly share code, notes, and snippets.

@dougallj
Created May 10, 2020 06:17
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save dougallj/a5ea8c25395dc59bcbb2b443ca7e29d1 to your computer and use it in GitHub Desktop.
#define CONSTEXPR constexpr static inline __attribute__((always_inline))
CONSTEXPR double p2(int x)
{
double powers_of_two[] = {0x1p-1074, 0x1p-1073, 0x1p-1072, 0x1p-1071, 0x1p-1070, 0x1p-1069, 0x1p-1068, 0x1p-1067, 0x1p-1066, 0x1p-1065, 0x1p-1064, 0x1p-1063, 0x1p-1062, 0x1p-1061, 0x1p-1060, 0x1p-1059, 0x1p-1058, 0x1p-1057, 0x1p-1056, 0x1p-1055, 0x1p-1054, 0x1p-1053, 0x1p-1052, 0x1p-1051, 0x1p-1050, 0x1p-1049, 0x1p-1048, 0x1p-1047, 0x1p-1046, 0x1p-1045, 0x1p-1044, 0x1p-1043, 0x1p-1042, 0x1p-1041, 0x1p-1040, 0x1p-1039, 0x1p-1038, 0x1p-1037, 0x1p-1036, 0x1p-1035, 0x1p-1034, 0x1p-1033, 0x1p-1032, 0x1p-1031, 0x1p-1030, 0x1p-1029, 0x1p-1028, 0x1p-1027, 0x1p-1026, 0x1p-1025, 0x1p-1024, 0x1p-1023, 0x1p-1022, 0x1p-1021, 0x1p-1020, 0x1p-1019, 0x1p-1018, 0x1p-1017, 0x1p-1016, 0x1p-1015, 0x1p-1014, 0x1p-1013, 0x1p-1012, 0x1p-1011, 0x1p-1010, 0x1p-1009, 0x1p-1008, 0x1p-1007, 0x1p-1006, 0x1p-1005, 0x1p-1004, 0x1p-1003, 0x1p-1002, 0x1p-1001, 0x1p-1000, 0x1p-999, 0x1p-998, 0x1p-997, 0x1p-996, 0x1p-995, 0x1p-994, 0x1p-993, 0x1p-992, 0x1p-991, 0x1p-990, 0x1p-989, 0x1p-988, 0x1p-987, 0x1p-986, 0x1p-985, 0x1p-984, 0x1p-983, 0x1p-982, 0x1p-981, 0x1p-980, 0x1p-979, 0x1p-978, 0x1p-977, 0x1p-976, 0x1p-975, 0x1p-974, 0x1p-973, 0x1p-972, 0x1p-971, 0x1p-970, 0x1p-969, 0x1p-968, 0x1p-967, 0x1p-966, 0x1p-965, 0x1p-964, 0x1p-963, 0x1p-962, 0x1p-961, 0x1p-960, 0x1p-959, 0x1p-958, 0x1p-957, 0x1p-956, 0x1p-955, 0x1p-954, 0x1p-953, 0x1p-952, 0x1p-951, 0x1p-950, 0x1p-949, 0x1p-948, 0x1p-947, 0x1p-946, 0x1p-945, 0x1p-944, 0x1p-943, 0x1p-942, 0x1p-941, 0x1p-940, 0x1p-939, 0x1p-938, 0x1p-937, 0x1p-936, 0x1p-935, 0x1p-934, 0x1p-933, 0x1p-932, 0x1p-931, 0x1p-930, 0x1p-929, 0x1p-928, 0x1p-927, 0x1p-926, 0x1p-925, 0x1p-924, 0x1p-923, 0x1p-922, 0x1p-921, 0x1p-920, 0x1p-919, 0x1p-918, 0x1p-917, 0x1p-916, 0x1p-915, 0x1p-914, 0x1p-913, 0x1p-912, 0x1p-911, 0x1p-910, 0x1p-909, 0x1p-908, 0x1p-907, 0x1p-906, 0x1p-905, 0x1p-904, 0x1p-903, 0x1p-902, 0x1p-901, 0x1p-900, 0x1p-899, 0x1p-898, 0x1p-897, 0x1p-896, 0x1p-895, 0x1p-894, 0x1p-893, 0x1p-892, 0x1p-891, 0x1p-890, 0x1p-889, 0x1p-888, 0x1p-887, 0x1p-886, 0x1p-885, 0x1p-884, 0x1p-883, 0x1p-882, 0x1p-881, 0x1p-880, 0x1p-879, 0x1p-878, 0x1p-877, 0x1p-876, 0x1p-875, 0x1p-874, 0x1p-873, 0x1p-872, 0x1p-871, 0x1p-870, 0x1p-869, 0x1p-868, 0x1p-867, 0x1p-866, 0x1p-865, 0x1p-864, 0x1p-863, 0x1p-862, 0x1p-861, 0x1p-860, 0x1p-859, 0x1p-858, 0x1p-857, 0x1p-856, 0x1p-855, 0x1p-854, 0x1p-853, 0x1p-852, 0x1p-851, 0x1p-850, 0x1p-849, 0x1p-848, 0x1p-847, 0x1p-846, 0x1p-845, 0x1p-844, 0x1p-843, 0x1p-842, 0x1p-841, 0x1p-840, 0x1p-839, 0x1p-838, 0x1p-837, 0x1p-836, 0x1p-835, 0x1p-834, 0x1p-833, 0x1p-832, 0x1p-831, 0x1p-830, 0x1p-829, 0x1p-828, 0x1p-827, 0x1p-826, 0x1p-825, 0x1p-824, 0x1p-823, 0x1p-822, 0x1p-821, 0x1p-820, 0x1p-819, 0x1p-818, 0x1p-817, 0x1p-816, 0x1p-815, 0x1p-814, 0x1p-813, 0x1p-812, 0x1p-811, 0x1p-810, 0x1p-809, 0x1p-808, 0x1p-807, 0x1p-806, 0x1p-805, 0x1p-804, 0x1p-803, 0x1p-802, 0x1p-801, 0x1p-800, 0x1p-799, 0x1p-798, 0x1p-797, 0x1p-796, 0x1p-795, 0x1p-794, 0x1p-793, 0x1p-792, 0x1p-791, 0x1p-790, 0x1p-789, 0x1p-788, 0x1p-787, 0x1p-786, 0x1p-785, 0x1p-784, 0x1p-783, 0x1p-782, 0x1p-781, 0x1p-780, 0x1p-779, 0x1p-778, 0x1p-777, 0x1p-776, 0x1p-775, 0x1p-774, 0x1p-773, 0x1p-772, 0x1p-771, 0x1p-770, 0x1p-769, 0x1p-768, 0x1p-767, 0x1p-766, 0x1p-765, 0x1p-764, 0x1p-763, 0x1p-762, 0x1p-761, 0x1p-760, 0x1p-759, 0x1p-758, 0x1p-757, 0x1p-756, 0x1p-755, 0x1p-754, 0x1p-753, 0x1p-752, 0x1p-751, 0x1p-750, 0x1p-749, 0x1p-748, 0x1p-747, 0x1p-746, 0x1p-745, 0x1p-744, 0x1p-743, 0x1p-742, 0x1p-741, 0x1p-740, 0x1p-739, 0x1p-738, 0x1p-737, 0x1p-736, 0x1p-735, 0x1p-734, 0x1p-733, 0x1p-732, 0x1p-731, 0x1p-730, 0x1p-729, 0x1p-728, 0x1p-727, 0x1p-726, 0x1p-725, 0x1p-724, 0x1p-723, 0x1p-722, 0x1p-721, 0x1p-720, 0x1p-719, 0x1p-718, 0x1p-717, 0x1p-716, 0x1p-715, 0x1p-714, 0x1p-713, 0x1p-712, 0x1p-711, 0x1p-710, 0x1p-709, 0x1p-708, 0x1p-707, 0x1p-706, 0x1p-705, 0x1p-704, 0x1p-703, 0x1p-702, 0x1p-701, 0x1p-700, 0x1p-699, 0x1p-698, 0x1p-697, 0x1p-696, 0x1p-695, 0x1p-694, 0x1p-693, 0x1p-692, 0x1p-691, 0x1p-690, 0x1p-689, 0x1p-688, 0x1p-687, 0x1p-686, 0x1p-685, 0x1p-684, 0x1p-683, 0x1p-682, 0x1p-681, 0x1p-680, 0x1p-679, 0x1p-678, 0x1p-677, 0x1p-676, 0x1p-675, 0x1p-674, 0x1p-673, 0x1p-672, 0x1p-671, 0x1p-670, 0x1p-669, 0x1p-668, 0x1p-667, 0x1p-666, 0x1p-665, 0x1p-664, 0x1p-663, 0x1p-662, 0x1p-661, 0x1p-660, 0x1p-659, 0x1p-658, 0x1p-657, 0x1p-656, 0x1p-655, 0x1p-654, 0x1p-653, 0x1p-652, 0x1p-651, 0x1p-650, 0x1p-649, 0x1p-648, 0x1p-647, 0x1p-646, 0x1p-645, 0x1p-644, 0x1p-643, 0x1p-642, 0x1p-641, 0x1p-640, 0x1p-639, 0x1p-638, 0x1p-637, 0x1p-636, 0x1p-635, 0x1p-634, 0x1p-633, 0x1p-632, 0x1p-631, 0x1p-630, 0x1p-629, 0x1p-628, 0x1p-627, 0x1p-626, 0x1p-625, 0x1p-624, 0x1p-623, 0x1p-622, 0x1p-621, 0x1p-620, 0x1p-619, 0x1p-618, 0x1p-617, 0x1p-616, 0x1p-615, 0x1p-614, 0x1p-613, 0x1p-612, 0x1p-611, 0x1p-610, 0x1p-609, 0x1p-608, 0x1p-607, 0x1p-606, 0x1p-605, 0x1p-604, 0x1p-603, 0x1p-602, 0x1p-601, 0x1p-600, 0x1p-599, 0x1p-598, 0x1p-597, 0x1p-596, 0x1p-595, 0x1p-594, 0x1p-593, 0x1p-592, 0x1p-591, 0x1p-590, 0x1p-589, 0x1p-588, 0x1p-587, 0x1p-586, 0x1p-585, 0x1p-584, 0x1p-583, 0x1p-582, 0x1p-581, 0x1p-580, 0x1p-579, 0x1p-578, 0x1p-577, 0x1p-576, 0x1p-575, 0x1p-574, 0x1p-573, 0x1p-572, 0x1p-571, 0x1p-570, 0x1p-569, 0x1p-568, 0x1p-567, 0x1p-566, 0x1p-565, 0x1p-564, 0x1p-563, 0x1p-562, 0x1p-561, 0x1p-560, 0x1p-559, 0x1p-558, 0x1p-557, 0x1p-556, 0x1p-555, 0x1p-554, 0x1p-553, 0x1p-552, 0x1p-551, 0x1p-550, 0x1p-549, 0x1p-548, 0x1p-547, 0x1p-546, 0x1p-545, 0x1p-544, 0x1p-543, 0x1p-542, 0x1p-541, 0x1p-540, 0x1p-539, 0x1p-538, 0x1p-537, 0x1p-536, 0x1p-535, 0x1p-534, 0x1p-533, 0x1p-532, 0x1p-531, 0x1p-530, 0x1p-529, 0x1p-528, 0x1p-527, 0x1p-526, 0x1p-525, 0x1p-524, 0x1p-523, 0x1p-522, 0x1p-521, 0x1p-520, 0x1p-519, 0x1p-518, 0x1p-517, 0x1p-516, 0x1p-515, 0x1p-514, 0x1p-513, 0x1p-512, 0x1p-511, 0x1p-510, 0x1p-509, 0x1p-508, 0x1p-507, 0x1p-506, 0x1p-505, 0x1p-504, 0x1p-503, 0x1p-502, 0x1p-501, 0x1p-500, 0x1p-499, 0x1p-498, 0x1p-497, 0x1p-496, 0x1p-495, 0x1p-494, 0x1p-493, 0x1p-492, 0x1p-491, 0x1p-490, 0x1p-489, 0x1p-488, 0x1p-487, 0x1p-486, 0x1p-485, 0x1p-484, 0x1p-483, 0x1p-482, 0x1p-481, 0x1p-480, 0x1p-479, 0x1p-478, 0x1p-477, 0x1p-476, 0x1p-475, 0x1p-474, 0x1p-473, 0x1p-472, 0x1p-471, 0x1p-470, 0x1p-469, 0x1p-468, 0x1p-467, 0x1p-466, 0x1p-465, 0x1p-464, 0x1p-463, 0x1p-462, 0x1p-461, 0x1p-460, 0x1p-459, 0x1p-458, 0x1p-457, 0x1p-456, 0x1p-455, 0x1p-454, 0x1p-453, 0x1p-452, 0x1p-451, 0x1p-450, 0x1p-449, 0x1p-448, 0x1p-447, 0x1p-446, 0x1p-445, 0x1p-444, 0x1p-443, 0x1p-442, 0x1p-441, 0x1p-440, 0x1p-439, 0x1p-438, 0x1p-437, 0x1p-436, 0x1p-435, 0x1p-434, 0x1p-433, 0x1p-432, 0x1p-431, 0x1p-430, 0x1p-429, 0x1p-428, 0x1p-427, 0x1p-426, 0x1p-425, 0x1p-424, 0x1p-423, 0x1p-422, 0x1p-421, 0x1p-420, 0x1p-419, 0x1p-418, 0x1p-417, 0x1p-416, 0x1p-415, 0x1p-414, 0x1p-413, 0x1p-412, 0x1p-411, 0x1p-410, 0x1p-409, 0x1p-408, 0x1p-407, 0x1p-406, 0x1p-405, 0x1p-404, 0x1p-403, 0x1p-402, 0x1p-401, 0x1p-400, 0x1p-399, 0x1p-398, 0x1p-397, 0x1p-396, 0x1p-395, 0x1p-394, 0x1p-393, 0x1p-392, 0x1p-391, 0x1p-390, 0x1p-389, 0x1p-388, 0x1p-387, 0x1p-386, 0x1p-385, 0x1p-384, 0x1p-383, 0x1p-382, 0x1p-381, 0x1p-380, 0x1p-379, 0x1p-378, 0x1p-377, 0x1p-376, 0x1p-375, 0x1p-374, 0x1p-373, 0x1p-372, 0x1p-371, 0x1p-370, 0x1p-369, 0x1p-368, 0x1p-367, 0x1p-366, 0x1p-365, 0x1p-364, 0x1p-363, 0x1p-362, 0x1p-361, 0x1p-360, 0x1p-359, 0x1p-358, 0x1p-357, 0x1p-356, 0x1p-355, 0x1p-354, 0x1p-353, 0x1p-352, 0x1p-351, 0x1p-350, 0x1p-349, 0x1p-348, 0x1p-347, 0x1p-346, 0x1p-345, 0x1p-344, 0x1p-343, 0x1p-342, 0x1p-341, 0x1p-340, 0x1p-339, 0x1p-338, 0x1p-337, 0x1p-336, 0x1p-335, 0x1p-334, 0x1p-333, 0x1p-332, 0x1p-331, 0x1p-330, 0x1p-329, 0x1p-328, 0x1p-327, 0x1p-326, 0x1p-325, 0x1p-324, 0x1p-323, 0x1p-322, 0x1p-321, 0x1p-320, 0x1p-319, 0x1p-318, 0x1p-317, 0x1p-316, 0x1p-315, 0x1p-314, 0x1p-313, 0x1p-312, 0x1p-311, 0x1p-310, 0x1p-309, 0x1p-308, 0x1p-307, 0x1p-306, 0x1p-305, 0x1p-304, 0x1p-303, 0x1p-302, 0x1p-301, 0x1p-300, 0x1p-299, 0x1p-298, 0x1p-297, 0x1p-296, 0x1p-295, 0x1p-294, 0x1p-293, 0x1p-292, 0x1p-291, 0x1p-290, 0x1p-289, 0x1p-288, 0x1p-287, 0x1p-286, 0x1p-285, 0x1p-284, 0x1p-283, 0x1p-282, 0x1p-281, 0x1p-280, 0x1p-279, 0x1p-278, 0x1p-277, 0x1p-276, 0x1p-275, 0x1p-274, 0x1p-273, 0x1p-272, 0x1p-271, 0x1p-270, 0x1p-269, 0x1p-268, 0x1p-267, 0x1p-266, 0x1p-265, 0x1p-264, 0x1p-263, 0x1p-262, 0x1p-261, 0x1p-260, 0x1p-259, 0x1p-258, 0x1p-257, 0x1p-256, 0x1p-255, 0x1p-254, 0x1p-253, 0x1p-252, 0x1p-251, 0x1p-250, 0x1p-249, 0x1p-248, 0x1p-247, 0x1p-246, 0x1p-245, 0x1p-244, 0x1p-243, 0x1p-242, 0x1p-241, 0x1p-240, 0x1p-239, 0x1p-238, 0x1p-237, 0x1p-236, 0x1p-235, 0x1p-234, 0x1p-233, 0x1p-232, 0x1p-231, 0x1p-230, 0x1p-229, 0x1p-228, 0x1p-227, 0x1p-226, 0x1p-225, 0x1p-224, 0x1p-223, 0x1p-222, 0x1p-221, 0x1p-220, 0x1p-219, 0x1p-218, 0x1p-217, 0x1p-216, 0x1p-215, 0x1p-214, 0x1p-213, 0x1p-212, 0x1p-211, 0x1p-210, 0x1p-209, 0x1p-208, 0x1p-207, 0x1p-206, 0x1p-205, 0x1p-204, 0x1p-203, 0x1p-202, 0x1p-201, 0x1p-200, 0x1p-199, 0x1p-198, 0x1p-197, 0x1p-196, 0x1p-195, 0x1p-194, 0x1p-193, 0x1p-192, 0x1p-191, 0x1p-190, 0x1p-189, 0x1p-188, 0x1p-187, 0x1p-186, 0x1p-185, 0x1p-184, 0x1p-183, 0x1p-182, 0x1p-181, 0x1p-180, 0x1p-179, 0x1p-178, 0x1p-177, 0x1p-176, 0x1p-175, 0x1p-174, 0x1p-173, 0x1p-172, 0x1p-171, 0x1p-170, 0x1p-169, 0x1p-168, 0x1p-167, 0x1p-166, 0x1p-165, 0x1p-164, 0x1p-163, 0x1p-162, 0x1p-161, 0x1p-160, 0x1p-159, 0x1p-158, 0x1p-157, 0x1p-156, 0x1p-155, 0x1p-154, 0x1p-153, 0x1p-152, 0x1p-151, 0x1p-150, 0x1p-149, 0x1p-148, 0x1p-147, 0x1p-146, 0x1p-145, 0x1p-144, 0x1p-143, 0x1p-142, 0x1p-141, 0x1p-140, 0x1p-139, 0x1p-138, 0x1p-137, 0x1p-136, 0x1p-135, 0x1p-134, 0x1p-133, 0x1p-132, 0x1p-131, 0x1p-130, 0x1p-129, 0x1p-128, 0x1p-127, 0x1p-126, 0x1p-125, 0x1p-124, 0x1p-123, 0x1p-122, 0x1p-121, 0x1p-120, 0x1p-119, 0x1p-118, 0x1p-117, 0x1p-116, 0x1p-115, 0x1p-114, 0x1p-113, 0x1p-112, 0x1p-111, 0x1p-110, 0x1p-109, 0x1p-108, 0x1p-107, 0x1p-106, 0x1p-105, 0x1p-104, 0x1p-103, 0x1p-102, 0x1p-101, 0x1p-100, 0x1p-99, 0x1p-98, 0x1p-97, 0x1p-96, 0x1p-95, 0x1p-94, 0x1p-93, 0x1p-92, 0x1p-91, 0x1p-90, 0x1p-89, 0x1p-88, 0x1p-87, 0x1p-86, 0x1p-85, 0x1p-84, 0x1p-83, 0x1p-82, 0x1p-81, 0x1p-80, 0x1p-79, 0x1p-78, 0x1p-77, 0x1p-76, 0x1p-75, 0x1p-74, 0x1p-73, 0x1p-72, 0x1p-71, 0x1p-70, 0x1p-69, 0x1p-68, 0x1p-67, 0x1p-66, 0x1p-65, 0x1p-64, 0x1p-63, 0x1p-62, 0x1p-61, 0x1p-60, 0x1p-59, 0x1p-58, 0x1p-57, 0x1p-56, 0x1p-55, 0x1p-54, 0x1p-53, 0x1p-52, 0x1p-51, 0x1p-50, 0x1p-49, 0x1p-48, 0x1p-47, 0x1p-46, 0x1p-45, 0x1p-44, 0x1p-43, 0x1p-42, 0x1p-41, 0x1p-40, 0x1p-39, 0x1p-38, 0x1p-37, 0x1p-36, 0x1p-35, 0x1p-34, 0x1p-33, 0x1p-32, 0x1p-31, 0x1p-30, 0x1p-29, 0x1p-28, 0x1p-27, 0x1p-26, 0x1p-25, 0x1p-24, 0x1p-23, 0x1p-22, 0x1p-21, 0x1p-20, 0x1p-19, 0x1p-18, 0x1p-17, 0x1p-16, 0x1p-15, 0x1p-14, 0x1p-13, 0x1p-12, 0x1p-11, 0x1p-10, 0x1p-9, 0x1p-8, 0x1p-7, 0x1p-6, 0x1p-5, 0x1p-4, 0x1p-3, 0x1p-2, 0x1p-1, 0x1p0, 0x1p1, 0x1p2, 0x1p3, 0x1p4, 0x1p5, 0x1p6, 0x1p7, 0x1p8, 0x1p9, 0x1p10, 0x1p11, 0x1p12, 0x1p13, 0x1p14, 0x1p15, 0x1p16, 0x1p17, 0x1p18, 0x1p19, 0x1p20, 0x1p21, 0x1p22, 0x1p23, 0x1p24, 0x1p25, 0x1p26, 0x1p27, 0x1p28, 0x1p29, 0x1p30, 0x1p31, 0x1p32, 0x1p33, 0x1p34, 0x1p35, 0x1p36, 0x1p37, 0x1p38, 0x1p39, 0x1p40, 0x1p41, 0x1p42, 0x1p43, 0x1p44, 0x1p45, 0x1p46, 0x1p47, 0x1p48, 0x1p49, 0x1p50, 0x1p51, 0x1p52, 0x1p53, 0x1p54, 0x1p55, 0x1p56, 0x1p57, 0x1p58, 0x1p59, 0x1p60, 0x1p61, 0x1p62, 0x1p63, 0x1p64, 0x1p65, 0x1p66, 0x1p67, 0x1p68, 0x1p69, 0x1p70, 0x1p71, 0x1p72, 0x1p73, 0x1p74, 0x1p75, 0x1p76, 0x1p77, 0x1p78, 0x1p79, 0x1p80, 0x1p81, 0x1p82, 0x1p83, 0x1p84, 0x1p85, 0x1p86, 0x1p87, 0x1p88, 0x1p89, 0x1p90, 0x1p91, 0x1p92, 0x1p93, 0x1p94, 0x1p95, 0x1p96, 0x1p97, 0x1p98, 0x1p99, 0x1p100, 0x1p101, 0x1p102, 0x1p103, 0x1p104, 0x1p105, 0x1p106, 0x1p107, 0x1p108, 0x1p109, 0x1p110, 0x1p111, 0x1p112, 0x1p113, 0x1p114, 0x1p115, 0x1p116, 0x1p117, 0x1p118, 0x1p119, 0x1p120, 0x1p121, 0x1p122, 0x1p123, 0x1p124, 0x1p125, 0x1p126, 0x1p127, 0x1p128, 0x1p129, 0x1p130, 0x1p131, 0x1p132, 0x1p133, 0x1p134, 0x1p135, 0x1p136, 0x1p137, 0x1p138, 0x1p139, 0x1p140, 0x1p141, 0x1p142, 0x1p143, 0x1p144, 0x1p145, 0x1p146, 0x1p147, 0x1p148, 0x1p149, 0x1p150, 0x1p151, 0x1p152, 0x1p153, 0x1p154, 0x1p155, 0x1p156, 0x1p157, 0x1p158, 0x1p159, 0x1p160, 0x1p161, 0x1p162, 0x1p163, 0x1p164, 0x1p165, 0x1p166, 0x1p167, 0x1p168, 0x1p169, 0x1p170, 0x1p171, 0x1p172, 0x1p173, 0x1p174, 0x1p175, 0x1p176, 0x1p177, 0x1p178, 0x1p179, 0x1p180, 0x1p181, 0x1p182, 0x1p183, 0x1p184, 0x1p185, 0x1p186, 0x1p187, 0x1p188, 0x1p189, 0x1p190, 0x1p191, 0x1p192, 0x1p193, 0x1p194, 0x1p195, 0x1p196, 0x1p197, 0x1p198, 0x1p199, 0x1p200, 0x1p201, 0x1p202, 0x1p203, 0x1p204, 0x1p205, 0x1p206, 0x1p207, 0x1p208, 0x1p209, 0x1p210, 0x1p211, 0x1p212, 0x1p213, 0x1p214, 0x1p215, 0x1p216, 0x1p217, 0x1p218, 0x1p219, 0x1p220, 0x1p221, 0x1p222, 0x1p223, 0x1p224, 0x1p225, 0x1p226, 0x1p227, 0x1p228, 0x1p229, 0x1p230, 0x1p231, 0x1p232, 0x1p233, 0x1p234, 0x1p235, 0x1p236, 0x1p237, 0x1p238, 0x1p239, 0x1p240, 0x1p241, 0x1p242, 0x1p243, 0x1p244, 0x1p245, 0x1p246, 0x1p247, 0x1p248, 0x1p249, 0x1p250, 0x1p251, 0x1p252, 0x1p253, 0x1p254, 0x1p255, 0x1p256, 0x1p257, 0x1p258, 0x1p259, 0x1p260, 0x1p261, 0x1p262, 0x1p263, 0x1p264, 0x1p265, 0x1p266, 0x1p267, 0x1p268, 0x1p269, 0x1p270, 0x1p271, 0x1p272, 0x1p273, 0x1p274, 0x1p275, 0x1p276, 0x1p277, 0x1p278, 0x1p279, 0x1p280, 0x1p281, 0x1p282, 0x1p283, 0x1p284, 0x1p285, 0x1p286, 0x1p287, 0x1p288, 0x1p289, 0x1p290, 0x1p291, 0x1p292, 0x1p293, 0x1p294, 0x1p295, 0x1p296, 0x1p297, 0x1p298, 0x1p299, 0x1p300, 0x1p301, 0x1p302, 0x1p303, 0x1p304, 0x1p305, 0x1p306, 0x1p307, 0x1p308, 0x1p309, 0x1p310, 0x1p311, 0x1p312, 0x1p313, 0x1p314, 0x1p315, 0x1p316, 0x1p317, 0x1p318, 0x1p319, 0x1p320, 0x1p321, 0x1p322, 0x1p323, 0x1p324, 0x1p325, 0x1p326, 0x1p327, 0x1p328, 0x1p329, 0x1p330, 0x1p331, 0x1p332, 0x1p333, 0x1p334, 0x1p335, 0x1p336, 0x1p337, 0x1p338, 0x1p339, 0x1p340, 0x1p341, 0x1p342, 0x1p343, 0x1p344, 0x1p345, 0x1p346, 0x1p347, 0x1p348, 0x1p349, 0x1p350, 0x1p351, 0x1p352, 0x1p353, 0x1p354, 0x1p355, 0x1p356, 0x1p357, 0x1p358, 0x1p359, 0x1p360, 0x1p361, 0x1p362, 0x1p363, 0x1p364, 0x1p365, 0x1p366, 0x1p367, 0x1p368, 0x1p369, 0x1p370, 0x1p371, 0x1p372, 0x1p373, 0x1p374, 0x1p375, 0x1p376, 0x1p377, 0x1p378, 0x1p379, 0x1p380, 0x1p381, 0x1p382, 0x1p383, 0x1p384, 0x1p385, 0x1p386, 0x1p387, 0x1p388, 0x1p389, 0x1p390, 0x1p391, 0x1p392, 0x1p393, 0x1p394, 0x1p395, 0x1p396, 0x1p397, 0x1p398, 0x1p399, 0x1p400, 0x1p401, 0x1p402, 0x1p403, 0x1p404, 0x1p405, 0x1p406, 0x1p407, 0x1p408, 0x1p409, 0x1p410, 0x1p411, 0x1p412, 0x1p413, 0x1p414, 0x1p415, 0x1p416, 0x1p417, 0x1p418, 0x1p419, 0x1p420, 0x1p421, 0x1p422, 0x1p423, 0x1p424, 0x1p425, 0x1p426, 0x1p427, 0x1p428, 0x1p429, 0x1p430, 0x1p431, 0x1p432, 0x1p433, 0x1p434, 0x1p435, 0x1p436, 0x1p437, 0x1p438, 0x1p439, 0x1p440, 0x1p441, 0x1p442, 0x1p443, 0x1p444, 0x1p445, 0x1p446, 0x1p447, 0x1p448, 0x1p449, 0x1p450, 0x1p451, 0x1p452, 0x1p453, 0x1p454, 0x1p455, 0x1p456, 0x1p457, 0x1p458, 0x1p459, 0x1p460, 0x1p461, 0x1p462, 0x1p463, 0x1p464, 0x1p465, 0x1p466, 0x1p467, 0x1p468, 0x1p469, 0x1p470, 0x1p471, 0x1p472, 0x1p473, 0x1p474, 0x1p475, 0x1p476, 0x1p477, 0x1p478, 0x1p479, 0x1p480, 0x1p481, 0x1p482, 0x1p483, 0x1p484, 0x1p485, 0x1p486, 0x1p487, 0x1p488, 0x1p489, 0x1p490, 0x1p491, 0x1p492, 0x1p493, 0x1p494, 0x1p495, 0x1p496, 0x1p497, 0x1p498, 0x1p499, 0x1p500, 0x1p501, 0x1p502, 0x1p503, 0x1p504, 0x1p505, 0x1p506, 0x1p507, 0x1p508, 0x1p509, 0x1p510, 0x1p511, 0x1p512, 0x1p513, 0x1p514, 0x1p515, 0x1p516, 0x1p517, 0x1p518, 0x1p519, 0x1p520, 0x1p521, 0x1p522, 0x1p523, 0x1p524, 0x1p525, 0x1p526, 0x1p527, 0x1p528, 0x1p529, 0x1p530, 0x1p531, 0x1p532, 0x1p533, 0x1p534, 0x1p535, 0x1p536, 0x1p537, 0x1p538, 0x1p539, 0x1p540, 0x1p541, 0x1p542, 0x1p543, 0x1p544, 0x1p545, 0x1p546, 0x1p547, 0x1p548, 0x1p549, 0x1p550, 0x1p551, 0x1p552, 0x1p553, 0x1p554, 0x1p555, 0x1p556, 0x1p557, 0x1p558, 0x1p559, 0x1p560, 0x1p561, 0x1p562, 0x1p563, 0x1p564, 0x1p565, 0x1p566, 0x1p567, 0x1p568, 0x1p569, 0x1p570, 0x1p571, 0x1p572, 0x1p573, 0x1p574, 0x1p575, 0x1p576, 0x1p577, 0x1p578, 0x1p579, 0x1p580, 0x1p581, 0x1p582, 0x1p583, 0x1p584, 0x1p585, 0x1p586, 0x1p587, 0x1p588, 0x1p589, 0x1p590, 0x1p591, 0x1p592, 0x1p593, 0x1p594, 0x1p595, 0x1p596, 0x1p597, 0x1p598, 0x1p599, 0x1p600, 0x1p601, 0x1p602, 0x1p603, 0x1p604, 0x1p605, 0x1p606, 0x1p607, 0x1p608, 0x1p609, 0x1p610, 0x1p611, 0x1p612, 0x1p613, 0x1p614, 0x1p615, 0x1p616, 0x1p617, 0x1p618, 0x1p619, 0x1p620, 0x1p621, 0x1p622, 0x1p623, 0x1p624, 0x1p625, 0x1p626, 0x1p627, 0x1p628, 0x1p629, 0x1p630, 0x1p631, 0x1p632, 0x1p633, 0x1p634, 0x1p635, 0x1p636, 0x1p637, 0x1p638, 0x1p639, 0x1p640, 0x1p641, 0x1p642, 0x1p643, 0x1p644, 0x1p645, 0x1p646, 0x1p647, 0x1p648, 0x1p649, 0x1p650, 0x1p651, 0x1p652, 0x1p653, 0x1p654, 0x1p655, 0x1p656, 0x1p657, 0x1p658, 0x1p659, 0x1p660, 0x1p661, 0x1p662, 0x1p663, 0x1p664, 0x1p665, 0x1p666, 0x1p667, 0x1p668, 0x1p669, 0x1p670, 0x1p671, 0x1p672, 0x1p673, 0x1p674, 0x1p675, 0x1p676, 0x1p677, 0x1p678, 0x1p679, 0x1p680, 0x1p681, 0x1p682, 0x1p683, 0x1p684, 0x1p685, 0x1p686, 0x1p687, 0x1p688, 0x1p689, 0x1p690, 0x1p691, 0x1p692, 0x1p693, 0x1p694, 0x1p695, 0x1p696, 0x1p697, 0x1p698, 0x1p699, 0x1p700, 0x1p701, 0x1p702, 0x1p703, 0x1p704, 0x1p705, 0x1p706, 0x1p707, 0x1p708, 0x1p709, 0x1p710, 0x1p711, 0x1p712, 0x1p713, 0x1p714, 0x1p715, 0x1p716, 0x1p717, 0x1p718, 0x1p719, 0x1p720, 0x1p721, 0x1p722, 0x1p723, 0x1p724, 0x1p725, 0x1p726, 0x1p727, 0x1p728, 0x1p729, 0x1p730, 0x1p731, 0x1p732, 0x1p733, 0x1p734, 0x1p735, 0x1p736, 0x1p737, 0x1p738, 0x1p739, 0x1p740, 0x1p741, 0x1p742, 0x1p743, 0x1p744, 0x1p745, 0x1p746, 0x1p747, 0x1p748, 0x1p749, 0x1p750, 0x1p751, 0x1p752, 0x1p753, 0x1p754, 0x1p755, 0x1p756, 0x1p757, 0x1p758, 0x1p759, 0x1p760, 0x1p761, 0x1p762, 0x1p763, 0x1p764, 0x1p765, 0x1p766, 0x1p767, 0x1p768, 0x1p769, 0x1p770, 0x1p771, 0x1p772, 0x1p773, 0x1p774, 0x1p775, 0x1p776, 0x1p777, 0x1p778, 0x1p779, 0x1p780, 0x1p781, 0x1p782, 0x1p783, 0x1p784, 0x1p785, 0x1p786, 0x1p787, 0x1p788, 0x1p789, 0x1p790, 0x1p791, 0x1p792, 0x1p793, 0x1p794, 0x1p795, 0x1p796, 0x1p797, 0x1p798, 0x1p799, 0x1p800, 0x1p801, 0x1p802, 0x1p803, 0x1p804, 0x1p805, 0x1p806, 0x1p807, 0x1p808, 0x1p809, 0x1p810, 0x1p811, 0x1p812, 0x1p813, 0x1p814, 0x1p815, 0x1p816, 0x1p817, 0x1p818, 0x1p819, 0x1p820, 0x1p821, 0x1p822, 0x1p823, 0x1p824, 0x1p825, 0x1p826, 0x1p827, 0x1p828, 0x1p829, 0x1p830, 0x1p831, 0x1p832, 0x1p833, 0x1p834, 0x1p835, 0x1p836, 0x1p837, 0x1p838, 0x1p839, 0x1p840, 0x1p841, 0x1p842, 0x1p843, 0x1p844, 0x1p845, 0x1p846, 0x1p847, 0x1p848, 0x1p849, 0x1p850, 0x1p851, 0x1p852, 0x1p853, 0x1p854, 0x1p855, 0x1p856, 0x1p857, 0x1p858, 0x1p859, 0x1p860, 0x1p861, 0x1p862, 0x1p863, 0x1p864, 0x1p865, 0x1p866, 0x1p867, 0x1p868, 0x1p869, 0x1p870, 0x1p871, 0x1p872, 0x1p873, 0x1p874, 0x1p875, 0x1p876, 0x1p877, 0x1p878, 0x1p879, 0x1p880, 0x1p881, 0x1p882, 0x1p883, 0x1p884, 0x1p885, 0x1p886, 0x1p887, 0x1p888, 0x1p889, 0x1p890, 0x1p891, 0x1p892, 0x1p893, 0x1p894, 0x1p895, 0x1p896, 0x1p897, 0x1p898, 0x1p899, 0x1p900, 0x1p901, 0x1p902, 0x1p903, 0x1p904, 0x1p905, 0x1p906, 0x1p907, 0x1p908, 0x1p909, 0x1p910, 0x1p911, 0x1p912, 0x1p913, 0x1p914, 0x1p915, 0x1p916, 0x1p917, 0x1p918, 0x1p919, 0x1p920, 0x1p921, 0x1p922, 0x1p923, 0x1p924, 0x1p925, 0x1p926, 0x1p927, 0x1p928, 0x1p929, 0x1p930, 0x1p931, 0x1p932, 0x1p933, 0x1p934, 0x1p935, 0x1p936, 0x1p937, 0x1p938, 0x1p939, 0x1p940, 0x1p941, 0x1p942, 0x1p943, 0x1p944, 0x1p945, 0x1p946, 0x1p947, 0x1p948, 0x1p949, 0x1p950, 0x1p951, 0x1p952, 0x1p953, 0x1p954, 0x1p955, 0x1p956, 0x1p957, 0x1p958, 0x1p959, 0x1p960, 0x1p961, 0x1p962, 0x1p963, 0x1p964, 0x1p965, 0x1p966, 0x1p967, 0x1p968, 0x1p969, 0x1p970, 0x1p971, 0x1p972, 0x1p973, 0x1p974, 0x1p975, 0x1p976, 0x1p977, 0x1p978, 0x1p979, 0x1p980, 0x1p981, 0x1p982, 0x1p983, 0x1p984, 0x1p985, 0x1p986, 0x1p987, 0x1p988, 0x1p989, 0x1p990, 0x1p991, 0x1p992, 0x1p993, 0x1p994, 0x1p995, 0x1p996, 0x1p997, 0x1p998, 0x1p999, 0x1p1000, 0x1p1001, 0x1p1002, 0x1p1003, 0x1p1004, 0x1p1005, 0x1p1006, 0x1p1007, 0x1p1008, 0x1p1009, 0x1p1010, 0x1p1011, 0x1p1012, 0x1p1013, 0x1p1014, 0x1p1015, 0x1p1016, 0x1p1017, 0x1p1018, 0x1p1019, 0x1p1020, 0x1p1021, 0x1p1022, 0x1p1023};
return powers_of_two[x+1074];
}
CONSTEXPR double double_and(double a, double b) {
return a * b;
}
CONSTEXPR double double_not(double a) {
return 1 - a;
}
CONSTEXPR double double_or(double a, double b) {
return a + b - a * b;
}
CONSTEXPR double select(double condition, double if_true, double if_false) {
return condition * if_true + double_not(condition) * if_false;
}
CONSTEXPR double adding_smallest_is_precise(double x) {
double add_error = x + p2(-1074) - x; // {0, p2(-1074), 2 * p2(-1074)}
add_error -= p2(-1074); // {-p2(-1074), 0, p2(-1074)}
add_error = add_error * p2(1074/2) * p2(1074/2); // {-1, 0, 1}
add_error *= add_error; // {1, 0, 1} - square to make a boolean
return double_not(add_error);
}
CONSTEXPR double is_exp_0_or_1(double x) {
double precise_add = adding_smallest_is_precise(x);
double precise_sub = adding_smallest_is_precise(-x);
return double_and(precise_add, precise_sub);
}
CONSTEXPR double is_zero(double x) {
double magic = p2(-1021) - p2(-1074);
return double_and(is_exp_0_or_1(x + magic), is_exp_0_or_1(x - magic));
}
CONSTEXPR double cmp_sub(double x, double y) {
// return a number with the same sign as x-y (or zero
// if x-y==0), while avoiding returning infinity.
double small = double_or(is_exp_0_or_1(x), is_exp_0_or_1(y));
double multiplier = (small + 1) * p2(-1);
return (x * multiplier) - (y * multiplier);
}
CONSTEXPR double is_equal(double x, double y) {
return is_zero(cmp_sub(x, y));
}
CONSTEXPR double is_exp_0(double x) {
return double_and(is_exp_0_or_1(x + p2(-1022)), is_exp_0_or_1(x - p2(-1022)));
}
CONSTEXPR double get_encoded_exponent(double v)
{
double tmp = v;
double exponent = 0;
#pragma unroll
for (int test = 1024; test >= 1; test /= 2) {
double trial = tmp * p2(-test);
double too_small = is_exp_0_or_1(trial);
tmp = select(too_small, tmp, trial);
exponent += select(too_small, 0, test);
}
return select(is_exp_0_or_1(v), double_not(is_exp_0(v)), exponent + 2);
}
CONSTEXPR double normalise(double v) {
// returns a value with the same sign and fraction bits, but the encoded
// exponent set to one, unless it was a subnormal, in which case it is
// left as zero. probably a bad name.
double res = v;
#pragma unroll
for (int test = 1024; test >= 1; test /= 2) {
double trial = res * p2(-test);
res = select(is_exp_0_or_1(trial), res, trial);
}
return select(is_exp_0_or_1(res), res, res * p2(-1));
}
CONSTEXPR double is_less_than_zero(double v) {
return is_exp_0(normalise(v) + p2(-1022));
}
CONSTEXPR double is_less_than(double a, double b) {
return is_less_than_zero(cmp_sub(a, b));
}
CONSTEXPR double small_positive_floor(double v) {
// WARNING: incorrect for negative numbers and some values over p2(52)
double r = v + p2(52) - p2(52);
return select(is_less_than(v, r), r - 1, r);
}
CONSTEXPR double get_fraction(double v) {
double result = normalise(v) * select(is_less_than_zero(v), -1.0, 1.0);
result -= select(is_exp_0(v), 0, p2(-1022));
result = result * p2(1074 / 2) * p2(1074 / 2);
return result;
}
struct low_high_doubles {
double low;
double high;
};
CONSTEXPR struct low_high_doubles constexpr_double_as_ints(double v)
{
double sign = is_less_than_zero(v);
double exponent = get_encoded_exponent(v);
double fraction = get_fraction(v);
double high_fraction = small_positive_floor(fraction * p2(-32));
double high = sign * p2(31) + exponent * p2(20) + high_fraction;
double low = fraction - high_fraction * p2(32);
return { low, high };
}
CONSTEXPR double double_from_sign_exponent_fraction(double sign, double exponent, double fraction)
{
double exp_is_non_zero = double_not(is_zero(exponent));
// scale fraction down to exponent 0
double v = fraction * p2(-1074);
// add the implicit leading one if needed (so exponent = 1)
v += select(exp_is_non_zero, p2(-1022), 0);
// compute how much we need to increment the exponent by
double e = select(exp_is_non_zero, exponent - 1, 0);
// shift it so that all but the first bit is after the point
e *= p2(-10);
#pragma unroll
for (int test = 1024; test >= 1; test >>= 1) {
// cond will be 1 if the relevant bit is set, otherwise zero
double cond = small_positive_floor(e);
// clear the current bit and shift the next bit into the ones place.
e = (e - cond) * 2;
if (test == 1024) {
// p2(1024) is out-of bounds, so multiply by its square
// root twice
v *= select(cond, p2(512), 1.0);
v *= select(cond, p2(512), 1.0);
} else {
v *= select(cond, p2(test), 1.0);
}
}
// generate a NaN value if one is expected.
double is_nan = double_and(is_equal(exponent, 2047), double_not(is_zero(fraction)));
// if it's a NaN, "v" will already be Infinity, so multiply by zero to make it NaN,
// otherwise multiply by one to leave it as-is.
v *= double_not(is_nan);
// set the sign bit
v *= select(sign, -1, 1);
return v;
}
CONSTEXPR double constexpr_ints_as_double(double l, double h)
{
double exp_and_sign = small_positive_floor(h * p2(-20));
double sign = small_positive_floor(h * p2(-31));
double exponent = exp_and_sign - sign * p2(11);
double fraction = (h - exp_and_sign * p2(20)) * p2(32) + l;
return double_from_sign_exponent_fraction(sign, exponent, fraction);
}
double ints_as_double(double low_u32, double high_u32)
{
return constexpr_ints_as_double(low_u32, high_u32);
}
low_high_doubles double_as_ints(double v)
{
return constexpr_double_as_ints(v);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment