diff --git a/gpt2.f90 b/gpt2.f90
index 6b3f9b5..bfcaa61 100644
--- a/gpt2.f90
+++ b/gpt2.f90
@@ -7,9 +7,28 @@ module gpt2_mod
 
 contains
 
+elemental real(sp) function fast_tanh(x) result(y)
+real(sp), intent(in) :: x
+real(sp) :: x2
+x2 = x*x
+y = x * (0.98569772605911309407 + x2 *(-0.2794500993392901382 &
+    + x2 * (6.8280504526399188164e-2 + x2 * (-1.0972014877337651823e-2 &
+    + x2 * (1.1132367134444316902e-3 + x2 * (-7.018851897305717565e-5 &
+    + x2 * (2.656616768082727089e-6 + x2 * (-5.5138381821615909058e-8 &
+    + x2 * 4.8162484477588665996e-10))))))))
+if (x > 5) then
+    y = 1
+    ! Does not work yet, but would eliminate the if condition
+    !y = 0.5_sp * (1-sign(1.0, x-5)) * (y-1) + 1
+end if
+if (x < -5) then
+    y = -1
+end if
+end function
+
 elemental real(sp) function gelu(x) result(y)
 real(sp), intent(in) :: x
-y = 0.5_sp * x * (1 + tanh(sqrt(2 / pi) * (x + 0.044715_sp * x**3)))
+y = 0.5_sp * x * (1 + fast_tanh(sqrt(2 / pi) * (x + 0.044715_sp * x**3)))
 end function
 
 function softmax(x) result(y)
@@ -56,9 +75,16 @@ function linear(x, w, b) result(y)
 function ffn(x, fc_w, fc_b, proj_w, proj_b) result(y)
 real(sp), intent(in) :: x(:,:), fc_w(:,:), fc_b(:), proj_w(:,:), proj_b(:)
 real(sp) :: y(size(x,1),size(x,2))
-!real(sp) :: a(4*size(x,1),size(x,2))
-!a = gelu(linear(x, fc_w, fc_b))
-y = linear(gelu(linear(x, fc_w, fc_b)), proj_w, proj_b)
+real(sp) :: a(4*size(x,1),size(x,2))
+integer :: i, j
+!a = linear(x, fc_w, fc_b)
+call matmul_2d(fc_w, x, a)
+do j = 1, size(a,1)
+do i = 1, size(a,2)
+    a(j,i) = gelu(a(j,i) + fc_b(j))
+end do
+end do
+y = linear(a, proj_w, proj_b)
 end function
 
 function attention(q, k, v, mask) result(y)