# Find the magic numbers to divide a number without division

An integer $$\x\in[0,2^{32}-1]\$$ divided by an integer $$\d\in{[1,2^{31}]}\$$ will produce an integral quotient $$\q\$$ and a remainder $$\r\$$, so that $$\x=d\times q+r\$$.

Any $$\q\$$, in fact, can be calculated by the following formula.

$$y=\left\lfloor\frac{mx}{2^{32}}\right\rfloor\\q=\left\lfloor\frac{y+f\lfloor(x-y)/2\rfloor}{2^s}\right\rfloor$$

where

$$m,s,f\in Integer\\m\in[2,2^{32}]\\s\in[0,31]\\f\in\{0,1\}$$

The formula can be rewritten with bitwise operators, to be more friendly with computers.

$$y=mx\gg32\\q=y+(x-y\gg1\ \&\ {-f})\gg s$$

Given $$\d\$$, produce a valid tuple $$\(m,s,f)\$$, so that applying the formula above will compute the correct value of $$\q\$$ for all $$\x\$$.

There may be multiple valid solutions, in which case you can output any number of any valid solution.

### Test Cases

(d -> m, s, f)
1 -> 4294967296, 0, 0
2 -> 2147483648, 0, 0
3 -> 2863311531, 1, 0
7 -> 613566757, 2, 1
11 -> 3123612579, 3, 0
111 -> 657787785, 6, 1
1111 -> 989659431, 8, 0
11111 -> 2038283013, 13, 1
1111111111 -> 2075258709, 29, 0
2147483647 -> 3, 30, 1
2147483648 -> 2, 0, 0

Each of the listed output is not the only possible solution. For example, both of these are correct.

1111111111 -> 2075258709, 29, 0
1111111111 -> 4150517417, 30, 0

# Python 3.8 (pre-release), 271 bytes

lambda d,P=2**32:[s:=min(range(32),key=lambda n:abs(ceil((P<<n)/d)*d-(P<<n))),m:=ceil((P<<s)/d),0]*(1<m<=P and m*d>>(32+s))or[s:=min([n for n in range(32)if 1<2**(33+n)/d-P<=P],key=lambda n:abs(ceil(2**(33+n)/d)*d-(2*P<<n))),m:=floor(2**(33+s)/d+1)-P,1]
from math import*

Try it online!

Probably severely golfable.

Returns a list of three elements in the order [s, m, f].

Here is a calculator that uses my function. You input a Python division expression (dividend/divisor), and it uses my function to return a quotient without directly dividing the dividend by the divisor.

Thanks @pxeger for suggesting the use of bitwise operators (2**(32+s) -> P<<s) in some places, which saves ten bytes.

## Explanation

The problem can be solved by comparing the "magic quotient" expression to the actual quotient. Since f is either zero or one, we can split the problem into two parts.

If f is zero, then the "magic quotient" expression becomes:

$$\left\lfloor \frac{y}{2^s} \right\rfloor = \left\lfloor \frac{\left\lfloor \frac{mx}{2^{32}} \right\rfloor}{2^s} \right\rfloor$$

We can ignore the floor functions so that the problem becomes one of approximation rather than complete accuracy. Then, the magic quotient becomes:

$$\frac{mx}{2^{32+s}}$$

What we want to do is find $$\m\$$ and $$\s\$$ such that the value above most closely approximates $$\\frac{x}{d}\$$. We want

$$\frac{mx}{2^{32+s}} \approx \frac{x}{d} \implies mdx \approx x2^{32+s} \implies md \approx 2^{32+s}$$

This problem is a lot easier to solve than the daunting task of dealing with floors and correcting errors, or, even worse, dealing with those pesky fractions. In fact, this part is very simple. We can just loop through the 32 possible values of $$\s\$$ and find the one with the lowest

$$\left\lceil \frac{2^{32+s}}{d} \right\rceil \times d - 2^{32+s}$$

We use the ceiling, and not the floor, because later on we use the floor to compute the quotient, and the ceiling allows us to get a more precise value.

The above expression is evaluated for each $$\s\$$ from 0 to 31, and the $$\s\$$ for which that expression is the lowest is the one that is to be outputted (or is at least considered for output). Next, we evaluate $$\m\$$, which is the ceiling expression above.

Then $$\f\$$ is zero (assumed above).

Now, of course, $$\f\$$ cannot always be zero, so some checks are done to determine whether or not $$\f = 0\$$ in fact works. $$\f = 0\$$ works if:

