def _apply_padding(lhs_arr, rhs_arr, offset, pad_mode, direction):
    """Apply padding to ``lhs_arr`` according to ``pad_mode``.

    This helper assigns the values in the excess parts (if existent)
    of ``lhs_arr`` according to the provided padding mode.

    This applies to the following values for ``pad_mode``:
    ``periodic``, ``symmetric``, ``order0``, ``order1``

    See `the online documentation
    <https://odlgroup.github.io/odl/math/resizing_ops.html>`_
    on resizing operators for details.
    """
    if pad_mode not in ('periodic', 'symmetric', 'order0', 'order1'):
        return

    full_slc = [slice(None)] * lhs_arr.ndim
    intersec_slc, _ = _intersection_slice_tuples(lhs_arr, rhs_arr, offset)

    if direction == 'forward':
        working_slc = list(intersec_slc)
    else:
        working_slc = list(full_slc)

    # TODO: order axes according to padding size for optimization (largest
    # last)? Axis strides could be important, too.
    for axis, (n_lhs, n_rhs) in enumerate(zip(lhs_arr.shape, rhs_arr.shape)):

        if n_lhs <= n_rhs:
            continue  # restriction, nothing to do

        n_pad_l = offset[axis]
        n_pad_r = n_lhs - n_rhs - n_pad_l

        # Error scenarios with illegal lengths
        if pad_mode == 'order0' and n_rhs == 0:
            raise ValueError('in axis {}: the smaller array must have size '
                             '>= 1 for order 0 padding, got 0'
                             ''.format(axis))

        if pad_mode == 'order1' and n_rhs < 2:
            raise ValueError('in axis {}: the smaller array must have size '
                             '>= 2 for order 1 padding, got {}'
                             ''.format(axis, n_rhs))

        for lr, pad_len in [('left', n_pad_l), ('right', n_pad_r)]:
            if pad_mode == 'periodic' and pad_len > n_rhs:
                raise ValueError('in axis {}: {} padding length {} exceeds '
                                 'the size {} of the smaller array; this is '
                                 'not allowed for periodic padding'
                                 ''.format(axis, lr, pad_len, n_rhs))

            elif pad_mode == 'symmetric' and pad_len >= n_rhs:
                raise ValueError('in axis {}: {} padding length {} is larger '
                                 'or equal to the size {} of the smaller '
                                 'array; this is not allowed for symmetric '
                                 'padding'
                                 ''.format(axis, lr, pad_len, n_rhs))

        # Slice tuples used to index LHS and RHS for left and right padding,
        # respectively; we make 4 copies of `working_slc` as lists
        lhs_slc_l, lhs_slc_r, rhs_slc_l, rhs_slc_r = map(
            list, [working_slc] * 4)

        # We're always using the outer (excess) parts involved in padding
        # on the LHS of the assignment, so we set them here.
        pad_slc_outer_l, pad_slc_outer_r = _padding_slices_outer(
            lhs_arr, rhs_arr, axis, offset)

        if direction == 'forward':
            lhs_slc_l[axis] = pad_slc_outer_l
            lhs_slc_r[axis] = pad_slc_outer_r
        else:
            rhs_slc_l[axis] = pad_slc_outer_l
            rhs_slc_r[axis] = pad_slc_outer_r

        if pad_mode in ('periodic', 'symmetric'):
            pad_slc_inner_l, pad_slc_inner_r = _padding_slices_inner(
                lhs_arr, rhs_arr, axis, offset, pad_mode)

            # Using `lhs_arr` on both sides of the assignment such that the
            # shapes match and the "corner" blocks are properly assigned
            # or used in the addition for the adjoint, respectively.
            if direction == 'forward':
                rhs_slc_l[axis] = pad_slc_inner_l
                rhs_slc_r[axis] = pad_slc_inner_r
                lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r = map(
                    tuple, [lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r])

                lhs_arr[lhs_slc_l] = lhs_arr[rhs_slc_l]
                lhs_arr[lhs_slc_r] = lhs_arr[rhs_slc_r]
            else:
                lhs_slc_l[axis] = pad_slc_inner_l
                lhs_slc_r[axis] = pad_slc_inner_r
                lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r = map(
                    tuple, [lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r])

                lhs_arr[lhs_slc_l] += lhs_arr[rhs_slc_l]
                lhs_arr[lhs_slc_r] += lhs_arr[rhs_slc_r]

        elif pad_mode == 'order0':
            # The `_padding_slices_inner` helper returns the slices for the
            # boundary values.
            left_slc, right_slc = _padding_slices_inner(
                lhs_arr, rhs_arr, axis, offset, pad_mode)

            if direction == 'forward':
                rhs_slc_l[axis] = left_slc
                rhs_slc_r[axis] = right_slc
                lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r = map(
                    tuple, [lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r])

                lhs_arr[lhs_slc_l] = lhs_arr[rhs_slc_l]
                lhs_arr[lhs_slc_r] = lhs_arr[rhs_slc_r]
            else:
                lhs_slc_l[axis] = left_slc
                lhs_slc_r[axis] = right_slc
                lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r = map(
                    tuple, [lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r])

                lhs_arr[lhs_slc_l] += np.sum(
                    lhs_arr[rhs_slc_l],
                    axis=axis, keepdims=True, dtype=lhs_arr.dtype)
                lhs_arr[lhs_slc_r] += np.sum(
                    lhs_arr[rhs_slc_r],
                    axis=axis, keepdims=True, dtype=lhs_arr.dtype)

        elif pad_mode == 'order1':
            # Some extra work necessary: need to compute the derivative at
            # the boundary and use that to continue with constant derivative.

            # Slice for broadcasting of a 1D array along `axis`
            bcast_slc = [None] * lhs_arr.ndim
            bcast_slc[axis] = slice(None)
            bcast_slc = tuple(bcast_slc)

            # Slices for the boundary in `axis`
            left_slc, right_slc = _padding_slices_inner(
                lhs_arr, rhs_arr, axis, offset, pad_mode)

            # Create slice tuples for indexing of the boundary values
            bdry_slc_l = list(working_slc)
            bdry_slc_l[axis] = left_slc
            bdry_slc_l = tuple(bdry_slc_l)
            bdry_slc_r = list(working_slc)
            bdry_slc_r[axis] = right_slc
            bdry_slc_r = tuple(bdry_slc_r)

            # For the slope at the boundary, we need two neighboring points.
            # We create the corresponding slices from the boundary slices.
            slope_slc_l = list(working_slc)
            slope_slc_l[axis] = slice(left_slc.start, left_slc.stop + 1)
            slope_slc_l = tuple(slope_slc_l)
            slope_slc_r = list(working_slc)
            slope_slc_r[axis] = slice(right_slc.start - 1, right_slc.stop)
            slope_slc_r = tuple(slope_slc_r)

            # The `np.arange`s, broadcast along `axis`, are used to create the
            # constant-slope continuation (forward) or to calculate the
            # first order moments (adjoint).
            arange_l = np.arange(-n_pad_l, 0,
                                 dtype=lhs_arr.dtype)[bcast_slc]
            arange_r = np.arange(1, n_pad_r + 1,
                                 dtype=lhs_arr.dtype)[bcast_slc]

            lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r = map(
                tuple, [lhs_slc_l, rhs_slc_l, lhs_slc_r, rhs_slc_r])

            if direction == 'forward':
                # Take first order difference to get the derivative
                # along `axis`.
                slope_l = np.diff(lhs_arr[slope_slc_l], n=1, axis=axis)
                slope_r = np.diff(lhs_arr[slope_slc_r], n=1, axis=axis)

                # Finally assign the constant slope values
                lhs_arr[lhs_slc_l] = lhs_arr[bdry_slc_l] + arange_l * slope_l
                lhs_arr[lhs_slc_r] = lhs_arr[bdry_slc_r] + arange_r * slope_r
            else:
                # Same as in 'order0'
                lhs_arr[bdry_slc_l] += np.sum(lhs_arr[rhs_slc_l],
                                              axis=axis, keepdims=True,
                                              dtype=lhs_arr.dtype)
                lhs_arr[bdry_slc_r] += np.sum(lhs_arr[rhs_slc_r],
                                              axis=axis, keepdims=True,
                                              dtype=lhs_arr.dtype)

                # Calculate the order 1 moments
                moment1_l = np.sum(arange_l * lhs_arr[rhs_slc_l],
                                   axis=axis, keepdims=True,
                                   dtype=lhs_arr.dtype)
                moment1_r = np.sum(arange_r * lhs_arr[rhs_slc_r],
                                   axis=axis, keepdims=True,
                                   dtype=lhs_arr.dtype)

                # Add moment1 at the "width-2 boundary layers", with the sign
                # corresponding to the sign in the derivative calculation
                # of the forward padding.
                sign = np.array([-1, 1])[bcast_slc]
                lhs_arr[slope_slc_l] += moment1_l * sign
                lhs_arr[slope_slc_r] += moment1_r * sign

        if direction == 'forward':
            working_slc[axis] = full_slc[axis]
        else:
            working_slc[axis] = intersec_slc[axis]
Apply padding to ``lhs_arr`` according to ``pad_mode``.

    This helper assigns the values in the excess parts (if existent)
    of ``lhs_arr`` according to the provided padding mode.

    This applies to the following values for ``pad_mode``:
    ``periodic``, ``symmetric``, ``order0``, ``order1``

    See `the online documentation
    <https://odlgroup.github.io/odl/math/resizing_ops.html>`_
    on resizing operators for details.