diff --git a/MANIFEST.in b/MANIFEST.in
index 6a9fa8a4..06b0e1d3 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -4,4 +4,5 @@ recursive-include m2cgen linear_algebra.*
 recursive-include m2cgen log1p.*
 recursive-include m2cgen tanh.*
 recursive-include m2cgen atan.*
+recursive-include m2cgen softmax.*
 global-exclude *.py[cod]
diff --git a/m2cgen/assemblers/boosting.py b/m2cgen/assemblers/boosting.py
index 99be8267..0469981e 100644
--- a/m2cgen/assemblers/boosting.py
+++ b/m2cgen/assemblers/boosting.py
@@ -72,8 +72,7 @@ def _assemble_multi_class_output(self, estimator_params):
             for i, e in enumerate(splits)
         ]
 
-        proba_exprs = self._multi_class_convert_output(exprs)
-        return ast.VectorVal(proba_exprs)
+        return self._multi_class_convert_output(exprs)
 
     def _assemble_bin_class_output(self, estimator_params):
         # Base score is calculated based on
@@ -97,7 +96,7 @@ def _final_transform(self, ast_to_transform):
         return ast_to_transform
 
     def _multi_class_convert_output(self, exprs):
-        return fallback_expressions.softmax(exprs)
+        return ast.SoftmaxExpr(exprs)
 
     def _bin_class_convert_output(self, expr, to_reuse=True):
         return fallback_expressions.sigmoid(expr, to_reuse=to_reuse)
@@ -250,8 +249,9 @@ def _multi_class_convert_output(self, exprs):
         return supported_objectives[self.objective_name](exprs)
 
     def _multi_class_sigmoid_transform(self, exprs):
-        return [self._bin_class_sigmoid_transform(expr, to_reuse=False)
-                for expr in exprs]
+        return ast.VectorVal([
+            self._bin_class_sigmoid_transform(expr, to_reuse=False)
+            for expr in exprs])
 
     def _bin_class_convert_output(self, expr, to_reuse=True):
         supported_objectives = {
diff --git a/m2cgen/assemblers/fallback_expressions.py b/m2cgen/assemblers/fallback_expressions.py
index e42329ed..56230484 100644
--- a/m2cgen/assemblers/fallback_expressions.py
+++ b/m2cgen/assemblers/fallback_expressions.py
@@ -190,7 +190,7 @@ def softmax(exprs):
     exp_exprs = [ast.ExpExpr(e, to_reuse=True) for e in exprs]
     exp_sum_expr = utils.apply_op_to_expressions(
         ast.BinNumOpType.ADD, *exp_exprs, to_reuse=True)
-    return [
+    return ast.VectorVal([
         ast.BinNumExpr(e, exp_sum_expr, ast.BinNumOpType.DIV)
         for e in exp_exprs
-    ]
+    ])
diff --git a/m2cgen/ast.py b/m2cgen/ast.py
index 4c167677..87f090e8 100644
--- a/m2cgen/ast.py
+++ b/m2cgen/ast.py
@@ -266,6 +266,29 @@ def __hash__(self):
         return hash(tuple(self.exprs))
 
 
+class SoftmaxExpr(VectorExpr):
+
+    def __init__(self, exprs, to_reuse=False):
+        assert all(e.output_size == 1 for e in exprs), (
+            "All expressions for SoftmaxExpr must be scalar")
+
+        self.exprs = exprs
+        self.to_reuse = to_reuse
+        self.output_size = len(exprs)
+
+    def __str__(self):
+        args = ",".join([str(e) for e in self.exprs])
+        return f"SoftmaxExpr({args},to_reuse={self.to_reuse})"
+
+    def __eq__(self, other):
+        return (type(other) is SoftmaxExpr and
+                self.output_size == other.output_size and
+                all(i == j for i, j in zip(self.exprs, other.exprs)))
+
+    def __hash__(self):
+        return hash(tuple(self.exprs))
+
+
 class BinVectorExpr(VectorExpr, BinExpr):
 
     def __init__(self, left, right, op):
@@ -384,9 +407,9 @@ def __hash__(self):
 
 NESTED_EXPRS_MAPPINGS = [
     ((BinExpr, CompExpr), lambda e: [e.left, e.right]),
-    (PowExpr, lambda e: [e.base_expr, e.exp_expr]),
-    (VectorVal, lambda e: e.exprs),
-    (IfExpr, lambda e: [e.test, e.body, e.orelse]),
+    ((PowExpr), lambda e: [e.base_expr, e.exp_expr]),
+    ((VectorVal, SoftmaxExpr), lambda e: e.exprs),
+    ((IfExpr), lambda e: [e.test, e.body, e.orelse]),
     ((AbsExpr, AtanExpr, ExpExpr, IdExpr, LogExpr, Log1pExpr,
       SqrtExpr, TanhExpr),
      lambda e: [e.expr]),
diff --git a/m2cgen/interpreters/c/interpreter.py b/m2cgen/interpreters/c/interpreter.py
index a8e55c8a..9f5aa00a 100644
--- a/m2cgen/interpreters/c/interpreter.py
+++ b/m2cgen/interpreters/c/interpreter.py
@@ -23,9 +23,12 @@ class CInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "log"
     log1p_function_name = "log1p"
     power_function_name = "pow"
+    softmax_function_name = "softmax"
     sqrt_function_name = "sqrt"
     tanh_function_name = "tanh"
 
+    with_softmax_expr = False
+
     def __init__(self, indent=4, function_name="score", *args, **kwargs):
         self.function_name = function_name
 
@@ -61,6 +64,11 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "linear_algebra.c")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.c")
+            self._cg.prepend_code_lines(utils.get_file_content(filename))
+
         if self.with_vectors:
             self._cg.add_dependency("<string.h>")
 
@@ -99,3 +107,18 @@ def interpret_bin_vector_num_expr(self, expr, **kwargs):
         self._cg.add_code_line(f"{func_inv};")
 
         return var_name
+
+    # Do the same things for softmax as for linear algebra.
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_vectors = True
+        self.with_softmax_expr = True
+
+        var_name = self._cg.add_var_declaration(expr.output_size)
+        nested = [self._do_interpret(expr, **kwargs) for expr in expr.exprs]
+        func_inv = self._cg.function_invocation(
+            self.softmax_function_name,
+            self._cg.vector_init(nested),
+            expr.output_size,
+            var_name)
+        self._cg.add_code_line(f"{func_inv};")
+        return var_name
diff --git a/m2cgen/interpreters/c_sharp/interpreter.py b/m2cgen/interpreters/c_sharp/interpreter.py
index 7a8b721c..17fc3363 100644
--- a/m2cgen/interpreters/c_sharp/interpreter.py
+++ b/m2cgen/interpreters/c_sharp/interpreter.py
@@ -24,10 +24,12 @@ class CSharpInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "Log"
     log1p_function_name = "Log1p"
     power_function_name = "Pow"
+    softmax_function_name = "Softmax"
     sqrt_function_name = "Sqrt"
     tanh_function_name = "Tanh"
 
     with_log1p_expr = False
+    with_softmax_expr = False
 
     def __init__(self, namespace="ML", class_name="Model", indent=4,
                  function_name="Score", *args, **kwargs):
@@ -66,6 +68,11 @@ def interpret(self, expr):
                         os.path.dirname(__file__), "log1p.cs")
                     self._cg.add_code_lines(utils.get_file_content(filename))
 
