diff --git a/.circleci/config.yml b/.circleci/config.yml
index d740dbb..bf4eb2b 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -14,7 +14,7 @@ jobs:
         default: "buster"
       rust-features:
         type: string
-        default: ""
+        default: "--all-targets"
     docker:
       - image: rust:<< parameters.rust-version >>-<< parameters.debian-version >>
     environment:
@@ -31,7 +31,7 @@ jobs:
             - bigdecimal-cargo-
       - run:
           name: Check
-          command: cargo check
+          command: cargo check << parameters.rust-features >>
       - save_cache:
           paths:
             - /usr/local/cargo
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index f938a8a..d67a2f9 100755
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -28,13 +28,13 @@
     before_script:
       - rustc --version && cargo --version
     script:
-      - cargo check --tests
+      - cargo check --all-targets
 
 .script-cargo-build: &script-cargo-build
   before_script:
     - rustc --version && cargo --version
   script:
-    - cargo build --tests
+    - cargo build --all-targets
 
 .script-cargo-test: &script-cargo-test
   before_script:
@@ -118,7 +118,7 @@ cargo:check:
   <<: *script-cargo-check
   script:
     # enable property tests for the stable 'pipeline'
-    - scripts/bigdecimal-property-tests cargo check --tests
+    - scripts/bigdecimal-property-tests run cargo check --all-targets
 
 cargo:clippy:
   stage: check
@@ -158,7 +158,7 @@ cargo:build-stable:
   <<: *script-cargo-build
   script:
     # enable property tests for the stable 'pipeline'
-    - scripts/bigdecimal-property-tests cargo build --tests
+    - scripts/bigdecimal-property-tests run cargo build --all-targets
 
 
 cargo:test-stable:
@@ -171,7 +171,7 @@ cargo:test-stable:
   <<: *script-cargo-test
   script:
     # enable property tests for the stable 'pipeline'
-    - scripts/bigdecimal-property-tests cargo test
+    - scripts/bigdecimal-property-tests run cargo test
 
 
 cargo:build:no-std:
@@ -197,6 +197,29 @@ cargo:test:no-std:
     - cargo test --no-default-features --lib
 
 
+cargo:build:serde:
+  stage: build
+  image: akubera/rust:stable
+  needs:
+    - cargo:check
+  variables:
+    RUST_CACHE_KEY: "stable+serde"
+  <<: *script-cargo-build
+  script:
+    - cargo build --features=serde --all-targets
+
+cargo:test:serde:
+  stage: test
+  image: akubera/rust:stable
+  needs:
+    - "cargo:build:serde"
+  variables:
+    RUST_CACHE_KEY: "stable+serde"
+  <<: *script-cargo-test
+  script:
+    - cargo test --features=serde --all-targets
+
+
 cargo:build-nightly:
   stage: build
   image: rustlang/rust:nightly
@@ -243,39 +266,9 @@ cargo:test-1.43:
   <<: *script-cargo-test
 
 
-cargo:check-1.54:
-  stage: check
-  image: "akubera/rust-kcov:1.54.0-bullseye"
-  rules:
-    *rules-always-master-otherwise-manual
-  variables:
-    RUST_CACHE_KEY: "1.54"
-  <<: *script-cargo-check
-
-cargo:build-1.54:
-  stage: build
-  image: "akubera/rust-kcov:1.54.0-bullseye"
-  needs:
-    - "cargo:check-1.54"
-  variables:
-    RUST_CACHE_KEY: "1.54"
-  <<: *script-cargo-build
-
-cargo:test-1.54:
-  stage: test
-  needs:
-    - "cargo:build-1.54"
-  image: "akubera/rust-kcov:1.54.0-bullseye"
-  variables:
-    RUST_CACHE_KEY: "1.54"
-  <<: *script-cargo-test
-
-
 cargo:check-1.70:
   stage: check
   image: "akubera/rust-grcov:1.70.0-bullseye"
-  rules:
-    *rules-always-master-otherwise-manual
   variables:
     RUST_CACHE_KEY: "1.70"
   <<: *script-cargo-check
diff --git a/Cargo.toml b/Cargo.toml
index 25c2935..49bd37e 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "bigdecimal"
-version = "0.4.1"
+version = "0.4.2"
 authors = ["Andrew Kubera"]
 description = "Arbitrary precision decimal numbers"
 documentation = "https://docs.rs/bigdecimal"
@@ -47,6 +47,6 @@ default = ["std"]
 string-only = []
 std = ["num-bigint/std", "num-integer/std", "num-traits/std"]
 
-[[bench]]
-name = "arithmetic"
-harness = false
+# BENCH: [[bench]]
+# BENCH: name = "arithmetic"
+# BENCH: harness = false
diff --git a/benches/arithmetic.rs b/benches/arithmetic.rs
index b6b3299..04ff995 100644
--- a/benches/arithmetic.rs
+++ b/benches/arithmetic.rs
@@ -1,4 +1,4 @@
-//! Benchmarks for arithmetic opertaion
+//! Benchmarks for arithmetic operation
 
 extern crate criterion;
 extern crate bigdecimal;
diff --git a/build.rs b/build.rs
index 7b35ce4..8e43e1a 100644
--- a/build.rs
+++ b/build.rs
@@ -1,50 +1,48 @@
-#![allow(clippy::style)]
-
-extern crate autocfg;
 
 use std::env;
