import math


def dbsimpson(f, limits: list, d: list):
    """Simpson's 1/3 rule for double integration

    int_{ay}^{by} int_{ax}^{bx} f(x,y) dxdy

    Args:
        f (func): two variable function, must return float or ndarray
        limits (list): limits of integration [ax, bx, ay, by]
        d (lsit): list of integral resolution [dx, dy]

    Returns:
        float: double integral of f(x,y) between the limits
    """
    ax, bx, ay, by = limits
    dx, dy = d
    nx = math.floor((bx - ax) / dx)
    ny = math.floor((by - ay) / dy)
    s = 0
    for i in range(ny + 1):  # loop of outer integral
        if i == 0 | i == ny:
            p = 1
        elif i % 2 != 0:
            p = 4
        else:
            p = 2

        for j in range(nx + 1):  # loop of inner integral
            if j == 0 | j == nx:
                q = 1
            elif j % 2 != 0:
                q = 4
            else:
                q = 2
            x = ax + j * dx
            y = ay + i * dy
            s += p * q * f(x, y)

    return dx * dy / 9 * s