+                if self.with_softmax_expr:
+                    filename = os.path.join(
+                        os.path.dirname(__file__), "softmax.cs")
+                    self._cg.add_code_lines(utils.get_file_content(filename))
+
         if self.with_math_module:
             self._cg.add_dependency("System.Math")
 
@@ -74,3 +81,7 @@ def interpret(self, expr):
     def interpret_log1p_expr(self, expr, **kwargs):
         self.with_log1p_expr = True
         return super().interpret_log1p_expr(expr, **kwargs)
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/dart/interpreter.py b/m2cgen/interpreters/dart/interpreter.py
index 1ed22067..e72d8f11 100644
--- a/m2cgen/interpreters/dart/interpreter.py
+++ b/m2cgen/interpreters/dart/interpreter.py
@@ -27,10 +27,12 @@ class DartInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "log"
     log1p_function_name = "log1p"
     power_function_name = "pow"
+    softmax_function_name = "softmax"
     sqrt_function_name = "sqrt"
     tanh_function_name = "tanh"
 
     with_log1p_expr = False
+    with_softmax_expr = False
     with_tanh_expr = False
 
     def __init__(self, indent=4, function_name="score", *args, **kwargs):
@@ -63,6 +65,11 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "log1p.dart")
             self._cg.add_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.dart")
+            self._cg.add_code_lines(utils.get_file_content(filename))
+
         if self.with_tanh_expr:
             filename = os.path.join(
                 os.path.dirname(__file__), "tanh.dart")
@@ -86,3 +93,7 @@ def interpret_log1p_expr(self, expr, **kwargs):
     def interpret_tanh_expr(self, expr, **kwargs):
         self.with_tanh_expr = True
         return super().interpret_tanh_expr(expr, **kwargs)
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/f_sharp/interpreter.py b/m2cgen/interpreters/f_sharp/interpreter.py
index 89300078..e9250d7b 100644
--- a/m2cgen/interpreters/f_sharp/interpreter.py
+++ b/m2cgen/interpreters/f_sharp/interpreter.py
@@ -30,10 +30,12 @@ class FSharpInterpreter(FunctionalToCodeInterpreter,
     exponent_function_name = "exp"
     logarithm_function_name = "log"
     log1p_function_name = "log1p"
+    softmax_function_name = "softmax"
     sqrt_function_name = "sqrt"
     tanh_function_name = "tanh"
 
     with_log1p_expr = False
+    with_softmax_expr = False
 
     def __init__(self, indent=4, function_name="score", *args, **kwargs):
         self.indent = indent
@@ -63,6 +65,11 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "log1p.fs")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.fs")
+            self._cg.prepend_code_lines(utils.get_file_content(filename))
+
         return self._cg.finalize_and_get_generated_code()
 
     def create_code_generator(self):
@@ -78,6 +85,10 @@ def interpret_log1p_expr(self, expr, **kwargs):
         self.with_log1p_expr = True
         return super().interpret_log1p_expr(expr, **kwargs)
 
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
+
     def _dump_cache(self):
         if self._cached_expr_results:
             for func_name, expr_result in self._cached_expr_results.values():
diff --git a/m2cgen/interpreters/go/interpreter.py b/m2cgen/interpreters/go/interpreter.py
index a325d4db..21663ee8 100644
--- a/m2cgen/interpreters/go/interpreter.py
+++ b/m2cgen/interpreters/go/interpreter.py
@@ -22,9 +22,12 @@ class GoInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "math.Log"
     log1p_function_name = "math.Log1p"
     power_function_name = "math.Pow"
+    softmax_function_name = "softmax"
     sqrt_function_name = "math.Sqrt"
     tanh_function_name = "math.Tanh"
 
+    with_softmax_expr = False
+
     def __init__(self, indent=4, function_name="score", *args, **kwargs):
         self.function_name = function_name
 
@@ -51,7 +54,16 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "linear_algebra.go")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.go")
+            self._cg.prepend_code_lines(utils.get_file_content(filename))
+
         if self.with_math_module:
             self._cg.add_dependency("math")
 
         return self._cg.finalize_and_get_generated_code()
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/haskell/interpreter.py b/m2cgen/interpreters/haskell/interpreter.py
index 6a435418..6c6abee4 100644
--- a/m2cgen/interpreters/haskell/interpreter.py
+++ b/m2cgen/interpreters/haskell/interpreter.py
@@ -21,10 +21,12 @@ class HaskellInterpreter(FunctionalToCodeInterpreter,
     exponent_function_name = "exp"
     logarithm_function_name = "log"
     log1p_function_name = "log1p"
+    softmax_function_name = "softmax"
     sqrt_function_name = "sqrt"
     tanh_function_name = "tanh"
 
     with_log1p_expr = False
+    with_softmax_expr = False
 
     def __init__(self,  module_name="Model", indent=4, function_name="score",
                  *args, **kwargs):
@@ -59,6 +61,11 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "log1p.hs")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.hs")
+            self._cg.prepend_code_lines(utils.get_file_content(filename))
+
         self._cg.prepend_code_line(self._cg.tpl_module_definition(
             module_name=self.module_name))
 
