diff --git a/numerical/bin_pow.go b/numerical/bin_pow.go
new file mode 100644
index 0000000..e0fd05d
--- /dev/null
+++ b/numerical/bin_pow.go
@@ -0,0 +1,14 @@
+package numerical
+
+// BinPow evaluates (base ^ deg) % rem
+func BinPow(base int, deg int, rem int) int {
+	var res = 1
+	for deg > 0 {
+		if (deg & 1) > 0 {
+			res = int(int64(res) * int64(base) % int64(rem))
+		}
+		base = int((int64(base) * int64(base)) % int64(rem))
+		deg >>= 1
+	}
+	return res
+}
diff --git a/numerical/bin_pow_test.go b/numerical/bin_pow_test.go
new file mode 100644
index 0000000..27fc934
--- /dev/null
+++ b/numerical/bin_pow_test.go
@@ -0,0 +1,19 @@
+package numerical
+
+import "testing"
+
+// TestBinPow tests binpow function
+func TestBinPow(t *testing.T) {
+
+	if BinPow(2, 10, 121323) != 1024 {
+		t.Error("[Error] BinPow(2, 10) is wrong")
+	}
+
+	if BinPow(1, 10, 121323) != 1 {
+		t.Error("[Error] BinPow(1, 10) is wrong")
+	}
+
+	if BinPow(0, 123123, 2) != 0 {
+		t.Error("[Error] BinPow(0, 123123) is wrong")
+	}
+}
diff --git a/numerical/factorial.go b/numerical/factorial.go
index 6536860..51cb6a7 100644
--- a/numerical/factorial.go
+++ b/numerical/factorial.go
@@ -1,15 +1,8 @@
-package main
-import "fmt"
+package numerical
 
 func factorial(num int) int {
-  if num == 0 {
-    return 1
-  }
-  return num * factorial(num - 1)
-}
-
-func main() {
-  num := 10
-  result := factorial(num)
-  fmt.Println(result)
+	if num == 0 {
+		return 1
+	}
+	return num * factorial(num-1)
 }
diff --git a/numerical/fibonacci.go b/numerical/fibonacci.go
index 047dc1a..3b61b76 100644
--- a/numerical/fibonacci.go
+++ b/numerical/fibonacci.go
@@ -1,16 +1,9 @@
-package main
-import "fmt"
+package numerical
 
 //using recursion
 func fibo(num int) int {
-  if num <= 1 {
-    return num
-  }
-  return fibo(num -1) + fibo(num - 2)
-}
-
-func main(){
-  num := 10
-  result := fibo(num)
-  fmt.Println(result)
+	if num <= 1 {
+		return num
+	}
+	return fibo(num-1) + fibo(num-2)
 }
diff --git a/numerical/gcd.go b/numerical/gcd.go
index ef69478..72af038 100644
--- a/numerical/gcd.go
+++ b/numerical/gcd.go
@@ -1,47 +1,48 @@
-package gcd
+package numerical
 
-func gcd(x uint, y uint) uint {
-    
-    var shift uint = 0
+// GCD returns gcd of x and y
+func GCD(x uint, y uint) uint {
 
-    if x == y {
-        return x
-    }
+	var shift uint = 0
 
-    if x == 0 {
-        return y
-    }
+	if x == y {
+		return x
+	}
 
-    if y == 0 {
-        return x
-    }
+	if x == 0 {
+		return y
+	}
 
-    for shift := 0; (x | y) & 1 == 0; shift++ {
-        x = x >> 1
-        y = y >> 1
-    }
+	if y == 0 {
+		return x
+	}
 
-    for ; (x & 1) == 0 ; {
-        x = x >> 1
-    }
+	for shift := 0; (x|y)&1 == 0; shift++ {
+		x = x >> 1
+		y = y >> 1
+	}
 
-    for ; y == 0 ; {
-        
-        for ; (y & 1) == 0 ; {
-            y = y >> 1
-        }
+	for (x & 1) == 0 {
+		x = x >> 1
+	}
 
-        if x > y {
-            t := x
-            x = y
-            y = t
-        }
+	for y == 0 {
 
-        y = y - x
+		for (y & 1) == 0 {
+			y = y >> 1
+		}
 
-    }
+		if x > y {
+			t := x
+			x = y
+			y = t
+		}
 
-    y = y << shift
+		y = y - x
 
-    return y
+	}
+
+	y = y << shift
+
+	return y
 }
diff --git a/numerical/gcd_test.go b/numerical/gcd_test.go
index 4cc4b59..b7883ee 100644
--- a/numerical/gcd_test.go
+++ b/numerical/gcd_test.go
@@ -1,18 +1,19 @@
-package gcd
+package numerical
 
 import "testing"
 
-func Test_gcd(t *testing.T) {
+// TestGcd tests gcd
+func TestGcd(t *testing.T) {
 
-	if gcd(100, 200) != 50 {
-		t.Error("[Error] gcd(100, 200) is wrong")
+	if GCD(100, 200) != 50 {
+		t.Error("[Error] GCD(100, 200) is wrong")
 	}
 
-	if gcd(4, 2) != 1 {
-		t.Error("[Error] gcd(4,2) is wrong")
+	if GCD(4, 2) != 1 {
+		t.Error("[Error] GCD(4,2) is wrong")
 	}
 
-	if gcd(6, 3) != 3 {
-		t.Error("[Error] gcd(6,3) is wrong")
+	if GCD(6, 3) != 3 {
+		t.Error("[Error] GCD(6,3) is wrong")
 	}
 }
diff --git a/numerical/prime_finder.go b/numerical/prime_finder.go
new file mode 100644
index 0000000..7594407
--- /dev/null
+++ b/numerical/prime_finder.go
@@ -0,0 +1,25 @@
+package numerical
+
+// PrimesUpTo finds all prime numbers from 1 to upperBound
+// It's implemented using eratosthenes sieve
+// Works in O(upperBound) time and space
+func PrimesUpTo(upperBound int) []int {
+	// lp array stores minimal prime divisor for every number from 2 to upperBound
+	var lp []int
+	lp = make([]int, upperBound+1)
+
+	// primes array stores primes
+	var primes []int
+
+	for i := 2; i <= upperBound; i++ {
+		if lp[i] == 0 {
+			lp[i] = i
+			primes = append(primes, i)
+		}
+		for j := 0; j < len(primes) && primes[j] <= lp[i] && i*primes[j] <= upperBound; j++ {
+			lp[i*primes[j]] = primes[j]
+		}
+	}
+
+	return primes
+}
diff --git a/numerical/prime_finder_test.go b/numerical/prime_finder_test.go
new file mode 100644
index 0000000..b7770bf
--- /dev/null
+++ b/numerical/prime_finder_test.go
@@ -0,0 +1,26 @@
+package numerical
+
+import "testing"
+
+func arrayEquals(a, b []int) bool {
+	if len(a) != len(b) {
+		return false
+	}
+	for i, v := range a {
+		if v != b[i] {
+			return false
+		}
+	}
+	return true
+}
+
+// TestPrimeFinder tests prime finding
+func TestPrimeFinder(t *testing.T) {
+
+	if !arrayEquals(PrimesUpTo(10), []int{2, 3, 5, 7}) {
+		t.Error("[Error] PrimesUpTo(10) is wrong")
+	}
+	if len(PrimesUpTo(100)) != 25 {
+		t.Error("[Error] PrimesUpTo(100) is wrong")
+	}
+}