diff --git a/python/_long.py b/python/_long.py
index e2ff938f..0e5aa3dc 100644
--- a/python/_long.py
+++ b/python/_long.py
@@ -99,11 +99,38 @@ def ulong_divmodi(a: list, b: int):
     ulong_unpad_(res)
     return res, carry
 
+
 def ulong_divmod(a: list, b: list):
-    q = [0]
-    while ulong_cmp(a, b) >= 0:
-        ulong_inc_(q)
-        a = ulong_sub(a, b)
+
+    if ulong_cmp(a, b) < 0:
+        return [0], a
+
+    if len(b) == 1:
+        q, r = ulong_divmodi(a, b[0])
+        r, _ = ulong_fromint(r)
+        return q, r
+
+    max = (len(a) - len(b)) * PyLong_SHIFT + \
+        (a[-1].bit_length() - b[-1].bit_length())
+
+    low = [0]
+
+    high = (max // PyLong_SHIFT) * [0] + \
+        [(2**(max % PyLong_SHIFT)) & PyLong_MASK]
+
+    while ulong_cmp(low, high) < 0:
+        ulong_inc_(high)
+        mid, r = ulong_divmodi(ulong_add(low, high), 2)
+        if ulong_cmp(a, ulong_mul(b, mid)) >= 0:
+            low = mid
+        else:
+            high = ulong_sub(mid, [1])
+
+    q = [0] * (len(a) - len(b) + 1)
+    while ulong_cmp(a, ulong_mul(b, low)) >= 0:
+        q = ulong_add(q, low)
+        a = ulong_sub(a, ulong_mul(b, low))
+    ulong_unpad_(q)
     return q, a
 
 def ulong_floordivi(a: list, b: int):
diff --git a/tests/09_long.py b/tests/09_long.py
index 412a2d38..9026bdf9 100644
--- a/tests/09_long.py
+++ b/tests/09_long.py
@@ -15,7 +15,6 @@ assert 1 + a == 3L
 assert 1 - a == -1L
 assert 2 * a == 4L
 
-
 # __lshift__ and __rshift__
 for i in range(29):
     assert 1L << i == 2 ** i
@@ -35,4 +34,17 @@ assert s == [4, 6, 7, 2, 3]
 
 assert 1 < 2L < 3 < 6.6
 assert 1L < 2 < 9.6 >= 7 > 2L
-assert 1L < 2 < 3 < 6.6
\ No newline at end of file
+assert 1L < 2 < 3 < 6.6
+
+assert 10000000000000000000000L // 3333L == 3000300030003000300L
+assert 10000000000000000000000L % 3333L == 100L
+assert 2L ** 100 // 3L ** 50 == 1765780L
+assert 2L ** 200 // 3L ** 100 == 3117982410207L
+assert 2L ** 500 // 3L ** 200 == 12323863745843010927046405923587284941366070573310012484L
+assert 2L ** 500 % 3L ** 200 == 242990057207501525999897657892105676264485903550870122812212566096970021710762636168532352280892L
+
+assert 3L**500 // 3L**400 == 3L**100
+assert 4562645248L // 3L == 1520881749L
+assert 3L**10 // 2L**200 == 0
+assert 2L**500 % 3L**200 == 242990057207501525999897657892105676264485903550870122812212566096970021710762636168532352280892L
+assert divmod(10L**115, 3L**47) == (376098003656605353510839433041026531338835641081336270795298342081499446459650747504311257564L, 10150482013473427429132L)