@@ -77,6 +84,10 @@ def interpret_log1p_expr(self, expr, **kwargs):
         self.with_log1p_expr = True
         return super().interpret_log1p_expr(expr, **kwargs)
 
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
+
     def _dump_cache(self):
         if self._cached_expr_results:
             self._cg.add_code_line("where")
diff --git a/m2cgen/interpreters/interpreter.py b/m2cgen/interpreters/interpreter.py
index fd684673..052263ba 100644
--- a/m2cgen/interpreters/interpreter.py
+++ b/m2cgen/interpreters/interpreter.py
@@ -87,6 +87,7 @@ class ToCodeInterpreter(BaseToCodeInterpreter):
     logarithm_function_name = NotImplemented
     log1p_function_name = NotImplemented
     power_function_name = NotImplemented
+    softmax_function_name = NotImplemented
     sqrt_function_name = NotImplemented
     tanh_function_name = NotImplemented
 
@@ -169,6 +170,17 @@ def interpret_log1p_expr(self, expr, **kwargs):
         return self._cg.function_invocation(
             self.log1p_function_name, nested_result)
 
+    def interpret_softmax_expr(self, expr, **kwargs):
+        if self.softmax_function_name is NotImplemented:
+            return self._do_interpret(
+                fallback_expressions.softmax(expr.exprs),
+                **kwargs)
+        self.with_vectors = True
+        self.with_math_module = True
+        nested = [self._do_interpret(expr, **kwargs) for expr in expr.exprs]
+        return self._cg.function_invocation(
+            self.softmax_function_name, self._cg.vector_init(nested))
+
     def interpret_sqrt_expr(self, expr, **kwargs):
         if self.sqrt_function_name is NotImplemented:
             return self._do_interpret(
diff --git a/m2cgen/interpreters/java/interpreter.py b/m2cgen/interpreters/java/interpreter.py
index c67532dd..45762949 100644
--- a/m2cgen/interpreters/java/interpreter.py
+++ b/m2cgen/interpreters/java/interpreter.py
@@ -30,9 +30,12 @@ class JavaInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "Math.log"
     log1p_function_name = "Math.log1p"
     power_function_name = "Math.pow"
+    softmax_function_name = "softmax"
     sqrt_function_name = "Math.sqrt"
     tanh_function_name = "Math.tanh"
 
+    with_softmax_expr = False
+
     def __init__(self, package_name=None, class_name="Model", indent=4,
                  function_name="score", *args, **kwargs):
         self.package_name = package_name
@@ -64,9 +67,18 @@ def interpret(self, expr):
                     os.path.dirname(__file__), "linear_algebra.java")
                 top_cg.add_code_lines(utils.get_file_content(filename))
 
+            if self.with_softmax_expr:
+                filename = os.path.join(
+                    os.path.dirname(__file__), "softmax.java")
+                top_cg.add_code_lines(utils.get_file_content(filename))
+
         return top_cg.finalize_and_get_generated_code()
 
     # Required by SubroutinesMixin to create new code generator for
     # each subroutine.
     def create_code_generator(self):
         return JavaCodeGenerator(indent=self.indent)
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/javascript/interpreter.py b/m2cgen/interpreters/javascript/interpreter.py
index 4ac4afcd..359c68c1 100644
--- a/m2cgen/interpreters/javascript/interpreter.py
+++ b/m2cgen/interpreters/javascript/interpreter.py
@@ -25,9 +25,12 @@ class JavascriptInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "Math.log"
     log1p_function_name = "Math.log1p"
     power_function_name = "Math.pow"
+    softmax_function_name = "softmax"
     sqrt_function_name = "Math.sqrt"
     tanh_function_name = "Math.tanh"
 
+    with_softmax_expr = False
+
     def __init__(self, indent=4, function_name="score",
                  *args, **kwargs):
         self.indent = indent
@@ -53,4 +56,13 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "linear_algebra.js")
             self._cg.add_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.js")
+            self._cg.add_code_lines(utils.get_file_content(filename))
+
         return self._cg.finalize_and_get_generated_code()
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/php/interpreter.py b/m2cgen/interpreters/php/interpreter.py
index d25fa5a6..3957827c 100644
--- a/m2cgen/interpreters/php/interpreter.py
+++ b/m2cgen/interpreters/php/interpreter.py
@@ -23,9 +23,12 @@ class PhpInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "log"
     log1p_function_name = "log1p"
     power_function_name = "pow"
+    softmax_function_name = "softmax"
     sqrt_function_name = "sqrt"
     tanh_function_name = "tanh"
 
+    with_softmax_expr = False
+
     def __init__(self, indent=4, function_name="score", *args, **kwargs):
         self.function_name = function_name
 
@@ -47,6 +50,15 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "linear_algebra.php")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.php")
+            self._cg.prepend_code_lines(utils.get_file_content(filename))
+
         self._cg.prepend_code_line("<?php")
 
         return self._cg.finalize_and_get_generated_code()
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/powershell/interpreter.py b/m2cgen/interpreters/powershell/interpreter.py
index 4fb63dbc..bc979c82 100644
--- a/m2cgen/interpreters/powershell/interpreter.py
+++ b/m2cgen/interpreters/powershell/interpreter.py
@@ -24,10 +24,12 @@ class PowershellInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "[math]::Log"
     log1p_function_name = "Log1p"
     power_function_name = "[math]::Pow"
+    softmax_function_name = "Softmax"
     sqrt_function_name = "[math]::Sqrt"
     tanh_function_name = "[math]::Tanh"
 
     with_log1p_expr = False
+    with_softmax_expr = False
 
     def __init__(self, indent=4, function_name="Score", *args, **kwargs):
         self.function_name = function_name