-use std::path::PathBuf;
+use std::path::{Path, PathBuf};
 
 
-fn main() -> std::io::Result<()> {
+fn main() {
     let ac = autocfg::new();
     ac.emit_rustc_version(1, 70);
 
-    let outdir = match std::env::var_os("OUT_DIR") {
-        None => return Ok(()),
-        Some(outdir) => outdir,
-    };
-    let outdir_path = PathBuf::from(outdir);
+    // Option::zip
+    ac.emit_rustc_version(1, 46);
+
+    // Remove this comment if enabled proptests
+    // ::PROPERTY-TESTS:: autocfg::emit("property_tests");
 
-    write_default_precision(&outdir_path, "default_precision.rs")?;
-    Ok(())
+    let outdir: PathBuf = std::env::var_os("OUT_DIR").unwrap().into();
+    write_default_precision_file(&outdir);
+    write_default_rounding_mode(&outdir);
 }
 
-/// Create default_precision.rs, containg definition of constant DEFAULT_PRECISION
-fn write_default_precision(outdir_path: &PathBuf, filename: &str) -> std::io::Result<()>
-{
 
-    let default_prec = env::var("RUST_BIGDECIMAL_DEFAULT_PRECISION")
-        .map(|s| s.parse::<std::num::NonZeroU32>().expect("$RUST_BIGDECIMAL_DEFAULT_PRECISION must be an integer > 0"))
-        .map(|nz_num| nz_num.into())
-        .unwrap_or(100u32);
+/// Create default_precision.rs, containing definition of constant DEFAULT_PRECISION loaded in src/lib.rs
+fn write_default_precision_file(outdir: &Path) {
+    let env_var = env::var("RUST_BIGDECIMAL_DEFAULT_PRECISION").unwrap_or_else(|_| "100".to_owned());
+    println!("cargo:rerun-if-env-changed=RUST_BIGDECIMAL_DEFAULT_PRECISION");
+
+    let rust_file_path = outdir.join("default_precision.rs");
 
-    let default_precision_rs_path = outdir_path.join(filename);
+    let default_prec: u32 = env_var
+        .parse::<std::num::NonZeroU32>()
+        .expect("$RUST_BIGDECIMAL_DEFAULT_PRECISION must be an integer > 0")
+        .into();
 
-    let default_precision = format!("const DEFAULT_PRECISION: u64 = {};", default_prec);
+    let rust_file_contents = format!("const DEFAULT_PRECISION: u64 = {};", default_prec);
+
+    std::fs::write(rust_file_path, rust_file_contents).unwrap();
+}
 
-    // Rewriting the file if it already exists with the same contents
-    // would force a rebuild.
-    match std::fs::read_to_string(&default_precision_rs_path) {
-        Ok(existing_contents) if existing_contents == default_precision => {},
-        _ => {
-            std::fs::write(&default_precision_rs_path, default_precision)
-                    .expect("Could not write big decimal default-precision file");
-        }
-    };
+/// Create default_rounding_mode.rs, using value of RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE environment variable
+fn write_default_rounding_mode(outdir: &Path) {
+    let rounding_mode_name = env::var("RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE").unwrap_or_else(|_| "HalfEven".to_owned());
+    println!("cargo:rerun-if-env-changed=RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE");
 
-    println!("cargo:rerun-if-changed={}", default_precision_rs_path.display());
-    println!("cargo:rerun-if-env-changed={}", "RUST_BIGDECIMAL_DEFAULT_PRECISION");
+    let rust_file_path = outdir.join("default_rounding_mode.rs");
+    let rust_file_contents = format!("const DEFAULT_ROUNDING_MODE: RoundingMode = RoundingMode::{};", rounding_mode_name);
 
-    Ok(())
+    std::fs::write(rust_file_path, rust_file_contents).unwrap();
 }
diff --git a/scripts/benchmark-bigdecimal b/scripts/benchmark-bigdecimal
index 04c08f2..bca871b 100755
--- a/scripts/benchmark-bigdecimal
+++ b/scripts/benchmark-bigdecimal
@@ -1,4 +1,4 @@
-#!/bin/sh
+#!/bin/bash
 #
 # Run Criterion Benchmarks
 #
@@ -14,7 +14,7 @@ mv Cargo.toml.bak Cargo.toml
 
 
 # store extra things for the benchmarking report
-if [ ! -z "$BENCHMARK_EXTRAS" ]; then
+if [ -n "$BENCHMARK_EXTRAS" ]; then
 cat <<EOF > index.html
 <!doctype html>
 <head>
@@ -36,12 +36,12 @@ EOF
 # Add svg plots to index html
 find target/criterion -name 'pdf.svg' -type f -print0 |
 sort -z |
-while read -d $'\0' svg_file
+while read -r -d $'\0' svg_file
 do
-  name=$(echo $svg_file | cut -d '/' -f 3)
+  name=$(echo "$svg_file" | cut -d '/' -f 3)
 
   sample_datafile=target/criterion/$name/new/sample.json
-  if [ -f $sample_datafile ]; then
+  if [ -f "$sample_datafile" ]; then
     echo "<p><a href='${sample_datafile}' class='sample-json'>$name</a>" >> index.html
   else
     echo "<p>$name</p>" >> index.html
diff --git a/scripts/bigdecimal-property-tests b/scripts/bigdecimal-property-tests
index bce2e79..8196fcb 100755
--- a/scripts/bigdecimal-property-tests
+++ b/scripts/bigdecimal-property-tests
@@ -5,15 +5,45 @@
 # Tests are defined in src/lib.tests.property-test.rs
 #
 
-# enable property-test dependencies in Cargo
-sed -i.bak -e 's|# PROPERTY-TESTS: ||' Cargo.toml
+enable_property_tests() {
+  # enable property-test dependencies in Cargo
+  sed -i.bak -e 's|# PROPERTY-TESTS: ||' Cargo.toml
 
-# include the property-test file in lib.rs
-sed -i.bak -e 's|// ::PROPERTY-TESTS:: ||' src/lib.rs
+  # add the property-test configuration in build.rs
+  sed -i.bak -e 's|// ::PROPERTY-TESTS:: ||' build.rs
+}
 
-# Run commands
-"$@"
 
-# Restore Cargo.toml with backup
-mv Cargo.toml.bak Cargo.toml
-mv src/lib.rs.bak src/lib.rs
+restore_disabled_property_tests() {
+  # Restore Cargo.toml with backup
+  mv Cargo.toml.bak Cargo.toml
+  mv build.rs.bak build.rs
+}
+
+
+DEFAULT_CMD=run
+CMD=${1:-$DEFAULT_CMD}
+shift
+
+case "${CMD}" in
+  run)
+    enable_property_tests
+    # Run commands
+    "$@"
+    restore_disabled_property_tests
+    ;;
+
+  test)
+    enable_property_tests
+    cargo test "$@"
+    restore_disabled_property_tests
+    ;;
+
+  enable)
+    enable_property_tests
+    ;;
+
+  disable)
+    restore_disabled_property_tests
+    ;;
+esac
diff --git a/scripts/fetch-benchmark-data.sh b/scripts/fetch-benchmark-data.sh
index 83af9b1..26d8d73 100755
--- a/scripts/fetch-benchmark-data.sh
+++ b/scripts/fetch-benchmark-data.sh
@@ -16,16 +16,16 @@ function fetch_benchmark_bigdecimal_file() {
 	local FILENAME="random-bigdecimals-$1.txt"
 	local FILEPATH=$TEST_DATA_DIR/$FILENAME
 
-	if [ -e $FILEPATH ]; then
+	if [ -e "$FILEPATH" ]; then
 		echo "exists: $FILEPATH"
 	else
 		local URL=${GITLAB_URL_PATTERN//<FILENAME>/$FILENAME}
 		echo "fetching: $FILEPATH from $URL"
 
-		if [ $CURL ]; then
-			$CURL -s --fail -L $URL -o "$FILEPATH"
-		elif [ $WGET ]; then
-			$WGET --quiet $URL -O "$FILEPATH"
+		if [ "$CURL" ]; then
+			$CURL -s --fail -L "$URL" -o "$FILEPATH"
+		elif [ "$WGET" ]; then
+			"$WGET" --quiet "$URL" -O "$FILEPATH"
 		else
 			echo "No supported fetching program"
 		fi
diff --git a/src/arithmetic/cbrt.rs b/src/arithmetic/cbrt.rs
new file mode 100644
index 0000000..8887fb9
--- /dev/null
+++ b/src/arithmetic/cbrt.rs
@@ -0,0 +1,151 @@
+//! Implementation of cube-root algorithm
+
+use crate::*;
+use num_bigint::BigUint;
+
+
+/// implementation of cuberoot - always positive
+pub(crate) fn impl_cbrt_uint_scale(n: &BigUint, scale: i64, ctx: &Context) -> BigDecimal {
+    // make guess based on number of bits in the number
+    let guess = make_cbrt_guess(n.bits(), scale);
+
+    let three = BigInt::from(3);
+
+    let n = BigInt::from_biguint(Sign::Plus, n.clone());
+
+    let max_precision = ctx.precision().get();
+
+    let next_iteration = move |r: BigDecimal| {
+        let sqrd = r.square();
+        let tmp = impl_division(
+            n.clone(),
+            &sqrd.int_val,
+            scale - sqrd.scale,
+            max_precision + 1,
+        );
+        let tmp = tmp + r.double();
+        impl_division(tmp.int_val, &three, tmp.scale, max_precision + 3)
+    };
+
+    // result initial
+    let mut running_result = next_iteration(guess);
+
+    let mut prev_result = BigDecimal::one();
+    let mut result = BigDecimal::zero();
+
+    // TODO: Prove that we don't need to arbitrarily limit iterations
+    // and that convergence can be calculated
+    while prev_result != result {
+        // store current result to test for convergence
+        prev_result = result;
+
+        running_result = next_iteration(running_result);
+
+        // result has clipped precision, running_result has full precision
+        result = if running_result.digits() > max_precision {
+            running_result.with_precision_round(ctx.precision(), ctx.rounding_mode())
+        } else {
+            running_result.clone()
+        };
+    }
+
+    return result;
+}
+
+/// Find initial cbrt guess based on number of bits in integer and the scale
+///
+/// ```math
+/// 2^bit_count * 10^-scale <= *n* < 2^(bit_count+1) * 10^-scale
+///
+/// cbrt(n2^bit_count * 10^-scale)
+/// cbrt(2^bit_count * 10^-scale)
+///    => Exp2[1/3 * Log2[2^bit_count * 10^-scale]]
+///    => Exp2[1/3 * (bit_count - scale * Log2[10]]
+/// ```
+///
+fn make_cbrt_guess(bit_count: u64, scale: i64) -> BigDecimal {
+    // weight of cube root average above minimum within range: 3/4*2^(4/3)
+    let magic_guess_scale = 1.1398815748423097_f64;
+
+    let bit_count = bit_count as f64;
+    let scale = scale as f64;
+
+    let initial_guess = (bit_count - scale * LOG2_10) / 3.0;
+    let res = magic_guess_scale * exp2(initial_guess);
+
+    match BigDecimal::try_from(res) {
+        Ok(res) => res,
+        Err(_) => {
+            // can't guess with float - just guess magnitude
+            let scale = (scale - bit_count / LOG2_10).round() as i64;
+            BigDecimal::new(BigInt::from(1), scale / 3)
+        }
+    }
+}
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use stdlib::num::NonZeroU64;
+
+    #[test]
+    fn test_cbrt() {
+        let vals = vec![
+            ("0.00", "0"),
+            ("1.00", "1"),
+            ("1.001", "1.000333222283909495175449559955220102010284758197360454054345461242739715702641939155238095670636841"),
+            ("10", "2.154434690031883721759293566519350495259344942192108582489235506346411106648340800185441503543243276"),
+            ("13409.179789484375", "23.7575"),
+            ("-59283293e25", "-84006090355.84281237113712383191213626687332139035750444925827809487776780721673264524620270275301685"),
+            ("94213372931e-127", "2.112049945275324414051072540210070583697242797173805198575907094646677475250362108901530353886613160E-39"),
+        ];
+        for &(x, y) in vals.iter() {
+            let a = BigDecimal::from_str(x).unwrap().cbrt();
+            let b = BigDecimal::from_str(y).unwrap();
+            assert_eq!(a, b);
+        }
+    }
+
+
+    #[test]
+    fn test_cbrt_prec15() {
+        let vals = vec![
+            ("0.00", "0"),
+            ("1.00", "1"),
+            ("1.001", "1.00033322228390"),
+            ("10", "2.15443469003188"),
+            ("13409.179789484375", "23.7575"),
+            ("-59283293e25", "-84006090355.8428"),
+            ("94213372931e-127", "2.11204994527532E-39"),
+        ];
+
+        let ctx = Context::new(NonZeroU64::new(15).unwrap(), RoundingMode::Down);
+
+        for &(x, y) in vals.iter() {
+            let a = BigDecimal::from_str(x).unwrap().cbrt_with_context(&ctx);
+            let b = y.parse().unwrap();
+            assert_eq!(a, b);
+        }
+    }
+
+    #[cfg(property_tests)]
+    mod prop {
+        use super::*;
+        use proptest::*;
+        use num_traits::FromPrimitive;
+
+        proptest! {
+            #[test]
+            fn cbrt_of_cube_is_self(f: f64, prec in 15..50u64) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+
+                let n = BigDecimal::from_f64(f).unwrap().with_prec(prec);
+                let n_cubed = n.cube();
+                let x = n_cubed.cbrt();
+                prop_assert_eq!(x, n);
+            }
+        }
+    }
+}
diff --git a/src/arithmetic/inverse.rs b/src/arithmetic/inverse.rs
new file mode 100644
index 0000000..398463e
--- /dev/null
+++ b/src/arithmetic/inverse.rs
@@ -0,0 +1,233 @@
+//! inverse implementation
+
+use crate::*;
+
+/// Implementation of inverse: (1/n)
+pub(crate) fn impl_inverse_uint_scale(n: &BigUint, scale: i64, ctx: &Context) -> BigDecimal {
+    let guess = make_inv_guess(n.bits(), scale);
+    let max_precision = ctx.precision().get();
+
+    let s = BigDecimal::new(BigInt::from_biguint(Sign::Plus, n.clone()), scale);
+    let two = BigDecimal::from(2);
+
+    let next_iteration = move |r: BigDecimal| {
+        let tmp = &two - &s * &r;
+        r * tmp
+    };
+
+    // calculate first iteration
+    let mut running_result = next_iteration(guess);
+    debug_assert!(!running_result.is_zero(), "Zero detected in inverse calculation of {}e{}", n, -scale);
+
+    let mut prev_result = BigDecimal::one();
+    let mut result = BigDecimal::zero();
+
+    // TODO: Prove that we don't need to arbitrarily limit iterations
+    // and that convergence can be calculated
+    while prev_result != result {
+        // store current result to test for convergence
+        prev_result = result;
+
+        // calculate next iteration
+        running_result = next_iteration(running_result).with_prec(max_precision + 2);
+
+        // 'result' has clipped precision, 'running_result' has full precision
+        result = if running_result.digits() > max_precision {
+            running_result.with_precision_round(ctx.precision(), ctx.rounding_mode())
+        } else {
+            running_result.clone()
+        };
+    }
+
+    return result;
+}
+
+
+/// guess inverse based on the number of bits in the integer and decimal's scale
+fn make_inv_guess(bit_count: u64, scale: i64) -> BigDecimal {
+    // scale by ln(2)
+    let magic_factor = stdlib::f64::consts::LN_2;
+
+    let bit_count = bit_count as f64;
+    let initial_guess = magic_factor * exp2(-bit_count);
+    if initial_guess.is_finite() && initial_guess != 0.0 {
+        if let Ok(mut result) = BigDecimal::try_from(initial_guess) {
+            result.scale -= scale;
+            return result;
+        }
+    }
+
+    // backup guess for out-of-range integers
+
+    let approx_scale = bit_count * stdlib::f64::consts::LOG10_2;
+    let approx_scale_int = approx_scale.trunc();
+    let approx_scale_frac = approx_scale - approx_scale_int;
+
+    let recip = libm::exp10(-approx_scale_frac);
+    let mut res = BigDecimal::from_f32((magic_factor * recip) as f32).unwrap();
+    res.scale += approx_scale_int as i64;
+    res.scale -= scale;
+    return res;
+}
+
+
+#[cfg(test)]
+mod test_make_inv_guess {
+    use super::*;
+    use paste::paste;
+
+    macro_rules! impl_case {
+        ( $bin_count:literal, -$scale:literal => $expected:literal ) => {
+            paste! { impl_case!( [< case_ $bin_count _n $scale >]: $bin_count, -$scale => $expected); }
+        };
+        ( $bin_count:literal, $scale:literal => $expected:literal ) => {
+            paste! { impl_case!( [< case_ $bin_count _ $scale >]: $bin_count, $scale => $expected); }
+        };
+        ( $name:ident: $bin_count:expr, $scale:expr => $expected:literal ) => {
+            impl_case!($name: $bin_count, $scale, prec=5 => $expected);
+        };
+        ( $name:ident: $bin_count:expr, $scale:expr, prec=$prec:literal => $expected:literal ) => {
+            #[test]
+            fn $name() {
+                let guess = make_inv_guess($bin_count, $scale);
+                let expected: BigDecimal = $expected.parse().unwrap();
+                assert_eq!(guess.with_prec($prec), expected.with_prec($prec));
+            }
+        };
+    }
+
+    impl_case!(0, 0 => "0.69315");
+    impl_case!(1, 0 => "0.34657");
+    impl_case!(2, 0 => "0.17329");
+    impl_case!(2, 1 => "1.7329");
+
+    // 1 / (2^3 * 10^5) ~
+    impl_case!(3, -5 => "8.6643e-07");
+
+    // 2^-20
+    impl_case!(20, 0 => "6.6104e-07");
+    impl_case!(20, -900 => "6.6104E-907");
+    impl_case!(20, 800 => "6.6104E+793");
+
+    impl_case!(40, 10000 => "6.3041E+9987");
+
+    impl_case!(70, -5 => "5.8712e-27");
+    impl_case!(70, 5 => "5.8712e-17");
+    impl_case!(70, 50 => "5.8712e+28");
+
+    impl_case!(888, -300 => "3.3588E-568");
+    impl_case!(888, -19 => "3.3588E-287");
+    impl_case!(888, 0 => "3.3588E-268");
+    impl_case!(888, 270 => "335.88");
+
+    impl_case!(1022, 10 => "1.5423e-298");
+    impl_case!(1022, 308 => "1.5423");
+
+    impl_case!(1038, 316 => "2353.4");
+
+    impl_case!(case_31028_n659: 31028, -659 => "3.0347E-10000");
+    impl_case!(case_31028_0: 31028, 0 => "3.0347E-9341");
+    impl_case!(case_31028_1: 31028, 1 => "3.0347E-9340");
+    impl_case!(case_31028_9340: 31028, 9340 => ".30347");
+    impl_case!(case_31028_10000: 31028, 10000 => "3.0347E+659");
+
+    // impl_case!(case_max: u64::MAX, 270 => "335.88");
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use stdlib::num::NonZeroU64;
+
+    #[test]
+    fn test_inverse_35543972957198043e291() {
+        let v = vec![
+            0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+            0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+            2324389888, 849200558
+        ];
+        let x = BigInt::new(Sign::Minus, v);
+        let d = BigDecimal::from(x);
+        let expected = "-2.813416500187520746852694701086705659180043761702417561798711758892800449936819185796527214192677476E-308".parse().unwrap();
+        assert_eq!(d.inverse(), expected);
+
+        assert_eq!(d.neg().inverse(), expected.neg());
+    }
+
+    macro_rules! impl_case {
+        ($name:ident: $prec:literal, $round:ident => $expected:literal) => {
+            #[test]
+            fn $name() {
+                let n = test_input();
+                let prec = NonZeroU64::new($prec).unwrap();
+                let rounding = RoundingMode::$round;
+                let ctx = Context::new(prec, rounding);
+
+                let result = n.inverse_with_context(&ctx);
+
+                let expected = $expected.parse().unwrap();
+                assert_eq!(&result, &expected);
+
+                let product = result * &n;
+                let epsilon = BigDecimal::new(BigInt::one(), $prec - 1);
+                let diff = (BigDecimal::one() - &product).abs();
+                assert!(diff < epsilon);
+            }
+        };
+    }
+
+    mod invert_seven {
+        use super::*;
+
+        fn test_input() -> BigDecimal {
+            BigDecimal::from(7u8)
+        }
+
+        impl_case!(case_prec10_round_down: 10, Down => "0.1428571428");
+        impl_case!(case_prec10_round_up: 10, Up => "0.1428571429");
+
+        impl_case!(case_prec11_round_ceiling: 11, Ceiling => "0.14285714286");
+    }
+
+
+    #[test]
+    fn inv_random_number() {
+       let n = BigDecimal::try_from(0.08121970592310568).unwrap();
+
+       let ctx = Context::new(NonZeroU64::new(40).unwrap(), RoundingMode::Down);
+       let i = n.inverse_with_context(&ctx);
+       assert_eq!(&i, &"12.31228294456944530942557443718279245563".parse().unwrap());
+
+       let product = i * &n;
+       assert!(BigDecimal::one() - &product < "1e-39".parse().unwrap());
+    }
+
+    #[cfg(property_tests)]
+    mod prop {
+        use super::*;
+        use proptest::*;
+        use num_traits::FromPrimitive;
+
+        proptest! {
+
+            #[test]
+            fn inverse_multiplies_to_one(f: f64, prec in 1..100u64) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+                prop_assume!(f != 0.0);
+
+                let n = BigDecimal::from_f64(f).unwrap();
+
+                let ctx = Context::new(NonZeroU64::new(prec).unwrap(), RoundingMode::Up);
+                let i = n.inverse_with_context(&ctx);
+                let product = &i * &n;
+
+                // accurate to precision minus one (due to rounding)
+                let epsilon = BigDecimal::new(1.into(), prec as i64 - 1);
+                let diff_from_one = BigDecimal::one() - &product;
+
+                prop_assert!(diff_from_one.abs() < epsilon, "{} >= {}", diff_from_one.abs(), epsilon);
+            }
+        }
+    }
+}
diff --git a/src/arithmetic/mod.rs b/src/arithmetic/mod.rs
new file mode 100644
index 0000000..614503d
--- /dev/null
+++ b/src/arithmetic/mod.rs
@@ -0,0 +1,76 @@
+//! arithmetic routines
+
+use crate::*;
+
+pub(crate) mod sqrt;
+pub(crate) mod cbrt;
+pub(crate) mod inverse;
+
+/// Return 10^pow
+///
+/// Try to calculate this with fewest number of allocations
+///
+pub(crate) fn ten_to_the(pow: u64) -> BigInt {
+    ten_to_the_uint(pow).into()
+}
+
+/// Return 10^pow
+pub(crate) fn ten_to_the_uint(pow: u64) -> BigUint {
+    if pow < 20 {
+        return BigUint::from(10u64.pow(pow as u32));
+    }
+
+    // linear case of 10^pow = 10^(19 * count + rem)
+    if pow < 590 {
+        let ten_to_nineteen = 10u64.pow(19);
+
+        // count factors of 19
+        let (count, rem) = pow.div_rem(&19);
+
+        let mut res = BigUint::from(ten_to_nineteen);
+        for _ in 1..count {
+            res *= ten_to_nineteen;
+        }
+        if rem != 0 {
+            res *= 10u64.pow(rem as u32);
+        }
+
+        return res;
+    }
+
+    // use recursive algorithm where linear case might be too slow
+    let (quotient, rem) = pow.div_rem(&16);
+    let x = ten_to_the_uint(quotient);
+
+    let x2 = &x * &x;
+    let x4 = &x2 * &x2;
+    let x8 = &x4 * &x4;
+    let res = &x8 * &x8;
+
+    if rem == 0 {
+        res
+    } else {
+        res * 10u64.pow(rem as u32)
+    }
+}
+
+
+/// Return number of decimal digits in integer
+pub(crate) fn count_decimal_digits(int: &BigInt) -> u64 {
+    count_decimal_digits_uint(int.magnitude())
+}
+
+/// Return number of decimal digits in unsigned integer
+pub(crate) fn count_decimal_digits_uint(uint: &BigUint) -> u64 {
+    if uint.is_zero() {
+        return 1;
+    }
+    let mut digits = (uint.bits() as f64 / LOG2_10) as u64;
+    // guess number of digits based on number of bits in UInt
+    let mut num = ten_to_the(digits).to_biguint().expect("Ten to power is negative");
+    while *uint >= num {
+        num *= 10u8;
+        digits += 1;
+    }
+    digits
+}
diff --git a/src/arithmetic/sqrt.rs b/src/arithmetic/sqrt.rs
new file mode 100644
index 0000000..d2a56df
--- /dev/null
+++ b/src/arithmetic/sqrt.rs
@@ -0,0 +1,258 @@
+//! square root implementation
+
+use crate::*;
+
+
+fn impl_division(mut num: BigUint, den: &BigUint, mut scale: i64, max_precision: u64) -> BigDecimal {
+    use Sign::Plus;
+
+    // quick zero check
+    if num.is_zero() {
+        return BigDecimal::new(BigInt::from_biguint(Plus, num), 0);
+    }
+
+    // shift digits until numerator is larger than denominator (set scale appropriately)
+    while num < *den {
+        scale += 1;
+        num *= 10u8;
+    }
+
+    // first division
+    let (mut quotient, mut remainder) = num.div_rem(den);
+
+    // division complete
+    if remainder.is_zero() {
+        return BigDecimal {
+            int_val: BigInt::from_biguint(Plus, quotient),
+            scale: scale,
+        };
+    }
+
+    let mut precision = count_decimal_digits_uint(&quotient);
+
+    // shift remainder by 1 decimal;
+    // quotient will be 1 digit upon next division
+    remainder *= 10u8;
+
+    while !remainder.is_zero() && precision < max_precision {
+        let (q, r) = remainder.div_rem(den);
+        quotient = quotient * 10u8 + q;
+        remainder = r * 10u8;
+
+        precision += 1;
+        scale += 1;
+    }
+
+    if !remainder.is_zero() {
+        // round final number with remainder
+        quotient += get_rounding_term_uint(&remainder.div(den));
+    }
+
+    BigDecimal::new(BigInt::from_biguint(Plus, quotient), scale)
+}
+
+fn get_rounding_term_uint(num: &BigUint) -> u8 {
+    if num.is_zero() {
+        return 0;
+    }
+
+    let digits = (num.bits() as f64 / LOG2_10) as u64;
+    let mut n = ten_to_the_uint(digits);
+
+    // loop-method
+    loop {
+        if *num < n {
+            return 1;
+        }
+        n *= 5u8;
+        if *num < n {
+            return 0;
+        }
+        n *= 2u8;
+    }
+
+    // string-method
+    // let s = format!("{}", num);
+    // let high_digit = u8::from_str(&s[0..1]).unwrap();
+    // if high_digit < 5 { 0 } else { 1 }
+}
+
+pub(crate) fn impl_sqrt(n: &BigUint, scale: i64, ctx: &Context) -> BigDecimal {
+    // make guess
+    let guess = {
+        let magic_guess_scale = 1.1951678538495576_f64;
+        let initial_guess = (n.bits() as f64 - scale as f64 * LOG2_10) / 2.0;
+        let res = magic_guess_scale * exp2(initial_guess);
+
+        if res.is_normal() {
+            BigDecimal::try_from(res).unwrap()
+        } else {
+            // can't guess with float - just guess magnitude
+            let scale = (n.bits() as f64 / -LOG2_10 + scale as f64).round() as i64;
+            BigDecimal::new(BigInt::from(1), scale / 2)
+        }
+    };
+
+    // // wikipedia example - use for testing the algorithm
+    // if self == &BigDecimal::from_str("125348").unwrap() {
+    //     running_result = BigDecimal::from(600)
+    // }
+
+    // TODO: Use context variable to set precision
+    let max_precision = ctx.precision().get();
+
+    let next_iteration = move |r: BigDecimal| {
+        // division needs to be precise to (at least) one extra digit
+        let tmp = impl_division(
+            n.clone(),
+            &r.int_val.magnitude(),
+            scale - r.scale,
+            max_precision + 1,
+        );
+
+        // half will increase precision on each iteration
+        (tmp + r).half()
+    };
+
+    // calculate first iteration
+    let mut running_result = next_iteration(guess);
+
+    let mut prev_result = BigDecimal::one();
+    let mut result = BigDecimal::zero();
+
+    // TODO: Prove that we don't need to arbitrarily limit iterations
+    // and that convergence can be calculated
+    while prev_result != result {
+        // store current result to test for convergence
+        prev_result = result;
+
+        // calculate next iteration
+        running_result = next_iteration(running_result);
+
+        // 'result' has clipped precision, 'running_result' has full precision
+        result = if running_result.digits() > max_precision {
+            running_result.with_prec(max_precision)
+        } else {
+            running_result.clone()
+        };
+    }
+
+    result
+}
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    #[test]
+    fn test_sqrt() {
+        let vals = vec![
+            ("1e-232", "1e-116"),
+            ("1.00", "1"),
+            ("1.001", "1.000499875062460964823258287700109753027590031219780479551442971840836093890879944856933288426795152"),
+            ("100", "10"),
+            ("49", "7"),
+            (".25", ".5"),
+            ("0.0152399025", ".12345"),
+            ("152399025", "12345"),
+            (".00400", "0.06324555320336758663997787088865437067439110278650433653715009705585188877278476442688496216758600590"),
+            (".1", "0.3162277660168379331998893544432718533719555139325216826857504852792594438639238221344248108379300295"),
+            ("2", "1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573"),
+            ("125348", "354.0451948551201563108487193176101314241016013304294520812832530590100407318465590778759640828114535"),
+            ("18446744073709551616.1099511", "4294967296.000000000012799992691725492477397918722952224079252026972356303360555051219312462698703293"),
+            ("3.141592653589793115997963468544185161590576171875", "1.772453850905515992751519103139248439290428205003682302442979619028063165921408635567477284443197875"),
+            (".000000000089793115997963468544185161590576171875", "0.000009475922962855041517561783740144225422359796851494316346796373337470068631250135521161989831460407155"),
+            ("0.7177700109762963922745342343167413624881759290454997218753321040760896053150388903350654937434826216803814031987652326749140535150336357405672040727695124057298138872112244784753994931999476811850580200000000000000000000000000000000", "0.8472130847527653667042980517799020703921106560594525833177762276594388966885185567535692987624493813"),
+            ("0.01234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901", "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111"),
+            ("0.1108890000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000444", "0.3330000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000667"),
+        ];
+        for &(x, y) in vals.iter() {
+            let a = BigDecimal::from_str(x).unwrap().sqrt().unwrap();
+            let b = BigDecimal::from_str(y).unwrap();
+            assert_eq!(a, b);
+        }
+    }
+
+    #[test]
+    fn test_big_sqrt() {
+        use num_bigint::BigInt;
+        let vals = vec![
+            (("2", -70), "141421356237309504880168872420969807.8569671875376948073176679737990732478462107038850387534327641573"),
+            (("3", -170), "17320508075688772935274463415058723669428052538103806280558069794519330169088000370811.46186757248576"),
+            (("9", -199), "9486832980505137995996680633298155601158665417975650480572514558377783315917714664032744325137900886"),
+            (("7", -200), "26457513110645905905016157536392604257102591830824501803683344592010688232302836277603928864745436110"),
+            (("777", -204), "2.787471972953270789531596912111625325974789615194854615319795902911796043681078997362635440358922503E+103"),
+            (("7", -600), "2.645751311064590590501615753639260425710259183082450180368334459201068823230283627760392886474543611E+300"),
+            (("2", -900), "1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573E+450"),
+            (("7", -999), "8.366600265340755479781720257851874893928153692986721998111915430804187725943170098308147119649515362E+499"),
+            (("74908163946345982392040522594123773796", -999), "2.736935584670307552030924971360722787091742391079630976117950955395149091570790266754718322365663909E+518"),
+            (("20", -1024), "4.472135954999579392818347337462552470881236719223051448541794490821041851275609798828828816757564550E512"),
+            (("3", 1025), "5.477225575051661134569697828008021339527446949979832542268944497324932771227227338008584361638706258E-513"),
+        ];
+        for &((s, scale), e) in vals.iter() {
+            let expected = BigDecimal::from_str(e).unwrap();
+
+            let sqrt = BigDecimal::new(BigInt::from_str(s).unwrap(), scale).sqrt().unwrap();
+            assert_eq!(sqrt, expected);
+        }
+    }
+
+    #[test]
+    fn case_sqrt_3242053850483855em13() {
+        let d: BigDecimal = "324.2053850483855".parse().unwrap();
+
+        let digitref = d.to_ref();
+        let (_, scale, uint) = digitref.as_parts();
+        let ctx = Context::default()
+                          .with_prec(11).unwrap()
+                          .with_rounding_mode(RoundingMode::Down);
+
+        let s = impl_sqrt(uint, scale, &ctx);
+        let expected: BigDecimal = "18.005704236".parse().unwrap();
+        assert_eq!(s, expected);
+
+        let ctx = Context::default()
+                          .with_prec(31).unwrap()
+                          .with_rounding_mode(RoundingMode::Up);
+
+        let s = impl_sqrt(uint, scale, &ctx);
+        let expected: BigDecimal = "18.00570423639090823994825477228".parse().unwrap();
+        assert_eq!(s, expected);
+    }
+
+    #[test]
+    fn case_sqrt_5085019992340351em25() {
+        let d: BigDecimal = "5.085019992340351e-10".parse().unwrap();
+
+        let digitref = d.to_ref();
+        let (_, scale, uint) = digitref.as_parts();
+        let ctx = Context::default()
+                          .with_prec(25).unwrap()
+                          .with_rounding_mode(RoundingMode::Down);
+
+        let s = impl_sqrt(uint, scale, &ctx);
+        let expected: BigDecimal = "0.00002254998889653906459324292".parse().unwrap();
+        assert_eq!(s, expected);
+    }
+
+    #[cfg(property_tests)]
+    mod prop {
+        use super::*;
+        use proptest::*;
+        use num_traits::FromPrimitive;
+
+        proptest! {
+            #[test]
+            fn sqrt_of_square_is_self(f: f64, prec in 15..50u64) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+
+                let n = BigDecimal::from_f64(f.abs()).unwrap().with_prec(prec);
+                let n_squared = n.square();
+                let x = n_squared.sqrt().unwrap();
+                prop_assert_eq!(x, n);
+            }
+        }
+    }
+}
diff --git a/src/context.rs b/src/context.rs
new file mode 100644
index 0000000..c48cc84
--- /dev/null
+++ b/src/context.rs
@@ -0,0 +1,175 @@
+//! Define arithmetical context
+//!
+
+use crate::*;
+use stdlib::num::NonZeroU64;
+
+
+// const DEFAULT_PRECISION: u64 = ${RUST_BIGDECIMAL_DEFAULT_PRECISION} or 100;
+include!(concat!(env!("OUT_DIR"), "/default_precision.rs"));
+// const DEFAULT_ROUNDING_MODE: RoundingMode = ${RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE} or HalfUp;
+include!(concat!(env!("OUT_DIR"), "/default_rounding_mode.rs"));
+
+
+/// Mathematical Context
+///
+/// Stores rules for numerical operations, such as how to round and
+/// number of digits to keep.
+///
+/// Defaults are defined at compile time, determined by environment
+/// variables:
+///
+/// | Variable                                |   Descripiton   | default  |
+/// |-----------------------------------------|-----------------|----------|
+/// | `RUST_BIGDECIMAL_DEFAULT_PRECISION`     | digit precision | 100      |
+/// | `RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE` | rounding-mode   | HalfEven |
+///
+/// It is recommended that the user set explicit values of a Context and *not*
+/// rely on compile time constants, but the option is there if necessary.
+///
+#[derive(Debug, Clone)]
+pub struct Context {
+    /// total number of digits
+    precision: NonZeroU64,
+    /// how to round
+    rounding: RoundingMode,
+}
+
+impl Context {
+    /// Create context with precision and rounding mode
+    pub fn new(precision: NonZeroU64, rounding: RoundingMode) -> Self {
+        Context {
+            precision: precision,
+            rounding: rounding,
+        }
+    }
+
+    /// Copy context with new precision value
+    pub fn with_precision(&self, precision: NonZeroU64) -> Self {
+        Self {
+            precision: precision,
+            ..*self
+        }
+    }
+
+    /// Copy context with new precision value
+    pub fn with_prec<T: ToPrimitive>(&self, precision: T) -> Option<Self> {
+        precision
+            .to_u64()
+            .and_then(NonZeroU64::new)
+            .map(|prec| {
+                Self {
+                    precision: prec,
+                    ..*self
+                }
+            })
+    }
+
+    /// Copy context with new rounding mode
+    pub fn with_rounding_mode(&self, mode: RoundingMode) -> Self {
+        Self {
+            rounding: mode,
+            ..*self
+        }
+    }
+
+    /// Return maximum precision
+    pub fn precision(&self) -> NonZeroU64 {
+        self.precision
+    }
+
+    /// Return rounding mode
+    pub fn rounding_mode(&self) -> RoundingMode {
+        self.rounding
+    }
+
+    /// Round digits x and y with the rounding mode
+    pub(crate) fn round_pair(&self, sign: Sign, x: u8, y: u8, trailing_zeros: bool) -> u8 {
+        self.rounding.round_pair(sign, (x, y), trailing_zeros)
+    }
+
+    /// Round digits x and y with the rounding mode
+    #[allow(dead_code)]
+    pub(crate) fn round_pair_with_carry(&self, sign: Sign, x: u8, y: u8, trailing_zeros: bool, carry: &mut u8) -> u8 {
+        let r = self.round_pair(sign, x, y, trailing_zeros);
+        if r == 10 {
+            *carry = 1;
+            0
+        } else {
+            r
+        }
+    }
+}
+
+impl stdlib::default::Default for Context {
+    fn default() -> Self {
+        Self {
+            precision: NonZeroU64::new(DEFAULT_PRECISION).unwrap(),
+            rounding: DEFAULT_ROUNDING_MODE,
+        }
+    }
+}
+
+impl Context {
+    /// Add two big digit references
+    pub fn add_refs<'a, 'b, A, B>(&self, a: A, b: B) -> BigDecimal
+    where
+        A: Into<BigDecimalRef<'a>>,
+        B: Into<BigDecimalRef<'b>>,
+    {
+        let mut sum = BigDecimal::zero();
+        self.add_refs_into(a, b, &mut sum);
+        sum
+    }
+
+    /// Add two decimal refs, storing value in dest
+    pub fn add_refs_into<'a, 'b, A, B>(&self, a: A, b: B, dest: &mut BigDecimal)
+    where
+        A: Into<BigDecimalRef<'a>>,
+        B: Into<BigDecimalRef<'b>>,
+    {
+        let a = a.into();
+        let b = b.into();
+        let sum = a.to_owned() + b.to_owned();
+        *dest = sum.with_precision_round(self.precision, self.rounding)
+    }
+}
+
+
+#[cfg(test)]
+mod test_context {
+    use super::*;
+
+    #[test]
+    fn contstructor_and_setters() {
+        let ctx = Context::default();
+        let c = ctx.with_prec(44).unwrap();
+        assert_eq!(c.precision.get(), 44);
+        assert_eq!(c.rounding, RoundingMode::HalfEven);
+
+        let c = c.with_rounding_mode(RoundingMode::Down);
+        assert_eq!(c.precision.get(), 44);
+        assert_eq!(c.rounding, RoundingMode::Down);
+    }
+
+    #[test]
+    fn sum_two_references() {
+        use stdlib::ops::Neg;
+
+        let ctx = Context::default();
+        let a: BigDecimal = "209682.134972197168613072130300".parse().unwrap();
+        let b: BigDecimal = "3.0782968222271332463325639E-12".parse().unwrap();
+
+        let sum = ctx.add_refs(&a, &b);
+        assert_eq!(sum, "209682.1349721971716913689525271332463325639".parse().unwrap());
+
+        // make negative copy of b without cloning values
+        let neg_b = b.to_ref().neg();
+
+        let sum = ctx.add_refs(&a, neg_b);
+        assert_eq!(sum, "209682.1349721971655347753080728667536674361".parse().unwrap());
+
+        let sum = ctx.with_prec(27).unwrap().with_rounding_mode(RoundingMode::Up).add_refs(&a, neg_b);
+        assert_eq!(sum, "209682.134972197165534775309".parse().unwrap());
+    }
+}
diff --git a/src/impl_cmp.rs b/src/impl_cmp.rs
new file mode 100644
index 0000000..32f7019
--- /dev/null
+++ b/src/impl_cmp.rs
@@ -0,0 +1,99 @@
+//! Implementation of comparison operations
+//!
+//! Comparisons between decimals and decimal refs
+//! are not directly supported as we lose some type
+//! inference features at the savings of a single
+//! '&' character.
+//!
+//! &BigDecimal and BigDecimalRef are comparable.
+//!
+
+use crate::{
+    BigDecimal,
+    BigDecimalRef,
+    Sign,
+};
+
+use stdlib::cmp::Ordering;
+use stdlib::iter;
+
+impl PartialEq for BigDecimal
+{
+    fn eq(&self, rhs: &BigDecimal) -> bool {
+        self.to_ref() == rhs.to_ref()
+    }
+}
+
+impl<'rhs, T> PartialEq<T> for BigDecimalRef<'_>
+where
+    T: Into<BigDecimalRef<'rhs>> + Copy
+{
+    fn eq(&self, rhs: &T) -> bool {
+        let rhs: BigDecimalRef<'rhs> = (*rhs).into();
+
+        match (self.sign(), rhs.sign()) {
+            // both zero
+            (Sign::NoSign, Sign::NoSign) => return true,
+            // signs are different
+            (a, b) if a != b => return false,
+            // signs are same, do nothing
+            _ => {}
+        }
+
+        let unscaled_int;
+        let scaled_int;
+        let trailing_zero_count;
+        match self.scale.cmp(&rhs.scale) {
+            Ordering::Greater => {
+                unscaled_int = self.digits;
+                scaled_int = rhs.digits;
+                trailing_zero_count = (self.scale - rhs.scale) as usize;
+            }
+            Ordering::Less => {
+                unscaled_int = rhs.digits;
+                scaled_int = self.digits;
+                trailing_zero_count = (rhs.scale - self.scale) as usize;
+            }
+            Ordering::Equal => return self.digits == rhs.digits,
+        }
+
+        if trailing_zero_count < 20 {
+            let scaled_int = scaled_int * crate::ten_to_the(trailing_zero_count as u64).magnitude();
+            return &scaled_int == unscaled_int;
+        }
+
+        let unscaled_digits = unscaled_int.to_radix_le(10);
+        let scaled_digits = scaled_int.to_radix_le(10);
+
+        // different lengths with trailing zeros
+        if unscaled_digits.len() != scaled_digits.len() + trailing_zero_count {
+            return false;
+        }
+
+        // add leading zero digits to digits that need scaled
+        let scaled = iter::repeat(&0u8).take(trailing_zero_count).chain(scaled_digits.iter());
+
+        // return true if all digits are the same
+        unscaled_digits.iter().zip(scaled).all(|(digit_a, digit_b)| { digit_a == digit_b })
+    }
+}
+
+
+#[cfg(test)]
+mod test_bigintref {
+    use super::*;
+    use stdlib::ops::Neg;
+
+    #[test]
+    fn test_borrow_neg_cmp() {
+        let x: BigDecimal = "1514932018891593.916341142773".parse().unwrap();
+        let y: BigDecimal = "1514932018891593916341142773e-12".parse().unwrap();
+
+        assert_eq!(x, y);
+
+        let x_ref = x.to_ref();
+        assert_eq!(x_ref, &y);
+        assert_ne!(x_ref.neg(), x_ref);
+        assert_eq!(x_ref.neg().neg(), x_ref);
+    }
+}
diff --git a/src/impl_num.rs b/src/impl_num.rs
index 4953cd1..8136973 100644
--- a/src/impl_num.rs
+++ b/src/impl_num.rs
@@ -85,7 +85,15 @@ impl Num for BigDecimal {
             }
         };
 
