-
-
Save dougallj/f5b799b718f6c4fa5785434989cd0539 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
constexpr static inline 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 static inline double double_and(double a, double b) | |
{ | |
return a * b; | |
} | |
constexpr static inline double double_not(double a) | |
{ | |
return 1 - a; | |
} | |
constexpr static inline double double_xor(double a, double b) | |
{ | |
double v = (a - b); | |
return v * v; | |
} | |
constexpr static inline double double_or(double a, double b) | |
{ | |
return a + b - a * b; | |
} | |
constexpr static inline double select(double condition, double if_true, double if_false) | |
{ | |
return condition * if_true + double_not(condition) * if_false; | |
} | |
constexpr static inline double adding_smallest_is_precise(double v) { | |
double add_error = v + p2(-1074) - v; // {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); | |
} | |
// assumes round-to-nearest | |
constexpr static inline double is_exp_0_or_1(double v) | |
{ | |
double precise_add = adding_smallest_is_precise(v); | |
double precise_sub = adding_smallest_is_precise(-v); | |
return double_and(precise_add, precise_sub); | |
} | |
// only works for non-huge numbers | |
constexpr static inline double is_exp_around_0_or_1(double v) | |
{ | |
double biased = v - p2(-1074); | |
return (biased + p2(-1074) - biased) * p2(1074 / 2) * p2(1074 / 2); | |
} | |
// only works for non-huge numbers | |
constexpr static inline double unsafe_is_exp_0_or_1(double v) | |
{ | |
return is_exp_around_0_or_1(v * (p2(0) + p2(-52))); | |
} | |
// if we know it's an integer, we can use a much simpler test | |
// (still assumes round-to-nearest) | |
constexpr static inline double is_integer_zero(double v) | |
{ | |
return (v + p2(-1022) - v) * p2(1022); | |
} | |
// if we know it's an integer, we can use a much simpler test | |
// (still assumes round-to-nearest) | |
constexpr static inline double is_small_integer_equal(double a, double b) | |
{ | |
return is_integer_zero(a - b); | |
} | |
constexpr static inline 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 static inline 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 static inline double normalize(double v) | |
{ | |
double res = v; | |
#pragma unroll | |
for (int test = 0x400; test >= 1; test >>= 1) { | |
double trial = res * p2(-test); | |
res = select(test == 0x400 ? is_exp_0_or_1(trial) : unsafe_is_exp_0_or_1(trial), res, trial); | |
} | |
return select(unsafe_is_exp_0_or_1(res), res, res * p2(-1)); | |
} | |
constexpr static inline double is_less_than_zero(double v) | |
{ | |
return is_exp_0_or_1(normalize(v) + p2(-1021)); | |
} | |
constexpr static inline double cmp_sub(double a, double b) | |
{ | |
// return a number with the same sign as a-b (double_or zero | |
// if a-b==0), while avoiding returning infinity. | |
double small = double_or(is_exp_around_0_or_1(a), is_exp_around_0_or_1(b)); | |
double multiplier = (small + 1) * p2(-1); | |
return (a * multiplier) - (b * multiplier); | |
} | |
constexpr static inline double is_equal(double a, double b) | |
{ | |
return is_zero(cmp_sub(a, b)); | |
} | |
constexpr static inline double is_less_than(double a, double b) | |
{ | |
return is_less_than_zero(cmp_sub(a, b)); | |
} | |
constexpr static inline double positive_floor(double v) | |
{ | |
double r = v + p2(52) - p2(52); | |
return r - is_less_than(v, r); | |
} | |
// exhaustively tested for the places we use it, not correct for all doubles | |
constexpr static inline double dirty_floor(double v) | |
{ | |
// for values between 0 and 0x100000 with up to 32 significant bits | |
// after the "decimal" point. | |
return v - (0.5 - p2(-33)) + (p2(52)+1) - (p2(52)+1); | |
} | |
struct sign_exp_frac_doubles { | |
double sign; | |
double exponent; | |
double fraction; | |
}; | |
constexpr static inline sign_exp_frac_doubles double_to_sign_exp_frac(double v) | |
{ | |
double tmp = v; | |
double exponent = 2; | |
double not_subnormal = double_not(is_exp_0(v)); | |
#pragma unroll | |
for (int test = 0x400; test >= 1; test >>= 1) { | |
double trial = tmp * p2(-test); | |
//double too_small = is_exp_0_or_1(trial); | |
double too_small = test == 0x400 ? is_exp_0_or_1(trial) : unsafe_is_exp_0_or_1(trial); | |
tmp = select(too_small, tmp, trial); | |
exponent += double_not(too_small) * test; | |
} | |
exponent = select(is_exp_0_or_1(v), not_subnormal, exponent); | |
tmp = select(unsafe_is_exp_0_or_1(tmp), tmp, tmp * p2(-1)); | |
double sign = is_exp_0_or_1(p2(-1021) + tmp); | |
tmp *= 1 - (sign * 2); | |
tmp = tmp * p2(1074 / 2) * p2(1074 / 2); | |
double fraction = tmp - not_subnormal * p2(52); | |
return { sign, exponent, fraction }; | |
} | |
struct low_high_doubles { | |
double low; | |
double high; | |
}; | |
constexpr static inline low_high_doubles constexpr_double_as_ints(double v) | |
{ | |
sign_exp_frac_doubles values = double_to_sign_exp_frac(v); | |
double sign = values.sign, exponent = values.exponent, fraction = values.fraction; | |
// seems "dirty floor" is just possible, because fraction has one less significant bit | |
// than full double precision provides (as it's missing the implicit leading "1") | |
double high_fraction = dirty_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 static inline double constexpr_ints_as_double(double l, double h) | |
{ | |
double exp_and_sign = dirty_floor(h * p2(-20)); | |
double fraction_high = h * p2(-20) - exp_and_sign; | |
double sign = dirty_floor(exp_and_sign * p2(-11)); | |
double exponent = exp_and_sign * p2(-11) - sign; | |
double exp_is_non_zero = double_not(is_integer_zero(exponent)); | |
double fraction = fraction_high + l * p2(-52); | |
// is_integer_zero works here too, but it's not a great name | |
double v = (fraction + exp_is_non_zero) * p2(-1022); | |
double e = exp_is_non_zero * (exponent - p2(-11)); | |
double cond = 0; | |
e *= 2; | |
#pragma unroll | |
for (int test = 0x400; test >= 2; test >>= 1) { | |
cond = dirty_floor(e); | |
e = (e - cond) * 2; | |
if (test == 0x400) { | |
v *= (p2(512) - 1.0) * cond + 1.0; | |
v *= (p2(512) - 1.0) * cond + 1.0; | |
} else { | |
// surprisingly valid, because | |
// adding 1.0 and subtracting 1.0 | |
// from big powers of two has no effect | |
v *= (p2(test) - 1.0) * cond + 1.0; | |
} | |
} | |
// simplified version of the final iteration | |
v *= e + 1; | |
// carefully generate a NaN value if one is expected. | |
double is_nan = double_and(is_small_integer_equal(exponent, 2047.0 * p2(-11)), double_not(is_integer_zero(fraction))); | |
is_nan = is_nan * p2(512) * p2(512); | |
is_nan -= is_nan; | |
v += is_nan; | |
v *= -sign * 2 + 1; | |
return v; | |
} | |
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