@@ -57,6 +59,11 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "log1p.ps1")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.ps1")
+            self._cg.prepend_code_lines(utils.get_file_content(filename))
+
         return self._cg.finalize_and_get_generated_code()
 
     def interpret_abs_expr(self, expr, **kwargs):
@@ -98,3 +105,7 @@ def interpret_pow_expr(self, expr, **kwargs):
         exp_result = self._do_interpret(expr.exp_expr, **kwargs)
         return self._cg.math_function_invocation(
             self.power_function_name, base_result, exp_result)
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/python/interpreter.py b/m2cgen/interpreters/python/interpreter.py
index 72c283e7..d5d5cc4c 100644
--- a/m2cgen/interpreters/python/interpreter.py
+++ b/m2cgen/interpreters/python/interpreter.py
@@ -27,9 +27,12 @@ class PythonInterpreter(ImperativeToCodeInterpreter,
     logarithm_function_name = "math.log"
     log1p_function_name = "math.log1p"
     power_function_name = "math.pow"
+    softmax_function_name = "softmax"
     sqrt_function_name = "math.sqrt"
     tanh_function_name = "math.tanh"
 
+    with_softmax_expr = False
+
     def __init__(self, indent=4, function_name="score", *args, **kwargs):
         self.function_name = function_name
 
@@ -51,6 +54,11 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "linear_algebra.py")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.py")
+            self._cg.prepend_code_lines(utils.get_file_content(filename))
+
         if self.with_math_module:
             self._cg.add_dependency("math")
 
@@ -60,3 +68,7 @@ def interpret_abs_expr(self, expr, **kwargs):
         nested_result = self._do_interpret(expr.expr, **kwargs)
         return self._cg.function_invocation(
             self.abs_function_name, nested_result)
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/r/interpreter.py b/m2cgen/interpreters/r/interpreter.py
index 5c70626a..86b28a49 100644
--- a/m2cgen/interpreters/r/interpreter.py
+++ b/m2cgen/interpreters/r/interpreter.py
@@ -1,4 +1,6 @@
-from m2cgen.interpreters import mixins
+import os
+
+from m2cgen.interpreters import utils, mixins
 from m2cgen.interpreters.interpreter import ImperativeToCodeInterpreter
 from m2cgen.interpreters.r.code_generator import RCodeGenerator
 
