I would agree with John's suggestion that the most "fundamental" definition of momentum is the quantity that is conserved when the Lagrangian is invariant under spatial translations $\mathbf{x}\rightarrow\mathbf{x}+\mathbf{s}$.
However, do terms like Lagrangian even make sense when we try to include massless objects?
Yes, absolutely. In fact the machinery used to describe such particles, quantum field theory, is expressed in terms of Lagrangians.
Lagrangian mechanics and generalized momentum
Let's start with the simpler classical dynamics notion of momentum. We start with a Lagrangian, expressed in terms of the generalized coordinates $q_i$ and their time derivatives,
$$L(q_i, \dot{q}_i)=T(\dot{q}_i)-V(q_i)$$
where $T$ and $V$ are the kinetic and the potential energies, respectively. The reason we care about the Lagrangian is that it allows us to easily derive the equations of motion for the system under consideration using the Euler-Lagrange equations,
$$\frac{\partial L}{\partial q_i}-\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{q}_i}\right)=0$$
which are a direct result of the principle of least action. There is a famous and beautiful theorem in mechanics by Emmy Noether, which states that for every continuous symmetry of the Lagrangian there exists some quantity that stays the same throughout time. More precisely, for some continuous transformation parameterised by $s$, $q_i(t)\rightarrow Q_i(s,t)$, if
$$\frac{\partial}{\partial s}L(Q_i(s,t),\dot{Q}_i(s,t),t)=0$$
then the quantity
$$\sum_i\frac{\partial L}{\partial\dot{q}_i}\frac{\partial Q_i}{\partial s}$$
is conserved though all time. These $\partial L/\partial\dot{q}_i$ are known as generalized momenta.
Defining momentum
Consider the specific case of a collection of free particles, which has a Lagrangian
$$L=\frac{1}{2}\sum_i m_i\dot{\mathbf{x}}_i^2-V(|\mathbf{x}_i-\mathbf{x}_j|)$$
Now consider making a translation $\mathbf{x}_i\rightarrow\mathbf{x}_i+s\mathbf{n}$ where $\mathbf{n}$ is an arbitrary vector. The above Lagrangian is clearly invariant under such a transformation. Therefore we know per Noether's theorem that the quantity
$$\sum_i\frac{\partial L}{\partial\dot{\mathbf{x}}_i}\frac{\partial(s\mathbf{n})}{\partial s}=\sum_i m_i\dot{\mathbf{x}}_i\cdot\mathbf{n}$$
is conserved. But this is just our usual conservation of the total momentum $\sum\mathbf{p}_i=\sum m_i\dot{\mathbf{x}}_i$ along the direction of $\mathbf{n}$! So we see that the quantity momentum naturally arises as the thing that is conserved under space-translation symmetry.
But this notion of momentum extends much further than free particles. It can be applied to any Lagrangian with space-translation symmetry, with the conserved quantities that arise being given the label momentum.
So now, how do massless particles like photons fit into this? For describing particles such as photons we use quantum field theory, however it suffices here to stick to the classical theory.
Lagrangian of fields
Like the name suggests, in field theory our dynamical variables are fields which are themselves functions of space and time (or, spacetime). As such we replace the usual Lagrangian with a Lagrangian density, $\mathcal{L}$ which is defined by
$$L=\int\mathrm{d}^3\mathbf{x}\,\mathcal{L}(\phi(x),\partial_\mu\phi(x),t)$$
where $\phi(x)$ is the field as a function of spacetime coordinate $x=(t,\mathbf{x})$ and $\partial_\mu=\partial/\partial x^\mu=(\partial/\partial t,\nabla)$ is the four-gradient.
Despite all these apparent differences Noether's theorem remains more-or-less unchanged, except now what matters is symmetry in the fields rather than the coordinates. You know how we considered a spatial translation before? Well there is in fact also a conserved quantity under time translational symmetry: energy. In field theory, we can go one step further and consider spacetime translations, $x^\mu\rightarrow x^\mu+\varepsilon^\mu$. I won't go into the details of the derivation here (you can find it in any decent field theory book), but the conserved quantity (called a conserved current now) that pops out is something called the energy-momentum tensor:
$$\mathcal{T}^{\mu\nu}=\frac{\partial\mathcal{L}}{\partial(\partial_\mu\phi)}\partial^\nu\phi-\eta^{\mu\nu}\mathcal{L}$$
where $\eta^{\mu\nu}$ is the Minkowski metric. The four-momentum is then identified with the $0\mu$ components:
$$p^\mu=\int\mathrm{d}^3\mathbf{x}\,\mathcal{T}^{0\mu}$$
From which the 3-momentum can be defined as $p^i$.
As an example, consider the Lagrangian for electrodynamics in the absence of sources:
$$\mathcal{L}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$$
where
$$F^{\mu\nu}=\partial^\mu A^\nu-\partial^\nu A^\mu=\begin{pmatrix}0 & -E_x & -E_y & -E_z \\
E_x & 0 & -B_z & B_y \\
E_y & B_z & 0 & -B_x \\
E_z & -B_y & B_x & 0
\end{pmatrix}$$
is the EM tensor, with $A^\mu=(\phi,\mathbf{A})$ being the four-potential. The energy-momentum tensor for this theory is (in SI units now, for clarity)
$$\mathcal{T}^{\mu\nu} = \frac{1}{\mu_0} \left[ F^{\mu \alpha}F^\nu{}_{\alpha} - \frac{1}{4} \eta^{\mu\nu}F_{\alpha\beta} F^{\alpha\beta}\right]=\begin{pmatrix} \frac{1}{2}\left(\varepsilon_0 E^2+\frac{1}{\mu_0}B^2\right) & S_\text{x}/c & S_\text{y}/c & S_\text{z}/c \\
S_\text{x}/c & -\sigma_\text{xx} & -\sigma_\text{xy} & -\sigma_\text{xz} \\
S_\text{y}/c & -\sigma_\text{yx} & -\sigma_\text{yy} & -\sigma_\text{yz} \\
S_\text{z}/c & -\sigma_\text{zx} & -\sigma_\text{zy} & -\sigma_\text{zz} \end{pmatrix}$$
where $\mathbf{S}=(\mathbf{E}\times\mathbf{B})/\mu_0$ is the Poynting vector and $\sigma_{ij}$ is the Maxwell stress tensor. I leave it as an exercise to verify that for electromagnetic waves we indeed have $E=|\mathbf{p}|c$ (i.e. $(\mathcal{T}^{00})^2=\sum_i(\mathcal{T}^{0i})^2$).
