diff --git a/numpy_financial/_financial.py b/numpy_financial/_financial.py
index 534230f..25865d5 100644
--- a/numpy_financial/_financial.py
+++ b/numpy_financial/_financial.py
@@ -12,6 +12,7 @@
 """
 
 from decimal import Decimal
+from typing import Optional
 
 import numba as nb
 import numpy as np
@@ -60,6 +61,35 @@ def _convert_when(when):
         return [_when_to_num[x] for x in when]
 
 
+def _validate_arrays(
+    *,
+    rate: Optional[np.ndarray],
+    pmt: Optional[np.ndarray],
+    pv: Optional[np.ndarray],
+    fv: Optional[np.ndarray],
+    when: Optional[np.ndarray],
+):
+    if rate is not None and rate.ndim != 1:
+        msg = "invalid shape for rates. Rate must be either a scalar or 1d array"
+        raise ValueError(msg)
+
+    if pmt is not None and pmt.ndim != 1:
+        msg = "invalid shape for pmt. Payments must be either a scalar or 1d array"
+        raise ValueError(msg)
+
+    if pv is not None and pv.ndim != 1:
+        msg = "invalid shape for pv. Present value must be either a scalar or 1d array"
+        raise ValueError(msg)
+
+    if fv is not None and fv.ndim != 1:
+        msg = "invalid shape for fv. Future value must be either a scalar or 1d array"
+        raise ValueError(msg)
+
+    if when is not None and when.ndim != 1:
+        msg = "invalid shape for when. When must be either a scalar or 1d array"
+        raise ValueError(msg)
+
+
 def fv(rate, nper, pmt, pv, when='end'):
     """Compute the future value.
 
@@ -261,6 +291,40 @@ def pmt(rate, nper, pv, fv=0, when='end'):
     return -(fv + pv * temp) / fact
 
 
+@nb.njit
+def _nper_inner_loop(rate, pmt, pv, fv, when):
+    if rate == 0.0:
+        if pmt == 0.0:
+            # If no repayments are made the payments will go on forever
+            return np.inf
+        else:
+            return -(fv + pv) / pmt
+    else:
+        # We know that rate != 0.0, so we are sure this won't cause a ZeroDivisionError
+        z = pmt * (1.0 + rate * when) / rate
+        try:
+            numer = np.log((-fv + z) / (pv + z))
+            denom = np.log(1.0 + rate)
+            return numer / denom
+        except Exception:  # As of March 24, numba only supports generic exceptions
+            # TODO: There are several ``ZeroDivisionError``s here.
+            #       We need to figure out exactly what's causing these
+            #       and return financially sensible values.
+            return np.nan
+
+
+@nb.njit
+def _nper_native(rates, pmts, pvs, fvs, whens, out):
+    for rate in range(rates.shape[0]):
+        for pmt in range(pmts.shape[0]):
+            for pv in range(pvs.shape[0]):
+                for fv in range(fvs.shape[0]):
+                    for when in range(whens.shape[0]):
+                        out[rate, pmt, pv, fv, when] = _nper_inner_loop(
+                            rates[rate], pmts[pmt], pvs[pv], fvs[fv], whens[when]
+                        )
+
+
 def nper(rate, pmt, pv, fv=0, when='end'):
     """Compute the number of periodic payments.
 
@@ -297,7 +361,7 @@ def nper(rate, pmt, pv, fv=0, when='end'):
     If you only had $150/month to pay towards the loan, how long would it take
     to pay-off a loan of $8,000 at 7% annual interest?
 
-    >>> print(np.round(npf.nper(0.07/12, -150, 8000), 5))
+    >>> round(npf.nper(0.07/12, -150, 8000), 5)
     64.07335
 
     So, over 64 months would be required to pay off the loan.
@@ -305,35 +369,46 @@ def nper(rate, pmt, pv, fv=0, when='end'):
     The same analysis could be done with several different interest rates
     and/or payments and/or total amounts to produce an entire table.
 
-    >>> npf.nper(*(np.ogrid[0.07/12: 0.08/12: 0.01/12,
-    ...                     -150   : -99    : 50    ,
-    ...                     8000   : 9001   : 1000]))
-    array([[[ 64.07334877,  74.06368256],
-            [108.07548412, 127.99022654]],
+    >>> rates = [0.05, 0.06, 0.07]
+    >>> payments = [100, 200, 300]
+    >>> amounts = [7_000, 8_000, 9_000]
+    >>> npf.nper(rates, payments, amounts).round(3)
+    array([[[-30.827, -32.987, -34.94 ],
+            [-20.734, -22.517, -24.158],
+            [-15.847, -17.366, -18.78 ]],
+    <BLANKLINE>
+           [[-28.294, -30.168, -31.857],
+            [-19.417, -21.002, -22.453],
+            [-15.025, -16.398, -17.67 ]],
     <BLANKLINE>
-           [[ 66.12443902,  76.87897353],
-            [114.70165583, 137.90124779]]])
+           [[-26.234, -27.891, -29.381],
+            [-18.303, -19.731, -21.034],
+            [-14.311, -15.566, -16.722]]])
+
     """
     when = _convert_when(when)