@@ -27,9 +29,12 @@ class RInterpreter(ImperativeToCodeInterpreter,
     exponent_function_name = "exp"
     logarithm_function_name = "log"
     log1p_function_name = "log1p"
+    softmax_function_name = "softmax"
     sqrt_function_name = "sqrt"
     tanh_function_name = "tanh"
 
+    with_softmax_expr = False
+
     def __init__(self, indent=4, function_name="score", *args, **kwargs):
         self.indent = indent
         self.function_name = function_name
@@ -42,6 +47,11 @@ def interpret(self, expr):
         self.enqueue_subroutine(self.function_name, expr)
         self.process_subroutine_queue(top_cg)
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.r")
+            top_cg.prepend_code_lines(utils.get_file_content(filename))
+
         return top_cg.finalize_and_get_generated_code()
 
     def create_code_generator(self):
@@ -53,6 +63,10 @@ def interpret_pow_expr(self, expr, **kwargs):
         return self._cg.infix_expression(
             left=base_result, right=exp_result, op="^")
 
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
+
     def interpret_bin_vector_expr(self, expr, **kwargs):
         self.with_linear_algebra = True
         return self._cg.infix_expression(
diff --git a/m2cgen/interpreters/ruby/interpreter.py b/m2cgen/interpreters/ruby/interpreter.py
index e5e7e786..818c62cd 100644
--- a/m2cgen/interpreters/ruby/interpreter.py
+++ b/m2cgen/interpreters/ruby/interpreter.py
@@ -22,10 +22,12 @@ class RubyInterpreter(ImperativeToCodeInterpreter,
     exponent_function_name = "Math.exp"
     logarithm_function_name = "Math.log"
     log1p_function_name = "log1p"
+    softmax_function_name = "softmax"
     sqrt_function_name = "Math.sqrt"
     tanh_function_name = "Math.tanh"
 
     with_log1p_expr = False
+    with_softmax_expr = False
 
     def __init__(self, indent=4, function_name="score", *args, **kwargs):
         self.function_name = function_name
@@ -46,13 +48,18 @@ def interpret(self, expr):
         if self.with_linear_algebra:
             filename = os.path.join(
                 os.path.dirname(__file__), "linear_algebra.rb")
-            self._cg.prepend_code_lines(utils.get_file_content(filename))
+            self._cg.add_code_lines(utils.get_file_content(filename))
 
         if self.with_log1p_expr:
             filename = os.path.join(
                 os.path.dirname(__file__), "log1p.rb")
             self._cg.add_code_lines(utils.get_file_content(filename))
 
+        if self.with_softmax_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "softmax.rb")
+            self._cg.add_code_lines(utils.get_file_content(filename))
+
         return self._cg.finalize_and_get_generated_code()
 
     def interpret_bin_num_expr(self, expr, **kwargs):
@@ -80,3 +87,7 @@ def interpret_pow_expr(self, expr, **kwargs):
     def interpret_log1p_expr(self, expr, **kwargs):
         self.with_log1p_expr = True
         return super().interpret_log1p_expr(expr, **kwargs)
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/m2cgen/interpreters/visual_basic/interpreter.py b/m2cgen/interpreters/visual_basic/interpreter.py
index 6e32552f..66fbea6c 100644
--- a/m2cgen/interpreters/visual_basic/interpreter.py
+++ b/m2cgen/interpreters/visual_basic/interpreter.py
@@ -22,10 +22,12 @@ class VisualBasicInterpreter(ImperativeToCodeInterpreter,
     exponent_function_name = "Math.Exp"
     logarithm_function_name = "Math.Log"
     log1p_function_name = "Log1p"
+    softmax_function_name = "Softmax"
     tanh_function_name = "Tanh"
 
     with_atan_expr = False
     with_log1p_expr = False
+    with_softmax_expr = False
     with_tanh_expr = False
 
     def __init__(self, module_name="Model", indent=4, function_name="Score",
@@ -55,10 +57,9 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "linear_algebra.bas")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
-        # Use own Tanh function in order to be compatible with both VB and VBA
-        if self.with_tanh_expr:
+        if self.with_atan_expr:
             filename = os.path.join(
-                os.path.dirname(__file__), "tanh.bas")
+                os.path.dirname(__file__), "atan.bas")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
         if self.with_log1p_expr:
@@ -66,9 +67,15 @@ def interpret(self, expr):
                 os.path.dirname(__file__), "log1p.bas")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
-        if self.with_atan_expr:
+        if self.with_softmax_expr:
             filename = os.path.join(
-                os.path.dirname(__file__), "atan.bas")
+                os.path.dirname(__file__), "softmax.bas")
+            self._cg.prepend_code_lines(utils.get_file_content(filename))
+
+        # Use own Tanh function in order to be compatible with both VB and VBA
+        if self.with_tanh_expr:
+            filename = os.path.join(
+                os.path.dirname(__file__), "tanh.bas")
             self._cg.prepend_code_lines(utils.get_file_content(filename))
 
         self._cg.prepend_code_line(self._cg.tpl_module_definition(
@@ -95,3 +102,7 @@ def interpret_tanh_expr(self, expr, **kwargs):
     def interpret_atan_expr(self, expr, **kwargs):
         self.with_atan_expr = True
         return super().interpret_atan_expr(expr, **kwargs)
+
+    def interpret_softmax_expr(self, expr, **kwargs):
+        self.with_softmax_expr = True
+        return super().interpret_softmax_expr(expr, **kwargs)
diff --git a/setup.cfg b/setup.cfg
index aa78dd4c..f5dc4654 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -4,4 +4,4 @@ license-file = LICENSE
 version = m2cgen/VERSION.txt
 
 [flake8]
-exclude = .git,generated_code_examples/,m2cgen/interpreters/python/linear_algebra.py
+exclude = .git,generated_code_examples/,m2cgen/interpreters/python/linear_algebra.py,m2cgen/interpreters/python/softmax.py
diff --git a/tests/assemblers/test_lightgbm.py b/tests/assemblers/test_lightgbm.py
index a910ed42..4501aa6a 100644
--- a/tests/assemblers/test_lightgbm.py
+++ b/tests/assemblers/test_lightgbm.py
@@ -56,19 +56,8 @@ def test_multi_class():
     assembler = assemblers.LightGBMModelAssembler(estimator)
     actual = assembler.assemble()
 
-    exponent = ast.ExpExpr(
-        ast.NumVal(-1.0986122886681098),
-        to_reuse=True)
-
-    exponent_sum = ast.BinNumExpr(
-        ast.BinNumExpr(exponent, exponent, ast.BinNumOpType.ADD),
-        exponent,
-        ast.BinNumOpType.ADD,
-        to_reuse=True)
-
-    softmax = ast.BinNumExpr(exponent, exponent_sum, ast.BinNumOpType.DIV)
-
-    expected = ast.VectorVal([softmax] * 3)
+    num_expr = ast.NumVal(-1.0986122886681098)
+    expected = ast.SoftmaxExpr([num_expr] * 3)
 
     assert utils.cmp_exprs(actual, expected)
 
diff --git a/tests/assemblers/test_xgboost.py b/tests/assemblers/test_xgboost.py
index 6727da14..ee949e1a 100644
--- a/tests/assemblers/test_xgboost.py
+++ b/tests/assemblers/test_xgboost.py
@@ -57,22 +57,11 @@ def test_multi_class():
     assembler = assemblers.XGBoostModelAssemblerSelector(estimator)
     actual = assembler.assemble()
 
-    exponent = ast.ExpExpr(
-        ast.BinNumExpr(
-            ast.NumVal(0.5),
-            ast.NumVal(0.0),
-            ast.BinNumOpType.ADD),
-        to_reuse=True)
-
-    exponent_sum = ast.BinNumExpr(
-        ast.BinNumExpr(exponent, exponent, ast.BinNumOpType.ADD),
-        exponent,
-        ast.BinNumOpType.ADD,
-        to_reuse=True)
-
-    softmax = ast.BinNumExpr(exponent, exponent_sum, ast.BinNumOpType.DIV)
-
-    expected = ast.VectorVal([softmax] * 3)
+    num_expr = ast.BinNumExpr(
+        ast.NumVal(0.5),
+        ast.NumVal(-7.663454759665456e-09),
+        ast.BinNumOpType.ADD)
+    expected = ast.SoftmaxExpr([num_expr] * 3)
 
     assert utils.cmp_exprs(actual, expected)
 
@@ -156,68 +145,43 @@ def test_multi_class_best_ntree_limit():
     assembler = assemblers.XGBoostModelAssemblerSelector(estimator)
     actual = assembler.assemble()
 
-    estimator_exp_class1 = ast.ExpExpr(
-        ast.BinNumExpr(
-            ast.NumVal(0.5),
-            ast.IfExpr(
-                ast.CompExpr(
-                    ast.FeatureRef(2),
-                    ast.NumVal(2.45000005),
-                    ast.CompOpType.GTE),
-                ast.NumVal(-0.219950154),
-                ast.NumVal(0.430243909)),
-            ast.BinNumOpType.ADD),
-        to_reuse=True)
-
-    estimator_exp_class2 = ast.ExpExpr(
-        ast.BinNumExpr(
-            ast.NumVal(0.5),
-            ast.IfExpr(
-                ast.CompExpr(
-                    ast.FeatureRef(2),
-                    ast.NumVal(2.45000005),
-                    ast.CompOpType.GTE),
-                ast.NumVal(0.103241883),
-                ast.NumVal(-0.215121984)),
-            ast.BinNumOpType.ADD),
-        to_reuse=True)
+    estimator_class1 = ast.BinNumExpr(
+        ast.NumVal(0.5),
+        ast.IfExpr(
+            ast.CompExpr(
+                ast.FeatureRef(2),
+                ast.NumVal(2.450000047683716),
+                ast.CompOpType.GTE),
+            ast.NumVal(-0.21995015442371368),
+            ast.NumVal(0.43024390935897827)),
+        ast.BinNumOpType.ADD)
 
-    estimator_exp_class3 = ast.ExpExpr(
-        ast.BinNumExpr(
-            ast.NumVal(0.5),
-            ast.IfExpr(
-                ast.CompExpr(
-                    ast.FeatureRef(3),
-                    ast.NumVal(1.6500001),
-                    ast.CompOpType.GTE),
-                ast.NumVal(0.402985066),
-                ast.NumVal(-0.193333372)),
-            ast.BinNumOpType.ADD),
-        to_reuse=True)
+    estimator_class2 = ast.BinNumExpr(
+        ast.NumVal(0.5),
+        ast.IfExpr(
+            ast.CompExpr(
+                ast.FeatureRef(2),
+                ast.NumVal(2.450000047683716),
+                ast.CompOpType.GTE),
+            ast.NumVal(0.10324188321828842),
+            ast.NumVal(-0.21512198448181152)),
+        ast.BinNumOpType.ADD)
 
-    exp_sum = ast.BinNumExpr(
-        ast.BinNumExpr(
-            estimator_exp_class1,
-            estimator_exp_class2,
-            ast.BinNumOpType.ADD),
-        estimator_exp_class3,
-        ast.BinNumOpType.ADD,
-        to_reuse=True)
+    estimator_class3 = ast.BinNumExpr(
+        ast.NumVal(0.5),
+        ast.IfExpr(
+            ast.CompExpr(
+                ast.FeatureRef(3),
+                ast.NumVal(1.6500000953674316),
+                ast.CompOpType.GTE),
+            ast.NumVal(0.4029850661754608),
+            ast.NumVal(-0.19333337247371674)),
+        ast.BinNumOpType.ADD)
 
-    expected = ast.VectorVal([
-        ast.BinNumExpr(
-            estimator_exp_class1,
-            exp_sum,
-            ast.BinNumOpType.DIV),
-        ast.BinNumExpr(
-            estimator_exp_class2,
-            exp_sum,
-            ast.BinNumOpType.DIV),
-        ast.BinNumExpr(
-            estimator_exp_class3,
-            exp_sum,
-            ast.BinNumOpType.DIV)
-    ])
+    expected = ast.SoftmaxExpr([
+        estimator_class1,
+        estimator_class2,
+        estimator_class3])
 
     assert utils.cmp_exprs(actual, expected)
 
diff --git a/tests/interpreters/test_c.py b/tests/interpreters/test_c.py
index 39d283d6..443fa76a 100644
--- a/tests/interpreters/test_c.py
+++ b/tests/interpreters/test_c.py
@@ -324,6 +324,36 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    interpreter = interpreters.CInterpreter()
+
+    expected_code = """
+#include <string.h>
+void softmax(double *x, int size, double *result) {
+    double max = x[0];
+    for (int i = 1; i < size; ++i) {
+        if (x[i] > max)
+            max = x[i];
+    }
+    double sum = 0.0;
+    for (int i = 0; i < size; ++i) {
+        result[i] = exp(x[i] - max);
+        sum += result[i];
+    }
+    for (int i = 0; i < size; ++i)
+        result[i] /= sum;
+}
+void score(double * input, double * output) {
+    double var0[2];
+    softmax((double[]){2.0, 3.0}, 2, var0);
+    memcpy(output, var0, 2 * sizeof(double));
+}"""
+
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_c_sharp.py b/tests/interpreters/test_c_sharp.py
index 77873dac..b89672c6 100644
--- a/tests/interpreters/test_c_sharp.py
+++ b/tests/interpreters/test_c_sharp.py
@@ -489,6 +489,41 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = """
+using static System.Math;
+namespace ML {
+    public static class Model {
+        public static double[] Score(double[] input) {
+            return Softmax(new double[2] {2.0, 3.0});
+        }
+        private static double[] Softmax(double[] x) {
+            int size = x.Length;
+            double[] result = new double[size];
+            double max = x[0];
+            for (int i = 1; i < size; ++i) {
+                if (x[i] > max)
+                    max = x[i];
+            }
+            double sum = 0.0;
+            for (int i = 0; i < size; ++i) {
+                result[i] = Exp(x[i] - max);
+                sum += result[i];
+            }
+            for (int i = 0; i < size; ++i)
+                result[i] /= sum;
+            return result;
+        }
+    }
+}
+"""
+
+    interpreter = CSharpInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_dart.py b/tests/interpreters/test_dart.py
index f5643fbe..721f8179 100644
--- a/tests/interpreters/test_dart.py
+++ b/tests/interpreters/test_dart.py
@@ -566,6 +566,32 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = """
+import 'dart:math';
+List<double> score(List<double> input) {
+    return softmax([2.0, 3.0]);
+}
+List<double> softmax(List<double> x) {
+    int size = x.length;
+    List<double> result = new List<double>(size);
+    double maxElem = x.reduce(max);
+    double sum = 0.0;
+    for (int i = 0; i < size; ++i) {
+        result[i] = exp(x[i] - maxElem);
+        sum += result[i];
+    }
+    for (int i = 0; i < size; ++i)
+        result[i] /= sum;
+    return result;
+}
+"""
+    interpreter = DartInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_f_sharp.py b/tests/interpreters/test_f_sharp.py
index b8f2c595..b9fca010 100644
--- a/tests/interpreters/test_f_sharp.py
+++ b/tests/interpreters/test_f_sharp.py
@@ -455,6 +455,23 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = """
+let private softmax x =
+    let maxElem = List.reduce max x
+    let exps = List.map (fun i -> exp (i - maxElem)) x
+    let sumExps = List.sum exps
+    List.map (fun i -> i / sumExps) exps
+let score (input : double list) =
+    softmax ([2.0; 3.0])
+"""
+
+    interpreter = FSharpInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_go.py b/tests/interpreters/test_go.py
index c5ad342e..e2281326 100644
--- a/tests/interpreters/test_go.py
+++ b/tests/interpreters/test_go.py
@@ -327,6 +327,39 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    interpreter = interpreters.GoInterpreter()
+
+    expected_code = """
+import "math"
+func softmax(x []float64) []float64 {
+    size := len(x)
+    result := make([]float64, size)
+    max := x[0]
+    for _, v := range x {
+        if (v > max) {
+            max = v
+        }
+    }
+    sum := 0.0
+    for i := 0; i < size; i++ {
+        result[i] = math.Exp(x[i] - max)
+        sum += result[i]
+    }
+    for i := 0; i < size; i++ {
+        result[i] /= sum
+    }
+    return result
+}
+func score(input []float64) []float64 {
+    return softmax([]float64{2.0, 3.0})
+}"""
+
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_haskell.py b/tests/interpreters/test_haskell.py
index 490bb5d6..b8013cd5 100644
--- a/tests/interpreters/test_haskell.py
+++ b/tests/interpreters/test_haskell.py
@@ -361,6 +361,27 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = r"""
+module Model where
+softmax :: [Double] -> [Double]
+softmax x =
+    let
+        m = maximum x
+        exps = map (\i -> exp (i - m)) x
+        sumExps = sum exps
+    in map (\i -> i / sumExps) exps
+score :: [Double] -> [Double]
+score input =
+    softmax ([2.0, 3.0])
+"""
+
+    interpreter = HaskellInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_java.py b/tests/interpreters/test_java.py
index 2ebf0afe..67c5bf08 100644
--- a/tests/interpreters/test_java.py
+++ b/tests/interpreters/test_java.py
@@ -557,6 +557,38 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    interpreter = interpreters.JavaInterpreter()
+
+    expected_code = """
+public class Model {
+    public static double[] score(double[] input) {
+        return softmax(new double[] {2.0, 3.0});
+    }
+    private static double[] softmax(double[] x) {
+        int size = x.length;
+        double[] result = new double[size];
+        double max = x[0];
+        for (int i = 1; i < size; ++i) {
+            if (x[i] > max)
+                max = x[i];
+        }
+        double sum = 0.0;
+        for (int i = 0; i < size; ++i) {
+            result[i] = Math.exp(x[i] - max);
+            sum += result[i];
+        }
+        for (int i = 0; i < size; ++i)
+            result[i] /= sum;
+        return result;
+    }
+}"""
+
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_javascript.py b/tests/interpreters/test_javascript.py
index 5dfc3861..9fe9afb0 100644
--- a/tests/interpreters/test_javascript.py
+++ b/tests/interpreters/test_javascript.py
@@ -342,6 +342,37 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    interpreter = interpreters.JavascriptInterpreter()
+
+    expected_code = """
+function score(input) {
+    return softmax([2.0, 3.0]);
+}
+function softmax(x) {
+    let size = x.length;
+    let result = new Array(size);
+    let max = x[0];
+    for (let i = 1; i < size; ++i) {
+        if (x[i] > max)
+            max = x[i];
+    }
+    let sum = 0.0;
+    for (let i = 0; i < size; ++i) {
+        result[i] = Math.exp(x[i] - max);
+        sum += result[i];
+    }
+    for (let i = 0; i < size; ++i)
+        result[i] /= sum;
+    return result;
+}
+"""
+
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_php.py b/tests/interpreters/test_php.py
index d9b00ed1..cc9c5f59 100644
--- a/tests/interpreters/test_php.py
+++ b/tests/interpreters/test_php.py
@@ -347,6 +347,34 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = """
+<?php
+function softmax(array $x) {
+    $size = count($x);
+    $result = array();
+    $m = max($x);
+    $sum = 0.0;
+    for ($i = 0; $i < $size; ++$i) {
+        $result[$i] = exp($x[$i] - $m);
+        $sum += $result[$i];
+    }
+    for ($i = 0; $i < $size; ++$i) {
+        $result[$i] /= $sum;
+    }
+    return $result;
+}
+function score(array $input) {
+    return softmax(array(2.0, 3.0));
+}
+"""
+
+    interpreter = PhpInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_powershell.py b/tests/interpreters/test_powershell.py
index 816ae8bb..2c9114b2 100644
--- a/tests/interpreters/test_powershell.py
+++ b/tests/interpreters/test_powershell.py
@@ -383,6 +383,38 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = """
+function Softmax([double[]] $x) {
+    [int] $size = $x.Length
+    [double[]] $result = @(0) * $size
+    [double] $max = $x[0]
+    for ([int] $i = 1; $i -lt $size; ++$i) {
+        if ($x[$i] -gt $max) {
+            $max = $x[$i]
+        }
+    }
+    [double] $sum = 0.0
+    for ([int] $i = 0; $i -lt $size; ++$i) {
+        $result[$i] = [math]::Exp($x[$i] - $max)
+        $sum += $result[$i]
+    }
+    for ([int] $i = 0; $i -lt $size; ++$i) {
+        $result[$i] /= $sum;
+    }
+    return $result
+}
+function Score([double[]] $InputVector) {
+    return Softmax $(@($(2.0), $(3.0)))
+}
+"""
+
+    interpreter = PowershellInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_python.py b/tests/interpreters/test_python.py
index 8bf1aa71..0e638a93 100644
--- a/tests/interpreters/test_python.py
+++ b/tests/interpreters/test_python.py
@@ -446,6 +446,27 @@ def score(input):
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    interpreter = interpreters.PythonInterpreter()
+
+    expected_code = """
+import math
+def softmax(x):
+    m = max(x)
+    exps = [math.exp(i - m) for i in x]
+    s = sum(exps)
+    for idx, _ in enumerate(exps):
+        exps[idx] = exps[idx] / s
+    return exps
+def score(input):
+    return softmax([2.0, 3.0])
+    """
+
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_r.py b/tests/interpreters/test_r.py
index 3684b810..7f622226 100644
--- a/tests/interpreters/test_r.py
+++ b/tests/interpreters/test_r.py
@@ -463,6 +463,25 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = """
+softmax <- function (x) {
+    m <- max(x)
+    exps <- exp(x - m)
+    s <- sum(exps)
+    return(exps / s)
+}
+score <- function(input) {
+    return(softmax(c(2.0, 3.0)))
+}
+"""
+
+    interpreter = RInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_ruby.py b/tests/interpreters/test_ruby.py
index 8dd7d427..ff9dea60 100644
--- a/tests/interpreters/test_ruby.py
+++ b/tests/interpreters/test_ruby.py
@@ -167,15 +167,15 @@ def test_bin_vector_expr():
         ast.BinNumOpType.ADD)
 
     expected_code = """
+def score(input)
+    add_vectors([1.0, 2.0], [3.0, 4.0])
+end
 def add_vectors(v1, v2)
     v1.zip(v2).map { |x, y| x + y }
 end
 def mul_vector_number(v1, num)
     v1.map { |i| i * num }
 end
-def score(input)
-    add_vectors([1.0, 2.0], [3.0, 4.0])
-end
 """
 
     interpreter = RubyInterpreter()
@@ -189,15 +189,15 @@ def test_bin_vector_num_expr():
         ast.BinNumOpType.MUL)
 
     expected_code = """
+def score(input)
+    mul_vector_number([1.0, 2.0], 1.0)
+end
 def add_vectors(v1, v2)
     v1.zip(v2).map { |x, y| x + y }
 end
 def mul_vector_number(v1, num)
     v1.map { |i| i * num }
 end
-def score(input)
-    mul_vector_number([1.0, 2.0], 1.0)
-end
 """
 
     interpreter = RubyInterpreter()
