diff --git a/hls4ml/backends/fpga/fpga_backend.py b/hls4ml/backends/fpga/fpga_backend.py
index db0f256172..eebb8d939f 100644
--- a/hls4ml/backends/fpga/fpga_backend.py
+++ b/hls4ml/backends/fpga/fpga_backend.py
@@ -129,22 +129,29 @@ def __init__(self, name):
             ConfigurableAttribute('skip', value_type=bool, default=False, description=descriptions.softmax_skip),
             TypeAttribute(
                 'exp_table',
-                default=FixedPrecisionType(18, 8, rounding_mode=RoundingMode.RND, saturation_mode=SaturationMode.SAT),
+                default=FixedPrecisionType(
+                    18, 8, signed=False, rounding_mode=RoundingMode.RND, saturation_mode=SaturationMode.SAT
+                ),
                 description=descriptions.table_type,
             ),
             TypeAttribute(
                 'inv_table',
-                default=FixedPrecisionType(18, 8, rounding_mode=RoundingMode.RND, saturation_mode=SaturationMode.SAT),
+                default=FixedPrecisionType(
+                    18, 8, signed=False, rounding_mode=RoundingMode.RND, saturation_mode=SaturationMode.SAT
+                ),
                 description=descriptions.table_type,
             ),
             TypeAttribute(
                 'inv_inp',
-                default=FixedPrecisionType(18, 8, rounding_mode=RoundingMode.RND, saturation_mode=SaturationMode.SAT),
+                default=UnspecifiedPrecisionType(),
+                description='What the accumulated value is cast to before accessing the inversion table (only in stable)',
             ),
             TypeAttribute(
-                'accum',
-                default=FixedPrecisionType(18, 8, rounding_mode=RoundingMode.RND, saturation_mode=SaturationMode.SAT),
+                'inp_norm',
+                default=UnspecifiedPrecisionType(),
+                description='The internal width used for the exp table lookup (only in stable)',
             ),
+            TypeAttribute('accum', description=descriptions.accum_type),
         ]
         self.attribute_map[Softmax] = softmax_attrs
 
diff --git a/hls4ml/backends/oneapi/oneapi_backend.py b/hls4ml/backends/oneapi/oneapi_backend.py
index 0c11c16d09..d448cd41cc 100644
--- a/hls4ml/backends/oneapi/oneapi_backend.py
+++ b/hls4ml/backends/oneapi/oneapi_backend.py
@@ -19,7 +19,6 @@
     Embedding,
     Layer,
     SimpleRNN,
-    Softmax,
 )
 from hls4ml.model.optimizer import get_backend_passes, layer_optimizer
 from hls4ml.model.types import FixedPrecisionType, IntegerPrecisionType, NamedType
@@ -210,7 +209,11 @@ def build(self, model, build_type='fpga_emu', run=False):
         try:
             subprocess.run('which icpx', shell=True, cwd=builddir, check=True)
         except subprocess.CalledProcessError:
-            raise RuntimeError('Could not find icpx. Please configure oneAPI appropriately')
+            print('icpx not found. Trying ahls')
+            try:
+                subprocess.run('which ahls', shell=True, cwd=builddir, check=True)
+            except subprocess.CalledProcessError:
+                raise RuntimeError('Could not find icpx or ahls. Please configure oneAPI appropriately')
         subprocess.run('cmake ..', shell=True, cwd=builddir, check=True)
         subprocess.run(f'make {build_type}', shell=True, cwd=builddir, check=True)
 
@@ -257,13 +260,6 @@ def init_activation(self, layer):
         if layer.get_attr('recurrent_activation') == 'tanh':
             layer.set_attr('recurrent_activation', 'dense_tanh')
 
-    @layer_optimizer(Softmax)
-    def init_softmax(self, layer):
-        if layer.model.config.get_config_value('IOType') == 'io_parallel':
-            assert len(layer.get_input_variable().shape) == 1, (
-                'Softmax with io_parallel strategy cannot be used on multidimensional tensors.'
-            )
-
     @layer_optimizer(Embedding)
     def init_embed(self, layer):
         if layer.attributes['n_in'] is None:
diff --git a/hls4ml/backends/oneapi/passes/core_templates.py b/hls4ml/backends/oneapi/passes/core_templates.py
index 9602b2d0fc..bab202a1f5 100644
--- a/hls4ml/backends/oneapi/passes/core_templates.py
+++ b/hls4ml/backends/oneapi/passes/core_templates.py
@@ -194,12 +194,42 @@ def format(self, node):
 
 softmax_config_template = """struct {type}_config{index} : nnet::activ_config {{
     static constexpr unsigned n_in = {n_in};
+
+    // For multi-dim softmax
+    static const unsigned n_slice = {n_slice};
+    static const unsigned n_outer = {n_outer};
+    static const unsigned n_inner = {n_inner};
+
+    // For legacy softmax
+    typedef {table_t.name} table_t;
     static constexpr unsigned table_size = {table_size};
+
+    static constexpr unsigned exp_table_size = {exp_table_size};
+    static constexpr unsigned inv_table_size = {inv_table_size};
     static constexpr unsigned io_type = nnet::{iotype};
     static constexpr unsigned reuse_factor = {reuse};
+
     static constexpr nnet::softmax_implementation implementation = nnet::softmax_implementation::{implementation};
+    typedef {smax_accum_t} accum_t;
     typedef {exp_table_t.name} exp_table_t;
-    typedef {inv_table_t.name} inv_table_t;
+    typedef {inv_table_t.name} inv_table_t;"""
+
+softmax_config_table_template = """
+
+    using {exp_table_name}_arr_t = nnet::array<exp_table_t, exp_table_size>;
+    using {inv_table_name}_arr_t = nnet::array<inv_table_t, inv_table_size>;
+    static constexpr const {exp_table_name}_arr_t exp_table = {exp_table_name};
+    static constexpr const {inv_table_name}_arr_t invert_table = {inv_table_name};
+}};\n"""
+
+softmax_config_table_template_stable = """
+    typedef {inv_inp_t.name} inv_inp_t;
+    typedef {inp_norm_t.name} inp_norm_t;
+
+    using {exp_table_name}_arr_t = nnet::array<exp_table_t, exp_table_size>;
+    using {inv_table_name}_arr_t = nnet::array<inv_table_t, inv_table_size>;
+    static constexpr const {exp_table_name}_arr_t exp_table = {exp_table_name};
+    static constexpr const {inv_table_name}_arr_t invert_table = {inv_table_name};
 }};\n"""
 
 activ_function_template = 'nnet::{activation}<{input_t}, {output_t}, {config}>({input}, {output});'
@@ -219,7 +249,34 @@ def __init__(self):
 
     def format(self, node):
         params = self._default_config_params(node)
-        params['type'] = node.get_attr('activation')
+        params['type'] = node.get_attr('activation').lower()
+
+        if (params['type'] == 'softmax') or (params['type'] == 'softmax_multidim'):
+            # If no table size is specified, assume default size of 1024
+            params.setdefault('exp_table_size', params['table_size'])
+            params.setdefault('inv_table_size', params['table_size'])
+            params['exp_table_size'] = min(params['exp_table_size'], params['table_size'])
+            params['inv_table_size'] = min(params['inv_table_size'], params['table_size'])
+
+            # This is for non-quantised layers where table size is not a layer attribute
+            if node.get_attr('exp_table_size', -1) == -1:
+                node.set_attr('exp_table_size', params['exp_table_size'])
+
+            if node.get_attr('inv_table_size', -1) == -1:
+                node.set_attr('inv_table_size', params['inv_table_size'])
+
+            params.setdefault('exp_scale', 1.0)
+            params.setdefault('parallelization_factor', -1)
+
+            n_slice = params['n_in'] // params['n_inner'] // params['n_outer']
+            params['n_slice'] = n_slice
+
+            params['exp_table_name'] = node.name + '_exp_table'
+            params['inv_table_name'] = node.name + '_inv_table'
+            params['smax_accum_t'] = params['accum_t'].name
+
+            if params['implementation'] == 'stable':
+                self.template = softmax_config_template + softmax_config_table_template_stable
 
         return self.template.format(**params)
 
@@ -251,7 +308,7 @@ def format(self, node):
 class SoftmaxConfigTemplate(ActivationConfigTemplate):
     def __init__(self):
         super(ActivationConfigTemplate, self).__init__(Softmax)  # Skip ActivationConfigTemplate's __init__
-        self.template = softmax_config_template
+        self.template = softmax_config_template + softmax_config_table_template
 
 
 class ActivationFunctionTemplate(FunctionCallTemplate):
@@ -304,6 +361,7 @@ def format(self, node):
         params = self._default_function_params(node)
         params['activation'] = node.get_attr('activation').lower()
         params['config'] = f'{node.get_attr("activation")}_config{node.index}'
+
         return self.template.format(**params)
 
 
diff --git a/hls4ml/converters/keras_v3/hgq2/softmax.py b/hls4ml/converters/keras_v3/hgq2/softmax.py
index 136662beba..324ded98be 100644
--- a/hls4ml/converters/keras_v3/hgq2/softmax.py
+++ b/hls4ml/converters/keras_v3/hgq2/softmax.py
@@ -12,10 +12,15 @@
     from keras import KerasTensor
 
 
-def fixed_quantizer_to_hls4ml_t(q: 'FixedPointQuantizerBase', take_max=False):
+def fixed_quantizer_to_hls4ml_t(q: 'FixedPointQuantizerBase', take_max=False, force_unsigned=False):
     from keras import ops
 
     k, i, f = q.kif
+
+    if force_unsigned:
+        k = 0.0
+        i += 1.0
+
     k = ops.convert_to_numpy(k)
     i = ops.convert_to_numpy(i)
     f = ops.convert_to_numpy(f)
