| 880 | #endif |
| 881 | |
| 882 | static uint64_t |
| 883 | mul64_by_fraction(uint64_t a, uint64_t b, uint64_t c) |
| 884 | { |
| 885 | uint64_t acc, bh, bl; |
| 886 | int i, s, sa, sb; |
| 887 | |
| 888 | /* |
| 889 | * Calculate (a * b) / c accurately enough without overflowing. c |
| 890 | * must be nonzero, and its top bit must be 0. a or b must be |
| 891 | * <= c, and the implementation is tuned for b <= c. |
| 892 | * |
| 893 | * The comments about times are for use in calcru1() with units of |
| 894 | * microseconds for 'a' and stathz ticks at 128 Hz for b and c. |
| 895 | * |
| 896 | * Let n be the number of top zero bits in c. Each iteration |
| 897 | * either returns, or reduces b by right shifting it by at least n. |
| 898 | * The number of iterations is at most 1 + 64 / n, and the error is |
| 899 | * at most the number of iterations. |
| 900 | * |
| 901 | * It is very unusual to need even 2 iterations. Previous |
| 902 | * implementations overflowed essentially by returning early in the |
| 903 | * first iteration, with n = 38 giving overflow at 105+ hours and |
| 904 | * n = 32 giving overlow at at 388+ days despite a more careful |
| 905 | * calculation. 388 days is a reasonable uptime, and the calculation |
| 906 | * needs to work for the uptime times the number of CPUs since 'a' |
| 907 | * is per-process. |
| 908 | */ |
| 909 | if (a >= (uint64_t)1 << 63) |
| 910 | return (0); /* Unsupported arg -- can't happen. */ |
| 911 | acc = 0; |
| 912 | for (i = 0; i < 128; i++) { |
| 913 | sa = flsll(a); |
| 914 | sb = flsll(b); |
| 915 | if (sa + sb <= 64) |
| 916 | /* Up to 105 hours on first iteration. */ |
| 917 | return (acc + (a * b) / c); |
| 918 | if (a >= c) { |
| 919 | /* |
| 920 | * This reduction is based on a = q * c + r, with the |
| 921 | * remainder r < c. 'a' may be large to start, and |
| 922 | * moving bits from b into 'a' at the end of the loop |
| 923 | * sets the top bit of 'a', so the reduction makes |
| 924 | * significant progress. |
| 925 | */ |
| 926 | acc += (a / c) * b; |
| 927 | a %= c; |
| 928 | sa = flsll(a); |
| 929 | if (sa + sb <= 64) |
| 930 | /* Up to 388 days on first iteration. */ |
| 931 | return (acc + (a * b) / c); |
| 932 | } |
| 933 | |
| 934 | /* |
| 935 | * This step writes a * b as a * ((bh << s) + bl) = |
| 936 | * a * (bh << s) + a * bl = (a << s) * bh + a * bl. The 2 |
| 937 | * additive terms are handled separately. Splitting in |
| 938 | * this way is linear except for rounding errors. |
| 939 | * |