-        let scale = decimal_offset - exponent_value;
+        let scale = match decimal_offset.checked_sub(exponent_value) {
+            Some(scale) => scale,
+            None => {
+                return Err(ParseBigDecimalError::Other(
+                    format!("Exponent overflow when parsing '{}'", s)
+                ))
+            }
+        };
+
         let big_int = BigInt::from_str_radix(&digits, radix)?;
 
         Ok(BigDecimal::new(big_int, scale))
@@ -163,3 +171,20 @@ impl ToBigInt for BigDecimal {
         Some(self.with_scale(0).int_val)
     }
 }
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    mod from_str_radix {
+        use super::*;
+
+        #[test]
+        fn out_of_bounds() {
+            let d = BigDecimal::from_str_radix("1e-9223372036854775808", 10);
+            assert_eq!(d.unwrap_err(), ParseBigDecimalError::Other("Exponent overflow when parsing '1e-9223372036854775808'".to_string()));
+
+        }
+    }
+}
diff --git a/src/impl_ops.rs b/src/impl_ops.rs
index 8e91e60..565d94b 100644
--- a/src/impl_ops.rs
+++ b/src/impl_ops.rs
@@ -1,18 +1,6 @@
 //! Implement math operations: Add,Sub, etc
 
-use crate::BigDecimal;
-use crate::stdlib::ops::{
-    Add, AddAssign,
-    Sub, SubAssign,
-    Mul, MulAssign,
-    Div, DivAssign,
-    Neg,
-};
-
-
-use crate::stdlib::convert::TryFrom;
-
-use num_traits::{Zero, One};
+use crate::*;
 
 
 macro_rules! impl_add_for_primitive {
@@ -35,6 +23,14 @@ macro_rules! impl_add_for_primitive {
         impl Add<$t> for &BigDecimal {
             type Output = BigDecimal;
 
+            fn add(self, rhs: $t) -> BigDecimal {
+                self.to_ref() + rhs
+            }
+        }
+
+        impl Add<$t> for BigDecimalRef<'_> {
+            type Output = BigDecimal;
+
             fn add(self, rhs: $t) -> BigDecimal {
                 BigDecimal::from(rhs) + self
             }
@@ -308,7 +304,7 @@ macro_rules! impl_div_for_primitive {
             type Output = BigDecimal;
 
             fn div(self, denom: &BigDecimal) -> BigDecimal {
-                BigDecimal::from(self) / denom
+                self / denom.clone()
             }
         }
     };
@@ -439,3 +435,35 @@ impl_div_for_primitive!(i128);
 
 impl_div_for_primitive!(f32);
 impl_div_for_primitive!(f64);
+
+
+impl Neg for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn neg(mut self) -> BigDecimal {
+        self.int_val = -self.int_val;
+        self
+    }
+}
+
+impl<'a> Neg for &'a BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn neg(self) -> BigDecimal {
+        -self.clone()
+    }
+}
+
+impl Neg for BigDecimalRef<'_> {
+    type Output = Self;
+
+    fn neg(self) -> Self::Output {
+        Self {
+            sign: self.sign.neg(),
+            digits: self.digits,
+            scale: self.scale,
+        }
+    }
+}
diff --git a/src/impl_ops_add.rs b/src/impl_ops_add.rs
new file mode 100644
index 0000000..8b80bfd
--- /dev/null
+++ b/src/impl_ops_add.rs
@@ -0,0 +1,391 @@
+//! Addition operator trait implementation
+//!
+
+use super::*;
+use stdlib::borrow::ToOwned;
+
+
+impl Add<BigDecimal> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigDecimal) -> BigDecimal {
+        if rhs.is_zero() {
+            return self;
+        }
+        if self.is_zero() {
+            return rhs;
+        }
+
+        let mut lhs = self;
+        match lhs.scale.cmp(&rhs.scale) {
+            Ordering::Equal => {
+                lhs.int_val += rhs.int_val;
+                lhs
+            }
+            Ordering::Less => lhs.take_and_scale(rhs.scale) + rhs,
+            Ordering::Greater => rhs.take_and_scale(lhs.scale) + lhs,
+        }
+    }
+}
+
+impl<'a, T: Into<BigDecimalRef<'a>>> Add<T> for BigDecimal {
+    type Output = BigDecimal;
+
+    fn add(mut self, rhs: T) -> BigDecimal {
+        self.add_assign(rhs);
+        self
+    }
+}
+
+impl Add<BigInt> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigInt) -> BigDecimal {
+        self + BigDecimal::from(rhs)
+    }
+}
+
+
+
+impl Add<BigDecimal> for &'_ BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigDecimal) -> BigDecimal {
+        rhs + self
+    }
+}
+
+impl<'a, T: Into<BigDecimalRef<'a>>> Add<T> for &'_ BigDecimal {
+    type Output = BigDecimal;
+    fn add(self, rhs: T) -> BigDecimal {
+        let rhs = rhs.into();
+        if rhs.is_zero() {
+            return self.clone();
+        }
+        if self.is_zero() {
+            return rhs.to_owned();
+        }
+        if self.scale >= rhs.scale {
+            self.to_owned() + rhs
+        } else {
+            rhs.to_owned() + self
+        }
+    }
+}
+
+impl Add<BigInt> for &'_ BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigInt) -> BigDecimal {
+        self.to_ref() + rhs
+    }
+}
+
+
+impl Add<BigDecimal> for BigDecimalRef<'_> {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigDecimal) -> BigDecimal {
+        rhs + self
+    }
+}
+
+impl<'a, T: Into<BigDecimalRef<'a>>> Add<T> for BigDecimalRef<'_> {
+    type Output = BigDecimal;
+    fn add(self, rhs: T) -> BigDecimal {
+        let rhs = rhs.into();
+        if self.scale >= rhs.scale {
+            self.to_owned() + rhs
+        } else {
+            rhs.to_owned() + self
+        }
+    }
+}
+
+impl Add<BigInt> for BigDecimalRef<'_> {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigInt) -> BigDecimal {
+        self + BigDecimal::from(rhs)
+    }
+}
+
+
+impl Add<BigDecimal> for BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigDecimal) -> BigDecimal {
+        rhs + self
+    }
+}
+
+impl<'a> Add<&'a BigDecimal> for BigInt {
+    type Output = BigDecimal;
+
+    fn add(self, rhs: &BigDecimal) -> BigDecimal {
+        rhs.to_ref().add(self)
+    }
+}
+
+impl<'a> Add<BigDecimalRef<'a>> for BigInt {
+    type Output = BigDecimal;
+
+    fn add(self, rhs: BigDecimalRef<'_>) -> BigDecimal {
+        rhs.add(self)
+    }
+}
+
+
+impl Add<BigDecimal> for &BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigDecimal) -> BigDecimal {
+        rhs + self
+    }
+}
+
+impl<'a> Add<&'a BigDecimal> for &BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: &BigDecimal) -> BigDecimal {
+        rhs + self
+    }
+}
+
+impl<'a> Add<BigDecimalRef<'a>> for &BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn add(self, rhs: BigDecimalRef<'_>) -> BigDecimal {
+        rhs + self
+    }
+}
+
+
+impl AddAssign<BigDecimal> for BigDecimal {
+    fn add_assign(&mut self, rhs: BigDecimal) {
+        if rhs.is_zero() {
+            return;
+        }
+        if self.is_zero() {
+            *self = rhs;
+            return;
+        }
+        self.add_assign(rhs.to_ref());
+    }
+}
+
+impl<'a, N: Into<BigDecimalRef<'a>>> AddAssign<N> for BigDecimal {
+    #[inline]
+    fn add_assign(&mut self, rhs: N) {
+        // TODO: Replace to_owned() with efficient addition algorithm
+        let rhs = rhs.into().to_owned();
+        match self.scale.cmp(&rhs.scale) {
+            Ordering::Less => {
+                let scaled = self.with_scale(rhs.scale);
+                self.int_val = scaled.int_val + &rhs.int_val;
+                self.scale = rhs.scale;
+            }
+            Ordering::Greater => {
+                let scaled = rhs.with_scale(self.scale);
+                self.int_val += scaled.int_val;
+            }
+            Ordering::Equal => {
+                self.int_val += &rhs.int_val;
+            }
+        }
+    }
+}
+
+impl AddAssign<BigInt> for BigDecimal {
+    #[inline]
+    fn add_assign(&mut self, rhs: BigInt) {
+        self.add_assign(&rhs);
+    }
+}
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use paste::paste;
+
+    macro_rules! impl_case {
+        ( $name:ident: $a:literal + $b:literal => $c:literal ) => {
+            #[test]
+            fn $name() {
+                let a: BigDecimal = $a.parse().unwrap();
+                let b: BigDecimal = $b.parse().unwrap();
+                let c: BigDecimal = $c.parse().unwrap();
+
+                assert_eq!(c, a.clone() + b.clone());
+                assert_eq!(c, a.clone() + b.to_ref());
+                assert_eq!(c, a.clone() + &b);
+
+                assert_eq!(c, &a + b.clone());
+                assert_eq!(c, &a + b.to_ref());
+                assert_eq!(c, &a + &b);
+
+                assert_eq!(c, a.to_ref() + b.clone());
+                assert_eq!(c, a.to_ref() + b.to_ref());
+                assert_eq!(c, a.to_ref() + &b);
+
+                // Reversed
+
+                assert_eq!(c, b.clone() + a.clone());
+                assert_eq!(c, b.clone() + a.to_ref());
+                assert_eq!(c, b.clone() + &a);
+
+                assert_eq!(c, &b + a.clone());
+                assert_eq!(c, &b + a.to_ref());
+                assert_eq!(c, &b + &a);
+
+                assert_eq!(c, b.to_ref() + a.clone());
+                assert_eq!(c, b.to_ref() + a.to_ref());
+                assert_eq!(c, b.to_ref() + &a);
+
+                let mut n = a.clone();
+                n += b.clone();
+                assert_eq!(c, n);
+
+                let mut n = a.clone();
+                n += &b;
+                assert_eq!(c, n);
+
+                let mut n = a.clone();
+                n += b.to_ref();
+                assert_eq!(c, n);
+
+                let mut n = b.clone();
+                n += a.clone();
+                assert_eq!(c, n);
+
+                let mut n = b.clone();
+                n += &a;
+                assert_eq!(c, n);
+
+                let mut n = b.clone();
+                n += a.to_ref();
+                assert_eq!(c, n);
+            }
+        };
+        ( $name:ident: $a:literal + (int) $b:literal => $c:literal ) => {
+            #[test]
+            fn $name() {
+                let a: BigDecimal = $a.parse().unwrap();
+                let b: BigInt = $b.parse().unwrap();
+                let c: BigDecimal = $c.parse().unwrap();
+
+                assert_eq!(c, a.clone() + b.clone());
+                assert_eq!(c, a.clone() + &b);
+                assert_eq!(c, &a + &b);
+                assert_eq!(c, &a + b.clone());
+                assert_eq!(c, a.to_ref() + &b);
+
+                assert_eq!(c, b.clone() + a.clone());
+                assert_eq!(c, b.clone() + a.to_ref());
+                assert_eq!(c, b.clone() + &a);
+
+                assert_eq!(c, &b + a.clone());
+                assert_eq!(c, &b + a.to_ref());
+                assert_eq!(c, &b + &a);
+
+                let mut n = a.clone();
+                n += b.clone();
+                assert_eq!(c, n);
+
+                let mut n = a.clone();
+                n += &b;
+                assert_eq!(c, n);
+            }
+        };
+    }
+
+    impl_case!(case_1234en2_1234en3: "12.34" + "1.234" => "13.574");
+    impl_case!(case_1234en2_n1234en3: "12.34" + "-1.234" => "11.106");
+    impl_case!(case_1234en2_n1234en2: "12.34" + "-12.34" => "0");
+    impl_case!(case_1234e6_1234en6: "1234e6" + "1234e-6" => "1234000000.001234");
+    impl_case!(case_1234en6_1234e6: "1234e6" + "1234e-6" => "1234000000.001234");
+    impl_case!(case_18446744073709551616_1: "18446744073709551616.0" + "1" => "18446744073709551617");
+    impl_case!(case_184467440737e3380_0: "184467440737e3380" + "0" => "184467440737e3380");
+    impl_case!(case_0_776en1: "0" + "77.6" => "77.6");
+
+    impl_case!(case_80802295e5_int0: "80802295e5" + (int)"0" => "80802295e5");
+    impl_case!(case_239200en4_intneg101: "23.9200" + (int)"-101" => "-77.0800");
+    impl_case!(case_46636423395767125en15_int0: "46.636423395767125" + (int)"123" => "169.636423395767125");
+
+
+    #[cfg(property_tests)]
+    mod prop {
+        use super::*;
+        use proptest::*;
+        use num_traits::FromPrimitive;
+
+        proptest! {
+            #[test]
+            fn add_refs_and_owners(f: f32, g: f32) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+                prop_assume!(g.is_normal());
+
+                let a = BigDecimal::from_f32(f).unwrap();
+                let b = BigDecimal::from_f32(g).unwrap();
+                let own_plus_ref = a.clone() + &b;
+                let ref_plus_own = &a + b.clone();
+
+                let mut c = a.clone();
+                c += &b;
+
+                let mut d = a.clone();
+                d += b;
+
+                prop_assert_eq!(&own_plus_ref, &ref_plus_own);
+                prop_assert_eq!(&c, &ref_plus_own);
+                prop_assert_eq!(&d, &ref_plus_own);
+            }
+
+            #[test]
+            fn addition_is_communative(f: f32, g: f32) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+                prop_assume!(g.is_normal());
+
+                let a = BigDecimal::from_f32(f).unwrap();
+                let b = BigDecimal::from_f32(g).unwrap();
+                let a_plus_b = &a + &b;
+                let b_plus_a = &b + &a;
+
+                prop_assert_eq!(a_plus_b, b_plus_a)
+            }
+
+            #[test]
+            fn addition_is_associative(f: f32, g: f32, h: f32) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+                prop_assume!(g.is_normal());
+                prop_assume!(h.is_normal());
+
+                let a = BigDecimal::from_f32(f).unwrap();
+                let b = BigDecimal::from_f32(g).unwrap();
+                let c = BigDecimal::from_f32(h).unwrap();
+
+                let ab = &a + &b;
+                let ab_c = ab + &c;
+
+                let bc = &b + &c;
+                let a_bc = a + bc;
+
+                prop_assert_eq!(ab_c, a_bc)
+            }
+        }
+    }
+}
diff --git a/src/impl_ops_div.rs b/src/impl_ops_div.rs
new file mode 100644
index 0000000..9f36b9f
--- /dev/null
+++ b/src/impl_ops_div.rs
@@ -0,0 +1,168 @@
+//! Implement division
+
+use super::*;
+
+impl Div<BigDecimal> for BigDecimal {
+    type Output = BigDecimal;
+    #[inline]
+    fn div(self, other: BigDecimal) -> BigDecimal {
+        if other.is_zero() {
+            panic!("Division by zero");
+        }
+        if self.is_zero() || other.is_one() {
+            return self;
+        }
+
+        let scale = self.scale - other.scale;
+
+        if self.int_val == other.int_val {
+            return BigDecimal {
+                int_val: 1.into(),
+                scale: scale,
+            };
+        }
+
+        let max_precision = DEFAULT_PRECISION;
+
+        return impl_division(self.int_val, &other.int_val, scale, max_precision);
+    }
+}
+
+impl<'a> Div<&'a BigDecimal> for BigDecimal {
+    type Output = BigDecimal;
+    #[inline]
+    fn div(self, other: &'a BigDecimal) -> BigDecimal {
+        if other.is_zero() {
+            panic!("Division by zero");
+        }
+        if self.is_zero() || other.is_one() {
+            return self;
+        }
+
+        let scale = self.scale - other.scale;
+
+        if self.int_val == other.int_val {
+            return BigDecimal {
+                int_val: 1.into(),
+                scale: scale,
+            };
+        }
+
+        let max_precision = DEFAULT_PRECISION;
+
+        return impl_division(self.int_val, &other.int_val, scale, max_precision);
+    }
+}
+
+forward_ref_val_binop!(impl Div for BigDecimal, div);
+
+impl<'a, 'b> Div<&'b BigDecimal> for &'a BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn div(self, other: &BigDecimal) -> BigDecimal {
+        if other.is_zero() {
+            panic!("Division by zero");
+        }
+        // TODO: Fix setting scale
+        if self.is_zero() || other.is_one() {
+            return self.clone();
+        }
+
+        let scale = self.scale - other.scale;
+
+        let num_int = &self.int_val;
+        let den_int = &other.int_val;
+
+        if num_int == den_int {
+            return BigDecimal {
+                int_val: 1.into(),
+                scale: scale,
+            };
+        }
+
+        let max_precision = DEFAULT_PRECISION;
+
+        return impl_division(num_int.clone(), den_int, scale, max_precision);
+    }
+}
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    #[test]
+    fn test_div() {
+        let vals = vec![
+            ("0", "1", "0"),
+            ("0", "10", "0"),
+            ("2", "1", "2"),
+            ("2e1", "1", "2e1"),
+            ("10", "10", "1"),
+            ("100", "10.0", "1e1"),
+            ("20.0", "200", ".1"),
+            ("4", "2", "2.0"),
+            ("15", "3", "5.0"),
+            ("1", "2", "0.5"),
+            ("1", "2e-2", "5e1"),
+            ("1", "0.2", "5"),
+            ("1.0", "0.02", "50"),
+            ("1", "0.020", "5e1"),
+            ("5.0", "4.00", "1.25"),
+            ("5.0", "4.000", "1.25"),
+            ("5", "4.000", "1.25"),
+            ("5", "4", "125e-2"),
+            ("100", "5", "20"),
+            ("-50", "5", "-10"),
+            ("200", "-5", "-40."),
+            ("1", "3", ".3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333"),
+            ("-2", "-3", ".6666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667"),
+            ("-12.34", "1.233", "-10.00811030008110300081103000811030008110300081103000811030008110300081103000811030008110300081103001"),
+            ("125348", "352.2283", "355.8714617763535752237966114591019517738921035021887792661748076460636467881768727839301952739175132"),
+        ];
+
+        for &(x, y, z) in vals.iter() {
+
+            let a = BigDecimal::from_str(x).unwrap();
+            let b = BigDecimal::from_str(y).unwrap();
+            let c = BigDecimal::from_str(z).unwrap();
+
+            assert_eq!(a.clone() / b.clone(), c);
+            assert_eq!(a.clone() / &b, c);
+            assert_eq!(&a / b.clone(), c);
+            assert_eq!(&a / &b, c);
+            // assert_eq!(q.scale, c.scale);
+
+            // let mut q = a;
+            // q /= b;
+            // assert_eq!(q, c);
+        }
+    }
+
+    #[test]
+    #[should_panic(expected = "Division by zero")]
+    fn test_division_by_zero_panics() {
+        let x = BigDecimal::from_str("3.14").unwrap();
+        let _r = x / 0;
+    }
+
+    #[test]
+    #[should_panic(expected = "Division by zero")]
+    fn test_division_by_zero_panics_v2() {
+        let x = BigDecimal::from_str("3.14").unwrap();
+        let _r = x / BigDecimal::zero();
+    }
+
+    #[test]
+    fn test_division_by_large_number() {
+        let n = 1u8;
+        let d: BigDecimal = "79437738588056219546528239237352667078".parse().unwrap();
+
+        let quotient_n_ref_d = n / &d;
+        let quotient_n_d = n / d.clone();
+        assert_eq!(quotient_n_ref_d, quotient_n_d);
+
+        assert_eq!(quotient_n_ref_d, "1.258847517281104957975270408416632052090243053529147458917576143852500316808428812104171430669001064E-38".parse().unwrap());
+    }
+}
diff --git a/src/impl_ops_mul.rs b/src/impl_ops_mul.rs
new file mode 100644
index 0000000..2b1df37
--- /dev/null
+++ b/src/impl_ops_mul.rs
@@ -0,0 +1,284 @@
+//! Multiplication operator trait implementation
+//!
+
+use super::*;
+use crate::stdlib::mem::swap;
+
+
+impl Mul<BigDecimal> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(mut self, rhs: BigDecimal) -> BigDecimal {
+        if self.is_one() {
+            rhs
+        } else if rhs.is_one() {
+            self
+        } else {
+            self.scale += rhs.scale;
+            self.int_val *= rhs.int_val;
+            self
+        }
+    }
+}
+
+impl<'a> Mul<&'a BigDecimal> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(mut self, rhs: &'a BigDecimal) -> BigDecimal {
+        if self.is_one() {
+            self.scale = rhs.scale;
+            self.int_val.set_zero();
+            self.int_val.add_assign(&rhs.int_val);
+            self
+        } else if rhs.is_one() {
+            self
+        } else {
+            self.scale += rhs.scale;
+            MulAssign::mul_assign(&mut self.int_val, &rhs.int_val);
+            self
+        }
+    }
+}
+
+impl<'a> Mul<BigDecimal> for &'a BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(self, rhs: BigDecimal) -> BigDecimal {
+        rhs * self
+    }
+}
+
+impl<'a, 'b> Mul<&'b BigDecimal> for &'a BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(self, rhs: &BigDecimal) -> BigDecimal {
+        if self.is_one() {
+            rhs.normalized()
+        } else if rhs.is_one() {
+            self.normalized()
+        } else {
+            let scale = self.scale + rhs.scale;
+            BigDecimal::new(&self.int_val * &rhs.int_val, scale)
+        }
+    }
+}
+
+impl Mul<BigInt> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(mut self, rhs: BigInt) -> BigDecimal {
+        self.int_val *= rhs;
+        self
+    }
+}
+
+impl<'a> Mul<&'a BigInt> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(mut self, rhs: &BigInt) -> BigDecimal {
+        self.int_val *= rhs;
+        self
+    }
+}
+
+impl<'a> Mul<BigInt> for &'a BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(self, mut rhs: BigInt) -> BigDecimal {
+        rhs *= &self.int_val;
+        BigDecimal::new(rhs, self.scale)
+    }
+}
+
+impl<'a, 'b> Mul<&'a BigInt> for &'b BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(self, rhs: &BigInt) -> BigDecimal {
+        if rhs.is_one() {
+            self.normalized()
+        } else if self.is_one() {
+            BigDecimal::new(rhs.clone(), 0)
+        } else {
+            let value = &self.int_val * rhs;
+            BigDecimal::new(value, self.scale)
+        }
+    }
+}
+
+impl Mul<BigDecimal> for BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(mut self, mut rhs: BigDecimal) -> BigDecimal {
+        if rhs.is_one() {
+            rhs.scale = 0;
+            swap(&mut rhs.int_val, &mut self);
+        } else if !self.is_one() {
+            rhs.int_val *= self;
+        }
+        rhs
+    }
+}
+
+impl<'a> Mul<BigDecimal> for &'a BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(self, mut rhs: BigDecimal) -> BigDecimal {
+        if self.is_one() {
+            rhs.normalized()
+        } else if rhs.is_one() {
+            rhs.int_val.set_zero();
+            rhs.int_val += self;
+            rhs.scale = 0;
+            rhs
+        } else {
+            rhs.int_val *= self;
+            rhs
+        }
+    }
+}
+
+impl<'a, 'b> Mul<&'a BigDecimal> for &'b BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(self, rhs: &BigDecimal) -> BigDecimal {
+        if self.is_one() {
+            rhs.normalized()
+        } else if rhs.is_one() {
+            BigDecimal::new(self.clone(), 0)
+        } else {
+            let value = &rhs.int_val * self;
+            BigDecimal::new(value, rhs.scale)
+        }
+    }
+}
+
+impl<'a> Mul<&'a BigDecimal> for BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn mul(mut self, rhs: &BigDecimal) -> BigDecimal {
+        if self.is_one() {
+            rhs.normalized()
+        } else if rhs.is_one() {
+            BigDecimal::new(self, 0)
+        } else {
+            self *= &rhs.int_val;
+            BigDecimal::new(self, rhs.scale)
+        }
+    }
+}
+
+forward_val_assignop!(impl MulAssign for BigDecimal, mul_assign);
+
+impl<'a> MulAssign<&'a BigDecimal> for BigDecimal {
+    #[inline]
+    fn mul_assign(&mut self, rhs: &BigDecimal) {
+        if rhs.is_one() {
+            return;
+        }
+        self.scale += rhs.scale;
+        self.int_val = &self.int_val * &rhs.int_val;
+    }
+}
+
+impl<'a> MulAssign<&'a BigInt> for BigDecimal {
+    #[inline]
+    fn mul_assign(&mut self, rhs: &BigInt) {
+        if rhs.is_one() {
+            return;
+        }
+        self.int_val *= rhs;
+    }
+}
+
+impl MulAssign<BigInt> for BigDecimal {
+    #[inline]
+    fn mul_assign(&mut self, rhs: BigInt) {
+        *self *= &rhs
+    }
+}
+
+
+#[cfg(test)]
+#[allow(non_snake_case)]
+mod bigdecimal_tests {
+    use crate::{stdlib, BigDecimal, ToString, FromStr, TryFrom};
+    use num_traits::{ToPrimitive, FromPrimitive, Signed, Zero, One};
+    use num_bigint;
+    use paste::paste;
+
+    /// Test multiplication of two bigdecimals
+    #[test]
+    fn test_mul() {
+
+        let vals = vec![
+            ("2", "1", "2"),
+            ("12.34", "1.234", "15.22756"),
+            ("2e1", "1", "20"),
+            ("3", ".333333", "0.999999"),
+            ("2389472934723", "209481029831", "500549251119075878721813"),
+            ("1e-450", "1e500", ".1e51"),
+            ("-995052931372975485719.533153137", "4.523087321", "-4500711297616988541501.836966993116075977"),
+            ("995052931372975485719.533153137", "-4.523087321", "-4500711297616988541501.836966993116075977"),
+            ("-8.37664968", "-1.9086963714056968482094712882596748", "15.988480848752691653730876239769592670324064"),
+            ("-8.37664968", "0", "0"),
+        ];
+
+        for &(x, y, z) in vals.iter() {
+
+            let mut a = BigDecimal::from_str(x).unwrap();
+            let b = BigDecimal::from_str(y).unwrap();
+            let c = BigDecimal::from_str(z).unwrap();
+
+            assert_eq!(a.clone() * b.clone(), c);
+            assert_eq!(a.clone() * &b, c);
+            assert_eq!(&a * b.clone(), c);
+            assert_eq!(&a * &b, c);
+
+            a *= b;
+            assert_eq!(a, c);
+        }
+    }
+
+    /// Test multiplication between big decimal and big integer
+    #[test]
+    fn test_mul_bigint() {
+        let vals = vec![
+            ("2", "1", "2"),
+            ("8.561", "10", "85.61"),
+            ("1.0000", "638655273892892437", "638655273892892437"),
+            ("10000", "638655273892892437", "6386552738928924370000"),
+            (".0005", "368408638655273892892437473", "184204319327636946446218.7365"),
+            ("9e-1", "368408638655273892892437473", "331567774789746503603193725.7"),
+            ("-1.175470587012343730098", "577575785", "-678923347.038065234601180476930"),
+            ("-1.175470587012343730098", "-76527768352678", "89956140788267.069799533723307502444"),
+            ("-1.175470587012343730098", "0", "0"),
+        ];
+
+        for &(x, y, z) in vals.iter() {
+            let a = BigDecimal::from_str(x).unwrap();
+            let b = num_bigint::BigInt::from_str(y).unwrap();
+            let c = BigDecimal::from_str(z).unwrap();
+
+            assert_eq!(a.clone() * b.clone(), c);
+            assert_eq!(b.clone() * a.clone(), c);
+            assert_eq!(a.clone() * &b, c);
+            assert_eq!(b.clone() * &a, c);
+            assert_eq!(&a * b.clone(), c);
+            assert_eq!(&b * a.clone(), c);
+            assert_eq!(&a * &b, c);
+            assert_eq!(&b * &a, c);
+        }
+    }
+}
diff --git a/src/impl_ops_rem.rs b/src/impl_ops_rem.rs
new file mode 100644
index 0000000..56bad43
--- /dev/null
+++ b/src/impl_ops_rem.rs
@@ -0,0 +1,167 @@
+//! Remainder implementations
+
+use super::*;
+
+
+impl Rem<BigDecimal> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn rem(self, other: BigDecimal) -> BigDecimal {
+        let scale = cmp::max(self.scale, other.scale);
+
+        let num = self.take_and_scale(scale).int_val;
+        let den = other.take_and_scale(scale).int_val;
+
+        BigDecimal::new(num % den, scale)
+    }
+}
+
+impl<'a> Rem<&'a BigDecimal> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn rem(self, other: &BigDecimal) -> BigDecimal {
+        let scale = cmp::max(self.scale, other.scale);
+        let num = self.take_and_scale(scale).int_val;
+        let den = &other.int_val;
+
+        let result = if scale == other.scale {
+            num % den
+        } else {
+            num % (den * ten_to_the((scale - other.scale) as u64))
+        };
+        BigDecimal::new(result, scale)
+    }
+}
+
+impl<'a> Rem<BigDecimal> for &'a BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn rem(self, other: BigDecimal) -> BigDecimal {
+        let scale = cmp::max(self.scale, other.scale);
+        let num = &self.int_val;
+        let den = other.take_and_scale(scale).int_val;
+
+        let result = if scale == self.scale {
+            num % den
+        } else {
+            let scaled_num = num * ten_to_the((scale - self.scale) as u64);
+            scaled_num % den
+        };
+
+        BigDecimal::new(result, scale)
+    }
+}
+
+impl<'a, 'b> Rem<&'b BigDecimal> for &'a BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn rem(self, other: &BigDecimal) -> BigDecimal {
+        let scale = cmp::max(self.scale, other.scale);
+        let num = &self.int_val;
+        let den = &other.int_val;
+
+        let result = match self.scale.cmp(&other.scale) {
+            Ordering::Equal => num % den,
+            Ordering::Less => {
+                let scaled_num = num * ten_to_the((scale - self.scale) as u64);
+                scaled_num % den
+            }
+            Ordering::Greater => {
+                let scaled_den = den * ten_to_the((scale - other.scale) as u64);
+                num % scaled_den
+            }
+        };
+        BigDecimal::new(result, scale)
+    }
+}
+
+impl RemAssign<&BigDecimal> for BigDecimal {
+    fn rem_assign(&mut self, other: &BigDecimal) {
+        let rem = (&*self).rem(other);
+        *self = rem;
+    }
+}
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use paste::paste;
+
+    macro_rules! impl_case {
+        ($a:literal % $b:literal => $c:literal ) => {
+            paste! {
+                impl_case!([< case_ $a _ $b >]: $a % $b => $c);
+            }
+        };
+        ($name:ident: $a:literal % $b:literal => $c:literal ) => {
+            #[test]
+            fn $name() {
+                let mut a: BigDecimal = $a.parse().unwrap();
+                let b: BigDecimal = $b.parse().unwrap();
+                let c: BigDecimal = $c.parse().unwrap();
+
+                assert_eq!(a.clone() % b.clone(), c);
+
+                assert_eq!(a.clone() % &b, c);
+                assert_eq!(&a % b.clone(), c);
+                assert_eq!(&a % &b, c);
+
+                a %= &b;
+                assert_eq!(a, c);
+            }
+        };
+    }
+
+    impl_case!("100" % "5" => "0");
+    impl_case!("2e1" % "1" => "0");
+    impl_case!("2" % "1" => "0");
+    impl_case!("1" % "3" => "1");
+    impl_case!("1" % "5e-1" => "0");
+    impl_case!("15e-1" % "1" => "0.5");
+    impl_case!("1" % "3e-2" => "1e-2");
+    impl_case!("10" % "3e-3" => "0.001");
+    impl_case!("3" % "2" => "1");
+    impl_case!("1234e-2" % "1233e-3" => "0.01");
+
+    impl_case!(case_neg3_2: "-3" % "2" => "-1");
+    impl_case!(case_3_neg2: "3" % "-2" => "1");
+    impl_case!(case_neg3_neg2: "3" % "-2" => "1");
+
+    impl_case!(case_neg95eneg1_515eneg2: "-9.5" % "5.15" => "-4.35");
+
+
+    #[cfg(property_tests)]
+    mod prop {
+        use super::*;
+        use proptest::*;
+        use num_traits::FromPrimitive;
+
+        proptest! {
+            #[test]
+            fn quotient_and_remainder(f: f32, g: f32) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+                prop_assume!(g.is_normal());
+                prop_assume!(!g.is_zero());
+
+                let (f, g) = if f.abs() > g.abs() {
+                    (f, g)
+                } else {
+                    (g, f)
+                };
+
+                let a = BigDecimal::from_f32(f).unwrap();
+                let b = BigDecimal::from_f32(g).unwrap();
+
+                let r = &a % &b;
+                let q = (&a / &b).with_scale(0);
+                assert_eq!(a, q * b + r);
+            }
+        }
+    }
+}
diff --git a/src/impl_ops_sub.rs b/src/impl_ops_sub.rs
new file mode 100644
index 0000000..5f5b305
--- /dev/null
+++ b/src/impl_ops_sub.rs
@@ -0,0 +1,345 @@
+//!
+//! Subtraction operator trait implementation
+//!
+
+use crate::*;
+
+
+impl Sub<BigDecimal> for BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigDecimal) -> BigDecimal {
+        if rhs.is_zero() {
+            return self;
+        }
+
+        if self.is_zero() {
+            return rhs.neg();
+        }
+
+        let mut lhs = self;
+        match lhs.scale.cmp(&rhs.scale) {
+            Ordering::Equal => {
+                lhs.int_val -= rhs.int_val;
+                lhs
+            }
+            Ordering::Less => {
+                lhs.take_and_scale(rhs.scale) - rhs
+            }
+            Ordering::Greater => {
+                let rhs = rhs.take_and_scale(lhs.scale);
+                lhs - rhs
+            },
+        }
+    }
+}
+
+impl Sub<BigDecimal> for &'_ BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigDecimal) -> BigDecimal {
+        self.to_ref() - rhs
+    }
+}
+
+impl Sub<BigDecimal> for BigDecimalRef<'_> {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigDecimal) -> BigDecimal {
+        (rhs - self).neg()
+    }
+}
+
+impl<'a, T: Into<BigDecimalRef<'a>>> Sub<T> for BigDecimal {
+    type Output = BigDecimal;
+
+    fn sub(mut self, rhs: T) -> BigDecimal {
+        self.sub_assign(rhs);
+        self
+    }
+}
+
+impl<'a, T: Into<BigDecimalRef<'a>>> Sub<T> for &'_ BigDecimal {
+    type Output = BigDecimal;
+
+    fn sub(self, rhs: T) -> BigDecimal {
+        let rhs = rhs.into();
+
+        match self.scale.cmp(&rhs.scale) {
+            Ordering::Equal => {
+                self.clone() - rhs
+            }
+            Ordering::Less => {
+                self.with_scale(rhs.scale) - rhs
+            }
+            Ordering::Greater => {
+                self - rhs.to_owned_with_scale(self.scale)
+            }
+        }
+    }
+}
+
+impl<'a, T: Into<BigDecimalRef<'a>>> Sub<T> for BigDecimalRef<'_> {
+    type Output = BigDecimal;
+
+    fn sub(self, rhs: T) -> BigDecimal {
+        let rhs = rhs.into();
+
+        match self.scale.cmp(&rhs.scale) {
+            Ordering::Equal => self.to_owned() - rhs,
+            Ordering::Less => self.to_owned_with_scale(rhs.scale) - rhs,
+            Ordering::Greater => self - rhs.to_owned_with_scale(self.scale),
+        }
+    }
+}
+
+impl Sub<BigInt> for BigDecimal {
+    type Output = BigDecimal;
+
+    fn sub(mut self, rhs: BigInt) -> BigDecimal {
+        self.sub_assign(rhs);
+        self
+    }
+}
+
+
+impl Sub<BigInt> for &'_ BigDecimal {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigInt) -> BigDecimal {
+        self.to_ref() - rhs
+    }
+}
+
+impl Sub<BigInt> for BigDecimalRef<'_> {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigInt) -> BigDecimal {
+        self - BigDecimal::from(rhs)
+    }
+}
+
+impl Sub<BigDecimal> for BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigDecimal) -> BigDecimal {
+        (rhs - self).neg()
+    }
+}
+
+impl Sub<BigDecimal> for &BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigDecimal) -> BigDecimal {
+        (rhs - self).neg()
+    }
+}
+
+impl<'a> Sub<BigDecimalRef<'a>> for BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigDecimalRef<'a>) -> BigDecimal {
+        (rhs - &self).neg()
+    }
+}
+
+
+impl<'a> Sub<BigDecimalRef<'a>> for &BigInt {
+    type Output = BigDecimal;
+
+    #[inline]
+    fn sub(self, rhs: BigDecimalRef<'a>) -> BigDecimal {
+        (rhs - self).neg()
+    }
+}
+
+
+impl SubAssign<BigDecimal> for BigDecimal {
+    #[inline]
+    fn sub_assign(&mut self, rhs: BigDecimal) {
+        if rhs.is_zero() {
+            return;
+        }
+        if self.is_zero() {
+            *self = rhs.neg();
+            return;
+        }
+        match self.scale.cmp(&rhs.scale) {
+            Ordering::Equal => {
+                self.int_val -= rhs.int_val;
+            }
+            Ordering::Less => {
+                self.int_val *= ten_to_the((rhs.scale - self.scale) as u64);
+                self.int_val -= rhs.int_val;
+                self.scale = rhs.scale;
+            }
+            Ordering::Greater => {
+                let mut rhs_int_val = rhs.int_val;
+                rhs_int_val *= ten_to_the((self.scale - rhs.scale) as u64);
+                self.int_val -= rhs_int_val;
+            }
+        }
+    }
+}
+
+impl<'rhs, T: Into<BigDecimalRef<'rhs>>> SubAssign<T> for BigDecimal {
+    #[inline]
+    fn sub_assign(&mut self, rhs: T) {
+        let rhs = rhs.into();
+        if rhs.is_zero() {
+            return;
+        }
+        if self.is_zero() {
+            *self = rhs.neg().to_owned();
+            return;
+        }
+
+        match self.scale.cmp(&rhs.scale) {
+            Ordering::Equal => {
+                self.int_val -= rhs.to_owned().int_val;
+            }
+            Ordering::Less => {
+                self.int_val *= ten_to_the((rhs.scale - self.scale) as u64);
+                self.int_val -= rhs.to_owned().int_val;
+                self.scale = rhs.scale;
+            }
+            Ordering::Greater => {
+                *self -= rhs.to_owned_with_scale(self.scale);
+            }
+        }
+    }
+}
+
+impl SubAssign<BigInt> for BigDecimal {
+    #[inline(always)]
+    fn sub_assign(&mut self, rhs: BigInt) {
+        *self -= BigDecimal::new(rhs, 0)
+    }
+}
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use paste::paste;
+
+    macro_rules! impl_case {
+        ($name:ident: $a:literal - $b:literal => $c:literal ) => {
+            #[test]
+            fn $name() {
+                let a: BigDecimal = $a.parse().unwrap();
+                let b: BigDecimal = $b.parse().unwrap();
+                let c: BigDecimal = $c.parse().unwrap();
+
+                assert_eq!(c, a.clone() - b.clone());
+
+                assert_eq!(c, a.clone() - &b);
+                assert_eq!(c, &a - b.clone());
+                assert_eq!(c, &a - &b);
+
+                assert_eq!(c, a.to_ref() - &b);
+                assert_eq!(c, &a - b.to_ref());
+                assert_eq!(c, a.to_ref() - b.to_ref());
+
+                let mut n = a.clone();
+                n -= b.to_ref();
+                assert_eq!(n, c);
+
+                let mut n = a.clone();
+                n -= &b;
+                assert_eq!(n, c);
+
+                let mut n = a.clone();
+                n -= b.clone();
+                assert_eq!(n, c);
+
+                let mut n = a.clone();
+                (&mut n).sub_assign(b.clone());
+                assert_eq!(n, c);
+            }
+        };
+        ($name:ident: $a:literal - (int) $b:literal => $c:literal ) => {
+            #[test]
+            fn $name() {
+                let a: BigDecimal = $a.parse().unwrap();
+                let b: BigInt = $b.parse().unwrap();
+                let expected: BigDecimal = $c.parse().unwrap();
+
+                assert_eq!(expected, a.clone() - b.clone());
+                assert_eq!(expected, a.clone() - &b);
+                assert_eq!(expected, &a - &b);
+                assert_eq!(expected, &a - b.clone());
+                assert_eq!(expected, a.to_ref() - &b);
+
+                let expected_neg = expected.clone().neg();
+                assert_eq!(expected_neg, b.clone() - a.clone());
+                assert_eq!(expected_neg, &b - a.to_ref());
+                assert_eq!(expected_neg, &b - a.clone());
+            }
+        };
+    }
+
+    impl_case!(case_1234en2_1234en3: "12.34" - "1.234" => "11.106");
+    impl_case!(case_1234en2_n1234en3: "12.34" - "-1.234" => "13.574");
+    impl_case!(case_1234e6_1234en6: "1234e6" - "1234e-6" => "1233999999.998766");
+    impl_case!(case_1234en6_1234e6: "1234e-6" - "1234e6" => "-1233999999.998766");
+    impl_case!(case_712911676en6_4856259269250829: "712911676e-6" - "4856259269250829" => "-4856259269250116.088324");
+    impl_case!(case_85616001e4_0: "85616001e4" - "0" => "85616001e4");
+    impl_case!(case_0_520707672en5: "0" - "5207.07672" => "-520707672e-5");
+    impl_case!(case_99291289e5_int0: "99291289e5" - (int)"0" => "99291289e5");
+    impl_case!(case_7051277471570131en16_int1: "0.7051277471570131" - (int)"1" => "-0.2948722528429869");
+    impl_case!(case_4068603022763836en8_intneg10: "40686030.22763836" - (int)"-10" => "40686040.22763836");
+
+    #[cfg(property_tests)]
+    mod prop {
+        use super::*;
+        use proptest::*;
+        use num_traits::FromPrimitive;
+
+        proptest! {
+            #[test]
+            fn sub_refs_and_owners(f: f32, g: f32) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+                prop_assume!(g.is_normal());
+
+                let a = BigDecimal::from_f32(f).unwrap();
+                let b = BigDecimal::from_f32(g).unwrap();
+                let own_minus_ref = a.clone() - &b;
+                let ref_minus_own = &a - b.clone();
+
+                let mut c = a.clone();
+                c -= &b;
+
+                let mut d = a.clone();
+                d -= b;
+
+                prop_assert_eq!(&own_minus_ref, &ref_minus_own);
+                prop_assert_eq!(&c, &ref_minus_own);
+                prop_assert_eq!(&d, &ref_minus_own);
+            }
+
+            #[test]
+            fn subtraction_is_anticommunative(f: f32, g: f32) {
+                // ignore non-normal numbers
+                prop_assume!(f.is_normal());
+                prop_assume!(g.is_normal());
+
+                let a = BigDecimal::from_f32(f).unwrap();
+                let b = BigDecimal::from_f32(g).unwrap();
+                let a_minus_b = &a - &b;
+                let b_minus_a = &b - &a;
+
+                prop_assert_eq!(a_minus_b, -b_minus_a)
+            }
+        }
+    }
+}
diff --git a/src/lib.rs b/src/lib.rs
index a77801b..e461a6a 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -46,12 +46,16 @@
 #![allow(clippy::suspicious_arithmetic_impl)]
 #![allow(clippy::suspicious_op_assign_impl)]
 #![allow(clippy::redundant_field_names)]
