diff --git a/.github/workflows/pythonpackage.yml b/.github/workflows/pythonpackage.yml
index 5d32bea..ff67b07 100644
--- a/.github/workflows/pythonpackage.yml
+++ b/.github/workflows/pythonpackage.yml
@@ -8,7 +8,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.9", "3.10", "3.11"]
     steps:
     - uses: actions/checkout@v3
     - name: Set up Python ${{ matrix.python-version }}
diff --git a/benchmarks/benchmarks.py b/benchmarks/benchmarks.py
index c026937..23250cd 100644
--- a/benchmarks/benchmarks.py
+++ b/benchmarks/benchmarks.py
@@ -1,40 +1,55 @@
+from decimal import Decimal
+
 import numpy as np
 
 import numpy_financial as npf
 
 
-class Npv1DCashflow:
-
-    param_names = ["cashflow_length"]
-    params = [
-        (1, 10, 100, 1000),
-    ]
-
-    def __init__(self):
-        self.cashflows = None
+def _to_decimal_array_1d(array):
+    return np.array([Decimal(x) for x in array.tolist()])
 
-    def setup(self, cashflow_length):
-        rng = np.random.default_rng(0)
-        self.cashflows = rng.standard_normal(cashflow_length)
 
-    def time_1d_cashflow(self, cashflow_length):
-        npf.npv(0.08, self.cashflows)
+def _to_decimal_array_2d(array):
+    decimals = [Decimal(x) for row in array.tolist() for x in row]
+    return np.array(decimals).reshape(array.shape)
 
 
-class Npv2DCashflows:
+class Npv2D:
 
-    param_names = ["n_cashflows", "cashflow_lengths"]
+    param_names = ["n_cashflows", "cashflow_lengths", "rates_lengths"]
     params = [
-        (1, 10, 100, 1000),
-        (1, 10, 100, 1000),
+        (1, 10, 100),
+        (1, 10, 100),
+        (1, 10, 100),
     ]
 
     def __init__(self):
+        self.rates_decimal = None
+        self.rates = None
+        self.cashflows_decimal = None
         self.cashflows = None
 
-    def setup(self, n_cashflows, cashflow_lengths):
+    def setup(self, n_cashflows, cashflow_lengths, rates_lengths):
         rng = np.random.default_rng(0)
-        self.cashflows = rng.standard_normal((n_cashflows, cashflow_lengths))
+        cf_shape = (n_cashflows, cashflow_lengths)
+        self.cashflows = rng.standard_normal(cf_shape)
+        self.rates = rng.standard_normal(rates_lengths)
+        self.cashflows_decimal = _to_decimal_array_2d(self.cashflows)
+        self.rates_decimal = _to_decimal_array_1d(self.rates)
+
+    def time_broadcast(self, n_cashflows, cashflow_lengths, rates_lengths):
+        npf.npv(self.rates, self.cashflows)
+
+    def time_for_loop(self, n_cashflows, cashflow_lengths, rates_lengths):
+        for rate in self.rates:
+            for cashflow in self.cashflows:
+                npf.npv(rate, cashflow)
+
+    def time_broadcast_decimal(self, n_cashflows, cashflow_lengths, rates_lengths):
+        npf.npv(self.rates_decimal, self.cashflows_decimal)
+
+    def time_for_loop_decimal(self, n_cashflows, cashflow_lengths, rates_lengths):
+        for rate in self.rates_decimal:
+            for cashflow in self.cashflows_decimal:
+                npf.npv(rate, cashflow)
 
-    def time_2d_cashflow(self, n_cashflows, cashflow_lengths):
-        npf.npv(0.08, self.cashflows)
diff --git a/numpy_financial/_financial.py b/numpy_financial/_financial.py
index 50728d8..26e31fe 100644
--- a/numpy_financial/_financial.py
+++ b/numpy_financial/_financial.py
@@ -13,6 +13,7 @@
 
 from decimal import Decimal
 
