Divide by an odd number, 2-adically

Question

Given \$a\$ and \$b\$, both odd \$n+1\$-bit integers, compute \$a/b\$ to a precision of \$n+1\$ bits in the 2-adic integers. That is, compute \$c\$ such that \$a = bc\, (\mathop{\rm mod} 2^{n+1})\$. \$n\$ should be your language's native integer size, or if native integers are bigints, take it as a parameter. If your language uses trits (and presumably is either Setun assembly, TriINTERCAL, or Malbolge), you may instead compute in the 3-adics, in which case \$a\$ and \$b\$ should be multiples of 3 plus 1.

Inputs should be \$(a-1)/2\$ and \$(b-1)/2\$ (trits: \$(x-1)/3\$).

This is code-golf, so shortest answer in bytes (per language) wins.

Test cases:

All test cases are truncatable; if the last \$n\$ bits match the inputs, the last \$n\$ bits of the outputs match.
Test cases (in hex, 32bit): (apologies for poor vinculum placement)

|  (a-1)/2  |  (b-1)/2  | (a/b-1)/2 |
|-----------+-----------+-----------|
| …00000000 | …00000001 | …55555555 | (1/3 = A̅B)
| …00000000 | …00000002 | …66666666 | (1/5 = C̅D)
| …00000001 | …00000002 | …33333333 | (3/5 = 6̅7)
| …00000000 | …00000003 | …DB6DB6DB | (1/7 = 6̅D̅B̅7)
| …FFFFFFFF | …00000003 | …24924924 | (-1/7 = 2̅4̅9̅)
| …4620BA27 | …1876DCBC | …48CAF903 | (random)

More test cases may be generated by multiplying random \$n+1\$-bit odd integers and taking the last \$n+1\$ bits of the result (then shifting right by 1).

A few test cases for ternary computers (nonary this time):

|  (a-1)/3  |  (b-1)/3  | (a/b-1)/3 |
|-----------+-----------+-----------|
| …00000000 | …00000002 | …51251251 | (1/7 = 3̅7̅6̅4)
| …23472148 | …12435871 | …65732854 | (random again)

Similarly, do the same with \$n+1\$-trit integers ending with a 1 trit.

Answer 1

Python, 45 bytes

lambda a,b,n:(a-~a)*pow(b-~b,-1,2<<n)//2%2**n

Attempt This Online!

Beginning with Python 3.8 the builtin pow can directly compute the modular inverse.

Answer 2

PARI/GP, 32 bytes

f(a,b,n)=(2*a+1)/(2*b+1)%2^n++\2

Attempt This Online!

PARI/GP has a built-in p-adic type, but using the % operator is shorter.

Answer 3

C (gcc), 48 bytes

c;i;f(a,b){for(i=32;i;--i)c=a-b-2*b*c;return c;}

Try it online!

Explanation: since \$(2b+1) (2c+1) = 2(2bc+b+c)+1\$ should end up equal to \$2a+1\$, we want to find a root of \$f(c) = 2bc+b+c-a\$. We do so using Newton's method (as in the proof of Hensel's lemma), with the optimization that \$f'(c) = 2b+1 \equiv 1\pmod{2}\$, so it suffices to divide by 1 instead of by \$f'(c)\$. So, we iterate to find a fixed point of \$c \mapsto c - \frac{f(c)}{1} = a-b-2bc\$. (And since each iteration increases the number of correct binary digits by at least one, and we are working with 32-bit integers, then 32 iterations suffices.)

Note that it's interesting we do not need to initialize c to any particular starting value, since whatever the starting value is, we converge to the unique fixed point anyway.

Answer 4

C (gcc), 41 bytes

f(a,b){return b?f(2*a*b+a+b,b*b+b<<1):a;}

Try it online!

This uses the same method as my previous answer: apply \$a/b = ab/b^2\$ repeatedly until \$b=1\$, and then \$a/1=a\$. Because \$(\mathbb{Z}/2^{33}\mathbb{Z})^{\times}\$ is a group of order \$2^{32}\$, this finishes after at most 32 iterations.

Answer 5

ARM T32 machine code, 16 bytes

004a 3201 4350 1840 3201 4351 d1f8 4770

In assembly:

f:  lsls r2, r1, #1 @ Set r2 to b-1 ((b-1)/2 is in r1).
    adds r2, #1     @ r2 becomes b.
    muls r0, r2     @ r0 (which was (a-1)/2) becomes b(a-1)/2.
    adds r0, r1     @ r0 becomes b(a-1)/2 + (b-1)/2 = (ab-1)/2.
    adds r2, #1     @ r2 becomes b+1.
    muls r1, r2     @ r1 becomes (b+1)(b-1)/2 = (b²-1)/2.
    bne.n f         @ Jump back if that is nonzero.
    bx lr           @ Return.