+#![allow(unused_imports)]
 
 
 pub extern crate num_bigint;
 pub extern crate num_traits;
 extern crate num_integer;
 
+#[cfg(test)]
+extern crate paste;
+
 #[cfg(feature = "serde")]
 extern crate serde;
 
@@ -67,18 +71,19 @@ use self::stdlib::convert::TryFrom;
 use self::stdlib::default::Default;
 use self::stdlib::hash::{Hash, Hasher};
 use self::stdlib::num::{ParseFloatError, ParseIntError};
-use self::stdlib::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Rem, Sub, SubAssign};
+use self::stdlib::ops::{
+    Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign, Rem, RemAssign,
+};
 use self::stdlib::iter::Sum;
 use self::stdlib::str::FromStr;
 use self::stdlib::string::{String, ToString};
 use self::stdlib::fmt;
 
-use num_bigint::{BigInt, ParseBigIntError, Sign};
+use num_bigint::{BigInt, BigUint, ParseBigIntError, Sign};
 use num_integer::Integer as IntegerTrait;
 pub use num_traits::{FromPrimitive, Num, One, Signed, ToPrimitive, Zero};
 
-#[allow(clippy::approx_constant)] // requires Rust 1.43.0
-const LOG2_10: f64 = 3.321928094887362_f64;
+use stdlib::f64::consts::LOG2_10;
 
 
 // const DEFAULT_PRECISION: u64 = ${RUST_BIGDECIMAL_DEFAULT_PRECISION} or 100;