-    rate, pmt, pv, fv, when = np.broadcast_arrays(rate, pmt, pv, fv, when)
-    nper_array = np.empty_like(rate, dtype=np.float64)
 
-    zero = rate == 0
-    nonzero = ~zero
-
-    with np.errstate(divide='ignore'):
-        # Infinite numbers of payments are okay, so ignore the
-        # potential divide by zero.
-        nper_array[zero] = -(fv[zero] + pv[zero]) / pmt[zero]
+    rate_inner = np.atleast_1d(rate)
+    pmt_inner = np.atleast_1d(pmt)
+    pv_inner = np.atleast_1d(pv)
+    fv_inner = np.atleast_1d(fv)
+    when_inner = np.atleast_1d(when)
+
+    _validate_arrays(
+        rate=rate_inner,
+        pmt=pmt_inner,
+        pv=pv_inner,
+        fv=fv_inner,
+        when=when_inner,
+    )
 
-    nonzero_rate = rate[nonzero]
-    z = pmt[nonzero] * (1 + nonzero_rate * when[nonzero]) / nonzero_rate
-    nper_array[nonzero] = (
-            np.log((-fv[nonzero] + z) / (pv[nonzero] + z))
-            / np.log(1 + nonzero_rate)
+    out_shape = _get_output_array_shape(
+        rate_inner, pmt_inner, pv_inner, fv_inner, when_inner
     )
+    out = np.empty(out_shape)
+    _nper_native(rate_inner, pmt_inner, pv_inner, fv_inner, when_inner, out)
 
-    return nper_array
+    return _ufunc_like(out)
 
 
 def _value_like(arr, value):
diff --git a/tests/strategies.py b/tests/strategies.py
new file mode 100644
index 0000000..2b27985
--- /dev/null
+++ b/tests/strategies.py
@@ -0,0 +1,48 @@
+from typing import Literal
+
+from hypothesis import strategies as st
+from hypothesis.extra import numpy as npst
+
+# numba only supports 32-bit or 64-bit little endian values
+NUMBA_ALLOWED_SIZES: list[Literal[32, 64]] = [32, 64]
+NUMBA_ALLOWED_ENDIANNESS = "<"
+
+
+def float_dtype():
+    return npst.floating_dtypes(
+        sizes=NUMBA_ALLOWED_SIZES,
+        endianness=NUMBA_ALLOWED_ENDIANNESS,
+    )
+
+
+def int_dtype():
+    return npst.integer_dtypes(
+        sizes=NUMBA_ALLOWED_SIZES,
+        endianness=NUMBA_ALLOWED_ENDIANNESS,
+    )
+
+
+def uint_dtype():
+    return npst.unsigned_integer_dtypes(
+        sizes=NUMBA_ALLOWED_SIZES,
+        endianness=NUMBA_ALLOWED_ENDIANNESS,
+    )
+
+
+real_scalar_dtypes = st.one_of(float_dtype(), int_dtype(), uint_dtype())
+cashflow_array_strategy = npst.arrays(
+    dtype=real_scalar_dtypes,
+    shape=npst.array_shapes(min_dims=1, max_dims=2, min_side=0, max_side=25),
+)
+cashflow_list_strategy = cashflow_array_strategy.map(lambda x: x.tolist())
+cashflow_array_like_strategy = st.one_of(
+    cashflow_array_strategy,
+    cashflow_list_strategy,
+)
+short_scalar_array_strategy = npst.arrays(
+    dtype=real_scalar_dtypes,
+    shape=npst.array_shapes(min_dims=0, max_dims=1, min_side=0, max_side=5),
+)
+when_strategy = st.sampled_from(
+    ['end', 'begin', 'e', 'b', 0, 1, 'beginning', 'start', 'finish']
+)
diff --git a/tests/test_financial.py b/tests/test_financial.py
index 2d64f0b..ff172f2 100644
--- a/tests/test_financial.py
+++ b/tests/test_financial.py
@@ -1,12 +1,7 @@
 import math
 from decimal import Decimal
 
-import hypothesis.extra.numpy as npst
-import hypothesis.strategies as st
-
-# Don't use 'import numpy as np', to avoid accidentally testing
-# the versions in numpy instead of numpy_financial.
-import numpy
+import numpy as np
 import pytest
 from hypothesis import given, settings
 from numpy.testing import (
@@ -17,37 +12,10 @@
 )
 
 import numpy_financial as npf
-
-
-def float_dtype():
-    return npst.floating_dtypes(sizes=[32, 64], endianness="<")
-
-
-def int_dtype():
-    return npst.integer_dtypes(sizes=[32, 64], endianness="<")
-
-
-def uint_dtype():
-    return npst.unsigned_integer_dtypes(sizes=[32, 64], endianness="<")
-
-
-real_scalar_dtypes = st.one_of(float_dtype(), int_dtype(), uint_dtype())
-
-
-cashflow_array_strategy = npst.arrays(
-    dtype=real_scalar_dtypes,
-    shape=npst.array_shapes(min_dims=1, max_dims=2, min_side=0, max_side=25),
-)
-cashflow_list_strategy = cashflow_array_strategy.map(lambda x: x.tolist())
-
-cashflow_array_like_strategy = st.one_of(
-    cashflow_array_strategy,
-    cashflow_list_strategy,
-)
-
-short_scalar_array = npst.arrays(
-    dtype=real_scalar_dtypes,
-    shape=npst.array_shapes(min_dims=0, max_dims=1, min_side=0, max_side=5),
+from tests.strategies import (
+    cashflow_array_like_strategy,
+    short_scalar_array_strategy,
+    when_strategy,
 )
 
 
@@ -219,7 +187,7 @@ def test_rate_with_infeasible_solution(self, number_type, when):
             number_type(5000.0),
             when=when,
         )
-        is_nan = Decimal.is_nan if number_type == Decimal else numpy.isnan
+        is_nan = Decimal.is_nan if number_type == Decimal else np.isnan
         assert is_nan(result)
 
     def test_rate_decimal(self):
@@ -231,7 +199,7 @@ def test_gh48(self):
         Test the correct result is returned with only infeasible solutions
         converted to nan.
         """
-        des = [-0.39920185, -0.02305873, -0.41818459, 0.26513414, numpy.nan]
+        des = [-0.39920185, -0.02305873, -0.41818459, 0.26513414, np.nan]
         nper = 2
         pmt = 0
         pv = [-593.06, -4725.38, -662.05, -428.78, -13.65]
@@ -280,18 +248,18 @@ def test_npv(self):
             rtol=1e-2,
         )
 
-    @given(rates=short_scalar_array, values=cashflow_array_strategy)
+    @given(rates=short_scalar_array_strategy, values=cashflow_array_like_strategy)
     @settings(deadline=None)
     def test_fuzz(self, rates, values):
         npf.npv(rates, values)
 
-    @pytest.mark.parametrize("rates", ([[1, 2, 3]], numpy.empty(shape=(1, 1, 1))))
+    @pytest.mark.parametrize("rates", ([[1, 2, 3]], np.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))))
+    @pytest.mark.parametrize("cf", ([[[1, 2, 3]]], np.empty(shape=(1, 1, 1))))
     def test_invalid_cashflows_shape(self, cf):
         rates = [1, 2, 3]
         with pytest.raises(ValueError):
@@ -299,8 +267,8 @@ def test_invalid_cashflows_shape(self, cf):
 
     @pytest.mark.parametrize("rate", (-1, -1.0))
     def test_rate_of_negative_one_returns_nan(self, rate):
-        cashflow = numpy.arange(5)
-        assert numpy.isnan(npf.npv(rate, cashflow))
+        cashflow = np.arange(5)
+        assert np.isnan(npf.npv(rate, cashflow))
 
 
 class TestPmt:
@@ -319,7 +287,7 @@ def test_pmt_broadcast(self):
         # Test the case where we use broadcast and
         # the arguments passed in are arrays.
         res = npf.pmt([[0.0, 0.8], [0.3, 0.8]], [12, 3], [2000, 20000])
-        tgt = numpy.array([[-166.66667, -19311.258], [-626.90814, -19311.258]])
+        tgt = np.array([[-166.66667, -19311.258], [-626.90814, -19311.258]])
         assert_allclose(res, tgt)
 
     def test_pmt_decimal_simple(self):
@@ -341,7 +309,7 @@ def test_pmt_decimal_broadcast(self):
             [Decimal("12"), Decimal("3")],
             [Decimal("2000"), Decimal("20000")],
         )
-        tgt = numpy.array(
+        tgt = np.array(
             [
                 [
                     Decimal("-166.6666666666666666666666667"),
@@ -386,7 +354,7 @@ def test_mirr(self, values, finance_rate, reinvest_rate, expected):
             difference = 10**-decimal_part_len
             assert_allclose(result, expected, atol=difference)
         else:
-            assert_(numpy.isnan(result))
+            assert_(np.isnan(result))
 
     def test_mirr_no_real_solution_exception(self):
         # Test that if there is no solution because all the cashflows
@@ -407,16 +375,16 @@ def test_basic_values(self):
         )
 
     def test_gh_18(self):
-        with numpy.errstate(divide="raise"):
+        with np.errstate(divide="raise"):
             assert_allclose(
                 npf.nper(0.1, 0, -500, 1500),
                 11.52670461,  # Computed using Google Sheet's NPER
             )
 
     def test_infinite_payments(self):
-        with numpy.errstate(divide="raise"):
+        with np.errstate(divide="raise"):
             result = npf.nper(0, -0.0, 1000)
-        assert_(result == numpy.inf)
+        assert_(result == np.inf)
 
     def test_no_interest(self):
         assert_(npf.nper(0, -100, 1000) == 10)
@@ -426,6 +394,33 @@ def test_broadcast(self):
             npf.nper(0.075, -2000, 0, 100000.0, [0, 1]), [21.5449442, 20.76156441], 4
         )
 
+    @given(
+        rates=short_scalar_array_strategy,
+        payments=short_scalar_array_strategy,
+        present_values=short_scalar_array_strategy,
+        future_values=short_scalar_array_strategy,
+        whens=when_strategy,
+    )
+    @settings(deadline=None)  # ignore jit compilation of a function
+    def test_fuzz(self, rates, payments, present_values, future_values, whens):
+        npf.nper(rates, payments, present_values, future_values, whens)
+
+    @pytest.mark.parametrize(
+        "rate_,pmt_,pv_,fv_,when_",
+        [
+            (np.empty((5, 2)), np.empty(5), np.empty(5), np.empty(5), np.empty(5)),
+            (np.empty(5), np.empty((5, 2)), np.empty(5), np.empty(5), np.empty(5)),
+            (np.empty(5), np.empty(5), np.empty((5, 2)), np.empty(5), np.empty(5)),
+            (np.empty(5), np.empty(5), np.empty(5), np.empty((5, 2)), np.empty(5)),
+            (np.empty(5), np.empty(5), np.empty(5), np.empty(5), np.empty((5, 2))),
+        ]
+    )
+    def test_2d_array_for_1d_shapes(self, rate_, pmt_, pv_, fv_, when_):
+        # Array values do not matter, only that they have invalid dimensions
+        with pytest.raises(ValueError, match="invalid shape for"):
+            npf.nper(rate_, pmt_, pv_, fv_, when_)
+
+
 
 class TestPpmt:
     def test_float(self):
@@ -518,7 +513,7 @@ def test_invalid_per(self, args):
         ],
     )
     def test_broadcast(self, when, desired):
-        args = (0.1 / 12, numpy.arange(1, 5), 24, 2000, 0)
+        args = (0.1 / 12, np.arange(1, 5), 24, 2000, 0)
         result = npf.ppmt(*args) if when is None else npf.ppmt(*args, when)
         assert_allclose(result, desired, rtol=1e-5)
 
@@ -548,7 +543,7 @@ def test_broadcast(self, when, desired):
     def test_broadcast_decimal(self, when, desired):
         args = (
             Decimal("0.1") / Decimal("12"),
-            numpy.arange(1, 5),
+            np.arange(1, 5),
             Decimal("24"),
             Decimal("2000"),
             Decimal("0"),
@@ -610,7 +605,7 @@ def test_when_is_end_decimal(self, when):
     @pytest.mark.parametrize(
         "per, desired",
         [
-            (0, numpy.nan),
+            (0, np.nan),
             (1, 0),
             (2, -594.107158),
             (3, -592.971592),
@@ -620,15 +615,15 @@ def test_gh_17(self, per, desired):
         # All desired results computed using Google Sheet's IPMT
         rate = 0.001988079518355057
         result = npf.ipmt(rate, per, 360, 300000, when="begin")
-        if numpy.isnan(desired):
-            assert numpy.isnan(result)
+        if np.isnan(desired):
+            assert np.isnan(result)
         else:
             assert_allclose(result, desired, rtol=1e-6)
 
     def test_broadcasting(self):
-        desired = [numpy.nan, -16.66666667, -16.03647345, -15.40102862, -14.76028842]
+        desired = [np.nan, -16.66666667, -16.03647345, -15.40102862, -14.76028842]
         assert_allclose(
-            npf.ipmt(0.1 / 12, numpy.arange(5), 24, 2000),
+            npf.ipmt(0.1 / 12, np.arange(5), 24, 2000),
             desired,
             rtol=1e-6,
         )
@@ -651,10 +646,10 @@ def test_decimal_broadcasting(self):
     def test_0d_inputs(self):
         args = (0.1 / 12, 1, 24, 2000)
         # Scalar inputs should return a scalar.
-        assert numpy.isscalar(npf.ipmt(*args))
-        args = (numpy.array(args[0]),) + args[1:]
+        assert np.isscalar(npf.ipmt(*args))
+        args = (np.array(args[0]),) + args[1:]
         # 0d array inputs should return a scalar.
-        assert numpy.isscalar(npf.ipmt(*args))
+        assert np.isscalar(npf.ipmt(*args))
 
 
 class TestFv:
@@ -733,7 +728,7 @@ def test_npv_irr_congruence(self):
         # a series of cashflows to be zero, so we should have
         #
         # NPV(IRR(x), x) = 0.
-        cashflows = numpy.array([-40000, 5000, 8000, 12000, 30000])
+        cashflows = np.array([-40000, 5000, 8000, 12000, 30000])
         assert_allclose(
             npf.npv(npf.irr(cashflows), cashflows),
             0,
@@ -771,7 +766,7 @@ def test_trailing_zeros(self):
     )
     def test_numpy_gh_6744(self, v):
         # Test that if there is no solution then npf.irr returns nan.
-        assert numpy.isnan(npf.irr(v))
+        assert np.isnan(npf.irr(v))
 
     def test_gh_15(self):
         v = [
@@ -808,13 +803,13 @@ def test_gh_15(self):
             1.3133070599585015e-313,
         ]
         result = npf.irr(v)
-        assert numpy.isfinite(result)
+        assert np.isfinite(result)
         # Very rough approximation taken from the issue.
         desired = -0.9999999990596069
         assert_allclose(result, desired, rtol=1e-9)
 
     def test_gh_39(self):
-        cashflows = numpy.array(
+        cashflows = np.array(
             [
                 -217500.0,
                 -217500.0,
@@ -856,7 +851,7 @@ def test_irr_no_real_solution_exception(self):
         # Test that if there is no solution because all the cashflows
         # have the same sign, then npf.irr returns NoRealSolutionException
         # when raise_exceptions is set to True.
-        cashflows = numpy.array([40000, 5000, 8000, 12000, 30000])
+        cashflows = np.array([40000, 5000, 8000, 12000, 30000])
 
         with pytest.raises(npf.NoRealSolutionError):
             npf.irr(cashflows, raise_exceptions=True)