@@ -363,6 +363,29 @@ def score(input)
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = """
+def score(input)
+    softmax([2.0, 3.0])
+end
+def softmax(x)
+    m = x.max
+    exps = []
+    s = 0.0
+    x.each_with_index do |v, i|
+        exps[i] = Math.exp(v - m)
+        s += exps[i]
+    end
+    exps.map { |i| i / s }
+end
+"""
+
+    interpreter = RubyInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/interpreters/test_visual_basic.py b/tests/interpreters/test_visual_basic.py
index e25d9705..0647f10b 100644
--- a/tests/interpreters/test_visual_basic.py
+++ b/tests/interpreters/test_visual_basic.py
@@ -571,6 +571,48 @@ def test_atan_expr():
     utils.assert_code_equal(interpreter.interpret(expr), expected_code)
 
 
+def test_softmax_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    expected_code = """
+Module Model
+Function Softmax(ByRef x() As Double) As Double()
+    Dim size As Integer
+    size = UBound(x) - LBound(x)
+    Dim result() As Double
+    ReDim result(size)
+    Dim max As Double
+    max = x(LBound(x))
+    Dim i As Integer
+    For i = LBound(x) + 1 To UBound(x)
+        If x(i) > max Then
+            max = x(i)
+        End If
+    Next i
+    Dim sum As Double
+    sum = 0.0
+    For i = LBound(x) To UBound(x)
+        result(i) = Math.Exp(x(i) - max)
+        sum = sum + result(i)
+    Next i
+    For i = LBound(x) To UBound(x)
+        result(i) = result(i) / sum
+    Next i
+    Softmax = result
+End Function
+Function Score(ByRef inputVector() As Double) As Double()
+    Dim var0(1) As Double
+    var0(0) = 2.0
+    var0(1) = 3.0
+    Score = Softmax(var0)
+End Function
+End Module
+"""
+
+    interpreter = VisualBasicInterpreter()
+    utils.assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
 def test_reused_expr():
     reused_expr = ast.ExpExpr(ast.NumVal(1.0), to_reuse=True)
     expr = ast.BinNumExpr(reused_expr, reused_expr, ast.BinNumOpType.DIV)
