diff --git a/storm/Big.cpp b/storm/Big.cpp
index a8d7804..3b6731c 100644
--- a/storm/Big.cpp
+++ b/storm/Big.cpp
@@ -49,6 +49,10 @@ void SBigMul(BigData* a, BigData* b, BigData* c) {
     Mul(a->Primary(), b->Primary(), c->Primary(), a->Stack());
 }
 
+void SBigSquare(BigData* a, BigData* b) {
+    Square(a->Primary(), b->Primary(), a->Stack());
+}
+
 void SBigSub(BigData* a, BigData* b, BigData* c) {
     Sub(a->Primary(), b->Primary(), c->Primary());
 }
diff --git a/storm/Big.hpp b/storm/Big.hpp
index 291eea2..c12c1a6 100644
--- a/storm/Big.hpp
+++ b/storm/Big.hpp
@@ -20,6 +20,8 @@ void SBigMul(BigData* a, BigData* b, BigData* c);
 
 void SBigNew(BigData** num);
 
+void SBigSquare(BigData* a, BigData* b);
+
 void SBigSub(BigData* a, BigData* b, BigData* c);
 
 void SBigToBinaryBuffer(BigData* num, uint8_t* data, uint32_t maxBytes, uint32_t* bytes);
diff --git a/storm/big/Ops.cpp b/storm/big/Ops.cpp
index 2affe61..725a7e3 100644
--- a/storm/big/Ops.cpp
+++ b/storm/big/Ops.cpp
@@ -170,6 +170,32 @@ void SetZero(BigBuffer& buffer) {
     buffer.Clear();
 }
 
+void Square(BigBuffer& a, const BigBuffer& b, BigStack& stack) {
+    auto& aa = stack.MakeDistinct(a, &a == &b);
+    aa.Clear();
+
+    uint32_t i = 0;
+    for (i = 0; b.IsUsed(i); i++) {
+        uint64_t carry = 0;
+
+        uint32_t j = 0;
+        for (j = 0; j <= i; j++) {
+            auto mul = b[i] * static_cast<uint64_t>(b[j]);
+            auto add = mul + aa[i + j];
+
+            if (j < i) {
+                carry += mul;
+            }
+
+            aa[i + j] = ExtractLowPartLargeSum(carry, add);
+        }
+
+        aa[i + j] = ExtractLowPart(carry);
+    }
+
+    stack.UnmakeDistinct(a, aa);
+}
+
 void Sub(BigBuffer& a, const BigBuffer& b, const BigBuffer& c) {
     uint64_t borrow = 0;
     uint32_t i = 0;
diff --git a/storm/big/Ops.hpp b/storm/big/Ops.hpp
index 94b502b..7d46b91 100644
--- a/storm/big/Ops.hpp
+++ b/storm/big/Ops.hpp
@@ -37,6 +37,8 @@ void SetOne(BigBuffer& buffer);
 
 void SetZero(BigBuffer& buffer);
 
+void Square(BigBuffer& a, const BigBuffer& b, BigStack& stack);
+
 void Sub(BigBuffer& a, const BigBuffer& b, const BigBuffer& c);
 
 void ToBinary(TSGrowableArray<uint8_t>& output, const BigBuffer& buffer);
diff --git a/test/Big.cpp b/test/Big.cpp
index e22c214..2e66e86 100644
--- a/test/Big.cpp
+++ b/test/Big.cpp
@@ -855,6 +855,28 @@ TEST_CASE("SBigMul", "[big]") {
     }
 }
 