This applies \$a/b = ab/b^2\$ repeatedly until \$b=1\$, and then \$a/1=a\$.
Because \$(\mathbb{Z}/2^{33}\mathbb{Z})^{\times}\$ is a group of order \$2^{32}\$, this finishes after at most 32 iterations.

Answer 6

R, 36 bytes

\(a,b){while((b*F)%%2^32-a)F=F+1;+F}

Attempt This Online! or try the first few test-cases with n set to 15

Input a and b, outputs c = a/b(modulo 2^32).

Uses floating-point calculations to try each multiple of b until it finds the first one that equals a mod 2^32 (corresponding to R's integer size). Note that this means that some inputs may be subject to floating-point rounding errors when c*b exceeds 2^53.

---

R, 40 bytes

\(a,b){while((F=(F+b)%%2^32)-a)T=T+1;+T}

Attempt This Online! or try the first 2 test-cases with n set to 15

This also uses floating-point calculations to successively increment by b until the total mod 2^32 equals a, but here they are limited to a maximum of 2^33, which is easily within the range that is representable without loss-of-precision, so output for all inputs is exact.
It's much slower, though.

Answer 7

Husk, 11 bytes

€Mo%`^2→⁰*N

Try it online! or try the first test-case with n set to 15

Input is n, b, a; output is c = a/b (modulo 2^(n+1)).

One byte less (remove the →) it it's acceptable to input the number of bits of the modulus (n+1) directly. 4 bytes less if it's acceptable to input the modulus itself directly.

Answer 8

Z80, 14 bytes (n=8)

Input is \$(a-1)/2\$ in a, \$(b-1)/2\$ in c; output is in d. Trashes a, b, f. (Well, everything trashes f.)

91 06 08 CB 38 1F 30 03
91 CB FA 10 F6 C9

Disassembled:

f:
    sub c       // 91       r = n - d
    ld b,8      // 06 08    do 8 times:
l:
    srl d       // CB 38    q >>= 1
    rra         // 1F       r >>= 1
    jr nc,ns    // 30 03    Was r even?
    sub c       // 91       N: r -= d
    set 7,d     // CB FA    N: q |= 1 << 7
ns:
    djnz l      // 10 F6    in any case: jump back first 7 times.
    ret         // C9       return

Answer 9

Charcoal, 29 bytes

ＮθＮηＦＮ«⊞υ﹪θ²≧⁻×η﹪θ²θ≦⊘θ»Ｉ↨⮌υ²

Try it online! Link is to verbose version of code. Takes a, b and n+1 as input and outputs c=a/b. Explanation: This is like long division in reverse; to get the i+1th bit you compute c to i bits and see whether you want to add a 0 or 1 bit to c. Rather than calculating b*c each time, a is reduced each time a 1 bit is added to c, and also a power of 2 is eliminated by dividing a by 2 each step.

ＮθＮηＦＮ«

Input a, b and n+1 and repeat that many times.

⊞υ﹪θ²

Save the LSB of a into c.

≧⁻×η﹪θ²θ

If a was odd then subtract b from it.

≦⊘θ

Divide a by 2.

»Ｉ↨⮌υ²

Convert c from binary.

Answer 10

MMIX, 44 bytes (11 instructions) (n = 64)

Call sequence: get \$a\$ in $(X+1), \$b\$ in $(X+2), PUSHJ $X,div2ad, quotient is now in $X, and as a bonus $(X+2) is unchanged and $(X+1) contains a sort of remainder that can be used again starting an instruction later to get more bits.

Hexdump (jxd -T):

00000000: 26000001 e3020000 e3030001 76ff0003  &¡¡¢ẉ£¡¡ẉ¤¡¢v”¡¤
00000010: 3b030301 c00202ff 76ff0001 3f000001  ;¤¤¢Ċ££”v”¡¢?¡¡¢
00000020: 260000ff 5b03fffa f8030000           &¡¡”[¤”«ẏ¤¡¡

Dissassembled:

div2ad  SUBU $0,$0,$1   // r = a - b
        SET  $2,0       // q = 0
        SET  $3,1       // holds 2^i
0H      ZSOD $255,$0,$3 // loop: tmp = (r&1)? 2^i : 0
        SLU  $3,$3,1    // increment i
        OR   $2,$2,$255 // q |= tmp
        ZSOD $255,$0,$2 // tmp = (r&1)? b : 0
        SRU  $0,$0,1    // r >>= 1
        SUBU $0,$0,$255 // r -= tmp
        PBNZ $3,0B      // if(i <= 64) goto loop (that 1 bit got shifted off the high end)
        POP  3,0        // return q,r,b in that order