@@ -87,57 +92,43 @@ include!(concat!(env!("OUT_DIR"), "/default_precision.rs"));
 #[macro_use]
 mod macros;
 
-#[cfg(test)]
-extern crate paste;
+// "low level" functions
+mod arithmetic;
 
 // From<T>, To<T>, TryFrom<T> impls
 mod impl_convert;
+
 // Add<T>, Sub<T>, etc...
 mod impl_ops;
+mod impl_ops_add;
+mod impl_ops_sub;
+mod impl_ops_mul;
+mod impl_ops_div;
+mod impl_ops_rem;
+
+// PartialEq
+mod impl_cmp;
 
 // Implementations of num_traits
 mod impl_num;
 
+// construct BigDecimals from strings and floats
 mod parsing;
+
+// Routines for rounding
 pub mod rounding;
 pub use rounding::RoundingMode;
 
-#[inline(always)]
-fn ten_to_the(pow: u64) -> BigInt {
-    if pow < 20 {
-        BigInt::from(10u64.pow(pow as u32))
-    } else {
-        let (half, rem) = pow.div_rem(&16);
-
-        let mut x = ten_to_the(half);
-
-        for _ in 0..4 {
-            x = &x * &x;
-        }
+// Mathematical context
+mod context;
+pub use context::Context;
 
-        if rem == 0 {
-            x
-        } else {
-            x * ten_to_the(rem)
-        }
-    }
-}
-
-#[inline(always)]
-fn count_decimal_digits(int: &BigInt) -> u64 {
-    if int.is_zero() {
-        return 1;
-    }
-    let uint = int.magnitude();
-    let mut digits = (uint.bits() as f64 / LOG2_10) as u64;
-    // guess number of digits based on number of bits in UInt
-    let mut num = ten_to_the(digits).to_biguint().expect("Ten to power is negative");
-    while *uint >= num {
-        num *= 10u8;
-        digits += 1;
-    }
-    digits
-}
+use arithmetic::{
+    ten_to_the,
+    ten_to_the_uint,
+    count_decimal_digits,
+    count_decimal_digits_uint,
+};
 
 /// Internal function used for rounding
 ///
@@ -204,6 +195,37 @@ impl BigDecimal {
         }
     }
 
+    /// Make a BigDecimalRef of this value
+    pub fn to_ref(&self) -> BigDecimalRef<'_> {
+        // search for "From<&'a BigDecimal> for BigDecimalRef<'a>"
+        self.into()
+    }
+
+    /// Returns the scale of the BigDecimal, the total number of
+    /// digits to the right of the decimal point (including insignificant
+    /// leading zeros)
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use bigdecimal::BigDecimal;
+    /// use std::str::FromStr;
+    ///
+    /// let a = BigDecimal::from(12345);  // No fractional part
+    /// let b = BigDecimal::from_str("123.45").unwrap();  // Fractional part
+    /// let c = BigDecimal::from_str("0.0000012345").unwrap();  // Completely fractional part
+    /// let d = BigDecimal::from_str("5e9").unwrap();  // Negative-fractional part
+    ///
+    /// assert_eq!(a.fractional_digit_count(), 0);
+    /// assert_eq!(b.fractional_digit_count(), 2);
+    /// assert_eq!(c.fractional_digit_count(), 10);
+    /// assert_eq!(d.fractional_digit_count(), -9);
+    /// ```
+    #[inline]
+    pub fn fractional_digit_count(&self) -> i64 {
+        self.scale
+    }
+
     /// Creates and initializes a `BigDecimal`.
     ///
     /// Decodes using `str::from_utf8` and forwards to `BigDecimal::from_str_radix`.
@@ -359,6 +381,21 @@ impl BigDecimal {
         }
     }
 
+    /// Take and return bigdecimal with the given sign
+    ///
+    /// The Sign value `NoSign` is ignored: only use Plus & Minus
+    ///
+    pub(crate) fn take_with_sign(self, sign: Sign) -> BigDecimal {
+        let BigDecimal { scale, mut int_val } = self;
+        if int_val.sign() != sign && sign != Sign::NoSign {
+            int_val = int_val.neg();
+        }
+        BigDecimal {
+            int_val: int_val,
+            scale: scale,
+        }
+    }
+
     /// Return a new BigDecimal object with precision set to new value
     ///
     /// ```
@@ -404,6 +441,32 @@ impl BigDecimal {
         }
     }
 
+    /// Return this BigDecimal with the given precision, rounding if needed
+    #[cfg(rustc_1_46)]  // Option::zip
+    pub fn with_precision_round(&self, prec: stdlib::num::NonZeroU64, round: RoundingMode) -> BigDecimal {
+        let digit_count = self.digits();
+        let new_prec = prec.get().to_i64();
+        let new_scale = new_prec
+                        .zip(digit_count.to_i64())
+                        .and_then(|(new_prec, old_prec)| new_prec.checked_sub(old_prec))
+                        .and_then(|prec_diff| self.scale.checked_add(prec_diff))
+                        .expect("precision overflow");
+
+        self.with_scale_round(new_scale, round)
+    }
+
+    #[cfg(not(rustc_1_46))]
+    pub fn with_precision_round(&self, prec: stdlib::num::NonZeroU64, round: RoundingMode) -> BigDecimal {
+        let new_scale = self.digits().to_i64().and_then(
+                            |old_prec| {
+                                prec.get().to_i64().and_then(
+                                    |new_prec| { new_prec.checked_sub(old_prec) })})
+                            .and_then(|prec_diff| self.scale.checked_add(prec_diff))
+                            .expect("precision overflow");
+
+        self.with_scale_round(new_scale, round)
+    }
+
     /// Return the sign of the `BigDecimal` as `num::bigint::Sign`.
     ///
     /// ```