+TEST_CASE("SBigSquare", "[square]") {
+    SECTION("squares 0xFFFFFFFF") {
+        BigData* a;
+        SBigNew(&a);
+
+        BigData* b;
+        SBigNew(&b);
+        SBigFromUnsigned(b, 0xFFFFFFFF);
+
+        SBigSquare(a, b);
+
+        a->Primary().Trim();
+
+        CHECK(a->Primary().Count() == 2);
+        CHECK(a->Primary()[0] == 0x1);
+        CHECK(a->Primary()[1] == 0xFFFFFFFE);
+
+        SBigDel(a);
+        SBigDel(b);
+    }
+}
+
 TEST_CASE("SBigSub", "[big]") {
     SECTION("subtracts 1 from 2") {
         BigData* a;
@@ -961,6 +983,126 @@ TEST_CASE("SetZero", "[big]") {
     }
 }
 
+TEST_CASE("Square", "[big]") {
+    SECTION("squares 0") {
+        BigData* a;
+        SBigNew(&a);
+
+        BigData* b;
+        SBigNew(&b);
+        SBigFromUnsigned(b, 0);
+
+        Square(a->Primary(), b->Primary(), a->Stack());
+
+        a->Primary().Trim();
+
+        CHECK(a->Primary().Count() == 0);
+
+        SBigDel(a);
+        SBigDel(b);
+    }
+
+    SECTION("squares 1") {
+        BigData* a;
+        SBigNew(&a);
+
+        BigData* b;
+        SBigNew(&b);
+        SBigFromUnsigned(b, 1);
+
+        Square(a->Primary(), b->Primary(), a->Stack());
+
+        a->Primary().Trim();
+
+        CHECK(a->Primary().Count() == 1);
+        CHECK(a->Primary()[0] == 1);
+
+        SBigDel(a);
+        SBigDel(b);
+    }
+
+    SECTION("squares 2") {
+        BigData* a;
+        SBigNew(&a);
+
+        BigData* b;
+        SBigNew(&b);
+        SBigFromUnsigned(b, 2);
+
+        Square(a->Primary(), b->Primary(), a->Stack());
+
+        a->Primary().Trim();
+
+        CHECK(a->Primary().Count() == 1);
+        CHECK(a->Primary()[0] == 4);
+
+        SBigDel(a);
+        SBigDel(b);
+    }
+
+    SECTION("squares 2") {
+        BigData* a;
+        SBigNew(&a);
+
+        BigData* b;
+        SBigNew(&b);
+        SBigFromUnsigned(b, 2);
+
+        Square(a->Primary(), b->Primary(), a->Stack());
+
+        a->Primary().Trim();
+
+        CHECK(a->Primary().Count() == 1);
+        CHECK(a->Primary()[0] == 4);
+
+        SBigDel(a);
+        SBigDel(b);
+    }
+
+    SECTION("squares 0xFFFFFFFF") {
+        BigData* a;
+        SBigNew(&a);
+
+        BigData* b;
+        SBigNew(&b);
+        SBigFromUnsigned(b, 0xFFFFFFFF);
+
+        Square(a->Primary(), b->Primary(), a->Stack());
+
+        a->Primary().Trim();
+
+        CHECK(a->Primary().Count() == 2);
+        CHECK(a->Primary()[0] == 0x1);
+        CHECK(a->Primary()[1] == 0xFFFFFFFE);
+
+        SBigDel(a);
+        SBigDel(b);
+    }
+
+    SECTION("squares 0x1111111111111111") {
+        BigData* a;
+        SBigNew(&a);
+
+        BigData* b;
+        SBigNew(&b);
+        uint64_t b_ = 0x1111111111111111;
+        SBigFromBinary(b, reinterpret_cast<uint8_t*>(&b_), sizeof(b_));
+
+        Square(a->Primary(), b->Primary(), a->Stack());
+
+        a->Primary().Trim();
+
+        CHECK(a->Primary().Count() == 4);
+        CHECK(a->Primary()[0] == 0x87654321);
+        CHECK(a->Primary()[1] == 0xfedcba9);
+        CHECK(a->Primary()[2] == 0x89abcdf0);
+        CHECK(a->Primary()[3] == 0x1234567);
+
+        SBigDel(a);
+        SBigDel(b);
+    }
+}
+
 TEST_CASE("Sub", "[big]") {
     SECTION("subtracts 0 from 1") {
         BigData* a;