diff --git a/fix16.tal b/fix16.tal
index 2a9fe42..cf2850b 100644
--- a/fix16.tal
+++ b/fix16.tal
@@ -258,6 +258,23 @@
     &odd #0080
     &rest ADD2 DUP EOR JMP2r
 
+@x16-sqrt ( x* -> sqrt(x)* )
+    LIT2r ff00
+    LIT2r 0200                  ( [c* 2*] )
+    DUP2 STH2kr x16-div         ( x* s=x/2* [c* 2*] )
+    &loop                       ( x* s0* [c* 2*] )
+        OVR2 OVR2 x16-div       ( x* s0* x/s0* [c* 2*] )
+        OVR2 x16-add            ( x* s0* (x/s0)+s0* [c* 2*] )
+        STH2kr x16-div          ( x* s0* s1=((x/s0)+s0)/2* [c* 2*] )
+        SWP2 OVR2 x16-ne        ( x* s1* go^ [c* 2*] )
+        SWP2r INC2r ORAkr STHr  ( x* s1* go^ ok^ [2* c+1*] )
+        SWP2r ?&continue        ( x* s1* go^ [c+1* 2*] )
+        POP !&done              ( x* s1* [c+1* 2*] )
+    &continue                   ( x* s1* go^ [c+1* 2*] )
+        ?&loop                  ( x* s1* [c+1* 2*] )
+    &done                       ( x* s1* [c* 2*] )
+        POP2r POP2r NIP2 JMP2r  ( s1* )
+
 @x16-cos ( x* -> cos(x)* )
     x16-pi/2 ADD2 ,x16-sin JMP
 
diff --git a/test-fix16.py b/test-fix16.py
index e76e25e..18be196 100644
--- a/test-fix16.py
+++ b/test-fix16.py
@@ -1,6 +1,6 @@
 #!/usr/bin/python
 
-from math import ceil, copysign, cos, floor, log, sin, tan
+from math import ceil, copysign, cos, floor, log, sin, sqrt, tan
 from os import environ
 from random import randint
 from subprocess import Popen, PIPE, run
@@ -113,6 +113,8 @@ def x16_lteq(x, y):
     return int(tosigned(x) <= tosigned(y))
 def x16_gteq(x, y):
     return int(tosigned(x) >= tosigned(y))
+def x16_sqrt(x):
+    return int(sqrt(x / 256) * 256)
 def x16_sin(x):
     z = round(sin(x / 256) * 256)
     return z if z >= 0 else 65536 + z
@@ -162,6 +164,7 @@ def main():
     test(p, trials, b'>', [('x', u16), ('y', u16)], u8, x16_gt)
     test(p, trials, b'{', [('x', u16), ('y', u16)], u8, x16_lteq)
     test(p, trials, b'}', [('x', u16), ('y', u16)], u8, x16_gteq)
+    test(p, trials, b'r', [('x', p16)], u16, x16_sqrt, eq=releq)
     test(p, trials, b's', [('x', p16)], u16, x16_sin, eq=abseq)
     test(p, trials, b'c', [('x', p16)], u16, x16_cos, eq=abseq)
     test(p, trials, b't', [('x', t16)], u16, x16_tan, eq=taneq)
diff --git a/test-fix16.tal b/test-fix16.tal
index d6ca8b8..536252e 100644
--- a/test-fix16.tal
+++ b/test-fix16.tal
@@ -64,6 +64,7 @@
     ;buf LDA LIT "> EQU ;test-x16-gt JCN2
     ;buf LDA LIT "{ EQU ;test-x16-lteq JCN2
     ;buf LDA LIT "} EQU ;test-x16-gteq JCN2
+    ;buf LDA LIT "r EQU ;test-x16-sqrt JCN2
     ;buf LDA LIT "s EQU ;test-x16-sin JCN2
     ;buf LDA LIT "c EQU ;test-x16-cos JCN2
     ;buf LDA LIT "t EQU ;test-x16-tan JCN2
@@ -113,6 +114,7 @@
 @test-x16-lteq       #04 #01  ;x16-lteq            ;test JMP2
 @test-x16-gt         #04 #01  ;x16-gt              ;test JMP2
 @test-x16-gteq       #04 #01  ;x16-gteq            ;test JMP2
+@test-x16-sqrt       #02 #02  ;x16-sqrt            ;test JMP2
 @test-x16-sin        #02 #02  ;x16-sin             ;test JMP2
 @test-x16-cos        #02 #02  ;x16-cos             ;test JMP2
 @test-x16-tan        #02 #02  ;x16-tan             ;test JMP2