diff --git a/tests/test_ast.py b/tests/test_ast.py
index cf37911a..af43d85a 100644
--- a/tests/test_ast.py
+++ b/tests/test_ast.py
@@ -60,7 +60,7 @@ def test_count_exprs_exclude_list():
                 ast.BinNumOpType.ADD)
         ]),
         ast.IdExpr(
-            ast.VectorVal([
+            ast.SoftmaxExpr([
                 ast.NumVal(1),
                 ast.NumVal(2),
                 ast.NumVal(3),
@@ -106,9 +106,9 @@ def test_exprs_str():
 PowExpr(NumVal(2.0),NumVal(3.0),to_reuse=False),
 TanhExpr(NumVal(1.0),to_reuse=False),
 BinNumExpr(NumVal(0.0),FeatureRef(0),ADD,to_reuse=False)]),
-IdExpr(VectorVal([
+IdExpr(SoftmaxExpr(
 NumVal(1.0),NumVal(2.0),NumVal(3.0),NumVal(4.0),NumVal(5.0),
-NumVal(6.0),NumVal(7.0),NumVal(8.0),FeatureRef(1)]),to_reuse=False),SUB),
+NumVal(6.0),NumVal(7.0),NumVal(8.0),FeatureRef(1),to_reuse=False),to_reuse=False),SUB),
 IfExpr(CompExpr(NumVal(2.0),NumVal(0.0),GT),NumVal(3.0),NumVal(4.0)),MUL)
 """.strip().replace("\n", "")
 
diff --git a/tests/test_fallback_expressions.py b/tests/test_fallback_expressions.py
index 8b58c85b..261ecf88 100644
--- a/tests/test_fallback_expressions.py
+++ b/tests/test_fallback_expressions.py
@@ -147,3 +147,27 @@ def score(input):
         """(var5)) + (var6)) * (var7)""")
 
     assert_code_equal(interpreter.interpret(expr), expected_code)
+
+
+def test_softmax_fallback_expr():
+    expr = ast.SoftmaxExpr([ast.NumVal(2.0), ast.NumVal(3.0)])
+
+    class InterpreterWithoutSoftmax(PythonInterpreter):
+        softmax_function_name = NotImplemented
+
+        def interpret_softmax_expr(self, expr, **kwargs):
+            return super(PythonInterpreter, self).interpret_softmax_expr(
+                expr, **kwargs)
+
+    interpreter = InterpreterWithoutSoftmax()
+
+    expected_code = """
+import math
+def score(input):
+    var0 = math.exp(2.0)
+    var1 = math.exp(3.0)
+    var2 = (var0) + (var1)
+    return [(var0) / (var2), (var1) / (var2)]
+"""
+
+    assert_code_equal(interpreter.interpret(expr), expected_code)