@@ -78,10 +83,12 @@ def handle(
         exp_oq = layer.exp_table.oq.quantizer
         inv_oq = layer.inv_table.oq.quantizer
         inv_iq = layer.inv_table.iq.quantizer
+
         assert isinstance(exp_oq, FixedPointQuantizerBase), 'Only fixed-point quantizer is supported for exp_table'
-        exp_table_t = fixed_quantizer_to_hls4ml_t(exp_oq)
-        inv_table_t = fixed_quantizer_to_hls4ml_t(inv_oq)
-        inv_inp_t = fixed_quantizer_to_hls4ml_t(inv_iq)
+        # Enforce unsigned quantisers on all three types
+        exp_table_t = fixed_quantizer_to_hls4ml_t(exp_oq, force_unsigned=True)
+        inv_table_t = fixed_quantizer_to_hls4ml_t(inv_oq, force_unsigned=True)
+        inv_inp_t = fixed_quantizer_to_hls4ml_t(inv_iq, force_unsigned=True)
         exp_scale = layer.input_scaler
 
         inv_table_size = 2**inv_inp_t.width
@@ -99,13 +106,19 @@ def handle(
         else:
             raise ValueError(f'Too many inputs for softmax layer {layer.name}: expected 1 or 2, got {len(in_tensors)}')
 
+        # For masked implementation assume first input is the tensor we are operating on
+        activation = 'softmax'
+        if len(in_tensors[0].shape[1:]) > 1:
+            if (1 not in in_tensors[0].shape[1:]) or (len(in_tensors[0].shape[1:]) > 2):
+                activation = 'softmax_multidim'
+
         config = {}
         config.update(self.default_config)
         config.update(
             {
                 'axis': ax,
                 'n_in': n_in,
-                'activation': 'softmax',
+                'activation': activation,
                 'n_outer': n_outer,
                 'n_inner': n_inner,
                 'implementation': impl,
@@ -122,7 +135,8 @@ def handle(
         )
 
         if layer.stable:
-            inp_norm_t = fixed_quantizer_to_hls4ml_t(layer.exp_table.iq.quantizer)
+            # Force unsigned since norm >= 0
+            inp_norm_t = fixed_quantizer_to_hls4ml_t(layer.exp_table.iq.quantizer, force_unsigned=True)
             config['inp_norm_t'] = inp_norm_t
 
         return (config,)
diff --git a/hls4ml/model/layers.py b/hls4ml/model/layers.py
index 42afbabd3c..625f62366d 100644
--- a/hls4ml/model/layers.py
+++ b/hls4ml/model/layers.py
@@ -965,6 +965,9 @@ def initialize(self):
         if 'n_in' not in self.attributes:
             self.set_attr('n_in', self.get_input_variable().size())
 
+        # set the needed types if needed
+        self._set_type_t('table')
+
 
 class ParametrizedActivation(Activation):
     _expected_attributes = [
@@ -1031,6 +1034,12 @@ class Softmax(Activation):
     def initialize(self):
         super().initialize()
 
+        # set the needed types if needed
+        self._set_type_t('exp_table')
+        self._set_type_t('inv_table')
+        self._set_type_t('inv_inp')
+        self._set_type_t('inv_norm')
+
 
 class TernaryTanh(Activation):
     def initialize(self):
diff --git a/hls4ml/model/optimizer/passes/infer_precision.py b/hls4ml/model/optimizer/passes/infer_precision.py
index 0bc29a2955..208224944d 100644
--- a/hls4ml/model/optimizer/passes/infer_precision.py
+++ b/hls4ml/model/optimizer/passes/infer_precision.py
@@ -90,6 +90,9 @@ def _infer_precision(self, node, types_to_infer):
 
         if node_class in ['PReLU']:
             return self._infer_prelu_act_precision(node, types_to_infer)
+
+        if node_class in ['Softmax']:
+            return self._infer_softmax_precision(node, types_to_infer)
         # What about quantized activation layer? Setting it to 'auto' manually will break it here. We should prevent
         # this in config_from_* functions
 
@@ -605,6 +608,46 @@ def _infer_prelu_act_precision(self, node, types_to_infer):
 
         return inferred_types
 
+    def _infer_softmax_precision(self, node, types_to_infer):
+        inferred_types = []
+
+        # for softmax, the table parameters have a default setting, so they don't need to be inferred
+        # here. We never expect them to be of type auto.
+
+        # For result, we leave it to be set externally (model default if not set). We expect it to
+        # likely be the output value, in which case the output format would determine it's precision.
+        # Therefore, only the accum is configured here
+
+        if 'accum_t' in types_to_infer:
+            exp_w = node.types['exp_table_t'].precision.width
+            exp_i = node.types['exp_table_t'].precision.integer
+            exp_s = node.types['exp_table_t'].precision.signed
+            ceillog = math.ceil(np.log2(node.get_attr('n_in')))
+            node.types['accum_t'].precision = FixedPrecisionType(exp_w + ceillog, exp_i + ceillog, signed=exp_s)
+            inferred_types.append('accum_t')
+
+        if 'inv_inp_t' in types_to_infer:
+            # if not set, just choose the accumulator type
+            node.types['inv_inp_t'].precision = node.types['accum_t'].precision
+            if node.model.config.backend.name == 'oneAPI':
+                # TODO - Can introduce this on onther backends since it benefits from only using required number of bits
+                # for the integer part. + 1 is for guarding
+                node.types['inv_inp_t'].precision.integer = 1 + math.ceil(math.log2(max(node.get_input_variable().shape)))
+            inferred_types.append('inv_inp_t')
+
+        if 'inp_norm_t' in types_to_infer:
+            in_type = node.get_input_variable().type.precision
+            if node.model.config.backend.name == 'oneAPI':
+                inp_norm_width = in_type.width
+                inp_norm_int = in_type.integer
+            else:
+                inp_norm_width = in_type.width - in_type.signed
+                inp_norm_int = in_type.integer - in_type.signed
+            node.types['inp_norm_t'].precision = FixedPrecisionType(inp_norm_width, inp_norm_int, signed=False)
+            inferred_types.append('inp_norm_t')
+
+        return inferred_types
+
 
 def _get_precision_from_constant(value: int | float, max_width=8):
     """A utility function to find a fixed type to store the constant
diff --git a/hls4ml/model/types.py b/hls4ml/model/types.py
index aa69130d9b..85ed1c8241 100644
--- a/hls4ml/model/types.py
+++ b/hls4ml/model/types.py
@@ -434,6 +434,9 @@ class UnspecifiedPrecisionType(PrecisionType):
     def __init__(self):
         super().__init__(width=0, signed=False)
 
+    def __str__(self):
+        return 'auto'
+
 
 def find_minimum_width(data, signed=True):
     """
diff --git a/hls4ml/templates/oneapi/firmware/nnet_utils/nnet_activation.h b/hls4ml/templates/oneapi/firmware/nnet_utils/nnet_activation.h
index f118ecb05c..476e49c6b0 100644
--- a/hls4ml/templates/oneapi/firmware/nnet_utils/nnet_activation.h
+++ b/hls4ml/templates/oneapi/firmware/nnet_utils/nnet_activation.h
@@ -99,21 +99,22 @@ template <class data_T, class res_T, typename CONFIG_T> void sigmoid(const data_
 
 enum class softmax_implementation { latency = 0, legacy = 1, stable = 2, argmax = 3 };
 
-template <class data_T, typename CONFIG_T> inline unsigned softmax_stable_idx_from_real_val(const data_T x) {
+template <class data_T, unsigned table_size> inline unsigned softmax_idx_from_real_val_negative(const data_T x) {
     // Number of address bits for table
-    static constexpr int N = ceillog2<CONFIG_T::table_size>::val;
+    static constexpr int N = ceillog2<table_size>::val;
 
     // Slice the top N bits of the input
     [[intel::fpga_register]] ac_int<N, false> y = x.template slc<N>(x.width - N - 1);
+
     // If x is the most negative value, the slice will be 0, so we need to set the 0-th bit to ensure correctness
     if (x != 0 && y == 0)
         y[0] = 1;
     return y.to_uint();
 }
 
-template <class data_T, typename CONFIG_T> inline unsigned softmax_latency_idx_from_real_val(const data_T x) {
+template <class data_T, unsigned table_size> inline unsigned softmax_idx_from_real_val(const data_T x) {
     // Number of address bits for table
-    static constexpr int N = ceillog2<CONFIG_T::table_size>::val;
+    static constexpr int N = ceillog2<table_size>::val;
 
     // Slice the top N bits of the input
     [[intel::fpga_register]] ac_int<N, false> y = x.template slc<N>(x.width - N);
@@ -121,9 +122,6 @@ template <class data_T, typename CONFIG_T> inline unsigned softmax_latency_idx_f
 }
 
 template <class data_T, class res_T, typename CONFIG_T> void softmax_stable(const data_T &data, res_T &res) {
-// Look-up tables
-#include "activation_tables/exp_table.tb"
-#include "activation_tables/invert_table.tb"
 
     // Find maximum
     Op_max<typename data_T::value_type> op_max;
@@ -131,28 +129,29 @@ template <class data_T, class res_T, typename CONFIG_T> void softmax_stable(cons
         reduce<typename data_T::value_type, CONFIG_T::n_in, Op_max<typename data_T::value_type>>(data.data(), op_max);
 
     // For the diffs, use the same type as the input but force rounding and saturation
-    [[intel::fpga_register]] ac_fixed<data_T::value_type::width, data_T::value_type::i_width, true, AC_RND, AC_SAT>
-        d_xi_xmax[CONFIG_T::n_in];
+    [[intel::fpga_register]] typename CONFIG_T::inp_norm_t d_xi_xmax[CONFIG_T::n_in];
     #pragma unroll
     for (unsigned i = 0; i < CONFIG_T::n_in; i++) {
-        d_xi_xmax[i] = data[i] - x_max;
+        d_xi_xmax[i] = x_max - data[i];
     }
 
     // Calculate all the e^x's
-    [[intel::fpga_register]] typename CONFIG_T::exp_table_t exp_res[CONFIG_T::n_in];
+    [[intel::fpga_register]] typename CONFIG_T::accum_t exp_res[CONFIG_T::n_in];
     #pragma unroll
     for (unsigned i = 0; i < CONFIG_T::n_in; i++) {
-        exp_res[i] = exp_table[softmax_stable_idx_from_real_val<typename data_T::value_type, CONFIG_T>(d_xi_xmax[i])];
+        exp_res[i] = CONFIG_T::exp_table[softmax_idx_from_real_val<typename CONFIG_T::inp_norm_t, CONFIG_T::exp_table_size>(
+            d_xi_xmax[i])]; // input_t, CONFIG_T
     }
 
     // Explicitly sum previously calculated exponentials with an adder tree
-    Op_add<typename CONFIG_T::exp_table_t> op_add;
-    [[intel::fpga_register]] typename CONFIG_T::exp_table_t exp_sum =
-        reduce<typename CONFIG_T::exp_table_t, CONFIG_T::n_in, Op_add<typename CONFIG_T::exp_table_t>>(exp_res, op_add);
+    Op_add<typename CONFIG_T::accum_t> op_add;
+    [[intel::fpga_register]] typename CONFIG_T::inv_inp_t exp_sum =
+        reduce<typename CONFIG_T::accum_t, CONFIG_T::n_in, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
 
     // Multiply previously calculated exponetials with the reciprocal of the sum
     [[intel::fpga_register]] typename CONFIG_T::inv_table_t inv_exp_sum =
-        invert_table[softmax_stable_idx_from_real_val<typename CONFIG_T::exp_table_t, CONFIG_T>(exp_sum)];
+        CONFIG_T::invert_table[softmax_idx_from_real_val<typename CONFIG_T::inv_inp_t, CONFIG_T::inv_table_size>(exp_sum)];
+
     #pragma unroll
     for (unsigned i = 0; i < CONFIG_T::n_in; i++) {
         res[i] = exp_res[i] * inv_exp_sum;
@@ -161,14 +160,13 @@ template <class data_T, class res_T, typename CONFIG_T> void softmax_stable(cons
 
 // TODO - Improve accuracy
 template <class data_T, class res_T, typename CONFIG_T> void softmax_latency(const data_T &data, res_T &res) {
-#include "activation_tables/exp_table_latency.tb"
-#include "activation_tables/invert_table_latency.tb"
 
     // Calculate all the e^x's
     [[intel::fpga_register]] typename CONFIG_T::exp_table_t exp_res[CONFIG_T::n_in];
     #pragma unroll
     for (unsigned i = 0; i < CONFIG_T::n_in; i++) {
-        exp_res[i] = exp_table_latency[softmax_latency_idx_from_real_val<typename data_T::value_type, CONFIG_T>(data[i])];
+        exp_res[i] =
+            CONFIG_T::exp_table[softmax_idx_from_real_val<typename data_T::value_type, CONFIG_T::exp_table_size>(data[i])];
     }
 
     // Explicitly sum the results with an adder tree.
@@ -178,7 +176,7 @@ template <class data_T, class res_T, typename CONFIG_T> void softmax_latency(con
 
     // Multiply previously calculated exponetials with the reciprocal of the sum
     [[intel::fpga_register]] typename CONFIG_T::inv_table_t inv_exp_sum =
-        invert_table_latency[softmax_latency_idx_from_real_val<typename CONFIG_T::exp_table_t, CONFIG_T>(exp_sum)];
+        CONFIG_T::invert_table[softmax_idx_from_real_val<typename CONFIG_T::exp_table_t, CONFIG_T::inv_table_size>(exp_sum)];
     #pragma unroll
     for (unsigned i = 0; i < CONFIG_T::n_in; i++) {
         res[i] = exp_res[i] * inv_exp_sum;
@@ -186,8 +184,6 @@ template <class data_T, class res_T, typename CONFIG_T> void softmax_latency(con
 }
 
 template <class data_T, class res_T, typename CONFIG_T> void softmax_legacy(const data_T &data, res_T &res) {
-#include "activation_tables/exp_table_legacy.tb"
-#include "activation_tables/invert_table_legacy.tb"
 
     [[intel::fpga_register]] int data_round[CONFIG_T::n_in];
 New_loop:
@@ -213,7 +209,7 @@ template <class data_T, class res_T, typename CONFIG_T> void softmax_legacy(cons
                 if (index > CONFIG_T::table_size - 1)
                     index = CONFIG_T::table_size - 1;
 
-                typename CONFIG_T::exp_table_t temp_exp = exp_table_legacy[index];
+                typename CONFIG_T::exp_table_t temp_exp = CONFIG_T::exp_table[index];
                 exp_res_temp += temp_exp;
             }
         }
@@ -222,7 +218,7 @@ template <class data_T, class res_T, typename CONFIG_T> void softmax_legacy(cons
             exp_res_index = 0;
         if (exp_res_index > CONFIG_T::table_size - 1)
             exp_res_index = CONFIG_T::table_size - 1;
-        res[ii] = invert_table_legacy[exp_res_index];
+        res[ii] = CONFIG_T::invert_table[exp_res_index];
     }
 }
 
@@ -246,25 +242,56 @@ template <class data_T, class res_T, typename CONFIG_T> void softmax_argmax(cons
 }
 
 template <class data_T, class res_T, typename CONFIG_T> inline void softmax(const data_T &data, res_T &res) {
-    switch (CONFIG_T::implementation) {
-    case softmax_implementation::stable:
-        softmax_stable<data_T, res_T, CONFIG_T>(data, res);
-        break;
-    case softmax_implementation::latency:
+    if constexpr (CONFIG_T::implementation == softmax_implementation::latency) {
         softmax_latency<data_T, res_T, CONFIG_T>(data, res);
-        break;
-    case softmax_implementation::legacy:
+    } else if constexpr (CONFIG_T::implementation == softmax_implementation::argmax) {
+        softmax_argmax<data_T, res_T, CONFIG_T>(data, res);
+    } else if constexpr (CONFIG_T::implementation == softmax_implementation::legacy) {
         softmax_legacy<data_T, res_T, CONFIG_T>(data, res);
-        break;
-    default:
+    } else {
         softmax_stable<data_T, res_T, CONFIG_T>(data, res);
-        break;
-    case softmax_implementation::argmax:
-        softmax_argmax<data_T, res_T, CONFIG_T>(data, res);
-        break;
     }
 }
 
+// *************************************************
+//       Multidimensional Softmax
+// *************************************************
+
+// Helper to remap the config for the core softmax function
+template <class CONFIG_T> struct softmax_multidim_slice_config : CONFIG_T {
+    static constexpr unsigned n_in = CONFIG_T::n_slice;
+};
+
+template <class data_T, class res_T, typename CONFIG_T> inline void softmax_multidim(const data_T &data, res_T &res) {
+    using buffer_data_t = std::array<typename data_T::value_type, CONFIG_T::n_slice>;
+    using buffer_res_t = std::array<typename res_T::value_type, CONFIG_T::n_slice>;
+    using slice_config = softmax_multidim_slice_config<CONFIG_T>;
+
+    #pragma unroll
+    for (unsigned i = 0; i < CONFIG_T::n_outer; i++) {
+        #pragma unroll
+        for (unsigned k = 0; k < CONFIG_T::n_inner; k++) {
+
+            [[intel::fpga_register]] buffer_data_t buffer_in;
+            [[intel::fpga_register]] buffer_res_t buffer_out;
+
+            // Gather Phase
+            #pragma unroll
+            for (unsigned j = 0; j < CONFIG_T::n_slice; j++) {
+                unsigned idx = (i * CONFIG_T::n_slice * CONFIG_T::n_inner) + (j * CONFIG_T::n_inner) + k;
+                buffer_in[j] = data[idx];
+            }
+
+            nnet::softmax<buffer_data_t, buffer_res_t, slice_config>(buffer_in, buffer_out);
+
+            #pragma unroll
+            for (unsigned j = 0; j < CONFIG_T::n_slice; j++) {
+                unsigned idx = (i * CONFIG_T::n_slice * CONFIG_T::n_inner) + (j * CONFIG_T::n_inner) + k;
+                res[idx] = buffer_out[j];
+            }
+        }
+    }
+}
 // *************************************************
 //       TanH Activation
 // *************************************************
diff --git a/hls4ml/templates/oneapi/firmware/nnet_utils/nnet_activation_stream.h b/hls4ml/templates/oneapi/firmware/nnet_utils/nnet_activation_stream.h
index e860c38988..cb65137927 100644
--- a/hls4ml/templates/oneapi/firmware/nnet_utils/nnet_activation_stream.h
+++ b/hls4ml/templates/oneapi/firmware/nnet_utils/nnet_activation_stream.h
@@ -271,73 +271,62 @@ template <class data_pipe, class res_pipe, typename CONFIG_T> void softsign_stre
 // *************************************************
 
 template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_stable_stream() {
-#include "activation_tables/exp_table.tb"
-#include "activation_tables/invert_table.tb"
 
-    constexpr unsigned multiplier_limit =
-        DIV_ROUNDUP(std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}, CONFIG_T::reuse_factor);
-    constexpr unsigned pipeline = std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{} / multiplier_limit;
+    using input_arr_t = typename ExtractPipeType<data_pipe>::value_type;
+    using input_t = typename ExtractPipeType<data_pipe>::value_type::value_type;
+    constexpr unsigned input_arr_size = std::tuple_size<input_arr_t>{};
 
-    [[intel::fpga_register]] typename ExtractPipeType<data_pipe>::value_type::value_type
-        data_array[std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}];
+    constexpr unsigned multiplier_limit = DIV_ROUNDUP(input_arr_size, CONFIG_T::reuse_factor);
+    constexpr unsigned pipeline = input_arr_size / multiplier_limit;
+
+    [[intel::fpga_register]] input_t data_array[input_arr_size];
 
 SoftmaxArrayLoop:
-    [[intel::initiation_interval(pipeline)]] for (unsigned i = 0;
-                                                  i < CONFIG_T::n_in /
-                                                          std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{};
-                                                  i++) {
+    [[intel::initiation_interval(pipeline)]] for (unsigned i = 0; i < CONFIG_T::n_in / input_arr_size; i++) {
         auto in_pack = data_pipe::read();
 
     SoftmaxArrayPackLoop:
         #pragma unroll
-        for (unsigned j = 0; j < std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}; j++) {
+        for (unsigned j = 0; j < input_arr_size; j++) {
             data_array[j] = in_pack[j];
         }
 
         // Find the max and compute all delta(x_i, x_max)
-        Op_max<typename ExtractPipeType<data_pipe>::value_type::value_type> op_max;
-        [[intel::fpga_register]] typename ExtractPipeType<data_pipe>::value_type::value_type x_max =
-            reduce<typename ExtractPipeType<data_pipe>::value_type::value_type,
-                   std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{},
-                   Op_max<typename ExtractPipeType<data_pipe>::value_type::value_type>>(data_array, op_max);
-
-        // For the diffs, use the same type as the input but force rounding and saturation
-        [[intel::fpga_register]] ac_fixed<ExtractPipeType<data_pipe>::value_type::value_type::width,
-                                          ExtractPipeType<data_pipe>::value_type::value_type::i_width, true, AC_RND, AC_SAT>
-            d_xi_xmax[std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}];
+        Op_max<input_t> op_max;
+        [[intel::fpga_register]] input_t x_max = reduce<input_t, input_arr_size, Op_max<input_t>>(data_array, op_max);
+
+        [[intel::fpga_register]] typename CONFIG_T::inp_norm_t d_xi_xmax[input_arr_size];
+
         #pragma unroll
-        for (unsigned j = 0; j < std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}; j++) {
-            d_xi_xmax[j] = data_array[j] - x_max;
+        for (unsigned j = 0; j < input_arr_size; j++) {
+            d_xi_xmax[j] = x_max - data_array[j];
         }
 
         // Calculate all the e^x's
-        [[intel::fpga_register]]
-        typename CONFIG_T::exp_table_t exp_res[std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}];
+        [[intel::fpga_register]] typename CONFIG_T::accum_t exp_res[input_arr_size];
+
         #pragma unroll
-        for (unsigned j = 0; j < std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}; j++) {
+        for (unsigned j = 0; j < input_arr_size; j++) {
             exp_res[j] =
-                exp_table[softmax_stable_idx_from_real_val<typename ExtractPipeType<data_pipe>::value_type::value_type,
-                                                           CONFIG_T>(d_xi_xmax[j])];
+                CONFIG_T::exp_table[softmax_idx_from_real_val<typename CONFIG_T::inp_norm_t, CONFIG_T::exp_table_size>(
+                    d_xi_xmax[j])];
         }
 
         // Explicitly sum the results with an adder tree.
         // Rounding & Saturation mode, which improve accuracy, prevent Vivado from expression balancing
-        Op_add<typename CONFIG_T::exp_table_t> op_add;
-        [[intel::fpga_register]] typename CONFIG_T::exp_table_t exp_sum =
-            reduce<typename CONFIG_T::exp_table_t, std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{},
-                   Op_add<typename CONFIG_T::exp_table_t>>(exp_res, op_add);
+        Op_add<typename CONFIG_T::accum_t> op_add;
+        [[intel::fpga_register]] typename CONFIG_T::inv_inp_t exp_sum =
+            reduce<typename CONFIG_T::accum_t, input_arr_size, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
 
         [[intel::fpga_register]] typename CONFIG_T::inv_table_t inv_exp_sum =
-            invert_table[softmax_stable_idx_from_real_val<typename CONFIG_T::exp_table_t, CONFIG_T>(exp_sum)];
+            CONFIG_T::invert_table[softmax_idx_from_real_val<typename CONFIG_T::inv_inp_t, CONFIG_T::inv_table_size>(
+                exp_sum)];
+
         typename ExtractPipeType<res_pipe>::value_type out_pack;
 
     SoftmaxInvPackLoop:
         #pragma unroll
         for (unsigned j = 0; j < std::tuple_size<typename ExtractPipeType<res_pipe>::value_type>{}; j++) {
-
-            // TODO - Find Quartus-equivalent pragma
-            // #pragma HLS ALLOCATION instances=mul limit=multiplier_limit operation
-
             out_pack[j] = exp_res[j] * inv_exp_sum;
         }
 
@@ -346,8 +335,6 @@ template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_stabl
 }
 
 template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_latency_stream() {
-#include "activation_tables/exp_table_latency.tb"
-#include "activation_tables/invert_table_latency.tb"
 
     constexpr unsigned multiplier_limit =
         DIV_ROUNDUP(std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}, CONFIG_T::reuse_factor);
@@ -367,8 +354,9 @@ template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_laten
     SoftmaxExpPackLoop:
         #pragma unroll
         for (unsigned j = 0; j < std::tuple_size<typename ExtractPipeType<data_pipe>::value_type>{}; j++) {
-            exp_res[j] = exp_table_latency[softmax_latency_idx_from_real_val<
-                typename ExtractPipeType<data_pipe>::value_type::value_type, CONFIG_T>(in_pack[j])];
+            exp_res[j] =
+                CONFIG_T::exp_table[softmax_idx_from_real_val<typename ExtractPipeType<data_pipe>::value_type::value_type,
+                                                              CONFIG_T::exp_table_size>(in_pack[j])];
         }
 
         // Explicitly sum the results with an adder tree.
@@ -379,7 +367,8 @@ template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_laten
 
         // Multiply previously calculated exponetials with the reciprocal of the sum
         [[intel::fpga_register]] typename CONFIG_T::inv_table_t inv_exp_sum =
-            invert_table_latency[softmax_latency_idx_from_real_val<typename CONFIG_T::exp_table_t, CONFIG_T>(exp_sum)];
+            CONFIG_T::invert_table[softmax_idx_from_real_val<typename CONFIG_T::exp_table_t, CONFIG_T::inv_table_size>(
+                exp_sum)];
 
         typename ExtractPipeType<res_pipe>::value_type out_pack;
     SoftmaxInvPackLoop:
@@ -394,8 +383,6 @@ template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_laten
 }
 
 template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_legacy_stream() {
-#include "activation_tables/exp_table_legacy.tb"
-#include "activation_tables/invert_table_legacy.tb"
 
     // Index into the lookup table based on data for exponentials
     [[intel::fpga_register]]
@@ -433,7 +420,7 @@ template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_legac
                         index = 0;
                     if (index > CONFIG_T::table_size - 1)
                         index = CONFIG_T::table_size - 1;
-                    exp_diff_res = exp_table_legacy[index];
+                    exp_diff_res = CONFIG_T::exp_table[index];
                 }
                 exp_res[i] += exp_diff_res;
             }
@@ -448,8 +435,8 @@ template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_legac
                 exp_res_index = 0;
             if (exp_res_index > CONFIG_T::table_size - 1)
                 exp_res_index = CONFIG_T::table_size - 1;
-            out_pack[j] =
-                static_cast<typename ExtractPipeType<res_pipe>::value_type::value_type>(invert_table_legacy[exp_res_index]);
+            out_pack[j] = static_cast<typename ExtractPipeType<res_pipe>::value_type::value_type>(
+                CONFIG_T::invert_table[exp_res_index]);
         }
 
         res_pipe::write(out_pack);
@@ -484,25 +471,64 @@ template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_argma
 }
 
 template <class data_pipe, class res_pipe, typename CONFIG_T> void softmax_stream() {
-    switch (CONFIG_T::implementation) {
-    case softmax_implementation::latency:
+    if constexpr (CONFIG_T::implementation == softmax_implementation::latency) {
         softmax_latency_stream<data_pipe, res_pipe, CONFIG_T>();
-        break;
-    case softmax_implementation::stable:
-        softmax_stable_stream<data_pipe, res_pipe, CONFIG_T>();
-        break;
-    case softmax_implementation::legacy:
-        softmax_legacy_stream<data_pipe, res_pipe, CONFIG_T>();
-        break;
-    case softmax_implementation::argmax:
+    } else if constexpr (CONFIG_T::implementation == softmax_implementation::argmax) {
+        softmax_argmax_stream<data_pipe, res_pipe, CONFIG_T>();
+    } else if constexpr (CONFIG_T::implementation == softmax_implementation::legacy) {
         softmax_argmax_stream<data_pipe, res_pipe, CONFIG_T>();
-        break;
-    default:
+    } else { // Default to stable
         softmax_stable_stream<data_pipe, res_pipe, CONFIG_T>();
-        break;
     }
 }
 
+// *************************************************
+//       Multidimensional Softmax
+// *************************************************
+template <class data_pipe, class res_pipe, typename CONFIG_T> inline void softmax_multidim_stream() {
+
+    using data_pipe_T = typename ExtractPipeType<data_pipe>::value_type;
+    using data_T = typename data_pipe_T::value_type;
+    using res_pipe_T = typename ExtractPipeType<res_pipe>::value_type;
+    using res_T = typename res_pipe_T::value_type;
+    using in_slice_arr_T = nnet::array<data_T, CONFIG_T::n_slice>;
+    using out_slice_arr_T = nnet::array<res_T, CONFIG_T::n_slice>;
+
+    using slice_config = softmax_multidim_slice_config<CONFIG_T>;
+
+    [[intel::fpga_register]] data_pipe_T buffer_in;
+    [[intel::fpga_register]] res_pipe_T buffer_out;
+    [[intel::fpga_register]] in_slice_arr_T smax_slice_in;
+    [[intel::fpga_register]] out_slice_arr_T smax_slice_out;
+
+    buffer_in = data_pipe::read();
+
+    #pragma unroll
+    for (unsigned i = 0; i < CONFIG_T::n_outer; i++) {
+        unsigned outer_offset = i * CONFIG_T::n_slice * CONFIG_T::n_inner;
+        #pragma unroll
+        for (unsigned k = 0; k < CONFIG_T::n_inner; k++) {
+
+            // TODO: Access might be inefficient consider data rearrangement
+            #pragma unroll
+            for (unsigned j = 0; j < CONFIG_T::n_slice; j++) {
+                unsigned idx = outer_offset + j * CONFIG_T::n_inner + k;
+                smax_slice_in[j] = buffer_in[idx];
+            }
+
+            nnet::softmax<in_slice_arr_T, out_slice_arr_T, slice_config>(smax_slice_in, smax_slice_out);
+
+            #pragma unroll
+            for (unsigned j = 0; j < CONFIG_T::n_slice; j++) {
+                unsigned idx = outer_offset + j * CONFIG_T::n_inner + k;
+                buffer_out[idx] = smax_slice_out[j];
+            }
+        }
+    }
+
+    res_pipe::write(buffer_out);
+}
+
 // *************************************************
 //       TanH Activation
 // *************************************************
diff --git a/hls4ml/templates/oneapi/firmware/parameters.h b/hls4ml/templates/oneapi/firmware/parameters.h
index 717059f1e8..ef4e5d26b9 100644
--- a/hls4ml/templates/oneapi/firmware/parameters.h
+++ b/hls4ml/templates/oneapi/firmware/parameters.h
@@ -6,6 +6,8 @@
 #include "nnet_utils/nnet_code_gen.h"
 #include "nnet_utils/nnet_helpers.h"
 
+// hls-fpga-machine-learning insert softmax tables
+
 // hls-fpga-machine-learning insert includes
 
 // hls-fpga-machine-learning insert layer-config
diff --git a/hls4ml/templates/vivado/nnet_utils/nnet_activation.h b/hls4ml/templates/vivado/nnet_utils/nnet_activation.h
index ac85e0b2cc..b262c86859 100644
--- a/hls4ml/templates/vivado/nnet_utils/nnet_activation.h
+++ b/hls4ml/templates/vivado/nnet_utils/nnet_activation.h
@@ -189,14 +189,13 @@ void softmax_latency(data_T data[CONFIG_T::n_slice], res_T res[CONFIG_T::n_slice
         // Note we are exponentiating the inputs, which have type data_T
         init_exp_table<data_T, CONFIG_T>(exp_table);
         // Note we are inverting the exponentials, which have type exp_table_t
-        init_invert_table<typename CONFIG_T::inv_inp_t, CONFIG_T>(invert_table);
+        init_invert_table<typename CONFIG_T::accum_t, CONFIG_T>(invert_table);
         initialized = true;
     }
 
     // Calculate all the e^x's
     typename CONFIG_T::accum_t exp_res[CONFIG_T::n_slice];
     #pragma HLS array_partition variable=exp_res complete
-    typename CONFIG_T::inv_inp_t exp_sum(0);
     for (unsigned i = 0; i < CONFIG_T::n_slice; i++) {
         #pragma HLS unroll
         unsigned x = softmax_idx_from_real_val<data_T, CONFIG_T::exp_table_size>(data[i]);
@@ -206,10 +205,11 @@ void softmax_latency(data_T data[CONFIG_T::n_slice], res_T res[CONFIG_T::n_slice
     // Explicitly sum the results with an adder tree.
     // Rounding & Saturation mode, which improve accuracy, prevent Vivado from expression balancing
     Op_add<typename CONFIG_T::accum_t> op_add;
-    exp_sum = reduce<typename CONFIG_T::accum_t, CONFIG_T::n_slice, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
+    typename CONFIG_T::accum_t exp_sum =
+        reduce<typename CONFIG_T::accum_t, CONFIG_T::n_slice, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
 
     typename CONFIG_T::inv_table_t inv_exp_sum =
-        invert_table[softmax_idx_from_real_val<typename CONFIG_T::inv_inp_t, CONFIG_T::inv_table_size>(exp_sum)];
+        invert_table[softmax_idx_from_real_val<typename CONFIG_T::accum_t, CONFIG_T::inv_table_size>(exp_sum)];
     for (unsigned i = 0; i < CONFIG_T::n_slice; i++) {
         #pragma HLS unroll
         res[i] = exp_res[i] * inv_exp_sum;
@@ -251,7 +251,6 @@ void softmax_stable(data_T data[CONFIG_T::n_slice], res_T res[CONFIG_T::n_slice]
     // Calculate all the e^x's
     typename CONFIG_T::accum_t exp_res[CONFIG_T::n_slice];
     #pragma HLS array_partition variable=exp_res complete
-    typename CONFIG_T::inv_inp_t exp_sum(0);
     for (unsigned i = 0; i < CONFIG_T::n_slice; i++) {
         #pragma HLS unroll
         unsigned x = softmax_idx_from_real_val<typename CONFIG_T::inp_norm_t, CONFIG_T::exp_table_size>(d_xi_xmax[i]);
@@ -261,7 +260,8 @@ void softmax_stable(data_T data[CONFIG_T::n_slice], res_T res[CONFIG_T::n_slice]
     // Explicitly sum the results with an adder tree.
     // Rounding & Saturation mode, which improve accuracy, prevent Vivado from expression balancing
     Op_add<typename CONFIG_T::accum_t> op_add;
-    exp_sum = reduce<typename CONFIG_T::accum_t, CONFIG_T::n_slice, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
+    typename CONFIG_T::inv_inp_t exp_sum =
+        reduce<typename CONFIG_T::accum_t, CONFIG_T::n_slice, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
 
     typename CONFIG_T::inv_table_t inv_exp_sum =
         invert_table[softmax_idx_from_real_val<typename CONFIG_T::inv_inp_t, CONFIG_T::inv_table_size>(exp_sum)];
@@ -271,18 +271,18 @@ void softmax_stable(data_T data[CONFIG_T::n_slice], res_T res[CONFIG_T::n_slice]
     }
 }
 
-template <typename CONFIG_T, int N_TABLE> void init_exp_table_legacy(typename CONFIG_T::table_t table_out[N_TABLE]) {
+template <typename CONFIG_T, int N_TABLE> void init_exp_table_legacy(typename CONFIG_T::exp_table_t table_out[N_TABLE]) {
     for (int ii = 0; ii < N_TABLE; ii++) {
         // First, convert from table index to X-value (signed 8-bit, range -8 to +8)
         float in_val = 2 * 8.0 * (ii - float(N_TABLE) / 2.0) / float(N_TABLE);
         // Next, compute lookup table function
-        typename CONFIG_T::table_t real_val = exp_fcn_float(in_val);
+        typename CONFIG_T::exp_table_t real_val = exp_fcn_float(in_val);
         // std::cout << "Lookup table In Value: " << in_val << " Result: " << real_val << std::endl;
         table_out[ii] = real_val;
     }
 }
 
-template <typename CONFIG_T, int N_TABLE> void init_invert_table_legacy(typename CONFIG_T::table_t table_out[N_TABLE]) {
+template <typename CONFIG_T, int N_TABLE> void init_invert_table_legacy(typename CONFIG_T::inv_table_t table_out[N_TABLE]) {
     // Inversion function:
     //   result = 1/x
     for (int ii = 0; ii < N_TABLE; ii++) {
@@ -301,12 +301,12 @@ void softmax_legacy(data_T data[CONFIG_T::n_slice], res_T res[CONFIG_T::n_slice]
     // Initialize the lookup table
 #ifdef __HLS_SYN__
     bool initialized = false;
-    typename CONFIG_T::table_t exp_table[CONFIG_T::exp_table_size];
-    typename CONFIG_T::table_t invert_table[CONFIG_T::inv_table_size];
+    typename CONFIG_T::exp_table_t exp_table[CONFIG_T::exp_table_size];
+    typename CONFIG_T::inv_table_t invert_table[CONFIG_T::inv_table_size];
 #else
     static bool initialized = false;
-    static typename CONFIG_T::table_t exp_table[CONFIG_T::exp_table_size];
-    static typename CONFIG_T::table_t invert_table[CONFIG_T::inv_table_size];
+    static typename CONFIG_T::exp_table_t exp_table[CONFIG_T::exp_table_size];
+    static typename CONFIG_T::inv_table_t invert_table[CONFIG_T::inv_table_size];
 #endif
     if (!initialized) {
         init_exp_table_legacy<CONFIG_T, CONFIG_T::exp_table_size>(exp_table);
@@ -317,22 +317,23 @@ void softmax_legacy(data_T data[CONFIG_T::n_slice], res_T res[CONFIG_T::n_slice]
     #pragma HLS PIPELINE
 
     // Index into the lookup table based on data for exponentials
-    typename CONFIG_T::table_t exp_res[CONFIG_T::n_slice]; // different, independent, fixed point precision
-    typename CONFIG_T::table_t exp_diff_res;               // different, independent, fixed point precision
+    typename CONFIG_T::accum_t exp_res[CONFIG_T::n_slice]; // different, independent, fixed point precision
+    typename CONFIG_T::exp_table_t exp_diff_res;           // different, independent, fixed point precision
     data_T data_cache[CONFIG_T::n_slice];
-    int data_round;
     int index;
+
     for (int ii = 0; ii < CONFIG_T::n_slice; ii++) {
         data_cache[ii] = data[ii];
         exp_res[ii] = 0;
     }
 
+    // first calculate 1/softmax as a sum over fractions.
     for (int ii = 0; ii < CONFIG_T::n_slice; ii++) {
         for (int jj = 0; jj < CONFIG_T::n_slice; jj++) {
             if (ii == jj)
                 exp_diff_res = 1;
             else {
-                data_round = (data_cache[jj] - data_cache[ii]) * CONFIG_T::exp_table_size / 16;
+                auto data_round = (data_cache[jj] - data_cache[ii]) * CONFIG_T::exp_table_size / 16;
                 index = data_round + 8 * CONFIG_T::exp_table_size / 16;
                 if (index < 0)
                     index = 0;
@@ -352,7 +353,7 @@ void softmax_legacy(data_T data[CONFIG_T::n_slice], res_T res[CONFIG_T::n_slice]
         if (exp_res_index > CONFIG_T::inv_table_size - 1)
             exp_res_index = CONFIG_T::inv_table_size - 1;
         // typename CONFIG_T::table_t exp_res_invert = invert_table[exp_res_index];
-        res[ii] = (res_T)invert_table[exp_res_index];
+        res[ii] = static_cast<res_T>(invert_table[exp_res_index]);
     }
 }
 
diff --git a/hls4ml/templates/vivado/nnet_utils/nnet_activation_stream.h b/hls4ml/templates/vivado/nnet_utils/nnet_activation_stream.h
index 50c6c4068c..814ed2f50a 100644
--- a/hls4ml/templates/vivado/nnet_utils/nnet_activation_stream.h
+++ b/hls4ml/templates/vivado/nnet_utils/nnet_activation_stream.h
@@ -120,8 +120,8 @@ void softmax_latency(hls::stream<data_T> &data, hls::stream<res_T> &res) {
     if (!initialized) {
         // Note we are exponentiating the inputs, which have type data_T
         init_exp_table<typename data_T::value_type, CONFIG_T>(exp_table);
-        // Note we are inverting the exponentials, which have type exp_table_t
-        init_invert_table<typename CONFIG_T::inv_inp_t, CONFIG_T>(invert_table);
+        // Note we are inverting the summed exponentials, which have type accum_t
+        init_invert_table<typename CONFIG_T::accum_t, CONFIG_T>(invert_table);
         initialized = true;
     }
 
@@ -150,7 +150,7 @@ void softmax_latency(hls::stream<data_T> &data, hls::stream<res_T> &res) {
         exp_sum = reduce<typename CONFIG_T::accum_t, data_T::size, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
 
         typename CONFIG_T::inv_table_t inv_exp_sum =
-            invert_table[softmax_idx_from_real_val<typename CONFIG_T::inv_inp_t, CONFIG_T::inv_table_size>(exp_sum)];
+            invert_table[softmax_idx_from_real_val<typename CONFIG_T::accum_t, CONFIG_T::inv_table_size>(exp_sum)];
 
         res_T out_pack;
         PRAGMA_DATA_PACK(out_pack)
@@ -216,7 +216,6 @@ void softmax_stable(hls::stream<data_T> &data, hls::stream<res_T> &res) {
         // Calculate all the e^x's
         typename CONFIG_T::accum_t exp_res[data_T::size];
         #pragma HLS ARRAY_PARTITION variable=exp_res complete
-        typename CONFIG_T::inv_inp_t exp_sum(0);
         for (unsigned j = 0; j < data_T::size; j++) {
             #pragma HLS UNROLL
             unsigned x = softmax_idx_from_real_val<typename CONFIG_T::inp_norm_t, CONFIG_T::exp_table_size>(d_xi_xmax[j]);
@@ -226,7 +225,8 @@ void softmax_stable(hls::stream<data_T> &data, hls::stream<res_T> &res) {
         // Explicitly sum the results with an adder tree.
         // Rounding & Saturation mode, which improve accuracy, prevent Vivado from expression balancing
         Op_add<typename CONFIG_T::accum_t> op_add;
-        exp_sum = reduce<typename CONFIG_T::accum_t, data_T::size, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
+        typename CONFIG_T::inv_inp_t exp_sum =
+            reduce<typename CONFIG_T::accum_t, data_T::size, Op_add<typename CONFIG_T::accum_t>>(exp_res, op_add);
 
         typename CONFIG_T::inv_table_t inv_exp_sum =
             invert_table[softmax_idx_from_real_val<typename CONFIG_T::inv_inp_t, CONFIG_T::inv_table_size>(exp_sum)];
@@ -249,22 +249,22 @@ void softmax_legacy(hls::stream<data_T> &data, hls::stream<res_T> &res) {
     // Initialize the lookup table
 #ifdef __HLS_SYN__
     bool initialized = false;
-    typename CONFIG_T::table_t exp_table[CONFIG_T::table_size];
-    typename CONFIG_T::table_t invert_table[CONFIG_T::table_size];
+    typename CONFIG_T::exp_table_t exp_table[CONFIG_T::exp_table_size];
+    typename CONFIG_T::inv_table_t invert_table[CONFIG_T::inv_table_size];
 #else
     static bool initialized = false;
-    static typename CONFIG_T::table_t exp_table[CONFIG_T::table_size];
-    static typename CONFIG_T::table_t invert_table[CONFIG_T::table_size];
+    static typename CONFIG_T::exp_table_t exp_table[CONFIG_T::exp_table_size];
+    static typename CONFIG_T::inv_table_t invert_table[CONFIG_T::inv_table_size];
 #endif
     if (!initialized) {
-        init_exp_table_legacy<CONFIG_T, CONFIG_T::table_size>(exp_table);
-        init_invert_table_legacy<CONFIG_T, CONFIG_T::table_size>(invert_table);
+        init_exp_table_legacy<CONFIG_T, CONFIG_T::exp_table_size>(exp_table);
+        init_invert_table_legacy<CONFIG_T, CONFIG_T::inv_table_size>(invert_table);
         initialized = true;
     }
 
     // Index into the lookup table based on data for exponentials
-    typename CONFIG_T::table_t exp_res[data_T::size];
-    typename CONFIG_T::table_t exp_diff_res;
+    typename CONFIG_T::accum_t exp_res[data_T::size];
+    typename CONFIG_T::exp_table_t exp_diff_res;
     typename data_T::value_type data_cache[data_T::size];
 
 SoftmaxInitLoop:
@@ -288,12 +288,12 @@ void softmax_legacy(hls::stream<data_T> &data, hls::stream<res_T> &res) {
                 if (i == j) {
                     exp_diff_res = 1;
                 } else {
-                    int data_round = (data_cache[j] - data_cache[i]) * CONFIG_T::table_size / 16;
+                    auto data_round = (data_cache[j] - data_cache[i]) * CONFIG_T::table_size / 16;
                     int index = data_round + 8 * CONFIG_T::table_size / 16;
                     if (index < 0)
                         index = 0;
-                    if (index > CONFIG_T::table_size - 1)
-                        index = CONFIG_T::table_size - 1;
+                    if (index > CONFIG_T::exp_table_size - 1)
+                        index = CONFIG_T::exp_table_size - 1;
                     exp_diff_res = exp_table[index];
                 }
 
@@ -311,10 +311,10 @@ void softmax_legacy(hls::stream<data_T> &data, hls::stream<res_T> &res) {
             int exp_res_index = exp_res[j] * CONFIG_T::table_size / 64;
             if (exp_res_index < 0)
                 exp_res_index = 0;
-            if (exp_res_index > CONFIG_T::table_size - 1)
-                exp_res_index = CONFIG_T::table_size - 1;
+            if (exp_res_index > CONFIG_T::inv_table_size - 1)
+                exp_res_index = CONFIG_T::inv_table_size - 1;
 
-            out_pack[j] = (typename res_T::value_type)invert_table[exp_res_index];
+            out_pack[j] = static_cast<typename res_T::value_type>(invert_table[exp_res_index]);
         }
         res.write(out_pack);
     }
diff --git a/hls4ml/writer/oneapi_writer.py b/hls4ml/writer/oneapi_writer.py
index 3c0a778c50..d5d769e22e 100644
--- a/hls4ml/writer/oneapi_writer.py
+++ b/hls4ml/writer/oneapi_writer.py
@@ -302,6 +302,19 @@ def write_parameters(self, model):
                         config = layer.get_attr('config_cpp', None)
                         if config:
                             newline += config + '\n'
+
+                elif '// hls-fpga-machine-learning insert softmax tables' in line:
+                    newline = line
+                    for layer in model.get_layers():
+                        if (
+                            layer.get_attr('activation') == 'softmax'
+                            or layer.get_attr('recurrent_activation') == 'softmax'
+                            or layer.get_attr('activation') == 'softmax_multidim'
+                            or layer.get_attr('recurrent_activation') == 'softmax_multidim'
+                        ) and 'implementation' in layer.attributes:
+                            newline += f'#include "nnet_utils/activation_tables/{layer.name}_exp_table.h"\n'
+                            newline += f'#include "nnet_utils/activation_tables/{layer.name}_inv_table.h"\n'
+
                 else:
                     newline = line
                 fout.write(newline)
@@ -549,16 +562,16 @@ def write_nnet_utils(self, model):
             dstpath = f'{model.config.get_output_dir()}/src/firmware/{dst}'
             copyfile(srcpath, dstpath)
 
-    def __get_table_size(self, model, activation):
+    def __get_table_size(self, model, activation, size_attr_name='table_size'):
         for layer in model.get_layers():
             if (
                 layer.get_attr('activation') == activation or layer.get_attr('recurrent_activation') == activation
-            ) and layer.get_attr('table_size') is not None:
-                return int(layer.get_attr('table_size'))
+            ) and layer.get_attr(size_attr_name) is not None:
+                return int(layer.get_attr(size_attr_name))
         return 1024
 
-    def __get_table_header(self, table_name, table_size):
-        table_header = f'static const typename CONFIG_T::table_t {table_name}[{table_size}] = {{'
+    def __get_table_header(self, table_name, table_size, table_type='table_t'):
+        table_header = f'static const typename CONFIG_T::{table_type} {table_name}[{table_size}] = {{'
         return table_header
 
     def __write_elu_table(self, model, path):
@@ -687,208 +700,283 @@ def __write_selu_table(self, model, path):
         h_file.write('};\n')
         h_file.close()
 
-    def __write_exp_table(self, model, path):
-        table_name = 'exp_table'
-        table_size = self.__get_table_size(model, 'softmax')
-
-        h_file = open(f'{path}/{table_name}.tb', 'w')
-        h_file.write(self.__get_table_header(table_name, table_size))
-
-        # Default fixed point precision
-        # 6 bits for integer part, 10 bits for decimal - total, 16
-        fp_bits = 16
-        fp_integer = 6
-        fp_signed = True
-
-        # Exp table should use the same precision as exp_table, as seen in Vivado code
-        # init_exp_table<data_T, CONFIG_T>(exp_table);
+    def __get_table_precision(self, model, activation, table_name='table_precision'):
         for layer in model.get_layers():
-            if layer.name == 'softmax':
-                ac_type = layer.get_input_variable().type
-                if ac_type is not None:
-                    try:
-                        fp_bits = ac_type.precision.integer + ac_type.precision.fractional
-                        fp_integer = ac_type.precision.integer
-                        fp_signed = ac_type.precision.signed
-                    except Exception:
-                        # FixedPrecisionType wasn't correctly stored in layer attributes, use default values
-                        pass
-                    if fp_signed is False:
-                        raise Exception('Softmax types need to be signed')
-
-        sep = ''
-        N = ceil_log2(table_size)
-        for i in range(table_size):
-            f = FixedPointEmulator(fp_bits, fp_integer, signed=fp_signed)
-            b = uint_to_binary(i, N)
-            if i == 0:
-                b.insert(0, 0)
-            else:
-                b.insert(0, 1)
-            f.set_msb_bits(b)
-            real_val = f.exp_float()
-            h_file.write(sep + str(real_val))
-            sep = ', '
+            if layer.get_attr('activation') == activation and layer.get_attr(table_name) is not None:
+                precision = layer.get_attr(table_name)
+                return precision.precision
 
-        h_file.write('};\n')
-        h_file.close()
+        return None  # fp_bits, fp_integer, fp_signed
 
-    def __write_invert_table(self, model, path):
-        table_name = 'invert_table'
-        table_size = self.__get_table_size(model, 'softmax')
+    def __write_exp_tables_stable(self, model, path):
 
-        h_file = open(f'{path}/{table_name}.tb', 'w')
-        h_file.write(self.__get_table_header(table_name, table_size))
-
-        # Default fixed point precision, in case values from layer attributes cannot be extracted
-        # 8 bits for integer part, 10 bits for decimal - total, 18
-        fp_bits = 18
-        fp_integer = 8
-        fp_signed = True
-
-        # Invert table should use the same precision as exp_table, as seen in Vivado code
-        # init_invert_table<typename CONFIG_T::exp_table_t, CONFIG_T>(invert_table);
         for layer in model.get_layers():
-            if layer.name == 'softmax':
-                ac_type = layer.get_attr('exp_table_t')
-                if ac_type is not None:
-                    try:
-                        fp_bits = ac_type.precision.integer + ac_type.precision.fractional
-                        fp_integer = ac_type.precision.integer
-                        fp_signed = ac_type.precision.signed
-                    except Exception:
-                        # FixedPrecisionType wasn't correctly stored in layer attributes, use default values
-                        pass
-                    if fp_signed is False:
-                        raise Exception('Softmax types need to be signed')
+            # Last property is essential since it seperates layer with activation property from actual activation layers
 
-        sep = ''
-        N = ceil_log2(table_size)
-        for i in range(table_size):
-            f = FixedPointEmulator(fp_bits, fp_integer, signed=fp_signed)
-            b = uint_to_binary(i, N)
-            b.insert(0, 0)
-            f.set_msb_bits(b)
-            real_val = f.inv_float()
-            h_file.write(sep + str(real_val))
-            sep = ', '
-
-        h_file.write('};\n')
-        h_file.close()
+            if (
+                (
+                    layer.get_attr('activation') == 'softmax'
+                    or layer.get_attr('activation') == 'softmax_multidim'
+                    or layer.get_attr('recurrent_activation') == 'softmax'
+                    or layer.get_attr('recurrent_activation') == 'softmax_multidim'
+                )
+                and 'implementation' in layer.attributes
+                and layer.get_attr('implementation') == 'stable'
+            ):
+                table_name = layer.name + '_exp_table'
+                table_size = min(int(layer.get_attr('table_size')), int(layer.get_attr('exp_table_size')))
+
+                with open(f'{path}/{table_name}.h', 'w') as h_file:
+                    header_name = table_name
+                    h_file.write(f'#ifndef {header_name.upper()}_H_\n')
+                    h_file.write(f'#define {header_name.upper()}_H_\n\n')
+
+                    h_file.write(
+                        f'static constexpr nnet::array<{layer.get_attr("exp_table_t").name},{table_size}> {table_name} = {{'
+                    )
 
-    def __write_exp_table_latency(self, model, path):
-        table_name = 'exp_table_latency'
-        table_size = self.__get_table_size(model, 'softmax')
+                    ac_type = layer.get_attr('inp_norm_t')
+                    fp_bits = ac_type.precision.integer + ac_type.precision.fractional
+                    fp_integer = ac_type.precision.integer
 
-        h_file = open(f'{path}/{table_name}.tb', 'w')
-        h_file.write(self.__get_table_header(table_name, table_size))
-
-        # Default fixed point precision
-        # 6 bits for integer part, 10 bits for decimal - total, 16
-        fp_bits = 16
-        fp_integer = 6
-        fp_signed = True
+                    # Copy scaling from attributes
+                    scale = (
+                        layer.attributes['exp_scale']
+                        if (('exp_scale' in layer.attributes) and (layer.attributes['exp_scale'] is not None))
+                        else 1.0
+                    )
 
-        # Exp table should use the same precision as exp_table, as seen in Vivado code
-        # init_exp_table<data_T, CONFIG_T>(exp_table);
+                    N = ceil_log2(table_size)
+                    if N > 2**fp_bits:
+                        raise Exception('Table size is bigger than what precision allows')
+                    maxval = 2**fp_integer - 1
+
+                    sep = ''
+                    # Use the top bits if table_size < 2**bit_width
+                    for i in range(table_size):
+                        # Norm type is always > 1 so if input quantiser is set to be signed for any reason,
+                        # force unsigned but keep the width
+                        f = FixedPointEmulator(fp_bits, fp_integer, signed=False)
+                        b = uint_to_binary(i, N)
+                        f.set_msb_bits(b)
+                        real_val = (1.0 / f.exp_float()) * scale
+                        if real_val > maxval:
+                            real_val = maxval
+                        h_file.write(sep + str(real_val))
+                        sep = ', '
+
+                    h_file.write('};\n\n')
+                    h_file.write('#endif')
+
+    def __write_invert_tables_stable(self, model, path):
         for layer in model.get_layers():
-            if layer.name == 'softmax':
-                ac_type = layer.get_input_variable().type
-                if ac_type is not None:
-                    try:
-                        fp_bits = ac_type.precision.integer + ac_type.precision.fractional
-                        fp_integer = ac_type.precision.integer
-                        fp_signed = ac_type.precision.signed
-                    except Exception:
-                        # FixedPrecisionType wasn't correctly stored in layer attributes, use default values
-                        pass
+            # Last property is essential since it seperates layer with activation property from actual activation layers
 
-        sep = ''
-        N = ceil_log2(table_size)
-        for i in range(table_size):
-            f = FixedPointEmulator(fp_bits, fp_integer, signed=fp_signed)
-            f.set_msb_bits(uint_to_binary(i, N))
-            real_val = f.exp_float()
-            h_file.write(sep + str(real_val))
-            sep = ', '
+            if (
+                (
+                    layer.get_attr('activation') == 'softmax'
+                    or layer.get_attr('activation') == 'softmax_multidim'
+                    or layer.get_attr('recurrent_activation') == 'softmax'
+                    or layer.get_attr('recurrent_activation') == 'softmax_multidim'
+                )
+                and 'implementation' in layer.attributes
+                and layer.get_attr('implementation') == 'stable'
+            ):
+                table_name = layer.name + '_inv_table'
+                table_size = min(int(layer.get_attr('table_size')), int(layer.get_attr('inv_table_size')))
+
+                with open(f'{path}/{table_name}.h', 'w') as h_file:
+                    header_name = table_name
+                    h_file.write(f'#ifndef {header_name.upper()}_H_\n')
+                    h_file.write(f'#define {header_name.upper()}_H_\n\n')
+
+                    h_file.write(
+                        f'static constexpr nnet::array<{layer.get_attr("inv_table_t").name},{table_size}> {table_name} = {{'
+                    )
 
-        h_file.write('};\n')
-        h_file.close()
+                    ac_type = layer.get_attr('inv_inp_t')
+                    fp_bits = ac_type.precision.integer + ac_type.precision.fractional
+                    fp_integer = ac_type.precision.integer
+
+                    N = ceil_log2(table_size)
+                    if N > 2**fp_bits:
+                        raise Exception('Table size is bigger than what precision allows')
+                    maxval = 2**fp_integer - 1
+
+                    # Use the top bits if table_size < 2**bit_width
+                    sep = ''
+                    for i in range(table_size):
+                        # Norm type is always > 1 so if input quantiser is set to be signed for any reason,
+                        # force unsigned but keep the width
+                        f = FixedPointEmulator(fp_bits, fp_integer, signed=False)
+                        b = uint_to_binary(i, N)
+                        f.set_msb_bits(b)
+                        real_val = f.inv_float()
+                        if real_val > maxval:
+                            real_val = maxval
+                        h_file.write(sep + str(real_val))
+                        sep = ', '
+
+                    h_file.write('};\n\n')
+                    h_file.write('#endif')
+
+    def __write_exp_tables_latency(self, model, path):
+        for layer in model.get_layers():
+            # Last property is essential since it seperates layer with activation property from actual activation layers
+            if (
+                (
+                    layer.get_attr('activation') == 'softmax'
+                    or layer.get_attr('activation') == 'softmax_multidim'
+                    or layer.get_attr('recurrent_activation') == 'softmax'
+                    or layer.get_attr('recurrent_activation') == 'softmax_multidim'
+                )
+                and 'implementation' in layer.attributes
+                and layer.get_attr('implementation') == 'latency'
+            ):
+                table_name = layer.name + '_exp_table'
+                table_size = int(layer.get_attr('exp_table_size'))
+
+                with open(f'{path}/{table_name}.h', 'w') as h_file:
+                    header_name = table_name
+                    h_file.write(f'#ifndef {header_name.upper()}_H_\n')
+                    h_file.write(f'#define {header_name.upper()}_H_\n\n')
+
+                    h_file.write(
+                        f'static constexpr nnet::array<{layer.get_attr("exp_table_t").name},{table_size}> {table_name} = {{'
+                    )
 
-    def __write_invert_table_latency(self, model, path):
-        table_name = 'invert_table_latency'
-        table_size = self.__get_table_size(model, 'softmax')
+                    ac_type = layer.get_input_variable().type
+                    fp_bits = ac_type.precision.integer + ac_type.precision.fractional
+                    fp_integer = ac_type.precision.integer
+                    fp_signed = ac_type.precision.signed
 
-        h_file = open(f'{path}/{table_name}.tb', 'w')
-        h_file.write(self.__get_table_header(table_name, table_size))
+                    sep = ''
+                    N = ceil_log2(table_size)
+                    for i in range(table_size):
+                        f = FixedPointEmulator(fp_bits, fp_integer, signed=fp_signed)
+                        f.set_msb_bits(uint_to_binary(i, N))
+                        real_val = f.exp_float()
+                        h_file.write(sep + str(real_val))
+                        sep = ', '
 
-        # Default fixed point precision, in case values from layer attributes cannot be extracted
-        # 8 bits for integer part, 10 bits for decimal - total, 18
-        fp_bits = 18
-        fp_integer = 8
-        fp_signed = True
+                    h_file.write('};\n\n')
+                    h_file.write('#endif')
 
-        # Invert table should use the same precision as exp_table, as seen in Vivado code
-        # init_invert_table<typename CONFIG_T::exp_table_t, CONFIG_T>(invert_table);
+    def __write_invert_tables_latency(self, model, path):
         for layer in model.get_layers():
-            if layer.name == 'softmax':
-                ac_type = layer.get_attr('exp_table_t')
-                if ac_type is not None:
-                    try:
-                        fp_bits = ac_type.precision.integer + ac_type.precision.fractional
-                        fp_integer = ac_type.precision.integer
-                        fp_signed = ac_type.precision.signed
-                    except Exception:
-                        # FixedPrecisionType wasn't correctly stored in layer attributes, use default values
-                        pass
+            # Last property is essential since it seperates layer with activation property from actual activation layers
+            if (
+                (
+                    layer.get_attr('activation') == 'softmax'
+                    or layer.get_attr('activation') == 'softmax_multidim'
+                    or layer.get_attr('recurrent_activation') == 'softmax'
+                    or layer.get_attr('recurrent_activation') == 'softmax_multidim'
+                )
+                and 'implementation' in layer.attributes
+                and layer.get_attr('implementation') == 'latency'
+            ):
+                table_name = layer.name + '_inv_table'
+                table_size = int(layer.get_attr('inv_table_size'))
+
+                with open(f'{path}/{table_name}.h', 'w') as h_file:
+                    header_name = table_name
+                    h_file.write(f'#ifndef {header_name.upper()}_H_\n')
+                    h_file.write(f'#define {header_name.upper()}_H_\n\n')
+
+                    h_file.write(
+                        f'static constexpr nnet::array<{layer.get_attr("inv_table_t").name},{table_size}> {table_name} = {{'
+                    )
 
-        sep = ''
-        N = ceil_log2(table_size)
-        for i in range(table_size):
-            f = FixedPointEmulator(fp_bits, fp_integer, signed=fp_signed)
-            f.set_msb_bits(uint_to_binary(i, N))
-            real_val = f.inv_float()
-            h_file.write(sep + str(real_val))
-            sep = ', '
+                    ac_type = layer.get_attr('exp_table_t')
+                    fp_bits = ac_type.precision.integer + ac_type.precision.fractional
+                    fp_integer = ac_type.precision.integer
+                    fp_signed = ac_type.precision.signed
 
-        h_file.write('};\n')
-        h_file.close()
+                    sep = ''
+                    N = ceil_log2(table_size)
+                    for i in range(table_size):
+                        f = FixedPointEmulator(fp_bits, fp_integer, signed=fp_signed)
+                        f.set_msb_bits(uint_to_binary(i, N))
+                        real_val = f.inv_float()
+                        h_file.write(sep + str(real_val))
+                        sep = ', '
+
+                    h_file.write('};\n\n')
+                    h_file.write('#endif')
 
     def __write_exp_table_legacy(self, model, path):
-        table_name = 'exp_table_legacy'
-        table_size = self.__get_table_size(model, 'softmax')
 
-        h_file = open(f'{path}/{table_name}.tb', 'w')
-        h_file.write(self.__get_table_header(table_name, table_size))
+        for layer in model.get_layers():
+            # Last property is essential since it seperates layer with activation property from actual activation layers
+            if (
+                (
+                    layer.get_attr('activation') == 'softmax'
+                    or layer.get_attr('activation') == 'softmax_multidim'
+                    or layer.get_attr('recurrent_activation') == 'softmax'
+                    or layer.get_attr('recurrent_activation') == 'softmax_multidim'
+                )
+                and 'implementation' in layer.attributes
+                and layer.get_attr('implementation') == 'legacy'
+            ):
+                table_name = layer.name + '_exp_table'
+                table_size = int(layer.get_attr('exp_table_size'))  # not sure if it works, have to test first
+
+                with open(f'{path}/{table_name}.h', 'w') as h_file:
+                    header_name = table_name
+                    h_file.write(f'#ifndef {header_name.upper()}_H_\n')
+                    h_file.write(f'#define {header_name.upper()}_H_\n\n')
+
+                    h_file.write(
+                        f'static constexpr nnet::array<{layer.get_attr("exp_table_t").name},{table_size}> {table_name} = {{'
+                    )
 
-        sep = ''
-        for i in range(table_size):
-            in_val = 2 * 8.0 * (i - float(table_size) / 2.0) / float(table_size)
-            real_val = np.exp(in_val)
-            h_file.write(sep + str(real_val))
-            sep = ', '
+                    sep = ''
+                    for i in range(table_size):
+                        in_val = 2 * 8.0 * (i - float(table_size) / 2.0) / float(table_size)
+                        real_val = np.exp(in_val)
+                        h_file.write(sep + str(real_val))
+                        sep = ', '
 
-        h_file.write('};\n')
-        h_file.close()
+                    h_file.write('};\n\n')
+                    h_file.write('#endif')
 
     def __write_invert_table_legacy(self, model, path):
-        table_name = 'invert_table_legacy'
-        table_size = self.__get_table_size(model, 'softmax')
 
-        h_file = open(f'{path}/{table_name}.tb', 'w')
-        h_file.write(self.__get_table_header(table_name, table_size))
+        for layer in model.get_layers():
+            # Last property is essential since it seperates layer with activation property from actual activation layers
+            if (
+                (
+                    layer.get_attr('activation') == 'softmax'
+                    or layer.get_attr('activation') == 'softmax_multidim'
+                    or layer.get_attr('recurrent_activation') == 'softmax'
+                    or layer.get_attr('recurrent_activation') == 'softmax_multidim'
+                )
+                and 'implementation' in layer.attributes
+                and layer.get_attr('implementation') == 'legacy'
+            ):
+                table_name = layer.name + '_inv_table'
+                table_size = int(layer.get_attr('inv_table_size'))
+
+                with open(f'{path}/{table_name}.h', 'w') as h_file:
+                    header_name = table_name
+                    h_file.write(f'#ifndef {header_name.upper()}_H_\n')
+                    h_file.write(f'#define {header_name.upper()}_H_\n\n')
+
+                    h_file.write(
+                        f'static constexpr nnet::array<{layer.get_attr("inv_table_t").name},{table_size}> {table_name} = {{'
+                    )
 
-        sep = ''
-        for i in range(table_size):
-            real_val = 0
-            in_val = 64.0 * i / float(table_size)
-            if in_val > 0.0:
-                real_val = 1.0 / in_val
-            h_file.write(sep + str(real_val))
-            sep = ', '
+                    sep = ''
+                    for i in range(table_size):
+                        real_val = 0
+                        in_val = 64.0 * i / float(table_size)
+                        if in_val > 0.0:
+                            real_val = 1.0 / in_val
+                        h_file.write(sep + str(real_val))
+                        sep = ', '
 
-        h_file.write('};\n')
-        h_file.close()
+                    h_file.write('};\n\n')
+                    h_file.write('#endif')
 
     def write_activation_tables(self, model):
         """Write the lookup tables for activation functions
@@ -909,10 +997,10 @@ def write_activation_tables(self, model):
         self.__write_softplus_table(model, dstpath)
         self.__write_softsign_table(model, dstpath)
         self.__write_selu_table(model, dstpath)
-        self.__write_exp_table(model, dstpath)
-        self.__write_invert_table(model, dstpath)
-        self.__write_exp_table_latency(model, dstpath)
-        self.__write_invert_table_latency(model, dstpath)
+        self.__write_exp_tables_stable(model, dstpath)
+        self.__write_invert_tables_stable(model, dstpath)
+        self.__write_exp_tables_latency(model, dstpath)
+        self.__write_invert_tables_latency(model, dstpath)
         self.__write_exp_table_legacy(model, dstpath)
         self.__write_invert_table_legacy(model, dstpath)
 
diff --git a/hls4ml/writer/quartus_writer.py b/hls4ml/writer/quartus_writer.py
index e0d6338ac3..f5b56d9f7c 100644
--- a/hls4ml/writer/quartus_writer.py
+++ b/hls4ml/writer/quartus_writer.py
@@ -1097,8 +1097,6 @@ def __write_exp_table(self, model, path):
                     except Exception:
                         # FixedPrecisionType wasn't correctly stored in layer attributes, use default values
                         pass
-                    if fp_signed is False:
-                        raise Exception('Softmax types need to be signed')
 
         sep = ''
         N = ceil_log2(table_size)
@@ -1143,8 +1141,6 @@ def __write_invert_table(self, model, path):
                     except Exception:
                         # FixedPrecisionType wasn't correctly stored in layer attributes, use default values
                         pass
-                    if fp_signed is False:
-                        raise Exception('Softmax types need to be signed')
 
         sep = ''
         N = ceil_log2(table_size)
diff --git a/test/pytest/test_auto_precision.py b/test/pytest/test_auto_precision.py
index d3738c8461..57ea268340 100644
--- a/test/pytest/test_auto_precision.py
+++ b/test/pytest/test_auto_precision.py
@@ -1,3 +1,4 @@
+import math
 from pathlib import Path
 
 import numpy as np
@@ -13,10 +14,12 @@
     ReLU,
     SeparableConv1D,
     SeparableConv2D,
+    Softmax,
 )
 from tensorflow.keras.models import Sequential
 
 import hls4ml
+import hls4ml.model.layers
 from hls4ml.model.optimizer.passes.infer_precision import _get_precision_from_constant
 
 test_root_path = Path(__file__).parent
@@ -285,3 +288,37 @@ def test_precision_from_constant_unit(val, expected_width):
     quantum = 2.0**-fp.fractional
     if expected_width < max_width:
         assert val % quantum == 0
+
+
+@pytest.mark.parametrize('n_in', [4, 8, 16])
+@pytest.mark.parametrize('backend', ['Vitis', 'oneAPI'])
+def test_auto_precision_softmax(test_case_id, n_in, backend):
+    """Test that auto accumulator precision is correctly inferred for softmax layers."""
+    model = Sequential()
+    model.add(Softmax(input_shape=(n_in,)))
+    model.compile()
+
+    config = hls4ml.utils.config_from_keras_model(model, backend=backend, granularity='name')
+
+    odir = str(test_root_path / test_case_id)
+    hls_model = hls4ml.converters.convert_from_keras_model(model, hls_config=config, output_dir=odir, backend=backend)
+
+    # Find the Softmax layer and verify accum_t precision
+    softmax_layer = next((layer for layer in hls_model.get_layers() if isinstance(layer, hls4ml.model.layers.Softmax)), None)
+    assert softmax_layer is not None, 'No Softmax layer found in converted model'
+
+    accum_t = softmax_layer.types['accum_t'].precision
+    exp_table_t = softmax_layer.types['exp_table_t'].precision
+
+    ceillog = math.ceil(math.log2(n_in))
+    expected_width = exp_table_t.width + ceillog
+    expected_integer = exp_table_t.integer + ceillog
+    expected_signed = exp_table_t.signed
+
+    assert accum_t.width == expected_width, f'Expected accum_t width {expected_width}, got {accum_t.width} (n_in={n_in})'
+    assert accum_t.integer == expected_integer, (
+        f'Expected accum_t integer {expected_integer}, got {accum_t.integer} (n_in={n_in})'
+    )
+    assert accum_t.signed == expected_signed, (
+        f'Expected accum_t signed={expected_signed}, got {accum_t.signed} (n_in={n_in})'
+    )
diff --git a/test/pytest/test_softmax.py b/test/pytest/test_softmax.py
index 418f64b558..84925c445c 100644
--- a/test/pytest/test_softmax.py
+++ b/test/pytest/test_softmax.py
@@ -14,13 +14,13 @@
 def generate_data(input_shape):
     shape = (5000, *input_shape)
     d = np.random.normal(0, 2, shape)
-    modify_entries = np.random.randint(0, 1, shape) < 0.05
+    modify_entries = np.random.rand(*shape) < 0.05
     d[modify_entries] = d[modify_entries] * 5 + 10
     return np.clip(d, -32, 31)
 
 
-@pytest.mark.parametrize('backend', ['Vivado', 'Vitis', 'Quartus', 'Catapult'])
-@pytest.mark.parametrize('strategy', ['stable', 'latency', 'argmax'])
+@pytest.mark.parametrize('backend', ['Vivado', 'Vitis', 'Quartus', 'Catapult', 'oneAPI'])
+@pytest.mark.parametrize('implementation', ['stable', 'latency', 'argmax'])
 @pytest.mark.parametrize(
     'input_bits,input_shape,table_bits,io_type,custom_accum',
     [
@@ -35,26 +35,31 @@ def generate_data(input_shape):
         ('16,6', (8, 8, 3), '18,8', 'io_stream', False),
     ],
 )
-def test_softmax(test_case_id, backend, strategy, generate_data, input_bits, input_shape, table_bits, io_type, custom_accum):
+def test_softmax(
+    test_case_id, backend, implementation, generate_data, input_bits, input_shape, table_bits, io_type, custom_accum
+):
+
+    if backend == 'Catapult' and implementation == 'argmax':
+        pytest.skip('Argmax is not supported in cataplut')
+
     X = generate_data
     model = tf.keras.models.Sequential()
     model.add(tf.keras.layers.Activation(input_shape=input_shape, activation='softmax', name='softmax'))
     model.compile()
 
-    table_type = f'fixed<{table_bits}, RND, SAT>'
+    table_type = f'ufixed<{table_bits}, RND, SAT>' if backend != 'oneAPI' else f'ac_fixed<{table_bits}, AC_RND, AC_SAT>'
 
     cfg = hls4ml.utils.config_from_keras_model(model, granularity='name', backend=backend)
-    cfg['LayerName']['softmax']['Strategy'] = strategy
-    cfg['LayerName']['softmax']['inv_table_t'] = table_type
-    cfg['LayerName']['softmax']['exp_table_t'] = table_type
-    cfg['LayerName']['softmax']['accum_t'] = table_type
-    cfg['LayerName']['softmax']['inv_inp_t'] = table_type
+    cfg['LayerName']['softmax']['Implementation'] = implementation
+    cfg['LayerName']['softmax']['Precision']['inv_table'] = table_type
+    cfg['LayerName']['softmax']['Precision']['exp_table'] = table_type
+    cfg['LayerName']['softmax']['Precision']['inv_inp'] = table_type
     if custom_accum:
         if backend not in ['Vivado', 'Vitis']:
             pytest.skip('Custom accumulators are only supported for Vivado and Vitis backends')
         W, I = map(int, input_bits.split(','))  # noqa: E741
-        cfg['LayerName']['softmax']['accum_t'] = f'fixed<{W + 3},{I + 3}>'
-        cfg['LayerName']['softmax']['inv_inp_t'] = f'fixed<{W + 2},{I + 2}>'
+        cfg['LayerName']['softmax']['Precision']['inv_inp'] = f'ufixed<{W + 2},{I + 2}>'
+        
     inp_layer_name = next(iter(cfg['LayerName'].keys()))
     cfg['LayerName'][inp_layer_name]['Precision']['result'] = f'fixed<{input_bits}>'
 
@@ -73,7 +78,7 @@ def test_softmax(test_case_id, backend, strategy, generate_data, input_bits, inp
     assert acc_hls4ml >= 0.98
 
 
-@pytest.mark.parametrize('backend', ['Vivado', 'Vitis', 'Quartus', 'Catapult'])
+@pytest.mark.parametrize('backend', ['Vivado', 'Vitis', 'Quartus', 'Catapult', 'oneAPI'])
 @pytest.mark.parametrize('io_type', ['io_parallel', 'io_stream'])
 def test_softmax_skipped(test_case_id, backend, io_type):
     X = np.random.rand(100, 10)