@@ -598,6 +661,12 @@ impl BigDecimal {
     /// ```
     #[inline]
     pub fn sqrt(&self) -> Option<BigDecimal> {
+        self.sqrt_with_context(&Context::default())
+    }
+
+    /// Take the square root of the number, using context for precision and rounding
+    ///
+    pub fn sqrt_with_context(&self, ctx: &Context) -> Option<BigDecimal> {
         if self.is_zero() || self.is_one() {
             return Some(self.clone());
         }
@@ -605,194 +674,49 @@ impl BigDecimal {
             return None;
         }
 
-        // make guess
-        let guess = {
-            let magic_guess_scale = 1.1951678538495576_f64;
-            let initial_guess = (self.int_val.bits() as f64 - self.scale as f64 * LOG2_10) / 2.0;
-            let res = magic_guess_scale * exp2(initial_guess);
-
-            if res.is_normal() {
-                BigDecimal::try_from(res).unwrap()
-            } else {
-                // can't guess with float - just guess magnitude
-                let scale = (self.int_val.bits() as f64 / -LOG2_10 + self.scale as f64).round() as i64;
-                BigDecimal::new(BigInt::from(1), scale / 2)
-            }
-        };
+        let uint = self.int_val.magnitude();
+        let result = arithmetic::sqrt::impl_sqrt(uint, self.scale, ctx);
 
-        // // wikipedia example - use for testing the algorithm
-        // if self == &BigDecimal::from_str("125348").unwrap() {
-        //     running_result = BigDecimal::from(600)
-        // }
-
-        // TODO: Use context variable to set precision
-        let max_precision = DEFAULT_PRECISION;
-
-        let next_iteration = move |r: BigDecimal| {
-            // division needs to be precise to (at least) one extra digit
-            let tmp = impl_division(
-                self.int_val.clone(),
-                &r.int_val,
-                self.scale - r.scale,
-                max_precision + 1,
-            );
-
-            // half will increase precision on each iteration
-            (tmp + r).half()
-        };
-
-        // calculate first iteration
-        let mut running_result = next_iteration(guess);
-
-        let mut prev_result = BigDecimal::one();
-        let mut result = BigDecimal::zero();
-
-        // TODO: Prove that we don't need to arbitrarily limit iterations
-        // and that convergence can be calculated
-        while prev_result != result {
-            // store current result to test for convergence
-            prev_result = result;
-
-            // calculate next iteration
-            running_result = next_iteration(running_result);
-
-            // 'result' has clipped precision, 'running_result' has full precision
-            result = if running_result.digits() > max_precision {
-                running_result.with_prec(max_precision)
-            } else {
-                running_result.clone()
-            };
-        }
-
-        return Some(result);
+        Some(result)
     }
 
-    /// Take the cube root of the number
+    /// Take the cube root of the number, using default context
     ///
     #[inline]
     pub fn cbrt(&self) -> BigDecimal {
+        self.cbrt_with_context(&Context::default())
+    }
+
+    /// Take cube root of self, using properties of context
+    pub fn cbrt_with_context(&self, ctx: &Context) -> BigDecimal {
         if self.is_zero() || self.is_one() {
             return self.clone();
         }
-        if self.is_negative() {
-            return -self.abs().cbrt();
-        }
-
-        // make guess
-        let guess = {
-            let magic_guess_scale = 1.124960491619939_f64;
-            let initial_guess = (self.int_val.bits() as f64 - self.scale as f64 * LOG2_10) / 3.0;
-            let res = magic_guess_scale * exp2(initial_guess);
-
-            if res.is_normal() {
-                BigDecimal::try_from(res).unwrap()
-            } else {
-                // can't guess with float - just guess magnitude
-                let scale = (self.int_val.bits() as f64 / LOG2_10 - self.scale as f64).round() as i64;
-                BigDecimal::new(BigInt::from(1), -scale / 3)
-            }
-        };
-
-        // TODO: Use context variable to set precision
-        let max_precision = DEFAULT_PRECISION;
-
-        let three = BigDecimal::from(3);
-
-        let next_iteration = move |r: BigDecimal| {
-            let sqrd = r.square();
-            let tmp = impl_division(
-                self.int_val.clone(),
-                &sqrd.int_val,
-                self.scale - sqrd.scale,
-                max_precision + 1,
-            );
-            let tmp = tmp + r.double();
-            impl_division(tmp.int_val, &three.int_val, tmp.scale - three.scale, max_precision + 1)
-        };
-
-        // result initial
-        let mut running_result = next_iteration(guess);
-
-        let mut prev_result = BigDecimal::one();
-        let mut result = BigDecimal::zero();
 
-        // TODO: Prove that we don't need to arbitrarily limit iterations
-        // and that convergence can be calculated
-        while prev_result != result {
-            // store current result to test for convergence
-            prev_result = result;
+        let uint = self.int_val.magnitude();
+        let result = arithmetic::cbrt::impl_cbrt_uint_scale(uint, self.scale, ctx);
 
-            running_result = next_iteration(running_result);
-
-            // result has clipped precision, running_result has full precision
-            result = if running_result.digits() > max_precision {
-                running_result.with_prec(max_precision)
-            } else {
-                running_result.clone()
-            };
-        }
-
-        return result;
+        // always copy sign
+        result.take_with_sign(self.sign())
     }
 
     /// Compute the reciprical of the number: x<sup>-1</sup>
     #[inline]
     pub fn inverse(&self) -> BigDecimal {
+        self.inverse_with_context(&Context::default())
+    }
+
+    /// Return inverse of self, rounding with ctx
+    pub fn inverse_with_context(&self, ctx: &Context) -> BigDecimal {
         if self.is_zero() || self.is_one() {
             return self.clone();
         }
-        if self.is_negative() {
-            return self.abs().inverse().neg();
-        }
-        let guess = {
-            let bits = self.int_val.bits() as f64;
-            let scale = self.scale as f64;
-
-            let magic_factor = 0.721507597259061_f64;
-            let initial_guess = scale * LOG2_10 - bits;
-            let res = magic_factor * exp2(initial_guess);
-
-            if res.is_normal() {
-                BigDecimal::try_from(res).unwrap()
-            } else {
-                // can't guess with float - just guess magnitude
-                let scale = (bits / LOG2_10 + scale).round() as i64;
-                BigDecimal::new(BigInt::from(1), -scale)
-            }
-        };
-
-        let max_precision = DEFAULT_PRECISION;
-        let next_iteration = move |r: BigDecimal| {
-            let two = BigDecimal::from(2);
-            let tmp = two - self * &r;
 
-            r * tmp
-        };
-
-        // calculate first iteration
-        let mut running_result = next_iteration(guess);
-
-        let mut prev_result = BigDecimal::one();
-        let mut result = BigDecimal::zero();
-
-        // TODO: Prove that we don't need to arbitrarily limit iterations
-        // and that convergence can be calculated
-        while prev_result != result {
-            // store current result to test for convergence
-            prev_result = result;
+        let uint = self.int_val.magnitude();
+        let result = arithmetic::inverse::impl_inverse_uint_scale(uint, self.scale, ctx);
 
-            // calculate next iteration
-            running_result = next_iteration(running_result).with_prec(max_precision);
-
-            // 'result' has clipped precision, 'running_result' has full precision
-            result = if running_result.digits() > max_precision {
-                running_result.with_prec(max_precision)
-            } else {
-                running_result.clone()
-            };
-        }
-
-        return result;
+        // always copy sign
+        result.take_with_sign(self.sign())
     }
 
     /// Return number rounded to round_digits precision after the decimal point
@@ -1072,821 +996,91 @@ impl Ord for BigDecimal {
     }
 }
 
-impl PartialEq for BigDecimal {
-    #[inline]
-    fn eq(&self, rhs: &BigDecimal) -> bool {
-        // fix scale and test equality
-        match self.scale.cmp(&rhs.scale) {
-            Ordering::Greater => {
-                let scaled_int_val = &rhs.int_val * ten_to_the((self.scale - rhs.scale) as u64);
-                self.int_val == scaled_int_val
-            }
-            Ordering::Less => {
-                let scaled_int_val = &self.int_val * ten_to_the((rhs.scale - self.scale) as u64);
-                scaled_int_val == rhs.int_val
-            }
-            Ordering::Equal => self.int_val == rhs.int_val,
-        }
-    }
-}
 
 impl Default for BigDecimal {
     #[inline]
     fn default() -> BigDecimal {
-        Zero::zero()
-    }
-}
-
-impl Zero for BigDecimal {
-    #[inline]
-    fn zero() -> BigDecimal {
-        BigDecimal::new(BigInt::zero(), 0)
-    }
-
-    #[inline]
-    fn is_zero(&self) -> bool {
-        self.int_val.is_zero()
-    }
-}
-
-impl One for BigDecimal {
-    #[inline]
-    fn one() -> BigDecimal {
-        BigDecimal::new(BigInt::one(), 0)
-    }
-}
-
-
-impl Add<BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn add(self, rhs: BigDecimal) -> BigDecimal {
-        let mut lhs = self;
-
-        match lhs.scale.cmp(&rhs.scale) {
-            Ordering::Equal => {
-                lhs.int_val += rhs.int_val;
-                lhs
-            }
-            Ordering::Less => lhs.take_and_scale(rhs.scale) + rhs,
-            Ordering::Greater => rhs.take_and_scale(lhs.scale) + lhs,
-        }
-    }
-}
-
-impl<'a> Add<&'a BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn add(self, rhs: &'a BigDecimal) -> BigDecimal {
-        let mut lhs = self;
-
-        match lhs.scale.cmp(&rhs.scale) {
-            Ordering::Equal => {
-                lhs.int_val += &rhs.int_val;
-                lhs
-            }
-            Ordering::Less => lhs.take_and_scale(rhs.scale) + rhs,
-            Ordering::Greater => rhs.with_scale(lhs.scale) + lhs,
-        }
-    }
-}
-
-impl<'a> Add<BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn add(self, rhs: BigDecimal) -> BigDecimal {
-        rhs + self
-    }
-}
-
-impl<'a, 'b> Add<&'b BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn add(self, rhs: &BigDecimal) -> BigDecimal {
-        let lhs = self;
-        match self.scale.cmp(&rhs.scale) {
-            Ordering::Less => lhs.with_scale(rhs.scale) + rhs,
-            Ordering::Greater => rhs.with_scale(lhs.scale) + lhs,
-            Ordering::Equal => BigDecimal::new(lhs.int_val.clone() + &rhs.int_val, lhs.scale),
-        }
-    }
-}
-
-impl Add<BigInt> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn add(self, rhs: BigInt) -> BigDecimal {
-        let mut lhs = self;
-
-        match lhs.scale.cmp(&0) {
-            Ordering::Equal => {
-                lhs.int_val += rhs;
-                lhs
-            }
-            Ordering::Greater => {
-                lhs.int_val += rhs * ten_to_the(lhs.scale as u64);
-                lhs
-            }
-            Ordering::Less => lhs.take_and_scale(0) + rhs,
-        }
-    }
-}
-
-impl<'a> Add<&'a BigInt> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn add(self, rhs: &BigInt) -> BigDecimal {
-        let mut lhs = self;
-
-        match lhs.scale.cmp(&0) {
-            Ordering::Equal => {
-                lhs.int_val += rhs;
-                lhs
-            }
-            Ordering::Greater => {
-                lhs.int_val += rhs * ten_to_the(lhs.scale as u64);
-                lhs
-            }
-            Ordering::Less => lhs.take_and_scale(0) + rhs,
-        }
-    }
-}
-
-impl<'a> Add<BigInt> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn add(self, rhs: BigInt) -> BigDecimal {
-        BigDecimal::new(rhs, 0) + self
-    }
-}
-
-impl<'a, 'b> Add<&'a BigInt> for &'b BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn add(self, rhs: &BigInt) -> BigDecimal {
-        self.with_scale(0) + rhs
-    }
-}
-
-forward_val_assignop!(impl AddAssign for BigDecimal, add_assign);
-
-impl<'a> AddAssign<&'a BigDecimal> for BigDecimal {
-    #[inline]
-    fn add_assign(&mut self, rhs: &BigDecimal) {
-        match self.scale.cmp(&rhs.scale) {
-            Ordering::Less => {
-                let scaled = self.with_scale(rhs.scale);
-                self.int_val = scaled.int_val + &rhs.int_val;
-                self.scale = rhs.scale;
-            }
-            Ordering::Greater => {
-                let scaled = rhs.with_scale(self.scale);
-                self.int_val += scaled.int_val;
-            }
-            Ordering::Equal => {
-                self.int_val += &rhs.int_val;
-            }
-        }
-    }
-}
-
-impl AddAssign<BigInt> for BigDecimal {
-    #[inline]
-    fn add_assign(&mut self, rhs: BigInt) {
-        *self += BigDecimal::new(rhs, 0)
-    }
-}
-
-impl<'a> AddAssign<&'a BigInt> for BigDecimal {
-    #[inline]
-    fn add_assign(&mut self, rhs: &BigInt) {
-        match self.scale.cmp(&0) {
-            Ordering::Equal => self.int_val += rhs,
-            Ordering::Greater => self.int_val += rhs * ten_to_the(self.scale as u64),
-            Ordering::Less => {
-                // *self += BigDecimal::new(rhs, 0)
-                self.int_val *= ten_to_the((-self.scale) as u64);
-                self.int_val += rhs;
-                self.scale = 0;
-            }
-        }
-    }
-}
-
-impl Sub<BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn sub(self, rhs: BigDecimal) -> BigDecimal {
-        let mut lhs = self;
-        let scale = cmp::max(lhs.scale, rhs.scale);
-
-        match lhs.scale.cmp(&rhs.scale) {
-            Ordering::Equal => {
-                lhs.int_val -= rhs.int_val;
-                lhs
-            }
-            Ordering::Less => lhs.take_and_scale(scale) - rhs,
-            Ordering::Greater => lhs - rhs.take_and_scale(scale),
-        }
-    }
-}
-
-impl<'a> Sub<&'a BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn sub(self, rhs: &BigDecimal) -> BigDecimal {
-        let mut lhs = self;
-        let scale = cmp::max(lhs.scale, rhs.scale);
-
-        match lhs.scale.cmp(&rhs.scale) {
-            Ordering::Equal => {
-                lhs.int_val -= &rhs.int_val;
-                lhs
-            }
-            Ordering::Less => lhs.take_and_scale(rhs.scale) - rhs,
-            Ordering::Greater => lhs - rhs.with_scale(scale),
-        }
-    }
-}
-
-impl<'a> Sub<BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn sub(self, rhs: BigDecimal) -> BigDecimal {
-        -(rhs - self)
-    }
-}
-
-impl<'a, 'b> Sub<&'b BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn sub(self, rhs: &BigDecimal) -> BigDecimal {
-        match self.scale.cmp(&rhs.scale) {
-            Ordering::Greater => {
-                let rhs = rhs.with_scale(self.scale);
-                self - rhs
-            }
-            Ordering::Less => self.with_scale(rhs.scale) - rhs,
-            Ordering::Equal => BigDecimal::new(&self.int_val - &rhs.int_val, self.scale),
-        }
-    }
-}
-
-impl Sub<BigInt> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn sub(self, rhs: BigInt) -> BigDecimal {
-        let mut lhs = self;
-
-        match lhs.scale.cmp(&0) {
-            Ordering::Equal => {
-                lhs.int_val -= rhs;
-                lhs
-            }
-            Ordering::Greater => {
-                lhs.int_val -= rhs * ten_to_the(lhs.scale as u64);
-                lhs
-            }
-            Ordering::Less => lhs.take_and_scale(0) - rhs,
-        }
-    }
-}
-
-impl<'a> Sub<&'a BigInt> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn sub(self, rhs: &BigInt) -> BigDecimal {
-        let mut lhs = self;
-
-        match lhs.scale.cmp(&0) {
-            Ordering::Equal => {
-                lhs.int_val -= rhs;
-                lhs
-            }
-            Ordering::Greater => {
-                lhs.int_val -= rhs * ten_to_the(lhs.scale as u64);
-                lhs
-            }
-            Ordering::Less => lhs.take_and_scale(0) - rhs,
-        }
-    }
-}
-
-impl<'a> Sub<BigInt> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn sub(self, rhs: BigInt) -> BigDecimal {
-        BigDecimal::new(rhs, 0) - self
-    }
-}
-
-impl<'a, 'b> Sub<&'a BigInt> for &'b BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn sub(self, rhs: &BigInt) -> BigDecimal {
-        self.with_scale(0) - rhs
-    }
-}
-
-forward_val_assignop!(impl SubAssign for BigDecimal, sub_assign);
-
-impl<'a> SubAssign<&'a BigDecimal> for BigDecimal {
-    #[inline]
-    fn sub_assign(&mut self, rhs: &BigDecimal) {
-        match self.scale.cmp(&rhs.scale) {
-            Ordering::Less => {
-                let lhs = self.with_scale(rhs.scale);
-                self.int_val = lhs.int_val - &rhs.int_val;
-                self.scale = rhs.scale;
-            }
-            Ordering::Greater => {
-                self.int_val -= rhs.with_scale(self.scale).int_val;
-            }
-            Ordering::Equal => {
-                self.int_val = &self.int_val - &rhs.int_val;
-            }
-        }
-    }
-}
-
-impl SubAssign<BigInt> for BigDecimal {
-    #[inline(always)]
-    fn sub_assign(&mut self, rhs: BigInt) {
-        *self -= BigDecimal::new(rhs, 0)
-    }
-}
-
-impl<'a> SubAssign<&'a BigInt> for BigDecimal {
-    #[inline(always)]
-    fn sub_assign(&mut self, rhs: &BigInt) {
-        match self.scale.cmp(&0) {
-            Ordering::Equal => SubAssign::sub_assign(&mut self.int_val, rhs),
-            Ordering::Greater => SubAssign::sub_assign(&mut self.int_val, rhs * ten_to_the(self.scale as u64)),
-            Ordering::Less => {
-                self.int_val *= ten_to_the((-self.scale) as u64);
-                SubAssign::sub_assign(&mut self.int_val, rhs);
-                self.scale = 0;
-            }
-        }
-    }
-}
-
-impl Mul<BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(mut self, rhs: BigDecimal) -> BigDecimal {
-        if self.is_one() {
-            rhs
-        } else if rhs.is_one() {
-            self
-        } else {
-            self.scale += rhs.scale;
-            self.int_val *= rhs.int_val;
-            self
-        }
-    }
-}
-
-impl<'a> Mul<&'a BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(mut self, rhs: &'a BigDecimal) -> BigDecimal {
-        if self.is_one() {
-            self.scale = rhs.scale;
-            self.int_val.set_zero();
-            self.int_val.add_assign(&rhs.int_val);
-            self
-        } else if rhs.is_one() {
-            self
-        } else {
-            self.scale += rhs.scale;
-            MulAssign::mul_assign(&mut self.int_val, &rhs.int_val);
-            self
-        }
-    }
-}
-
-impl<'a> Mul<BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(self, rhs: BigDecimal) -> BigDecimal {
-        rhs * self
-    }
-}
-
-impl<'a, 'b> Mul<&'b BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(self, rhs: &BigDecimal) -> BigDecimal {
-        if self.is_one() {
-            rhs.normalized()
-        } else if rhs.is_one() {
-            self.normalized()
-        } else {
-            let scale = self.scale + rhs.scale;
-            BigDecimal::new(&self.int_val * &rhs.int_val, scale)
-        }
-    }
-}
-
-impl Mul<BigInt> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(mut self, rhs: BigInt) -> BigDecimal {
-        self.int_val *= rhs;
-        self
-    }
-}
-
-impl<'a> Mul<&'a BigInt> for BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(mut self, rhs: &BigInt) -> BigDecimal {
-        self.int_val *= rhs;
-        self
-    }
-}
-
-impl<'a> Mul<BigInt> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(self, mut rhs: BigInt) -> BigDecimal {
-        rhs *= &self.int_val;
-        BigDecimal::new(rhs, self.scale)
-    }
-}
-
-impl<'a, 'b> Mul<&'a BigInt> for &'b BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(self, rhs: &BigInt) -> BigDecimal {
-        if rhs.is_one() {
-            self.normalized()
-        } else if self.is_one() {
-            BigDecimal::new(rhs.clone(), 0)
-        } else {
-            let value = &self.int_val * rhs;
-            BigDecimal::new(value, self.scale)
-        }
-    }
-}
-
-impl Mul<BigDecimal> for BigInt {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(mut self, mut rhs: BigDecimal) -> BigDecimal {
-        if rhs.is_one() {
-            rhs.scale = 0;
-            stdlib::mem::swap(&mut rhs.int_val, &mut self);
-        } else if !self.is_one() {
-            rhs.int_val *= self;
-        }
-        rhs
-    }
-}
-
-impl<'a> Mul<BigDecimal> for &'a BigInt {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(self, mut rhs: BigDecimal) -> BigDecimal {
-        if self.is_one() {
-            rhs.normalized()
-        } else if rhs.is_one() {
-            rhs.int_val.set_zero();
-            rhs.int_val += self;
-            rhs.scale = 0;
-            rhs
-        } else {
-            rhs.int_val *= self;
-            rhs
-        }
-    }
-}
-
-impl<'a, 'b> Mul<&'a BigDecimal> for &'b BigInt {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(self, rhs: &BigDecimal) -> BigDecimal {
-        if self.is_one() {
-            rhs.normalized()
-        } else if rhs.is_one() {
-            BigDecimal::new(self.clone(), 0)
-        } else {
-            let value = &rhs.int_val * self;
-            BigDecimal::new(value, rhs.scale)
-        }
-    }
-}
-
-impl<'a> Mul<&'a BigDecimal> for BigInt {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn mul(mut self, rhs: &BigDecimal) -> BigDecimal {
-        if self.is_one() {
-            rhs.normalized()
-        } else if rhs.is_one() {
-            BigDecimal::new(self, 0)
-        } else {
-            self *= &rhs.int_val;
-            BigDecimal::new(self, rhs.scale)
-        }
-    }
-}
-
-forward_val_assignop!(impl MulAssign for BigDecimal, mul_assign);
-
-impl<'a> MulAssign<&'a BigDecimal> for BigDecimal {
-    #[inline]
-    fn mul_assign(&mut self, rhs: &BigDecimal) {
-        if rhs.is_one() {
-            return;
-        }
-        self.scale += rhs.scale;
-        self.int_val = &self.int_val * &rhs.int_val;
-    }
-}
-
-impl<'a> MulAssign<&'a BigInt> for BigDecimal {
-    #[inline]
-    fn mul_assign(&mut self, rhs: &BigInt) {
-        if rhs.is_one() {
-            return;
-        }
-        self.int_val *= rhs;
-    }
-}
-
-impl MulAssign<BigInt> for BigDecimal {
-    #[inline]
-    fn mul_assign(&mut self, rhs: BigInt) {
-        *self *= &rhs
-    }
-}
-
-
-#[inline(always)]
-fn impl_division(mut num: BigInt, den: &BigInt, mut scale: i64, max_precision: u64) -> BigDecimal {
-    // quick zero check
-    if num.is_zero() {
-        return BigDecimal::new(num, 0);
-    }
-
-    match (num.is_negative(), den.is_negative()) {
-        (true, true) => return impl_division(num.neg(), &den.neg(), scale, max_precision),
-        (true, false) => return -impl_division(num.neg(), den, scale, max_precision),
-        (false, true) => return -impl_division(num, &den.neg(), scale, max_precision),
-        (false, false) => (),
-    }
-
-    // shift digits until numerator is larger than denominator (set scale appropriately)
-    while num < *den {
-        scale += 1;
-        num *= 10;
-    }
-
-    // first division
-    let (mut quotient, mut remainder) = num.div_rem(den);
-
-    // division complete
-    if remainder.is_zero() {
-        return BigDecimal {
-            int_val: quotient,
-            scale: scale,
-        };
-    }
-
-    let mut precision = count_decimal_digits(&quotient);
-
-    // shift remainder by 1 decimal;
-    // quotient will be 1 digit upon next division
-    remainder *= 10;
-
-    while !remainder.is_zero() && precision < max_precision {
-        let (q, r) = remainder.div_rem(den);
-        quotient = quotient * 10 + q;
-        remainder = r * 10;
-
-        precision += 1;
-        scale += 1;
-    }
-
-    if !remainder.is_zero() {
-        // round final number with remainder
-        quotient += get_rounding_term(&remainder.div(den));
-    }
-
-    let result = BigDecimal::new(quotient, scale);
-    // println!(" {} / {}\n = {}\n", self, other, result);
-    return result;
-}
-
-impl Div<BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-    #[inline]
-    fn div(self, other: BigDecimal) -> BigDecimal {
-        if other.is_zero() {
-            panic!("Division by zero");
-        }
-        if self.is_zero() || other.is_one() {
-            return self;
-        }
-
-        let scale = self.scale - other.scale;
-
-        if self.int_val == other.int_val {
-            return BigDecimal {
-                int_val: 1.into(),
-                scale: scale,
-            };
-        }
-
-        let max_precision = DEFAULT_PRECISION;
-
-        return impl_division(self.int_val, &other.int_val, scale, max_precision);
-    }
-}
-
-impl<'a> Div<&'a BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-    #[inline]
-    fn div(self, other: &'a BigDecimal) -> BigDecimal {
-        if other.is_zero() {
-            panic!("Division by zero");
-        }
-        if self.is_zero() || other.is_one() {
-            return self;
-        }
-
-        let scale = self.scale - other.scale;
-
-        if self.int_val == other.int_val {
-            return BigDecimal {
-                int_val: 1.into(),
-                scale: scale,
-            };
-        }
-
-        let max_precision = DEFAULT_PRECISION;
-
-        return impl_division(self.int_val, &other.int_val, scale, max_precision);
-    }
-}
-
-forward_ref_val_binop!(impl Div for BigDecimal, div);
-
-impl<'a, 'b> Div<&'b BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
-
-    #[inline]
-    fn div(self, other: &BigDecimal) -> BigDecimal {
-        if other.is_zero() {
-            panic!("Division by zero");
-        }
-        // TODO: Fix setting scale
-        if self.is_zero() || other.is_one() {
-            return self.clone();
-        }
-
-        let scale = self.scale - other.scale;
-
-        let num_int = &self.int_val;
-        let den_int = &other.int_val;
-
-        if num_int == den_int {
-            return BigDecimal {
-                int_val: 1.into(),
-                scale: scale,
-            };
-        }
-
-        let max_precision = DEFAULT_PRECISION;
-
-        return impl_division(num_int.clone(), den_int, scale, max_precision);
+        Zero::zero()
     }
 }
 