• $$\m\$$ is within the required range; and
• the magic quotient of $$\d\$$ by $$\d\$$ is equal to one (if it isn't, then we need some additional value in the numerator, which is where $$\f = 1\$$ comes in).

The second condition is sufficient because, considering some nonnegative integer $$\ad+b\$$ where $$\a\$$ and $$\b\$$ are nonnegative integers, $$\b\$$ being less than $$\d\$$, the value of

$$\left\lfloor \frac{\left\lfloor \frac{mad+mb}{2^{32}} \right\rfloor}{2^s} \right\rfloor = \left\lfloor \frac{\left\lfloor \frac{amd}{2^{32}} \right\rfloor}{2^s} \right\rfloor$$

This is because $$\mb\$$ is less than $$\md\$$, so $$\\frac{mb}{2^{32}}\$$ is less than one.

If the magic quotient of $$\d\$$ by $$\d\$$ is one, then

$$\left\lfloor \frac{\left\lfloor \frac{md}{2^{32}} \right\rfloor}{2^s} \right\rfloor = 1$$

Consequently,

$$\left\lfloor \frac{\left\lfloor \frac{amd}{2^{32}} \right\rfloor}{2^s} \right\rfloor = a \left\lfloor \frac{\left\lfloor \frac{md}{2^{32}} \right\rfloor}{2^s} \right\rfloor = a$$

which is indeed the integer part of $$\\frac{ad+b}{d}\$$. Since this works for every $$\ad+b\$$, the $$\s\$$, $$\m\$$ and $$\f\$$ values are satisfactory.

Otherwise, we need to consider $$\f = 1\$$. This value of $$\f\$$ will definitely work.

The problem is more complicated if $$\f = 1\$$, but it is still a problem of approximation. I'll leave the work of finding the simple equation to you, but here it is for those who don't have the time / are too lazy to do the work themselves:

$$d(m+2^{32}) \approx 2^{33+s}$$

From here, the work is very similar to before, except that since we must get a good value, we restrict the search space to the values of $$\s\$$ such that the $$\m\$$ corresponding to $$\s\$$ is within the required range. Another difference is that $$\m\$$ is found using a Python equivalent of -~ rather than ceil, the difference between them being that -~ increments integers whereas ceil leaves them unchanged. I'm still unsure as to why this is needed, but it works so I'm not going to question it.

I don't do these long mathematical explanations very often, so feel free to critique my mathematics. If the LaTEX looks weird, please do tell me so I can fix it. My reasoning might also be wrong, so if it is, please tell me.

• The 2**(32+x)s can be replaced with P<<x (in a few places you need brackets though) Commented Mar 8, 2022 at 12:34
• @pxeger Yes ;) I'll keep golfing as I go Commented Mar 8, 2022 at 12:35
• This might be interesting to you. See libdivide_internal_u32_gen and libdivide_u32_do Commented Mar 8, 2022 at 14:55
• 1<m<=P and can be P>=m>1and for -1 Commented Mar 8, 2022 at 14:56
• I've fixed the floors/ceils in the MathJax so it uses the full vertical height. If you feel the previous looked better, feel free to revert it back. Commented Mar 8, 2022 at 15:09

# Charcoal, 64 61 bytes

Ｎθ≔Ｘ²¦³²η≔⌕ＥＸ²…³²¦⁶⁴÷×⁻⊖η﹪ηθ﹪±ιθι⁰ζ≔÷×ηＸ²ζθε≔›εηδＩ⟦⁻ε×δ⊖η⁻ζδδ

Try it online! Link is to verbose version of code. Explanation:

Ｎθ

Input d.

≔Ｘ²¦³²η

Precompute 2³² as it gets used several times.

≔⌕ＥＸ²…³²¦⁶⁴÷×⁻⊖η﹪ηθ﹪±ιθι⁰ζ

Find the first power of 2 between 2³² and 2⁶⁴ which satisfies the formula in A346496, even if d is a power of 2. Note that the formula given there calculates the first b where 2**b<=(2**32-(2**32%d)-1)*(d-1-(2**b-1)%d), but I actually calculate (2**32-(2**32%d)-1)*(-2**b%d)//(2**b) and find the first b where that is 0. For now, assume f=0 and s=b-32.

≔÷×ηＸ²ζθε

Assume m=A346495(d) (but using the value of b computed above even if d is a power of 2). This value is correct if it lies within the desired range [2..2³²].

≔›εηδ

Compute f=m>2³².

Ｉ⟦⁻ε×δ⊖η⁻ζδδ

If f=1, subtract 2³²-1 from m and 1 from s. (The OEIS algorithm for m>2³² seems to erroneously hard-code the value of s-1 as 2, which happens to work for its example of d=7.) Output the final values of m, s and f.