diff --git a/std/math/big/bits.jule b/std/math/big/bits.jule
index a0107a531..f0597f116 100644
--- a/std/math/big/bits.jule
+++ b/std/math/big/bits.jule
@@ -327,11 +327,14 @@ fn karatsuba(mut x: bits, mut y: bits): bits {
     ret karatsuba_add(karatsuba_add(c1, c2), p2)
 }
 
+// Recursion division algorithm. It will update left operand if necessary.
+// Uses bit shifting strategy.
+// Returns quotient.
 fn recursive_div(mut &x: bits, mut &y: bits): bits {
-    if x.len < y.len {
+    match cmp(x, y) {
+    | -1:
         ret nil
-    }
-    if cmp(x, y) == 0 {
+    | 0:
         ret [1]
     }
     let mut yq = make(bits, y.len, x.len)
@@ -363,3 +366,27 @@ fn recursive_div(mut &x: bits, mut &y: bits): bits {
     add_res(yq, k)
     ret yq
 }
+
+// Recursion modulo algorithm. It will update left operand if necessary.
+// Uses bit shifting strategy.
+// Returns remainder.
+fn recursive_mod(mut &x: bits, &y: bits): bits {
+    match cmp(x, y) {
+    | 0:
+        ret nil
+    | -1:
+        ret x
+    }
+    let mut yq = make(bits, y.len, x.len)
+    _ = copy(yq, y)
+    for cmp(yq, x) == -1 {
+        yq = append(yq[:1], yq...)
+        yq[0] = 0b0
+    }
+    if yq.len == y.len {
+        ret [1]
+    }
+    sub_res(x, yq[1:])
+    fit(x)
+    ret recursive_mod(x, y)
+}
diff --git a/std/math/big/int.jule b/std/math/big/int.jule
index dc00a32a1..a9eeb8449 100644
--- a/std/math/big/int.jule
+++ b/std/math/big/int.jule
@@ -157,6 +157,29 @@ impl Int {
         }
     }
 
+    // Modulo two Int and returns result.
+    pub fn mod(self, y: Int): Int {
+        let mut r = self
+        r %= y
+        ret r
+    }
+
+    // Modulo Int.
+    pub fn mod_assign(mut self, y: Int) {
+        if self.neg == y.neg {
+            self.nat %= y.nat
+            self.neg = self.len() > 0 && y.neg
+            ret
+        }
+        self.nat %= y.nat
+        if self.len() == 0 {
+            self.neg = false
+            ret
+        }
+        self.nat -= y.nat
+        self.neg = self.len() > 0 && y.neg
+    }
+
     // Bitwise left shift.
     pub fn shl(self, y: uint): Int {
         let mut r = self
diff --git a/std/math/big/int_test.jule b/std/math/big/int_test.jule
index 5d14b2595..8278c32ad 100644
--- a/std/math/big/int_test.jule
+++ b/std/math/big/int_test.jule
@@ -39,7 +39,7 @@ static cases_int_mul = [
 static cases_int_div = [
     ["1010", "+", "1", "-", "1010", "-"],
     ["1010", "+", "10", "+", "101", "+"],
-    ["1000000", "-", "110", "-", "1011", "+"],
+    ["1000000", "-", "110", "-", "1010", "+"],
     ["1001100", "+", "11", "+", "11001", "+"],
     ["1000000", "+", "100", "+", "10000", "+"],
     ["11111011", "-", "1101", "+", "10011", "-"],
@@ -49,6 +49,22 @@ static cases_int_div = [
     ["100111001110111001111001111111010101001001011000000000110101110110111110110000101101001", "+", "10", "+", "10011100111011100111100111111101010100100101100000000011010111011011111011000010110100", "+"],
 ]
 
+static cases_int_mod = [
+    ["1010", "-", "11", "+", "10", "+"],
+    ["1011110110110101010101010101010011001111010101101011100", "-", "1101101", "+", "111011", "+"],
+    ["1010", "+", "11", "-", "10", "-"],
+    ["1010", "+", "11", "+", "1", "+"],
+    ["1011110110110101010101010101010011001111010101101011100", "+", "1101101", "-", "111011", "-"],
+    ["1010", "+", "10000000000", "+", "1010", "+"],
+    ["10111101", "+", "10111111", "+", "10111101", "+"],
+    ["10001110000010", "+", "11000011010100000", "-", "10110001100011110", "-"],
+    ["1110101011011101", "+", "1000000000000000000000", "-", "111110001010100100011", "-"],
+    ["1110101011011101", "-", "1000000000000000000000", "-", "1110101011011101", "-"],
+    ["1000000000", "-", "1101", "-", "101", "-"],
+    ["1010", "-", "10000000000", "+", "1111110110", "+"],
+    ["1010", "+", "10000000000", "-", "1111110110", "-"],
+]
+
 #test
 fn test_int_add(mut t: &T) {
     for _, c in cases_int_add {
@@ -169,6 +185,36 @@ fn test_int_div(mut t: &T) {
     }
 }
 
+#test
+fn test_int_mod(mut t: &T) {
+    for _, c in cases_int_mod {
+        let n1 = Int.from_bits(c[0], c[1] == "-") else {
+            t.errorf("exception occurs: {}", error)
+            continue
+        }
+        let n2 = Int.from_bits(c[2], c[3] == "-") else {
+            t.errorf("exception occurs: {}", error)
+            continue
+        }
+        let mut r = n1.mod(n2)
+        let cr = c[4]
+        if r.len() != cr.len ||
+            r.neg && c[5] != "-" ||
+            !r.neg && c[5] != "+" {
+            t.errorf("{}{} % {}{} != {}{}", c[1], c[0], c[3], c[2], c[5], c[4])
+            continue
+        }
+        for i, b in r.nat.bits {
+            let cb = cr[cr.len - 1 - i]
+            if b == 0b1 && cb != '1' ||
+                b == 0b0 && cb != '0' {
+                t.errorf("{}{} % {}{} != {}{}", c[1], c[0], c[3], c[2], c[5], c[4])
+                break
+            }
+        }
+    }
+}
+
 #test
 fn test_int_lt(mut t: &T) {
     t.assert(!Int.from_bits("1011010", false)!.lt(Int.from_bits("00001011010", false)!), "1) 1011010 < 00001011010")
diff --git a/std/math/big/nat.jule b/std/math/big/nat.jule
index 9a991738b..60e8f005c 100644
--- a/std/math/big/nat.jule
+++ b/std/math/big/nat.jule
@@ -159,6 +159,12 @@ impl Nat {
         | 1:
             ret
         | 0:
+            // Right operand is zero, result is always zero.
+            self.bits = nil
+            ret
+        }
+        if self.len() == 0 {
+            // Right operand is zero, result is always zero.
             self.bits = nil
             ret
         }
@@ -171,6 +177,7 @@ impl Nat {
             ret
         }
 
+        // Make size of operands are equal if not.
         let mut xb = self.bits
         let mut yb = y.bits
         let n = max(xb.len, yb.len)
@@ -213,23 +220,80 @@ impl Nat {
 
     // Divides Nat.
     pub fn div_assign(mut self, y: Nat) {
+        if self.len() == 0 {
+            // Left operans is zero, remainder is always zero.
+            ret
+        }
+        match y.len() {
+        | 1:
+            // Right operand is 1, quotient is always equals to left operand.
+            ret
+        | 0:
+            panic("std::math::big: division by zero")
+        }
+        match self.cmp(y) {
+        | -1:
+            // Left operand is less than right oprand.
+            // Quotient is always equals to zero.
+            self.bits = nil
+            ret
+        | 0:
+            // Left oprand and right operand are equal.
+            // Quotient is always zero.
+            self.bits = [1]
+            ret
+        }
+        // Use clone because of recursive division can change left operand.
+        let mut xb = clone(self.bits)
+        self.bits = recursive_div(xb, y.bits)
+        self.fit()
+    }
+
+    // Modulo two Nat and returns result.
+    pub fn mod(self, y: Nat): Nat {
+        let mut r = self
+        r %= y
+        ret r
+    }
+
+    // Modulo Nat.
+    pub fn mod_assign(mut self, y: Nat) {
+        if self.len() == 0 {
+            // Left operans is zero, remainder is always zero.
+            ret
+        }
         match {
         | y.len() == 1:
+            // Right operand is 1, remainder is always zero.
+            self.bits = nil
+            ret
+        | y.len() == 2 && y.bits[0] == 0b0:
+            // Right operand is 2.
+            // If left opeand is even, remainder is always zero.
+            // If left operand is odd, remainder is always one.
+            if self.even() {
+               self.bits = nil
+            } else {
+               self.bits = [1]
+            }
             ret
         | y.len() == 0:
             panic("std::math::big: division by zero")
         }
         match self.cmp(y) {
         | -1:
-            self.bits = nil
+            // Left operand less than right operand.
+            // Remainder always equals to left operand.
             ret
         | 0:
-            self = Nat.one()
+            // Left oprand and right operand are equal.
+            // Remainder is always zero.
+            self.bits = nil
             ret
         }
         // Use clone because of recursive division can change left operand.
         let mut xb = clone(self.bits)
-        self.bits = recursive_div(xb, y.bits)
+        self.bits = recursive_mod(xb, y.bits)
         self.fit()
     }
 
diff --git a/std/math/big/nat_test.jule b/std/math/big/nat_test.jule
index 461a6a464..bfe59453f 100644
--- a/std/math/big/nat_test.jule
+++ b/std/math/big/nat_test.jule
@@ -31,6 +31,7 @@ static cases_nat_sub = [
     ["10001100001110", "111101101", "10000100100001"],
     ["1010", "10100", "1010"],
     ["101011", "11000", "10011"],
+    ["11101110111011111111111011110101010010101101111010101101010101010101011110011010010110010101110110010101101011001011", "11101110111011111111111011110101010010101101111010101101010101010101011110011010010110010101110110010101101011001010", "1"],
 ]
 
 static cases_nat_mul = [
@@ -51,7 +52,7 @@ static cases_nat_mul = [
 static cases_nat_div = [
     ["1010", "1", "1010"],
     ["1010", "10", "101"],
-    ["1000000", "110", "1011"],
+    ["1000000", "110", "1010"],
     ["1001100", "11", "11001"],
     ["1000000", "100", "10000"],
     ["11111011", "1101", "10011"],
@@ -60,6 +61,20 @@ static cases_nat_div = [
     ["10011001101010", "10", "1001100110101"],
     ["1010101010101010101010101010101010101010101010", "10", "101010101010101010101010101010101010101010101"],
     ["100111001110111001111001111111010101001001011000000000110101110110111110110000101101001", "10", "10011100111011100111100111111101010100100101100000000011010111011011111011000010110100"],
+    ["110101101011110101010100111", "110100", "1000010000100101101111"],
+    ["110101101011110101010100111", "11011001", "1111110101010101010"],
+]
+
+static cases_nat_mod = [
+    ["1010", "1010", ""],
+    ["1010", "1", ""],
+    ["111001010111011001", "1110110", "11111"],
+    ["11101110111011111111111011110101010010101101111010101101010101010101011110011010010110010101110110010101101011001011", "1110101", "1011001"],
+    ["11110110110101010101010101111111111011010101010101010101001111111001010101001010101010110000000011111101101001010111001111", "1110101", "1100001"],
+    ["0", "11100111", ""],
+    ["0", "0", ""],
+    ["110101101011110101010100111", "110100", "11011"],
+    ["1010", "11", "1"],
 ]
 
 static cases_nat_shl: [][]u64 = [
@@ -215,6 +230,34 @@ fn test_nat_div(mut t: &T) {
     }
 }
 
+#test
+fn test_nat_mod(mut t: &T) {
+    for _, c in cases_nat_mod {
+        let n1 = Nat.from_bits(c[0]) else {
+            t.errorf("exception occurs: {}", error)
+            continue
+        }
+        let n2 = Nat.from_bits(c[1]) else {
+            t.errorf("exception occurs: {}", error)
+            continue
+        }
+        let mut r = n1.mod(n2)
+        let cr = c[2]
+        if r.len() != cr.len {
+            t.errorf("{} % {} != {}", c[0], c[1], cr)
+            continue
+        }
+        for i, b in r.bits {
+            let cb = cr[cr.len - 1 - i]
+            if b == 0b1 && cb != '1' ||
+                b == 0b0 && cb != '0' {
+                t.errorf("{} % {} != {}", c[0], c[1], cr)
+                break
+            }
+        }
+    }
+}
+
 #test
 fn test_nat_shl(mut t: &T) {
     for _, c in cases_nat_shl {