-impl Rem<BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-
+impl Zero for BigDecimal {
     #[inline]
-    fn rem(self, other: BigDecimal) -> BigDecimal {
-        let scale = cmp::max(self.scale, other.scale);
-
-        let num = self.take_and_scale(scale).int_val;
-        let den = other.take_and_scale(scale).int_val;
+    fn zero() -> BigDecimal {
+        BigDecimal::new(BigInt::zero(), 0)
+    }
 
-        BigDecimal::new(num % den, scale)
+    #[inline]
+    fn is_zero(&self) -> bool {
+        self.int_val.is_zero()
     }
 }
 
-impl<'a> Rem<&'a BigDecimal> for BigDecimal {
-    type Output = BigDecimal;
-
+impl One for BigDecimal {
     #[inline]
-    fn rem(self, other: &BigDecimal) -> BigDecimal {
-        let scale = cmp::max(self.scale, other.scale);
-        let num = self.take_and_scale(scale).int_val;
-        let den = &other.int_val;
-
-        let result = if scale == other.scale {
-            num % den
-        } else {
-            num % (den * ten_to_the((scale - other.scale) as u64))
-        };
-        BigDecimal::new(result, scale)
+    fn one() -> BigDecimal {
+        BigDecimal::new(BigInt::one(), 0)
     }
 }
-impl<'a> Rem<BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
 
-    #[inline]
-    fn rem(self, other: BigDecimal) -> BigDecimal {
-        let scale = cmp::max(self.scale, other.scale);
-        let num = &self.int_val;
-        let den = other.take_and_scale(scale).int_val;
 
-        let result = if scale == self.scale {
-            num % den
-        } else {
-            let scaled_num = num * ten_to_the((scale - self.scale) as u64);
-            scaled_num % den
-        };
+fn impl_division(mut num: BigInt, den: &BigInt, mut scale: i64, max_precision: u64) -> BigDecimal {
+    // quick zero check
+    if num.is_zero() {
+        return BigDecimal::new(num, 0);
+    }
 
-        BigDecimal::new(result, scale)
+    match (num.is_negative(), den.is_negative()) {
+        (true, true) => return impl_division(num.neg(), &den.neg(), scale, max_precision),
+        (true, false) => return -impl_division(num.neg(), den, scale, max_precision),
+        (false, true) => return -impl_division(num, &den.neg(), scale, max_precision),
+        (false, false) => (),
     }
-}
 
-impl<'a, 'b> Rem<&'b BigDecimal> for &'a BigDecimal {
-    type Output = BigDecimal;
+    // shift digits until numerator is larger than denominator (set scale appropriately)
+    while num < *den {
+        scale += 1;
+        num *= 10;
+    }
 
-    #[inline]
-    fn rem(self, other: &BigDecimal) -> BigDecimal {
-        let scale = cmp::max(self.scale, other.scale);
-        let num = &self.int_val;
-        let den = &other.int_val;
+    // first division
+    let (mut quotient, mut remainder) = num.div_rem(den);
 
-        let result = match self.scale.cmp(&other.scale) {
-            Ordering::Equal => num % den,
-            Ordering::Less => {
-                let scaled_num = num * ten_to_the((scale - self.scale) as u64);
-                scaled_num % den
-            }
-            Ordering::Greater => {
-                let scaled_den = den * ten_to_the((scale - other.scale) as u64);
-                num % scaled_den
-            }
+    // division complete
+    if remainder.is_zero() {
+        return BigDecimal {
+            int_val: quotient,
+            scale: scale,
         };
-        BigDecimal::new(result, scale)
     }
-}
 
-impl Neg for BigDecimal {
-    type Output = BigDecimal;
+    let mut precision = count_decimal_digits(&quotient);
 
-    #[inline]
-    fn neg(mut self) -> BigDecimal {
-        self.int_val = -self.int_val;
-        self
-    }
-}
+    // shift remainder by 1 decimal;
+    // quotient will be 1 digit upon next division
+    remainder *= 10;
 
-impl<'a> Neg for &'a BigDecimal {
-    type Output = BigDecimal;
+    while !remainder.is_zero() && precision < max_precision {
+        let (q, r) = remainder.div_rem(den);
+        quotient = quotient * 10 + q;
+        remainder = r * 10;
 
-    #[inline]
-    fn neg(self) -> BigDecimal {
-        -self.clone()
+        precision += 1;
+        scale += 1;
+    }
+
+    if !remainder.is_zero() {
+        // round final number with remainder
+        quotient += get_rounding_term(&remainder.div(den));
     }
+
+    let result = BigDecimal::new(quotient, scale);
+    // println!(" {} / {}\n = {}\n", self, other, result);
+    return result;
 }
 
+
+
 impl Signed for BigDecimal {
     #[inline]
     fn abs(&self) -> BigDecimal {
@@ -2002,6 +1196,163 @@ impl fmt::Debug for BigDecimal {
 }
 
 
+/// Immutable big-decimal, referencing a borrowed buffer of digits
+///
+/// The non-digit information like `scale` and `sign` may be changed
+/// on these objects, which otherwise would require cloning the full
+/// digit buffer in the BigDecimal.
+///
+/// Built from full `BigDecimal` object using the `to_ref()` method.
+/// `BigDecimal` not implement `AsRef`, so we will reserve the method
+/// `as_ref()` for a later time.
+///
+/// May be transformed into full BigDecimal object using the `to_owned()`
+/// method.
+/// This clones the bigdecimal digits.
+///
+/// BigDecimalRef (or `Into<BigDecimalRef>`) should be preferred over
+/// using `&BigDecimal` for library functions that need an immutable
+/// reference to a bigdecimal, as it may be much more efficient.
+///
+/// NOTE: Using `&BigDecimalRef` is redundant, and not recommended.
+///
+/// ## Examples
+///
+/// ```
+/// # use bigdecimal::*; use std::ops::Neg;
+/// fn add_one<'a, N: Into<BigDecimalRef<'a>>>(n: N) -> BigDecimal {
+///     n.into() + 1
+/// }
+///
+/// let n: BigDecimal = "123.456".parse().unwrap();
+///
+/// // call via "standard" reference (implements Into)
+/// let m = add_one(&n);
+/// assert_eq!(m, "124.456".parse().unwrap());
+///
+/// // call by negating the reference (fast: no-digit cloning involved)
+/// let m = add_one(n.to_ref().neg());
+/// assert_eq!(m, "-122.456".parse().unwrap());
+/// ```
+///
+#[derive(Clone, Copy, Debug)]
+pub struct BigDecimalRef<'a> {
+    sign: Sign,
+    digits: &'a BigUint,
+    scale: i64,
+}
+
+impl BigDecimalRef<'_> {
+    /// Clone digits to make this reference a full BigDecimal object
+    pub fn to_owned(&self) -> BigDecimal {
+        BigDecimal {
+            scale: self.scale,
+            int_val: BigInt::from_biguint(self.sign, self.digits.clone()),
+        }
+    }
+
+    /// Clone digits, returning BigDecimal with given scale
+    ///
+    /// ```
+    /// # use bigdecimal::*;
+    ///
+    /// let n: BigDecimal = "123.45678".parse().unwrap();
+    /// let r = n.to_ref();
+    /// assert_eq!(r.to_owned_with_scale(5), n.clone());
+    /// assert_eq!(r.to_owned_with_scale(0), "123".parse().unwrap());
+    /// assert_eq!(r.to_owned_with_scale(-1), "12e1".parse().unwrap());
+    ///
+    /// let x = r.to_owned_with_scale(8);
+    /// assert_eq!(&x, &n);
+    /// assert_eq!(x.fractional_digit_count(), 8);
+    /// ```
+    pub fn to_owned_with_scale(&self, scale: i64) -> BigDecimal {
+        use stdlib::cmp::Ordering::*;
+
+        let digits = match scale.cmp(&self.scale) {
+            Equal => self.digits.clone(),
+            Greater => self.digits * ten_to_the_uint((scale - self.scale) as u64),
+            Less => self.digits / ten_to_the_uint((self.scale - scale) as u64)
+        };
+
+        BigDecimal {
+            scale: scale,
+            int_val: BigInt::from_biguint(self.sign, digits),
+        }
+    }
+
+    /// Sign of decimal
+    pub fn sign(&self) -> Sign {
+        self.sign
+    }
+
+    /// Return number of digits 'right' of the decimal point
+    /// (including leading zeros)
+    pub fn fractional_digit_count(&self) -> i64 {
+        self.scale
+    }
+
+    /// Count total number of decimal digits
+    pub fn count_digits(&self) -> u64 {
+        count_decimal_digits_uint(self.digits)
+    }
+
+    /// Split into components
+    pub(crate) fn as_parts(&self) -> (Sign, i64, &BigUint) {
+        (self.sign, self.scale, self.digits)
+    }
+
+    /// Take absolute value of the decimal (non-negative sign)
+    pub fn abs(&self) -> Self {
+        Self {
+            sign: self.sign * self.sign,
+            digits: self.digits,
+            scale: self.scale,
+        }
+    }
+
+    /// Take square root of this number
+    pub fn sqrt_with_context(&self, ctx: &Context) -> Option<BigDecimal> {
+        use Sign::*;
+
+        let (sign, scale, uint) = self.as_parts();
+
+        match sign {
+            Minus => None,
+            NoSign => Some(Zero::zero()),
+            Plus => Some(arithmetic::sqrt::impl_sqrt(uint, scale, ctx)),
+        }
+    }
+
+    /// Return if the referenced decimal is zero
+    pub fn is_zero(&self) -> bool {
+        self.digits.is_zero()
+    }
+}
+
+impl<'a> From<&'a BigDecimal> for BigDecimalRef<'a> {
+    fn from(n: &'a BigDecimal) -> Self {
+        let sign = n.int_val.sign();
+        let mag = n.int_val.magnitude();
+        Self {
+            sign: sign,
+            digits: mag,
+            scale: n.scale,
+        }
+    }
+}
+
+impl<'a> From<&'a BigInt> for BigDecimalRef<'a> {
+    fn from(n: &'a BigInt) -> Self {
+        Self {
+            sign: n.sign(),
+            digits: n.magnitude(),
+            scale: 0,
+        }
+    }
+}
+
+
 /// Tools to help serializing/deserializing `BigDecimal`s
 #[cfg(feature = "serde")]
 mod bigdecimal_serde {
@@ -2175,6 +1526,29 @@ mod bigdecimal_tests {
     use num_bigint;
     use paste::paste;
 
+    #[test]
+    fn test_fractional_digit_count() {
+        // Zero value
+        let vals = BigDecimal::from(0);
+        assert_eq!(vals.fractional_digit_count(), 0);
+        assert_eq!(vals.to_ref().fractional_digit_count(), 0);
+
+        // Fractional part with trailing zeros
+        let vals = BigDecimal::from_str("1.0").unwrap();
+        assert_eq!(vals.fractional_digit_count(), 1);
+        assert_eq!(vals.to_ref().fractional_digit_count(), 1);
+
+        // Fractional part
+        let vals = BigDecimal::from_str("1.23").unwrap();
+        assert_eq!(vals.fractional_digit_count(), 2);
+        assert_eq!(vals.to_ref().fractional_digit_count(), 2);
+
+        // shifted to 'left' has negative scale
+        let vals = BigDecimal::from_str("123e5").unwrap();
+        assert_eq!(vals.fractional_digit_count(), -5);
+        assert_eq!(vals.to_ref().fractional_digit_count(), -5);
+    }
+
     #[test]
     fn test_sum() {
         let vals = vec![
@@ -2342,234 +1716,6 @@ mod bigdecimal_tests {
         assert!(BigDecimal::try_from(f64::NAN).is_err());
     }
 
-    #[test]
-    fn test_add() {
-        let vals = vec![
-            ("12.34", "1.234", "13.574"),
-            ("12.34", "-1.234", "11.106"),
-            ("1234e6", "1234e-6", "1234000000.001234"),
-            ("1234e-6", "1234e6", "1234000000.001234"),
-            ("18446744073709551616.0", "1", "18446744073709551617"),
-            ("184467440737e3380", "0", "184467440737e3380"),
-        ];
-
-        for &(x, y, z) in vals.iter() {
-
-            let mut a = BigDecimal::from_str(x).unwrap();
-            let b = BigDecimal::from_str(y).unwrap();
-            let c = BigDecimal::from_str(z).unwrap();
-
-            assert_eq!(a.clone() + b.clone(), c);
-
-            assert_eq!(a.clone() + &b, c);
-            assert_eq!(&a + b.clone(), c);
-            assert_eq!(&a + &b, c);
-
-            a += b;
-            assert_eq!(a, c);
-        }
-    }
-
-    #[test]
-    fn test_sub() {
-        let vals = vec![
-            ("12.34", "1.234", "11.106"),
-            ("12.34", "-1.234", "13.574"),
-            ("1234e6", "1234e-6", "1233999999.998766"),
-        ];
-
-        for &(x, y, z) in vals.iter() {
-
-            let mut a = BigDecimal::from_str(x).unwrap();
-            let b = BigDecimal::from_str(y).unwrap();
-            let c = BigDecimal::from_str(z).unwrap();
-
-            assert_eq!(a.clone() - b.clone(), c);
-
-            assert_eq!(a.clone() - &b, c);
-            assert_eq!(&a - b.clone(), c);
-            assert_eq!(&a - &b, c);
-
-            a -= b;
-            assert_eq!(a, c);
-        }
-    }
-
-    /// Test multiplication of two bigdecimals
-    #[test]
-    fn test_mul() {
-
-        let vals = vec![
-            ("2", "1", "2"),
-            ("12.34", "1.234", "15.22756"),
-            ("2e1", "1", "20"),
-            ("3", ".333333", "0.999999"),
-            ("2389472934723", "209481029831", "500549251119075878721813"),
-            ("1e-450", "1e500", ".1e51"),
-            ("-995052931372975485719.533153137", "4.523087321", "-4500711297616988541501.836966993116075977"),
-            ("995052931372975485719.533153137", "-4.523087321", "-4500711297616988541501.836966993116075977"),
-            ("-8.37664968", "-1.9086963714056968482094712882596748", "15.988480848752691653730876239769592670324064"),
-            ("-8.37664968", "0", "0"),
-        ];
-
-        for &(x, y, z) in vals.iter() {
-
-            let mut a = BigDecimal::from_str(x).unwrap();
-            let b = BigDecimal::from_str(y).unwrap();
-            let c = BigDecimal::from_str(z).unwrap();
-
-            assert_eq!(a.clone() * b.clone(), c);
-            assert_eq!(a.clone() * &b, c);
-            assert_eq!(&a * b.clone(), c);
-            assert_eq!(&a * &b, c);
-
-            a *= b;
-            assert_eq!(a, c);
-        }
-    }
-
-    /// Test multiplication between big decimal and big integer
-    #[test]
-    fn test_mul_bigint() {
-        let vals = vec![
-            ("2", "1", "2"),
-            ("8.561", "10", "85.61"),
-            ("1.0000", "638655273892892437", "638655273892892437"),
-            ("10000", "638655273892892437", "6386552738928924370000"),
-            (".0005", "368408638655273892892437473", "184204319327636946446218.7365"),
-            ("9e-1", "368408638655273892892437473", "331567774789746503603193725.7"),
-            ("-1.175470587012343730098", "577575785", "-678923347.038065234601180476930"),
-            ("-1.175470587012343730098", "-76527768352678", "89956140788267.069799533723307502444"),
-            ("-1.175470587012343730098", "0", "0"),
-        ];
-
-        for &(x, y, z) in vals.iter() {
-            let a = BigDecimal::from_str(x).unwrap();
-            let b = num_bigint::BigInt::from_str(y).unwrap();
-            let c = BigDecimal::from_str(z).unwrap();
-
-            assert_eq!(a.clone() * b.clone(), c);
-            assert_eq!(b.clone() * a.clone(), c);
-            assert_eq!(a.clone() * &b, c);
-            assert_eq!(b.clone() * &a, c);
-            assert_eq!(&a * b.clone(), c);
-            assert_eq!(&b * a.clone(), c);
-            assert_eq!(&a * &b, c);
-            assert_eq!(&b * &a, c);
-        }
-    }
-
-    #[test]
-    fn test_div() {
-        let vals = vec![
-            ("0", "1", "0"),
-            ("0", "10", "0"),
-            ("2", "1", "2"),
-            ("2e1", "1", "2e1"),
-            ("10", "10", "1"),
-            ("100", "10.0", "1e1"),
-            ("20.0", "200", ".1"),
-            ("4", "2", "2.0"),
-            ("15", "3", "5.0"),
-            ("1", "2", "0.5"),
-            ("1", "2e-2", "5e1"),
-            ("1", "0.2", "5"),
-            ("1.0", "0.02", "50"),
-            ("1", "0.020", "5e1"),
-            ("5.0", "4.00", "1.25"),
-            ("5.0", "4.000", "1.25"),
-            ("5", "4.000", "1.25"),
-            ("5", "4", "125e-2"),
-            ("100", "5", "20"),
-            ("-50", "5", "-10"),
-            ("200", "-5", "-40."),
-            ("1", "3", ".3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333"),
-            ("-2", "-3", ".6666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667"),
-            ("-12.34", "1.233", "-10.00811030008110300081103000811030008110300081103000811030008110300081103000811030008110300081103001"),
-            ("125348", "352.2283", "355.8714617763535752237966114591019517738921035021887792661748076460636467881768727839301952739175132"),
-        ];
-
-        for &(x, y, z) in vals.iter() {
-
-            let a = BigDecimal::from_str(x).unwrap();
-            let b = BigDecimal::from_str(y).unwrap();
-            let c = BigDecimal::from_str(z).unwrap();
-
-            assert_eq!(a.clone() / b.clone(), c);
-            assert_eq!(a.clone() / &b, c);
-            assert_eq!(&a / b.clone(), c);
-            assert_eq!(&a / &b, c);
-            // assert_eq!(q.scale, c.scale);
-
-            // let mut q = a;
-            // q /= b;
-            // assert_eq!(q, c);
-        }
-    }
-
-    #[test]
-    #[should_panic(expected = "Division by zero")]
-    fn test_division_by_zero_panics() {
-        let x = BigDecimal::from_str("3.14").unwrap();
-        let _r = x / 0;
-    }
-
-    #[test]
-    #[should_panic(expected = "Division by zero")]
-    fn test_division_by_zero_panics_v2() {
-        let x = BigDecimal::from_str("3.14").unwrap();
-        let _r = x / BigDecimal::zero();
-    }
-
-    #[test]
-    fn test_rem() {
-        let vals = vec![
-            ("100", "5", "0"),
-            ("2e1", "1", "0"),
-            ("2", "1", "0"),
-            ("1", "3", "1"),
-            ("1", "0.5", "0"),
-            ("1.5", "1", "0.5"),
-            ("1", "3e-2", "1e-2"),
-            ("10", "0.003", "0.001"),
-            ("3", "2", "1"),
-            ("-3", "2", "-1"),
-            ("3", "-2", "1"),
-            ("-3", "-2", "-1"),
-            ("12.34", "1.233", "0.01"),
-        ];
-        for &(x, y, z) in vals.iter() {
-            let a = BigDecimal::from_str(x).unwrap();
-            let b = BigDecimal::from_str(y).unwrap();
-            let c = BigDecimal::from_str(z).unwrap();
-
-            let rem = &a % &b;
-            assert_eq!(rem, c, "{} [&{} % &{}] == {}", rem, a, b, c);
-
-            let rem = a.clone() % &b;
-            assert_eq!(rem, c, "{} [{} % &{}] == {}", rem, a, b, c);
-
-            let rem = &a % b.clone();
-            assert_eq!(rem, c, "{} [&{} % {}] == {}", rem, a, b, c);
-
-            let rem = a.clone() % b.clone();
-            assert_eq!(rem, c, "{} [{} % {}] == {}", rem, a, b, c);
-        }
-        let vals = vec![
-            (("100", -2), ("50", -1), "0"),
-            (("100", 0), ("50", -1), "0"),
-            (("100", -2), ("30", 0), "10"),
-            (("100", 0), ("30", -1), "10"),
-        ];
-        for &((x, xs), (y, ys), z) in vals.iter() {
-            let a = BigDecimal::from_str(x).unwrap().with_scale(xs);
-            let b = BigDecimal::from_str(y).unwrap().with_scale(ys);
-            let c = BigDecimal::from_str(z).unwrap();
-            let rem = &a % &b;
-            assert_eq!(rem, c, "{} [{} % {}] == {}", rem, a, b, c);
-        }
-    }
-
     #[test]
     fn test_equal() {
         let vals = vec![
@@ -2938,76 +2084,6 @@ mod bigdecimal_tests {
         }
     }
 
-    #[test]
-    fn test_sqrt() {
-        let vals = vec![
-            ("1e-232", "1e-116"),
-            ("1.00", "1"),
-            ("1.001", "1.000499875062460964823258287700109753027590031219780479551442971840836093890879944856933288426795152"),
-            ("100", "10"),
-            ("49", "7"),
-            (".25", ".5"),
-            ("0.0152399025", ".12345"),
-            ("152399025", "12345"),
-            (".00400", "0.06324555320336758663997787088865437067439110278650433653715009705585188877278476442688496216758600590"),
-            (".1", "0.3162277660168379331998893544432718533719555139325216826857504852792594438639238221344248108379300295"),
-            ("2", "1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573"),
-            ("125348", "354.0451948551201563108487193176101314241016013304294520812832530590100407318465590778759640828114535"),
-            ("18446744073709551616.1099511", "4294967296.000000000012799992691725492477397918722952224079252026972356303360555051219312462698703293"),
-            ("3.141592653589793115997963468544185161590576171875", "1.772453850905515992751519103139248439290428205003682302442979619028063165921408635567477284443197875"),
-            (".000000000089793115997963468544185161590576171875", "0.000009475922962855041517561783740144225422359796851494316346796373337470068631250135521161989831460407155"),
-            ("0.7177700109762963922745342343167413624881759290454997218753321040760896053150388903350654937434826216803814031987652326749140535150336357405672040727695124057298138872112244784753994931999476811850580200000000000000000000000000000000", "0.8472130847527653667042980517799020703921106560594525833177762276594388966885185567535692987624493813"),
-            ("0.01234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901", "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111"),
-            ("0.1108890000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000444", "0.3330000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000667"),
-        ];
-        for &(x, y) in vals.iter() {
-            let a = BigDecimal::from_str(x).unwrap().sqrt().unwrap();
-            let b = BigDecimal::from_str(y).unwrap();
-            assert_eq!(a, b);
-        }
-    }
-
-    #[test]
-    fn test_big_sqrt() {
-        use num_bigint::BigInt;
-        let vals = vec![
-            (("2", -70), "141421356237309504880168872420969807.8569671875376948073176679737990732478462107038850387534327641573"),
-            (("3", -170), "17320508075688772935274463415058723669428052538103806280558069794519330169088000370811.46186757248576"),
-            (("9", -199), "9486832980505137995996680633298155601158665417975650480572514558377783315917714664032744325137900886"),
-            (("7", -200), "26457513110645905905016157536392604257102591830824501803683344592010688232302836277603928864745436110"),
-            (("777", -204), "2.787471972953270789531596912111625325974789615194854615319795902911796043681078997362635440358922503E+103"),
-            (("7", -600), "2.645751311064590590501615753639260425710259183082450180368334459201068823230283627760392886474543611E+300"),
-            (("2", -900), "1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573E+450"),
-            (("7", -999), "8.366600265340755479781720257851874893928153692986721998111915430804187725943170098308147119649515362E+499"),
-            (("74908163946345982392040522594123773796", -999), "2.736935584670307552030924971360722787091742391079630976117950955395149091570790266754718322365663909E+518"),
-            (("20", -1024), "4.472135954999579392818347337462552470881236719223051448541794490821041851275609798828828816757564550E512"),
-            (("3", 1025), "5.477225575051661134569697828008021339527446949979832542268944497324932771227227338008584361638706258E-513"),
-        ];
-        for &((s, scale), e) in vals.iter() {
-            let expected = BigDecimal::from_str(e).unwrap();
-
-            let sqrt = BigDecimal::new(BigInt::from_str(s).unwrap(), scale).sqrt().unwrap();
-            assert_eq!(sqrt, expected);
-        }
-    }
-
-    #[test]
-    fn test_cbrt() {
-        let vals = vec![
-            ("0.00", "0"),
-            ("1.00", "1"),
-            ("1.001", "1.000333222283909495175449559955220102010284758197360454054345461242739715702641939155238095670636841"),
-            ("10", "2.154434690031883721759293566519350495259344942192108582489235506346411106648340800185441503543243276"),
-            ("-59283293e25", "-84006090355.84281237113712383191213626687332139035750444925827809487776780721673264524620270275301685"),
-            ("94213372931e-127", "2.112049945275324414051072540210070583697242797173805198575907094646677475250362108901530353886613160E-39"),
-        ];
-        for &(x, y) in vals.iter() {
-            let a = BigDecimal::from_str(x).unwrap().cbrt();
-            let b = BigDecimal::from_str(y).unwrap();
-            assert_eq!(a, b);
-        }
-    }
-
     mod double {
         use super::*;
 
@@ -3256,8 +2332,15 @@ mod test_with_scale_round {
     include!("lib.tests.with_scale_round.rs");
 }
 
-// enable these tests with scripts/bigdecimal-property-tests
-// ::PROPERTY-TESTS:: #[cfg(test)] #[macro_use] extern crate proptest;
-// ::PROPERTY-TESTS:: #[cfg(test)] mod property_tests {
-// ::PROPERTY-TESTS::       use super::*; use paste::paste;
-// ::PROPERTY-TESTS::       include!("lib.tests.property-tests.rs"); }
+
+#[cfg(all(test, property_tests))]
+extern crate proptest;
+
+#[cfg(all(test, property_tests))]
+mod proptests {
+    use super::*;
+    use paste::paste;
+    use proptest::*;
+
+    include!("lib.tests.property-tests.rs");
+}
diff --git a/src/parsing.rs b/src/parsing.rs
index 8bd917b..777c158 100644
--- a/src/parsing.rs
+++ b/src/parsing.rs
@@ -1,6 +1,7 @@
 //! Routines for parsing values into BigDecimals
 
 use super::{BigDecimal, ParseBigDecimalError};
+use stdlib::num::FpCategory;
 
 use stdlib::cmp::{self, Ordering};
 
@@ -17,7 +18,7 @@ pub(crate) fn try_parse_from_f32(n: f32) -> Result<BigDecimal, ParseBigDecimalEr
     match n.classify() {
         Nan => Err(ParseBigDecimalError::Other("NAN".into())),
         Infinite => Err(ParseBigDecimalError::Other("Infinite".into())),
-        Subnormal => Err(ParseBigDecimalError::Other("Subnormal".into())),
+        Subnormal => Ok(parse_from_f32_subnormal(n)),
         Normal | Zero => Ok(parse_from_f32(n)),
     }
 }
@@ -49,9 +50,10 @@ fn split_f32_into_parts(f: f32) -> (u32, i64, Sign) {
 
 /// Create bigdecimal from f32
 ///
-/// Non "normal" values is undefined behavior
-///
 pub(crate) fn parse_from_f32(n: f32) -> BigDecimal {
+    if n.classify() == FpCategory::Subnormal {
+        return parse_from_f32_subnormal(n);
+    }
     let bits = n.to_bits();
 
     if (bits << 1) == 0 {
@@ -94,6 +96,30 @@ pub(crate) fn parse_from_f32(n: f32) -> BigDecimal {
     }
 }
 
+/// Create bigdecimal from subnormal f32
+pub(crate) fn parse_from_f32_subnormal(n: f32) -> BigDecimal {
+    debug_assert_eq!(n.classify(), FpCategory::Subnormal);
+    let bits = n.to_bits();
+
+    let sign_bit = bits >> 31;
+    debug_assert_eq!(bits >> 24, sign_bit << 7);
+
+    let frac = bits - (sign_bit << 31);
+
+    // 5^149 = 5^126 + 5^23  (f32-bit-bias=126, fraction-bits=23)
+    let five_to_149 = BigUint::from_slice(&[
+        1466336501, 2126633373, 2856417274, 1232167559, 2512314040, 1644054862,
+        3843013918, 3873995871, 858643596, 3706384338, 65604258
+    ]);
+
+    let sign = if sign_bit == 0 { Sign::Plus } else { Sign::Minus };
+    let magnitude = BigUint::from(frac) * five_to_149;
+    let scale = 149;
+    let result = BigDecimal::new(BigInt::from_biguint(sign, magnitude), scale);
+    return result;
+}
+
+
 #[cfg(test)]
 #[allow(non_snake_case)]
 mod test_parse_from_f32 {
@@ -112,7 +138,7 @@ pub(crate) fn try_parse_from_f64(n: f64) -> Result<BigDecimal, ParseBigDecimalEr
     match n.classify() {
         Nan => Err(ParseBigDecimalError::Other("NAN".into())),
         Infinite => Err(ParseBigDecimalError::Other("Infinite".into())),
-        Subnormal => Err(ParseBigDecimalError::Other("Subnormal".into())),
+        Subnormal => Ok(parse_from_f64_subnormal(n)),
         Normal | Zero => Ok(parse_from_f64(n)),
     }
 }
@@ -141,12 +167,48 @@ fn split_f64_into_parts(f: f64) -> (u64, i64, Sign) {
     (frac, pow, sign)
 }
 
+/// Create bigdecimal from subnormal f64
+pub(crate) fn parse_from_f64_subnormal(n: f64) -> BigDecimal {
+    debug_assert_eq!(n.classify(), FpCategory::Subnormal);
+    let bits = n.to_bits();
+
+    let sign_bit = bits >> 63;
+    debug_assert_eq!(bits >> 52, sign_bit << 11);
+
+    // 5^1074 = 5^1022 + 5^52  (f64-bit-bias=1022, fraction-bits=52)
+    let five_to_1074 = BigUint::from_slice(&[
+        2993937753, 2678407619, 3969251600, 2340035423,  635686544, 3544357150, 2618749834,
+        3195461310, 2593678749, 4014721034, 2512738537, 1379014958, 2606506302, 1209795638,
+        3422246832, 2235398534, 2765471138, 3453720203, 3699786234, 1752628667, 3832472493,
+        2479745915, 4210941784, 2088904316, 4137646701, 3840319652, 3815898978, 2202136831,
+        1022273801, 1470939580, 2032173740, 4063736241, 2069243191, 4077145663, 4033014231,
+        1920904652, 4195885152, 3551517817, 4246423481, 2447790869, 1797774111,   11284306,
+         195273359, 3811183395, 4065514955, 3382133286, 1078447835, 2100087074, 3915378083,
+        1127077286, 1409634978, 2331452623, 1301118814, 3692061923, 2506161869, 4270519152,
+        1066095370,  212429084, 3729063602, 3175008277, 2075072468, 2136773221, 4247151843,
+        2395660055,  449096848, 2439918400, 1564416362, 3638689409, 3054795416, 1803373736,
+        1506581328, 2791252870, 3391180271, 1768177410, 3891987426, 3655546435, 3881223940,
+         903390128
+    ]);
+
+    let frac = bits - (sign_bit << 63);
+
+    let sign = if sign_bit == 0 { Sign::Plus } else { Sign::Minus };
+    let magnitude = BigUint::from(frac) * five_to_1074;
+    let scale = 1074;
+
+    return BigDecimal::new(BigInt::from_biguint(sign, magnitude), scale);
+}
 
 /// Create bigdecimal from f64
 ///
 /// Non "normal" values is undefined behavior
 ///
 pub(crate) fn parse_from_f64(n: f64) -> BigDecimal {
+    if n.classify() == FpCategory::Subnormal {
+        return parse_from_f64_subnormal(n);
+    }
+
     let bits = n.to_bits();
 
     // shift right by 1 bit to handle -0.0
diff --git a/src/parsing.tests.parse_from_f32.rs b/src/parsing.tests.parse_from_f32.rs
index ab9fd29..7c67b9f 100644
--- a/src/parsing.tests.parse_from_f32.rs
+++ b/src/parsing.tests.parse_from_f32.rs
@@ -58,6 +58,12 @@ impl_test!(_4294967295 : 4294967295. == "4294967296");
 impl_test!(_158456325029e18 : 1.58456325029e+29 == "158456325028528675187087900672");
 
 
+
+impl_test!(_1_40129846432e_45 : 1.40129846432e-45 == "1.40129846432481707092372958328991613128026194187651577175706828388979108268586060148663818836212158203125E-45");
+impl_test!(_1_e42 : 1e-42 == "1.0005271035279193886395429224690001177341070264998322610345467546973108330377044694614596664905548095703125e-42");
+impl_test!(_3_92E_n39 : 3.92E-39 == "3.91999933059456489828739575494312783522406115751507460249208160269472102366083987590172910131514072418212890625E-39");
+impl_test!(_2_81341650018752E_n308 : 2.81341650018752E-308 == "0");
+
 #[test]
 fn case_f32_min() {
     let n = f32::MIN;
diff --git a/src/parsing.tests.parse_from_f64.rs b/src/parsing.tests.parse_from_f64.rs
index 5f8c195..d3b9aef 100644
--- a/src/parsing.tests.parse_from_f64.rs
+++ b/src/parsing.tests.parse_from_f64.rs
@@ -50,6 +50,8 @@ impl_test!(_50 : 50. == "50");
 impl_test!(_nanbits : bits:0b_0_11111111111_1000000000000000000000000000000000000000000000000001 => "269653970229347426076201969312749943170150807578117307259170330445749843759196293443300553362892619730839480672521111823337121537071529813188030913831084401350087805833926634314566788423582671529934053315387252306324360914392174188827078768228648633522131134987762597502339006422840407304422939101316534763520");
 impl_test!(_3105036184601418e246 : bits:0b_0_11100000000_0000000000000000000000000000000000000000000000000000 => "3105036184601417870297958976925005110513772034233393222278104076052101905372753772661756817657292955900975461394262146412343160088229628782888574550082362278408909952041699811100530571263196889650525998387432937501785693707632115712");
 
+impl_test!(_2_81341650018752E_308 : 2.81341650018752E-308 == "2.8134165001875198278759275525943498067505063001967969175506480744152639496835355462897889950138699429916690515722729976876607247658891051736045520063301219592298855232146428654590713004216312194773871772185068366206180596731958890086634117134422695105490626598276746331472433159429067991016548063113298957324839879447939977012897422163463450947345510093578791948321798481101505330952230105511530048812659083481787407026258844307461890753626327683153826358878159001221539330872743255707112001100520519610144879206546597846231715071742093092641158571855689231930930474890818690333095288369471228217443460522531282790309374378111440076317827545086535792316428407651758951233693496387904508572484340169054222573303301594335791590596740352481219815672375261783599853515625E-308");
+
 
 #[test]
 fn case_f64_min() {
diff --git a/src/with_std.rs b/src/with_std.rs
index 4cddf7c..ace3082 100644
--- a/src/with_std.rs
+++ b/src/with_std.rs
@@ -23,4 +23,7 @@ mod stdlib {
 
     #[cfg(test)]
     pub use std::collections::hash_map::DefaultHasher;
+
+    pub use std::vec::Vec;
+    pub use std::borrow;
 }
diff --git a/src/without_std.rs b/src/without_std.rs
index 8e35874..deac488 100644
--- a/src/without_std.rs
+++ b/src/without_std.rs
@@ -33,5 +33,7 @@ mod stdlib {
     #[cfg(test)]
     pub use siphasher::sip::SipHasher as DefaultHasher;
 
+    pub use alloc::borrow;
     pub use alloc::string;
+    pub use alloc::vec::Vec;
 }