+import numba as nb
 import numpy as np
 
 __all__ = ['fv', 'pmt', 'nper', 'ipmt', 'ppmt', 'pv', 'rate',
@@ -46,6 +47,36 @@ def _convert_when(when):
         return [_when_to_num[x] for x in when]
 
 
+def _return_ufunc_like(array):
+    try:
+        # If size of array is one, return scalar
+        return array.item()
+    except ValueError:
+        # Otherwise, return entire array
+        return array
+
+
+def _is_object_array(array):
+    return array.dtype == np.dtype("O")
+
+
+def _use_decimal_dtype(*arrays):
+    return any(_is_object_array(array) for array in arrays)
+
+
+def _to_decimal_array_1d(array):
+    return np.array([Decimal(x) for x in array.tolist()])
+
+
+def _to_decimal_array_2d(array):
+    decimals = [Decimal(x) for row in array.tolist() for x in row]
+    return np.array(decimals).reshape(array.shape)
+
+
+def _get_output_array_shape(*arrays):
+    return tuple(array.shape[0] for array in arrays)
+
+
 def fv(rate, nper, pmt, pv, when='end'):
     """Compute the future value.
 
@@ -825,14 +856,35 @@ def irr(values, *, guess=None, tol=1e-12, maxiter=100, raise_exceptions=False):
     return np.nan
 
 
+@nb.njit(parallel=True)
+def _npv_native(rates, values, out):
+    for i in nb.prange(rates.shape[0]):
+        for j in nb.prange(values.shape[0]):
+            acc = 0.0
+            for t in range(values.shape[1]):
+                acc += values[j, t] / ((1.0 + rates[i]) ** t)
+            out[i, j] = acc
+
+
+# We require ``forceobj=True`` here to support decimal.Decimal types
+@nb.jit(forceobj=True)
+def _npv_decimal(rates, values, out):
+    for i in range(rates.shape[0]):
+        for j in range(values.shape[0]):
+            acc = Decimal("0.0")
+            for t in range(values.shape[1]):
+                acc += values[j, t] / ((Decimal("1.0") + rates[i]) ** t)
+            out[i, j] = acc
+
+
 def npv(rate, values):
     r"""Return the NPV (Net Present Value) of a cash flow series.
 
     Parameters
     ----------
-    rate : scalar
+    rate : scalar or array_like shape(K, )
         The discount rate.
-    values : array_like, shape(M, )
+    values : array_like, shape(M, ) or shape(M, N)
         The values of the time series of cash flows.  The (fixed) time
         interval between cash flow "events" must be the same as that for
         which `rate` is given (i.e., if `rate` is per year, then precisely
@@ -843,9 +895,10 @@ def npv(rate, values):
 
     Returns
     -------
-    out : float
+    out : float or array shape(K, M)
         The NPV of the input cash flow series `values` at the discount
-        `rate`.
+        `rate`. `out` follows the ufunc convention of returning scalars
+        instead of single element arrays.
 
     Warnings
     --------
@@ -891,16 +944,58 @@ def npv(rate, values):
     >>> np.round(npf.npv(rate, cashflows) + initial_cashflow, 5)
     3065.22267
 
+    The NPV calculation may be applied to several ``rates`` and ``cashflows``
+    simulatneously. This produces an array of shape
+    ``(len(rates), len(cashflows))``.
+
+    >>> rates = [0.00, 0.05, 0.10]
+    >>> cashflows = [[-4_000, 500, 800], [-5_000, 600, 900]]
+    >>> npf.npv(rates, cashflows).round(2)
+    array([[-2700.  , -3500.  ],
+           [-2798.19, -3612.24],
+           [-2884.3 , -3710.74]])
+
+    The NPV calculation also supports `decimal.Decimal` types, for example
+    if using Decimal ``rates``:
+
+    >>> rates = [Decimal("0.00"), Decimal("0.05"), Decimal("0.10")]
+    >>> cashflows = [[-4_000, 500, 800], [-5_000, 600, 900]]
+    >>> npf.npv(rates, cashflows)
+    array([[Decimal('-2700.0'), Decimal('-3500.0')],
+           [Decimal('-2798.185941043083900226757370'),
+            Decimal('-3612.244897959183673469387756')],
+           [Decimal('-2884.297520661157024793388430'),
+            Decimal('-3710.743801652892561983471074')]], dtype=object)
+
+    This also works for Decimal cashflows.
+
     """
+    rates = np.atleast_1d(rate)
     values = np.atleast_2d(values)
-    timestep_array = np.arange(0, values.shape[1])
-    npv = (values / (1 + rate) ** timestep_array).sum(axis=1)
-    try:
-        # If size of array is one, return scalar
-        return npv.item()
-    except ValueError:
-        # Otherwise, return entire array
-        return npv
+
+    if rates.ndim != 1:
+        msg = "invalid shape for rates. Rate must be either a scalar or 1d array"
+        raise ValueError(msg)
+
+    if values.ndim != 2:
+        msg = "invalid shape for values. Values must be either a 1d or 2d array"
+        raise ValueError(msg)
+
+    dtype = Decimal if _use_decimal_dtype(rates, values) else np.float64
+
+    if dtype == Decimal:
+        rates = _to_decimal_array_1d(rates)
+        values = _to_decimal_array_2d(values)
+
+    shape = _get_output_array_shape(rates, values)
+    out = np.empty(shape=shape, dtype=dtype)
+
+    if dtype == Decimal:
+        _npv_decimal(rates, values, out)
+    else:
+        _npv_native(rates, values, out)
+
+    return _return_ufunc_like(out)
 
 
 def mirr(values, finance_rate, reinvest_rate, *, raise_exceptions=False):
diff --git a/pyproject.toml b/pyproject.toml
index 62b758e..8fe3076 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -24,7 +24,6 @@ classifiers = [
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
-    "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3 :: Only",
     "Topic :: Software Development",
     "Topic :: Office/Business :: Financial :: Accounting",
@@ -38,8 +37,9 @@ classifiers = [
 packages = [{include = "numpy_financial"}]
 
 [tool.poetry.dependencies]
-python = "^3.9"
+python = "^3.9,<3.12"
 numpy = "^1.23"
+numba = "^0.58.1"
 
 
 [tool.poetry.group.test.dependencies]
diff --git a/tests/test_financial.py b/tests/test_financial.py
index ad01952..2f8f63d 100644
--- a/tests/test_financial.py
+++ b/tests/test_financial.py
@@ -164,7 +164,7 @@ def test_rate_maximum_iterations_exception_array(self):
 class TestNpv:
     def test_npv(self):
         assert_almost_equal(
-            npf.npv(0.05, [-15000, 1500, 2500, 3500, 4500, 6000]),
+            npf.npv(0.05, [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0]),
             122.89, 2)
 
     def test_npv_decimal(self):
@@ -174,17 +174,50 @@ def test_npv_decimal(self):
 
     def test_npv_broadcast(self):
         cashflows = [
-            [-15000, 1500, 2500, 3500, 4500, 6000],
-            [-15000, 1500, 2500, 3500, 4500, 6000],
-            [-15000, 1500, 2500, 3500, 4500, 6000],
-            [-15000, 1500, 2500, 3500, 4500, 6000],
+            [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0],
+            [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0],
+            [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0],
+            [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0],
         ]
         expected_npvs = [
-            122.8948549, 122.8948549, 122.8948549, 122.8948549
+            [122.8948549, 122.8948549, 122.8948549, 122.8948549]
         ]
         actual_npvs = npf.npv(0.05, cashflows)
         assert_allclose(actual_npvs, expected_npvs)
 
+    @pytest.mark.parametrize("dtype", [Decimal, float])
+    def test_npv_broadcast_equals_for_loop(self, dtype):
+        cashflows_str = [
+            ["-15000.0", "1500.0", "2500.0", "3500.0", "4500.0", "6000.0"],
+            ["-25000.0", "1500.0", "2500.0", "3500.0", "4500.0", "6000.0"],
+            ["-35000.0", "1500.0", "2500.0", "3500.0", "4500.0", "6000.0"],
+            ["-45000.0", "1500.0", "2500.0", "3500.0", "4500.0", "6000.0"],
+        ]
+        rates_str = ["-0.05", "0.00", "0.05", "0.10", "0.15"]
+
+        cashflows = numpy.array([[dtype(x) for x in cf] for cf in cashflows_str])
+        rates = numpy.array([dtype(x) for x in rates_str])
+
+        expected = numpy.empty((len(rates), len(cashflows)), dtype=dtype)
+        for i, r in enumerate(rates):
+            for j, cf in enumerate(cashflows):
+                expected[i, j] = npf.npv(r, cf)
+
+        actual = npf.npv(rates, cashflows)
+        assert_equal(actual, expected)
+
+    @pytest.mark.parametrize("rates", ([[1, 2, 3]], numpy.empty(shape=(1,1,1))))
+    def test_invalid_rates_shape(self, rates):
+        cashflows = [1, 2, 3]
+        with pytest.raises(ValueError):
+            npf.npv(rates, cashflows)
+
+    @pytest.mark.parametrize("cf", ([[[1, 2, 3]]], numpy.empty(shape=(1, 1, 1))))
+    def test_invalid_cashflows_shape(self, cf):
+        rates = [1, 2, 3]
+        with pytest.raises(ValueError):
+            npf.npv(rates, cf)
+
 
 class TestPmt:
     def test_pmt_simple(self):