シフト演算と割り算(python3)

# はじめに

前回の<a href="https://qiita.com/b1ueskydragon/items/0b8e0c382d782423c6d3">高速累乗計算</a>に続いて、シンプルなアルゴリズムに関するものです.
python3 で実装します.

# ゴール
割り算(/)・掛け算(*)・剰余(%)の演算子を使わずに二つの正の整数の割り算を実装する.
ただし、整数の商のみ求め、余りは捨てる.

求まった商が 

```.py
# 切り捨て除算
a // b 
```

と同じ結果になればいいかと思います.

# まず残された選択肢を検討する
単純に、四則演算をもとに考えてみましょう.
「割り算・掛け算・剰余の演算子を使わずに」ということで、
「足し算」と「引き算」が残っていることがわかります.

## たとえば引き算
割り算なので被除数をどんどん小さくしていく必要があります. そういった点を考えると、引き算は都合がいいかもしれません.

以下は、除数で被除数をどんどん引いていくことを実装しています.


```.py
def division_brute_force(a, b):
    """
    a // b

    :param a: dividend
    :param b: divisor
    :return: quotient
    """
    n = 0
    while a >= b:
        a -= b
        n += 1
    return n

```

ただし、このような力任せ方法だと、データが大きくなればなるほど、計算が遅くなってしまいます.
一例で、

```.math
2^{50} \div  1
```

このように被除数と除数の差が大きすぎるとき、
除数1は被除数2^50を2^50回引くことになります.

試しで値を入れてみましたが、私の環境だと、長々と返ってこないです...😉

そこで、他の演算子をも検討する必要性が出てきます.
ビット演算子です. なかでも、`シフト演算子`を利用します.

まずはシフト演算について、次にどう利用できるかについて考えてみます.

# シフト演算
二進数のビット(桁)、つまり0と1の列をずらします.
論理シフト(符号なし)と算術シフト(符号あり)でさらにわけられます.


```.math
nビット左にずらす = 2^{n} 倍
\\
nビット右にずらす = 2^{-n} 倍
```

## 論理シフト(符号なし)
### 例
10進数の`8`は2進数で表すと`1000`です.

####`1000`を3bit論理左シフト
3桁左へずらして`1000_ _ _`空いたところには`0`を埋めます.
`1000000` (10進数の`64`)

#### `1000`を3bit論理右シフト
右から数えて3桁捨てます. ぽん.
`1` (10進数の`1`)

## 算術シフト(符号あり)

算術シフトの場合は、一番左のビットが「符号」となり、`0`なら正の数、`1`なら負の数を意味します.
シフトを行う際は、符号ビットは位置を固定したまま、符号より右にあるビット列を左右にずらして計算します.

### 例
10進数の`-16`を8bitの補数表現の2進数で表すと`11110000`です. 
これは9bitの`100000000`から`10000`を引いた結果です.

#### `11110000`を2bit算術左シフト
一番左の符号の`1`を固定して、残りの`1110000`を左に2bitずらすと、`111000000`になりますが、
このままだと桁があふれるので、左側から1番目と2番目を捨てて`1000000`にします.
あとは固定しておいた符号`1`とくっつけると`11000000`が求まります(10進数の`-64`)

#### `11110000`を2bit算術右シフト
一番左の符号の`1`を固定して、残りの`1110000`を右に2bitずらすと、`_ _ 1110000`になりますが、
このままだと桁があふれるので、右側から1番目と2番目を捨てて`_ _ 11100`にします.
空いた桁には符号と同じ数字を埋めるので、`1`を埋めると`1111100`になります.
あとは固定しておいた符号`1`とくっつけると`11111100`が求まります(10進数の`-4`)

ただ、今回は正の整数のみを限定して扱うので、論理シフトをベースに考えていきます.

## 演算子

```.py
# shift Left
8 << 3 # 64

# shift Right
8 >> 3 # 1
```


# もう一度割り算

たとえばこれ.

```.math
100 \div  3
```

商は33、余りは1です.
余りは捨てるので、求める結果は33になります.

割り算を紙に書いて計算する時、こういう風に計算すると思います.

タイピングでは、うまく表現できず......

