OFFSET
1,3
COMMENTS
This sequence is used for division by multiplication with an approximation of the inverse of the divisor on computers if they support multiplying two 32-bit numbers for a 64-bit result. This is usually much faster than a division instruction because integer division is a very slow operation on most computers.
If x and n are unsigned 32-bit quantities, x/n is replaced by (in C notation) ((uint64_t) a(n) * (uint64_t) x) >> b where b is A346496(n).
If n = 2^k, then m=1 and b=k (where the multiplication does not have to be performed).
All a(n) are smaller than 2^33. Those a(n) larger than 2^32, such as a(7), cannot be used directly if the arithmetic is restricted to 32-bit. In this case, the following sequence can be used, in C code again, where all quantities (also intermediate) are unsigned 32-bit integers:
m = a(n) - 2^32; /* Precomputed constant */
mb = b - 33; /* Precomputed constant, b is A346496(n) */
q = ((uint64_t) m * (uint64_t) x) >> 32;
t = (x - q)/2 + q;
res = t >> mb; [Corrected by Baard Nossum, Jul 15 2024]
REFERENCES
Henry S. Warren, Hacker's Delight, 2nd edition, Addison-Wesley, 2013, pp. 230-234, "Unsigned Division by Divisors >= 1"
FORMULA
If n is a power of two, a(n) = 1.
Otherwise, let n_c = 2^32 - (2^32 mod n) - 1 and search for the lowest b so that 2^b > n_c * (n-1- ((2^b-1) mod n)). Then, a(n) = (2^b + n - 1 - ((2^b - 1) mod n))/n, where b is A346496(n).
EXAMPLE
For n=3, m = a(3) = 2863311531 = (2^33 + 1)/3 = AAAAAAAB in base 16, b = A346496(3) = 33, so for 0 <= x < 2^32, x/3 = floor (x*2863311531 / 2^33). For x=11111, a(n)*x = 31814254420941, and floor(31814254420941/2^33) = 3703 = floor(11111/3).
For n=7, on a system where 64-bit arithmetic is not available:
m = 4908534053 - 2^32 = 61356675, so for x=12345, m*x = 757448152875, q=(m*x)>>32 = 1763, x-q = 10582, (x-q)/2 + q = 7054, res = 7054/4 = 1763 = floor(12345/7).
PROG
(PARI) for(n=1, 25, my(k, j=ispower(n, , &k), n_c=2^32-2^32%n-1, b=1); if(j&&k==2, print1(1, ", "), while(b<=n_c*(n-1-(b-1)%n), b+=b); print1((b+n-1-((b-1)%n))/n, ", "))) \\ Hugo Pfoertner, Aug 24 2021
(Python)
def power2(n): return n == 2**(n.bit_length()-1)
def a(n):
if power2(n): return 1
n_c, b = 2**32 - (2**32%n) - 1, 1
while 2**b <= n_c * (n - 1 - ((2**b - 1)%n)): b += 1
return (2**b + n - 1 - ((2**b - 1)%n))//n
print([a(n) for n in range(1, 26)]) # Michael S. Branicky, Aug 24 2021
CROSSREFS
KEYWORD
nonn,easy,fini
AUTHOR
Thomas König, Jul 20 2021
STATUS
approved