--- a/rust/fpnum/src/lib.rs	Thu Jul 04 21:40:50 2019 +0300
+++ b/rust/fpnum/src/lib.rs	Fri Jul 05 21:16:33 2019 +0300
@@ -73,24 +73,9 @@
     pub fn sqrt(&self) -> Self {
         debug_assert!(self.is_positive());
 
-        let mut t: u64 = 0x4000_0000_0000_0000;
-        let mut r: u64 = 0;
-        let mut q = self.value;
-
-        for _ in 0..32 {
-            let s = r + t;
-            r >>= 1;
-
-            if s <= q {
-                q -= s;
-                r += t;
-            }
-            t >>= 2;
-        }
-
         Self {
             sign_mask: POSITIVE_MASK,
-            value: r << 16,
+            value: integral_sqrt(self.value) << 16,
         }
     }
 
@@ -354,25 +339,11 @@
         if r < LINEARIZE_TRESHOLD {
             FPNum::from(r as u32)
         } else {
-            let mut sqr: u128 = (self.x_value as u128).pow(2) + (self.y_value as u128).pow(2);
-
-            let mut t: u128 = 0x40000000_00000000_00000000_00000000;
-            let mut r: u128 = 0;
-
-            for _ in 0..64 {
-                let s = r + t;
-                r >>= 1;
-
-                if s <= sqr {
-                    sqr -= s;
-                    r += t;
-                }
-                t >>= 2;
-            }
+            let sqr: u128 = (self.x_value as u128).pow(2) + (self.y_value as u128).pow(2);
 
             FPNum {
                 sign_mask: POSITIVE_MASK,
-                value: r as u64,
+                value: integral_sqrt_ext(sqr) as u64,
             }
         }
     }
@@ -500,29 +471,50 @@
 bin_assign_op_impl!(FPPoint, ops::MulAssign<FPNum>, mul_assign, *);
 bin_assign_op_impl!(FPPoint, ops::DivAssign<FPNum>, div_assign, /);
 
+pub fn integral_sqrt(mut value: u64) -> u64 {
+    let mut digit_sqr = 0x4000_0000_0000_0000u64.wrapping_shr(value.leading_zeros() & 0xFFFF_FFFE);
+    let mut result = 0;
+
+    while digit_sqr != 0 {
+        let approx = digit_sqr + result;
+        result >>= 1;
+
+        if approx <= value {
+            value -= approx;
+            result += digit_sqr;
+        }
+        digit_sqr >>= 2;
+    }
+    result
+}
+
+pub fn integral_sqrt_ext(mut value: u128) -> u128 {
+    let mut digit_sqr =
+        0x40000000_00000000_00000000_00000000u128.wrapping_shr(value.leading_zeros() & 0xFFFF_FFFE);
+    let mut result = 0u128;
+
+    while digit_sqr != 0 {
+        let approx = result + digit_sqr;
+        result >>= 1;
+
+        if approx <= value {
+            value -= approx;
+            result += digit_sqr;
+        }
+        digit_sqr >>= 2;
+    }
+    result
+}
+
 pub fn distance<T>(x: T, y: T) -> FPNum
 where
     T: Into<i64> + std::fmt::Debug,
 {
-    let mut sqr: u128 = (x.into().pow(2) as u128).shl(64) + (y.into().pow(2) as u128).shl(64);
-
-    let mut t: u128 = 0x40000000_00000000_00000000_00000000;
-    let mut r: u128 = 0;
-
-    for _ in 0..64 {
-        let s = r + t;
-        r >>= 1;
-
-        if s <= sqr {
-            sqr -= s;
-            r += t;
-        }
-        t >>= 2;
-    }
+    let sqr: u128 = (x.into().pow(2) as u128).shl(64) + (y.into().pow(2) as u128).shl(64);
 
     FPNum {
         sign_mask: POSITIVE_MASK,
-        value: r as u64,
+        value: integral_sqrt_ext(sqr) as u64,
     }
 }
 

--- a/rust/integral-geometry/src/lib.rs	Thu Jul 04 21:40:50 2019 +0300
+++ b/rust/integral-geometry/src/lib.rs	Fri Jul 05 21:16:33 2019 +0300
@@ -1,4 +1,4 @@
-use fpnum::{distance, fp, FPNum, FPPoint};
+use fpnum::{fp, integral_sqrt, FPNum, FPPoint};
 use std::{
     cmp::{max, min},
     ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, RangeInclusive, Sub, SubAssign},
@@ -45,7 +45,8 @@
 
     #[inline]
     pub fn integral_norm(self) -> u32 {
-        distance(self.x, self.y).abs_round()
+        let sqr = (self.x as u64).pow(2) + (self.y as u64).pow(2);
+        integral_sqrt(sqr) as u32
     }
 
     #[inline]