![00.png](https://qiita-image-store.s3.amazonaws.com/0/194231/f9720511-830e-3955-33f2-b43158989ccc.png)

割り算は、被除数を除数で引いていくことで実装できる、という点を考慮すると、

被除数を

```.math
除数 \times k \times 10^{n}
```

で引いて小さくしていけばいいことがわかります.
10進数なので`k`は1から9の間の自然数になるでしょうし、自然数`n`は桁になります.
割った余りが次の被除数となり、上記の式が被除数よりは大きくならない、ぎりぎりの大きさになるように`n`を決めていきます.

この過程を、余り(次の被除数)が除数より小さくなるまで繰り返すことになります.

![02ss.png](https://qiita-image-store.s3.amazonaws.com/0/194231/a0d3cd9e-13e6-7bad-d3c1-62bbd7b51e1f.png)

途中から出てきた商を足し合わせると、最終的な商が求まります.

普通の割り算です.


## もう一度割り算 (二進数)

まったく同じ式を２進数に変換します.

```.math
1100100(2) \div  11(2)
```

![03s.png](https://qiita-image-store.s3.amazonaws.com/0/194231/11d83814-79cb-9a4e-894f-5de6770d427b.png)

被除数を

```.math
除数 \times 2^{n}
```

で引いて小さくしていけばいいことがわかります. 
2進数なので、先ほどの`k`は1となります.

下記のことを考慮すると、掛け算の演算子の代わりにシフト演算子を使うことでこの割り算が実装できることがわかります.

```.math
nビット左にずらす = 2^{n} 倍
```

余り(次の被除数)が除数より小さくなるまで繰り返します.

# シフト演算で割り算

ゼロ除算例外処理と、除数が被除数より大きい場合の処理も含んでいます.

### 最初に扱うビット数を決めて、除数のビットを削っていく実装

```.py
def division_binary_shift_m(a, b, m):
    """
    :param a: dividend
    :param b: divisor
    :param m: maximum bit capacity
    :return: quotient
    """
    if a < b:
        return 0
    if b == 0:
        raise ZeroDivisionError

    q = 0
    sb = b << m
    rem = a  # remainder (init is dividend)
    while rem >= b:
        while sb > rem:
            sb = sb >> 1
            m = m - 1

        q = q + (1 << m)
        rem = rem - sb

    return q
```

最初に扱うビット数`m`を決めて、除数のビット数を一旦最大ビットまで増やして(左シフトして)おきます.
余り`rem`の初期値は、まだ何も引いていないので、被除数になります.

被除数で割れるようになるまで、一つずつ除数のビットを削って(右にシフトして)いきます.
除数と共にビット数`m`も一つずつ減らしていきます.

内側の while loop により「除数より大なりイコールな被除数」な状態(割り算できる状態)になったら、桁数`m`が決まります.

商`q`には`2^m`を足します. 2のm乗を`1 << m`で表現しています.
余りからは除数を引いて、次の被除数にします.

> シフト演算子より加法(+ と -)の方が優先されるので、括弧で囲みます.
>
```.py
q = q + (1 << m)
```

> 複合代入(+= など)はシフト演算子より順位が低いのでまとめて書くと括弧要らず.
>
```.py
 q += 1 << m
```

### 最大桁を決めずに、除数のビットを増やしていく実装

```.py
def division_binary_shift(a, b):
    """
    a // b

    :param a: dividend
    :param b: divisor
    :return: quotient
    """
    if a < b:
        return 0
    if b == 0:
        raise ZeroDivisionError

    q = 0
    sq = 1
    ori_b = b
    while a >= b:
        b <<= 1
        if a >= b:
            sq <<= 1
        else:
            b >>= 1
            a -= b
            q += sq
            sq = 1
            b = ori_b

    return q

```

被除数のサイズのぎりぎりまで除数(`b`)を左シフトしていきます. 
各段階における商(`sq`)も、同じく桁を増やしていきます.

除数が被除数より大きくなったら除数を一つ少ない桁に戻して、被除数から除数を引きます.

各段階から出てきた商は、最終的な商(`q`)に足していきます.
各段階における一連の処理が終わったら(被除数から引いて、最終的な商に貯める)商を1に、除数をもともとの除数に戻します.



# 実行例

```.py
a, b = 14285699999999857143, 7
m = 64
expected = a // b

assert division_binary_shift_m(a, b, m) == expected
assert division_binary_shift(a, b) == expected

# 商は 2040814285714265306
```

### メモリ使用量
せっかく(?)だから、並べて比較してみました.
一応、一番最初の brute_force 版もかけてみましたけど、案の定まったく返ってこない 😉

```
$ python -m memory_profiler my_division.py

Filename: my_division.py

Line #    Mem usage    Increment   Line Contents
================================================
     4   26.395 MiB   26.395 MiB   @profile
     5                             def division_binary_shift(a, b):
     6   26.395 MiB    0.000 MiB       if a < b:
     7                                     return 0
     8   26.395 MiB    0.000 MiB       if b == 0:
     9                                     raise ZeroDivisionError
    10
    11   26.395 MiB    0.000 MiB       q = 0
    12   26.395 MiB    0.000 MiB       sq = 1
    13   26.395 MiB    0.000 MiB       ori_b = b
    14   26.398 MiB    0.000 MiB       while a >= b:
    15   26.398 MiB    0.000 MiB           b <<= 1
    16   26.398 MiB    0.004 MiB           if a >= b:
    17   26.398 MiB    0.000 MiB               sq <<= 1
    18                                     else:
    19   26.398 MiB    0.000 MiB               b >>= 1
    20   26.398 MiB    0.000 MiB               a -= b
    21   26.398 MiB    0.000 MiB               q += sq
    22   26.398 MiB    0.000 MiB               sq = 1
    23   26.398 MiB    0.000 MiB               b = ori_b
    24
    25   26.398 MiB    0.000 MiB       return q

Line #    Mem usage    Increment   Line Contents
================================================
    28   26.402 MiB   26.402 MiB   @profile
    29                             def division_binary_shift_m(a, b, m):
    30   26.402 MiB    0.000 MiB       if a < b:
    31                                     return 0
    32   26.402 MiB    0.000 MiB       if b == 0:
    33                                     raise ZeroDivisionError
    34
    35   26.402 MiB    0.000 MiB       q = 0
    36   26.402 MiB    0.000 MiB       sb = b << m
    37   26.402 MiB    0.000 MiB       rem = a
    38   26.402 MiB    0.000 MiB       while rem >= b:
    39   26.402 MiB    0.000 MiB           while sb > rem:
    40   26.402 MiB    0.000 MiB               sb = sb >> 1
    41   26.402 MiB    0.000 MiB               m = m - 1
    42
    43   26.402 MiB    0.000 MiB           q = q + (1 << m)
    44   26.402 MiB    0.000 MiB           rem = rem - sb
    45
    46   26.402 MiB    0.000 MiB       return q
```


# おわりまして
コンピューターが行う掛け算と割り算は、ビットを左右にずらすことがベースになっていると思いますが、改めてそういった演算や実装に関して考えるきっかけになって、興味深かったです.

ご指摘や補足などありましたら、随時教えてくださいませ.