diff --git a/README.md b/README.md
index 2b60397..f91ccc7 100644
--- a/README.md
+++ b/README.md
@@ -86,7 +86,7 @@ fn main() raises:
     print(a * b)  # 152415787654.32099750190521
     print(a.true_divide(b, precision=80))  # 100000.0001
 
-    # === Mathematical Functions === #
+    # === Exponential Functions === #
     print(a.sqrt(precision=80))
     # 11111.111066111110969430554981749302328338130654689094538188579359566416821203641
     print(a.cbrt(precision=80))
@@ -102,6 +102,24 @@ fn main() raises:
     print(a.ln(precision=80))
     # 18.631401767168018032693933348296537542797015174553735308351756611901741276655161
 
+    # === Trigonometric Functions === #
+    print(a.sin(precision=200))
+    # 0.99985093087193092464780008002600992896256609588456
+    #   91036188395766389946401881352599352354527727927177
+    #   79589259132243649550891532070326452232864052771477
+    #   31418817041042336608522984511928095747763538486886
+    print(b.cos(precision=1000))
+    # -0.9969577603867772005841841569997528013669868536239849713029893885930748434064450375775817720425329394
+    #    9756020177557431933434791661179643984869397089102223199519409695771607230176923201147218218258755323
+    #    7563476302904118661729889931783126826250691820526961290122532541861737355873869924820906724540889765
+    #    5940445990824482174517106016800118438405307801022739336016834311018727787337447844118359555063575166
+    #    5092352912854884589824773945355279792977596081915868398143592738704592059567683083454055626123436523
+    #    6998108941189617922049864138929932713499431655377552668020889456390832876383147018828166124313166286
+    #    6004871998201597316078894718748251490628361253685772937806895692619597915005978762245497623003811386
+    #    0913693867838452088431084666963414694032898497700907783878500297536425463212578556546527017688874265
+    #    0785862902484462361413598747384083001036443681873292719322642381945064144026145428927304407689433744
+    #    5821277763016669042385158254006302666602333649775547203560187716156055524418512492782302125286330865
+
     # === Internal representation of the number === #
     (
         BDec(
diff --git a/docs/examples_on_bdec.mojo b/docs/examples_on_bdec.mojo
index e1b4a06..c3c2417 100644
--- a/docs/examples_on_bdec.mojo
+++ b/docs/examples_on_bdec.mojo
@@ -11,7 +11,7 @@ fn main() raises:
     print(a * b)  # 152415787654.32099750190521
     print(a.true_divide(b, precision=80))  # 100000.0001
 
-    # === Mathematical Functions === #
+    # === Exponential Functions === #
     print(a.sqrt(precision=80))
     # 11111.111066111110969430554981749302328338130654689094538188579359566416821203641
     print(a.cbrt(precision=80))
@@ -27,6 +27,24 @@ fn main() raises:
     print(a.ln(precision=80))
     # 18.631401767168018032693933348296537542797015174553735308351756611901741276655161
 
+    # === Trigonometric Functions === #
+    print(a.sin(precision=200))
+    # 0.99985093087193092464780008002600992896256609588456
+    #   91036188395766389946401881352599352354527727927177
+    #   79589259132243649550891532070326452232864052771477
+    #   31418817041042336608522984511928095747763538486886
+    print(b.cos(precision=1000))
+    # -0.9969577603867772005841841569997528013669868536239849713029893885930748434064450375775817720425329394
+    #    9756020177557431933434791661179643984869397089102223199519409695771607230176923201147218218258755323
+    #    7563476302904118661729889931783126826250691820526961290122532541861737355873869924820906724540889765
+    #    5940445990824482174517106016800118438405307801022739336016834311018727787337447844118359555063575166
+    #    5092352912854884589824773945355279792977596081915868398143592738704592059567683083454055626123436523
+    #    6998108941189617922049864138929932713499431655377552668020889456390832876383147018828166124313166286
+    #    6004871998201597316078894718748251490628361253685772937806895692619597915005978762245497623003811386
+    #    0913693867838452088431084666963414694032898497700907783878500297536425463212578556546527017688874265
+    #    0785862902484462361413598747384083001036443681873292719322642381945064144026145428927304407689433744
+    #    5821277763016669042385158254006302666602333649775547203560187716156055524418512492782302125286330865
+
     # === Internal representation of the number === #
     (
         BDec(
diff --git a/docs/internal_notes.md b/docs/internal_notes.md
index 671f4ea..13d95c4 100644
--- a/docs/internal_notes.md
+++ b/docs/internal_notes.md
@@ -2,7 +2,18 @@
 
 ## Values and results
 
-- For power functionality: `BigDecimal.power()`, Python's decimal, WolframAlpha give the same result, but `mpmath` gives a different result. Eample: `0.123456789 ** 1000`, `1234523894766789 ** 1098.1209848`.
+- For power functionality: `BigDecimal.power()`, Python's decimal, WolframAlpha give the same result, but `mpmath` gives a different result. Eamples: 
+  - `0.123456789 ** 1000`
+  - `1234523894766789 ** 1098.1209848`
+- For sin functionality: `BigDecimal.sin()` and WolframAlpha give the same results, but `mpmath` gives a different result. This occurs mainly for pi-related values. Examples:
+  - `sin(3.1415926535897932384626433833)`, precision 50:
+    - Decimojo:     -2.0497115802830600624894179025055407692183593713791E-29
+    - WolframAlpha: -2.0497115802830600624894179025055407692183593713791 x 10-29
+    - mpmath:       -2.049711580283060062489453928920860542175349360102e-29
+  - `sin(6.2831853071795864769252867666)`, precision 50:
+    - Decimojo:     4.4.0994231605661201249788358050110815384367187427582E-29
+    - WolframAlpha: 4.4.0994231605661201249788358050110815384367187427582 x 10-29
+    - mpmath:       4.0994231605661201249789078578417210843506987202039e-29
 
 ## Time complexity
 
diff --git a/src/decimojo/bigdecimal/bigdecimal.mojo b/src/decimojo/bigdecimal/bigdecimal.mojo
index 74885e2..7e0388b 100644
--- a/src/decimojo/bigdecimal/bigdecimal.mojo
+++ b/src/decimojo/bigdecimal/bigdecimal.mojo
@@ -523,6 +523,15 @@ struct BigDecimal(
         """
         return decimojo.bigdecimal.arithmetics.truncate_divide(self, other)
 
+    @always_inline
+    fn __mod__(self, other: Self) raises -> Self:
+        """Returns the result of modulo operation.
+        See `arithmetics.truncate_modulo()` for more information.
+        """
+        return decimojo.bigdecimal.arithmetics.truncate_modulo(
+            self, other, precision=28
+        )
+
     @always_inline
     fn __pow__(self, exponent: Self) raises -> Self:
         """Returns the result of exponentiation."""
@@ -753,6 +762,16 @@ struct BigDecimal(
         return decimojo.bigdecimal.exponential.power(self, exponent, precision)
 
     # === Trigonometric operations === #
+    @always_inline
+    fn sin(self, precision: Int = 28) raises -> Self:
+        """Returns the sine of the BigDecimal number."""
+        return decimojo.bigdecimal.trigonometric.sin(self, precision)
+
+    @always_inline
+    fn cos(self, precision: Int = 28) raises -> Self:
+        """Returns the cosine of the BigDecimal number."""
+        return decimojo.bigdecimal.trigonometric.cos(self, precision)
+
     @always_inline
     fn arctan(self, precision: Int = 28) raises -> Self:
         """Returns the arctangent of the BigDecimal number."""
@@ -826,6 +845,13 @@ struct BigDecimal(
         fill_zeros_to_precision: Bool,
     ) raises:
         """Rounds the number to the specified precision in-place.
+
+        Notes:
+
+        Note that precision is the number of significant digits,
+        not the number of decimal places. If you want to round to a
+        specific number of decimal places, use `round()` instead.
+
         See `rounding.round_to_precision()` for more information.
         """
         decimojo.bigdecimal.rounding.round_to_precision(
diff --git a/src/decimojo/bigdecimal/constants.mojo b/src/decimojo/bigdecimal/constants.mojo
index 849f963..7b2c265 100644
--- a/src/decimojo/bigdecimal/constants.mojo
+++ b/src/decimojo/bigdecimal/constants.mojo
@@ -166,7 +166,14 @@ fn pi(precision: Int) raises -> BigDecimal:
 
     # Use precomputed value for precision ≤ 1024
     if precision <= 1024:
-        return PI_1024.round(precision, RoundingMode.ROUND_HALF_EVEN)
+        var result = PI_1024
+        result.round_to_precision(
+            precision,
+            RoundingMode.ROUND_HALF_EVEN,
+            remove_extra_digit_due_to_rounding=True,
+            fill_zeros_to_precision=False,
+        )
+        return result^
 
     # Use Chudnovsky with binary splitting for maximum speed
     return pi_chudnovsky_binary_split(precision)
@@ -212,7 +219,13 @@ fn pi_chudnovsky_binary_split(precision: Int) raises -> BigDecimal:
     # Final formula: π = 426880 * √10005 / sum_series
     var result = bdec_426880 * bdec_10005.sqrt(working_precision) * sum_series
 
-    return result.round(precision, RoundingMode.ROUND_HALF_EVEN)
+    result.round_to_precision(
+        precision,
+        RoundingMode.ROUND_HALF_EVEN,
+        remove_extra_digit_due_to_rounding=True,
+        fill_zeros_to_precision=False,
+    )
+    return result^
 
 
 fn chudnovsky_split(a: Int, b: Int, precision: Int) raises -> Rational:
@@ -309,4 +322,10 @@ fn pi_machin(precision: Int) raises -> BigDecimal:
     var pi_over_4 = term1 - term2
     var result = bdec_4 * pi_over_4
 
-    return result.round(precision, RoundingMode.ROUND_HALF_EVEN)
+    result.round_to_precision(
+        precision,
+        RoundingMode.ROUND_HALF_EVEN,
+        remove_extra_digit_due_to_rounding=True,
+        fill_zeros_to_precision=False,
+    )
+    return result^
diff --git a/src/decimojo/bigdecimal/rounding.mojo b/src/decimojo/bigdecimal/rounding.mojo
index cdf3b36..c0e7bfe 100644
--- a/src/decimojo/bigdecimal/rounding.mojo
+++ b/src/decimojo/bigdecimal/rounding.mojo
@@ -99,9 +99,9 @@ fn round_to_precision(
             RoundingMode.ROUND_UP: Round up.
             RoundingMode.ROUND_HALF_UP: Round half up.
             RoundingMode.ROUND_HALF_EVEN: Round half even.
-        remove_extra_digit_due_to_rounding: If True, remove an trailing.
-            digit if the rounding mode result in an extra digit.
-        fill_zeros_to_precision: If True, fill zeros to the precision.
+        remove_extra_digit_due_to_rounding: If True, remove a trailing digit if
+            the rounding mode result in an extra leading digit.
+        fill_zeros_to_precision: If True, fill trailing zeros to the precision.
     """
 
     var ndigits_coefficient = number.coefficient.number_of_digits()
diff --git a/src/decimojo/bigdecimal/trigonometric.mojo b/src/decimojo/bigdecimal/trigonometric.mojo
index 56be886..19af310 100644
--- a/src/decimojo/bigdecimal/trigonometric.mojo
+++ b/src/decimojo/bigdecimal/trigonometric.mojo
@@ -30,8 +30,274 @@ import decimojo.utility
 # ===----------------------------------------------------------------------=== #
 
 
-fn sin(x: BigDecimal, precision: Int) raises:
-    ...
+fn sin(x: BigDecimal, precision: Int) raises -> BigDecimal:
+    """Calculates sine (sin) of the number.
+
+    Args:
+        x: The input number in radians.
+        precision: The desired precision of the result.
+
+    Returns:
+        The sine of x with the specified precision.
+
+    Notes:
+    This function adopts range reduction for optimal convergence.
+    """
+
+    # Yuhao Zhu's notes:
+    # I use a very comservative number of buffer digits because we need to have
+    # a very high precision to calculate the pi so that we can conduct range
+    # reduction accurately.
+    # Otherwise, the result will be inaccurate when x is close to π-related
+    # values, e.g., π/2, π, 3π/2, 2π, etc.
+    alias BUFFER_DIGITS = 99
+    var working_precision = precision + BUFFER_DIGITS
+
+    var result: BigDecimal
+
+    if x.is_zero():
+        return BigDecimal.from_raw_components(
+            UInt32(0), scale=precision, sign=x.sign
+        )
+
+    var bdec_2 = BigDecimal.from_raw_components(UInt32(2), scale=0, sign=False)
+    var bdec_4 = BigDecimal.from_raw_components(UInt32(4), scale=0, sign=False)
+    var bdec_6 = BigDecimal.from_raw_components(UInt32(6), scale=0, sign=False)
+    var bdec_pi = decimojo.bigdecimal.constants.pi(precision=working_precision)
+    var bdec_2pi = bdec_2 * bdec_pi
+    var bdec_pi_div_2 = bdec_pi.true_divide(bdec_2, precision=working_precision)
+    var bdec_1d6 = BigDecimal.from_raw_components(
+        UInt32(16), scale=1, sign=False
+    )
+    var bdec_pi_div_4 = bdec_pi.true_divide(bdec_4, precision=working_precision)
+
+    # Step 1: Reduce to (-2π, 2π) using modulo and symmetry
+    # sin(x) = sin(x mod 2π)
+    var x_reduced: BigDecimal
+    if x.compare_absolute(bdec_2pi) >= 0:
+        # x_reduced = x mod 2π
+        x_reduced = x % bdec_2pi
+    else:
+        x_reduced = x
+
+    # Step 2: Reduce [-2π, -6] or [6, 2π] to [6-2π, 2π-6]
+    # sin(x) = sin(x - 2π)
+    # This is because 2π is an instable point for comparison.
+    # To avoid infinite recursion in the final step,
+    # we reduce it to [6-2π, 2π-6].
+    if x_reduced.compare_absolute(bdec_6) >= 0:
+        if x_reduced.sign:
+            # x in [-2π, -6], reduce to [0, 2π-6]
+            x_reduced += bdec_2pi
+        else:
+            # x in [6, 2π], reduce to [0, 2π-6]
+            x_reduced -= bdec_2pi
+
+    # Step 2: Reduce to [0, 2π) using symmetry
+    # At this stage, the value should be in the range [0, 6].
+    var is_negative: Bool
+    if x_reduced.sign:
+        is_negative = True
+        x_reduced = -x_reduced
+    else:
+        is_negative = False
+
+    # Step 3: Reduce to [0, π/4] with different cases
+
+    # |x| ≤ π/4: Use Taylor series directly
+    if x_reduced.compare_absolute(bdec_pi_div_4) <= 0:
+        result = sin_taylor_series(
+            x_reduced, minimum_precision=working_precision
+        )
+
+    # π/4 < |x| ≤ π/2: Use identity sin(x) = cos(π/2 - x)
+    # 0 ≤ (π/2 - x) < π/4
+    # Use 1.6 because π/4 is an instable point for the next case.
+    # To avoid infinite recursion, we use 1.6 as a threshold.
+    # π/4 < |x| ≤ 1.6
+    elif x_reduced.compare_absolute(bdec_1d6) <= 0:
+        x_reduced = bdec_pi_div_2 - x_reduced
+        result = cos_taylor_series(
+            x_reduced, minimum_precision=working_precision
+        )
+
+    # π/2 < |x| ≤ π: Use identity sin(x) = sin(π - x)
+    # 0 ≤ (π - x) < π/2
+    # Because 1.6 is used as a threshold before
+    # π/2 < 1.6 < |x| ≤ π
+    # 0 ≤ (π - x) < π - 1.6 < π/2
+    elif x_reduced.compare_absolute(bdec_pi) <= 0:
+        x_reduced = bdec_pi - x_reduced
+        result = sin(x_reduced, precision=precision)
+
+    # π < |x| < 2π: Use identity sin(x) = -sin(x - π)
+    # 0 < (x - π) < π
+    # Note tha the acutal range is (π, 6), so it is reduced to (0, 6 - π).
+    else:
+        x_reduced = x_reduced - bdec_pi
+        result = -sin(x_reduced, precision=precision)
+
+    if is_negative:
+        result = -result
+
+    result.round_to_precision(
+        precision,
+        RoundingMode.ROUND_HALF_EVEN,
+        remove_extra_digit_due_to_rounding=True,
+        fill_zeros_to_precision=False,
+    )
+
+    return result^
+
+
+fn sin_taylor_series(
+    x: BigDecimal, minimum_precision: Int
+) raises -> BigDecimal:
+    """Calculates sine of a number with Taylor series.
+
+    Args:
+        x: The input number in radians.
+        minimum_precision: The minimum precision of the result.
+
+    Returns:
+        The sine of the input number with the specified precision plus
+        some extra digits to ensure accuracy.
+
+    Notes:
+
+    Using Taylor series.
+    sin(x) = x - x³/3! + x⁵/5! - x⁷/7! + ...
+    """
+
+    alias BUFFER_DIGITS = 9  # word-length, easy to append and trim
+    var working_precision = minimum_precision + BUFFER_DIGITS
+
+    if x.is_zero():
+        return BigDecimal.from_raw_components(
+            UInt32(0), scale=minimum_precision, sign=x.sign
+        )
+
+    var term = x  # x^n / n!
+    var result = x
+    var x_squared = x * x
+    var n = 1
+    var sign = -1
+
+    # Continue until term is smaller than desired precision
+    var epsilon = BigDecimal(BigUInt.ONE, scale=working_precision, sign=False)
+
+    while term.compare_absolute(epsilon) > 0:
+        # x^n = x^(n-2) * x^2 / ((n-1)(n))
+        n += 2
+        term = term * x_squared
+        term = term.true_divide(
+            BigDecimal(n) * BigDecimal(n - 1), precision=working_precision
+        )
+        if sign == 1:
+            result += term
+        else:
+            result -= term
+        sign *= -1
+
+        # Ensure that the result will not explode in size
+        result.round_to_precision(
+            working_precision,
+            rounding_mode=RoundingMode.ROUND_DOWN,
+            remove_extra_digit_due_to_rounding=False,
+            fill_zeros_to_precision=False,
+        )
+
+    return result^
+
+
+fn cos(x: BigDecimal, precision: Int) raises -> BigDecimal:
+    """Calculates cosine (cos) of the number.
+
+    Args:
+        x: The input number in radians.
+        precision: The desired precision of the result.
+
+    Returns:
+        The cosine of x with the specified precision.
+
+    Notes:
+    This function adopts range reduction for optimal convergence.
+    """
+
+    alias BUFFER_DIGITS = 99
+    var working_precision = precision + BUFFER_DIGITS
+
+    if x.is_zero():
+        return BigDecimal.from_raw_components(
+            UInt32(1), scale=precision, sign=x.sign
+        )
+
+    # cos(x) = sin(π/2 - x)
+    var pi = decimojo.bigdecimal.constants.pi(precision=working_precision)
+    var pi_div_2 = pi.true_divide(2, precision=working_precision)
+    var result = sin(pi_div_2 - x, precision=precision)
+    return result^
+
+
+fn cos_taylor_series(
+    x: BigDecimal, minimum_precision: Int
+) raises -> BigDecimal:
+    """Calculates cosine using Taylor series.
+
+    Args:
+        x: The input number in radians.
+        minimum_precision: The minimum precision of the result.
+
+    Returns:
+        The cosine of the input number with the specified precision plus
+        some extra digits to ensure accuracy.
+
+    Notes:
+
+    Using Taylor series.
+    cos(x) = 1 - x²/2! + x⁴/4! - x⁶/6! + ...
+    """
+
+    alias BUFFER_DIGITS = 9
+    var working_precision = minimum_precision + BUFFER_DIGITS
+
+    if x.is_zero():
+        return BigDecimal.from_raw_components(
+            UInt32(1), scale=minimum_precision, sign=x.sign
+        )
+
+    var bdec_1 = BigDecimal.from_raw_components(UInt32(1), scale=0, sign=False)
+    var term = bdec_1  # Current term: x^n / n!
+    var result = bdec_1  # Start with 1
+    var x_squared = x * x
+    var n = 0  # Current power (0, 2, 4, 6, ...)
+    var sign = -1  # Alternating sign
+
+    var epsilon = BigDecimal(BigUInt.ONE, scale=working_precision, sign=False)
+
+    while term.compare_absolute(epsilon) > 0:
+        n += 2  # Next even power: 2, 4, 6, 8, ...
+        term = term * x_squared
+        term = term.true_divide(
+            BigDecimal(n) * BigDecimal(n - 1), precision=working_precision
+        )
+
+        if sign == 1:
+            result += term
+        else:
+            result -= term
+
+        sign *= -1
+
+        # # Prevent size explosion
+        result.round_to_precision(
+            working_precision,
+            rounding_mode=RoundingMode.ROUND_DOWN,
+            remove_extra_digit_due_to_rounding=False,
+            fill_zeros_to_precision=False,
+        )
+
+    return result^
 
 
 # ===----------------------------------------------------------------------=== #
@@ -56,12 +322,13 @@ fn arctan(x: BigDecimal, precision: Int) raises -> BigDecimal:
     bdec_2 = BigDecimal.from_raw_components(UInt32(2), scale=0, sign=False)
     bdec_0d5 = BigDecimal.from_raw_components(UInt32(5), scale=1, sign=False)
 
+    var result: BigDecimal
+
     if x.compare_absolute(bdec_0d5) <= 0:
         # |x| <= 0.5, use Taylor series:
         # print("Using Taylor series for arctan with |x| <= 0.5")
-        return arctan_taylor_series(x, minimum_precision=precision).round(
-            ndigits=precision, rounding_mode=RoundingMode.ROUND_HALF_EVEN
-        )
+        result = arctan_taylor_series(x, minimum_precision=precision)
+
     elif x.compare_absolute(bdec_2) <= 0:
         # |x| <= 2, use the identity:
         # arctan(x) = 2 * arctan(x / (1 + sqrt(1 + x²)))
@@ -71,12 +338,10 @@ fn arctan(x: BigDecimal, precision: Int) raises -> BigDecimal:
         var x_divided = x.true_divide(
             bdec_1 + sqrt_term, precision=working_precision
         )
-        var result = bdec_2 * arctan_taylor_series(
+        result = bdec_2 * arctan_taylor_series(
             x_divided, minimum_precision=precision
         )
-        return result.round(
-            ndigits=precision, rounding_mode=RoundingMode.ROUND_HALF_EVEN
-        )
+
     else:  # x.compare_absolute(bdec_1) > 0
         # |x| > 2, use the identity:
         # For x > 2: arctan(x) = π/2 - arctan(1/x)
@@ -91,15 +356,18 @@ fn arctan(x: BigDecimal, precision: Int) raises -> BigDecimal:
             reciprocal_x^, minimum_precision=precision
         )
 
-        var result: BigDecimal
         if x.sign:
             result = -half_pi - arctan_reciprocal
         else:
             result = half_pi - arctan_reciprocal
 
-        return result.round(
-            ndigits=precision, rounding_mode=RoundingMode.ROUND_HALF_EVEN
-        )
+    result.round_to_precision(
+        precision,
+        RoundingMode.ROUND_HALF_EVEN,
+        remove_extra_digit_due_to_rounding=True,
+        fill_zeros_to_precision=False,
+    )
+    return result^
 
 
 fn arctan_taylor_series(
@@ -120,16 +388,15 @@ fn arctan_taylor_series(
     Using Taylor series.
     arctan(x) = x - x³/3 + x⁵/5 - x⁷/7 + ...
     The input x must be in the range (-0.5, 0.5) for convergence.
-    Time complexity is O(n^4) for precision n.
-    Every time you double the precision, the time taken increases by a
-    factor of 16.
     """
 
     alias BUFFER_DIGITS = 9  # word-length, easy to append and trim
     var working_precision = minimum_precision + BUFFER_DIGITS
 
     if x.is_zero():
-        return BigDecimal(0)
+        return BigDecimal.from_raw_components(
+            UInt32(0), scale=minimum_precision, sign=x.sign
+        )
 
     var term = x  # x^n
     var term_divided = x  # x^n / n
diff --git a/tests/bigdecimal/test_bigdecimal_trigonometric.mojo b/tests/bigdecimal/test_bigdecimal_trigonometric.mojo
index d94f0a8..0704cdc 100644
--- a/tests/bigdecimal/test_bigdecimal_trigonometric.mojo
+++ b/tests/bigdecimal/test_bigdecimal_trigonometric.mojo
@@ -11,18 +11,15 @@ from decimojo.tests import TestCase, parse_file, load_test_cases
 alias file_path = "tests/bigdecimal/test_data/bigdecimal_trigonometric.toml"
 
 
-fn test_bigdecimal_trignometric() raises:
-    # Load test cases from TOML file
-    var toml = parse_file(file_path)
-    var test_cases: List[TestCase]
-
-    print("------------------------------------------------------")
-    print("Testing BigDecimal arctan...")
+fn run_test[
+    func: fn (BDec, Int) raises -> BDec
+](toml: tomlmojo.parser.TOMLDocument, table_name: String, msg: String) raises:
+    """Run a specific test case from the TOML document."""
     print("------------------------------------------------------")
-
-    test_cases = load_test_cases(toml, "arctan_tests")
+    print("Testing BigDecimal ", msg, "...", sep="")
+    var test_cases = load_test_cases(toml, table_name)
     for test_case in test_cases:
-        var result = BDec(test_case.a).arctan(precision=50)
+        var result = func(BDec(test_case.a), 50)
         testing.assert_equal(
             lhs=String(result),
             rhs=test_case.expected,
@@ -30,6 +27,22 @@ fn test_bigdecimal_trignometric() raises:
         )
 
 
+fn test_bigdecimal_trignometric() raises:
+    # Load test cases from TOML file
+    var toml = parse_file(file_path)
+
+    run_test[func = decimojo.bigdecimal.trigonometric.arctan](
+        toml,
+        "arctan_tests",
+        "arctan",
+    )
+    run_test[func = decimojo.bigdecimal.trigonometric.sin](
+        toml,
+        "sin_tests",
+        "sin",
+    )
+
+
 fn main() raises:
     print("Running BigDecimal trigonometric tests")
 
diff --git a/tests/bigdecimal/test_data/bigdecimal_trigonometric.toml b/tests/bigdecimal/test_data/bigdecimal_trigonometric.toml
index b77da5f..a600fa8 100644
--- a/tests/bigdecimal/test_data/bigdecimal_trigonometric.toml
+++ b/tests/bigdecimal/test_data/bigdecimal_trigonometric.toml
@@ -1,3 +1,127 @@
+# ===----------------------------------------------------------------------=== #
+# Test cases for sin
+# ===----------------------------------------------------------------------=== #
+
+[[sin_tests]]
+a = "0"
+b = ""
+expected = "0E-50"
+description = "sin(0) = 0"
+
+[[sin_tests]]
+a = "0.5235987755982988730771072305"
+b = ""
+expected = "0.49999999999999999999999999995965723364237186253345"
+description = "sin(π/6) = 0.5"
+
+[[sin_tests]]
+a = "0.7853981633974483096156608458"
+b = ""
+expected = "0.70710678118654752440084436209079478214990911530509"
+description = "sin(π/4) = √2/2"
+
+[[sin_tests]]
+a = "1.0471975511965977461542144611"
+b = ""
+expected = "0.86602540378443864676372317075635236943854106034267"
+description = "sin(π/3) = √3/2"
+
+[[sin_tests]]
+a = "1.5707963267948966192313216916"
+b = ""
+expected = "1.0000000000000000000000000000000000000000000000000"
+description = "sin(π/2) = 1"
+
+[[sin_tests]]
+a = "2.0943951023931954923084289221"
+b = ""
+expected = "0.86602540378443864676372317079610381153712576003023"
+description = "sin(2π/3) = √3/2"
+
+[[sin_tests]]
+a = "2.3561944901923449288469825375"
+b = ""
+expected = "0.70710678118654752440084436207630113257096165239855"
+description = "sin(3π/4) = √2/2"
+
+[[sin_tests]]
+a = "2.6179938779914943653855361833"
+b = ""
+expected = "0.49999999999999999999999997352813139522743036701006"
+description = "sin(5π/6) = 0.5"
+
+[[sin_tests]]
+a = "3.1415926535897932384626433833"
+b = ""
+expected = "-2.0497115802830600624894179025055407692183593713791E-29"
+description = "sin(π) = 0"
+
+[[sin_tests]]
+a = "4.7123889803846898576939650749"
+b = ""
+expected = "-1.0000000000000000000000000000000000000000000000000"
+description = "sin(3π/2) = -1"
+
+[[sin_tests]]
+a = "6.2831853071795864769252867666"
+b = ""
+expected = "4.0994231605661201249788358050110815384367187427582E-29"
+description = "sin(2π) = 0"
+
+[[sin_tests]]
+a = "-0.5235987755982988730771072305"
+b = ""
+expected = "-0.49999999999999999999999999995965723364237186253345"
+description = "sin(-π/6) = -0.5"
+
+[[sin_tests]]
+a = "-1.5707963267948966192313216916"
+b = ""
+expected = "-1.0000000000000000000000000000000000000000000000000"
+description = "sin(-π/2) = -1"
+
+[[sin_tests]]
+a = "-3.1415926535897932384626433833"
+b = ""
+expected = "2.0497115802830600624894179025055407692183593713791E-29"
+description = "sin(-π) = 0"
+
+[[sin_tests]]
+a = "0.1"
+b = ""
+expected = "0.099833416646828152306814198410622026989915388017982"
+description = "sin(0.1) - small positive value"
+
+[[sin_tests]]
+a = "0.3"
+b = ""
+expected = "0.29552020666133957510532074568502737367783211174262"
+description = "sin(0.3) - moderate small value"
+
+[[sin_tests]]
+a = "1.0"
+b = ""
+expected = "0.84147098480789650665250232163029899962256306079837"
+description = "sin(1) - one radian"
+
+[[sin_tests]]
+a = "2.0"
+b = ""
+expected = "0.90929742682568169539601986591174484270225497144789"
+description = "sin(2) - two radians"
+
+[[sin_tests]]
+a = "10.0"
+b = ""
+expected = "-0.54402111088936981340474766185137728168364301291622"
+description = "sin(10) - large value test"
+
+[[sin_tests]]
+a = "0.0001"
+b = ""
+expected = "0.000099999999833333333416666666646825396828152557318973"
+description = "sin(0.0001) - very small value (linear approximation region)"
+
 # ===----------------------------------------------------------------------=== #
 # Test cases for arctan
 # ===----------------------------------------------------------------------=== #
@@ -5,7 +129,7 @@
 [[arctan_tests]]
 a = "0"
 b = ""
-expected = "0.000000000000000000000000000000000000000000000E-6"
+expected = "0E-50"
 description = "Arctan of 0"
 
 [[arctan_tests]]
@@ -41,7 +165,7 @@ description = "Arctan of 0.2 (small positive value)"
 [[arctan_tests]]
 a = "0.1"
 b = ""
-expected = "0.09966865249116202737844611987802059024327832250431"
+expected = "0.099668652491162027378446119878020590243278322504315"
 description = "Arctan of 0.1 (very small value for good convergence)"
 
 [[arctan_tests]]
@@ -71,31 +195,31 @@ description = "Arctan of 0.9"
 [[arctan_tests]]
 a = "2"
 b = ""
-expected = "1.10714871779409050301706546017853704007004764540143"
+expected = "1.1071487177940905030170654601785370400700476454014"
 description = "Arctan of 2 (>1, tests convergence for larger values)"
 
 [[arctan_tests]]
 a = "3"
 b = ""
-expected = "1.24904577239825442582991707728109012307782940412990"
+expected = "1.2490457723982544258299170772810901230778294041299"
 description = "Arctan of 3"
 
 [[arctan_tests]]
 a = "5"
 b = ""
-expected = "1.37340076694501586086127192644496114865099959589970"
+expected = "1.3734007669450158608612719264449611486509995958997"
 description = "Arctan of 5 (large value)"
 
 [[arctan_tests]]
 a = "-10"
 b = ""
-expected = "-1.47112767430373459185287557176173085185530637718324"
+expected = "-1.4711276743037345918528755717617308518553063771832"
 description = "Arctan of -10 (negative large value)"
 
 [[arctan_tests]]
 a = "-20"
 b = ""
-expected = "-1.52083793107295385782131540460490656060730761926405"
+expected = "-1.5208379310729538578213154046049065606073076192640"
 description = "Arctan of -20 (negative large value)"
 
 [[arctan_tests]]
@@ -113,17 +237,17 @@ description = "Arctan of tan(π/6) = √3/3 should equal π/6"
 [[arctan_tests]]
 a = "-1.7320508075688772935274463415058723669428052538103806280558069794519330169088000"
 b = ""
-expected = "-1.04719755119659774615421446109316762806572313312504"
+expected = "-1.0471975511965977461542144610931676280657231331250"
 description = "Arctan of -√3 should equal -π/3"
 
 [[arctan_tests]]
 a = "0.00001"
 b = ""
-expected = "9.99999999966666666668666666666523809523820635E-6"
+expected = "9.9999999996666666666866666666652380952382063492063E-6"
 description = "Arctan of very small value (near linear approximation)"
 
 [[arctan_tests]]
 a = "100"
 b = ""
-expected = "1.56079666010823138102498157543047189353721534714318"
+expected = "1.5607966601082313810249815754304718935372153471432"
 description = "Arctan of 100 (very large value, should approach π/2)"