I take the question to be "How did someone come up with the particular choice of $x^2$ and $x$ to plug into $g$ and discover discontinuity?" Of course, not being a mind-reader, I don't know how it was actually done, but here's how it might have been done. The key observation is that, in the fraction that defines $g(x,y)$, the numerator $xy^2$ is the geometric mean of the two terms in the denominator, $x^2$ and $y^4$. (That might sound complicated, but it really just means that the exponents of $x$ and $y$ in the numerator are the averages of their exponents in the two terms of the denominator: For the exponents of $x$, $1$ is the average of $2$ and $0$, and for the exponents of $y$, $2$ is the average of $0$ and $4$.) So if we give $x$ and $y$ values that make the two denominator terms equal, then the numerator will automatically be equal to them also. So the numerator will match each of the terms in the denominator, and the fraction will be $1/2$. If we can find such values for $x$ and $y$ arbitrarily close to $0$, then that will make $g$ discontinuous at $(0,0)$. The easiest way to achieve this, meaning to make $x^2=y^4$, is to give $y$ an arbitrary value, say $t$, and to set $x=t^2$. So you set $(x,y)=(t^2,t)$ to find points, as close to $(0,0)$ as you like if $t$ is very small, where $g$ takes the value $1/2$.
Finally, if you're in a nasty mood, you re-name the variable $t$ as $x$, even though it serves as the value to substitute for $y$, just to confuse readers.
P.S. If you know how to use a program like Mathematica, you can have it plot the graph of $g$, and you'll probably be able to see a sort of a ridge in the graph, over the parabola $x=y^2$, at height $z=1/2$. So this might be another way to "guess" the substitution that proves discontinuity of $g$ at $(0,0)$.
I think you also forgot to normalize your vector $\,\vec p\,$ . Taking it from where you were, and assuming we already have $\,||\vec p||=1\,$ , we have by definition with $\,\bf x:=(a_1,a_2,...,a_n)\;$
$$D_{\vec p}f(\bf x):=\text{The directional derivative of $\,f\,$ at point $\,\bf x\,$ in the direction of $\,\vec p$}:=$$
$$:= \lim_{t\to 0}\frac{f(x+t\vec p)-f(x)}t\stackrel{\text{Chain Rule}}=\sum_{k=1}^n \vec p\cdot \frac{\partial f}{\partial x_k}(a_1,...,a_n)=\nabla f(\bf{x})\cdot\bf \vec p$$
So if the function is differentiable at any point (say, when the partial derivatives exist and are continuous at any point), its directional derivative in any direction exists and it's pretty easy to calculate by means of the gradient of $\,f\,$ .
Added: We want the directional derivative of $\;f\;,\;f:\Bbb R^n\to\Bbb R\;$ at the point $\,x:=(x_1,...,x_n)\,$ and in the direction $\,p:=(p_1,...,p_n)\,$ (no little arrows, no nothing) .
We also assume $\,||p||:=\sqrt{\sum_{k=1}^n p_i^2}=1\;$ (note then that you will have to normalize the vector in which direction you want the derivative in case it is not normalized!).
Now we define a new function
$$g:\Bbb R\to\Bbb R^n\;,\;\;g(t):=x+tp=\left(x_1+tp_1\,,\,x_2+tp_2\,,\ldots,x_n+tp_n\right):=(g_1(t),...,g_n(t))$$
Note that $\;\forall\,t\;,\;\,g'(t)=p\implies g'(0)=p\,$ , and we also denote the derivative of some function $\,k\,$ by $\,Dk\,$ and by $\,D_uk\,$ the directional derivative of $\,k\,$ in the direction of $\,u\,$, for simplicity.
Thus, by definition we get that
$$D_pf(x):=\lim_{t\to 0}\frac{f(x+tp)-f(x)}t=\left.\frac d{dt}(f(x+tp))\right|_{t=0}=\left.D(f\circ g)(t)\right|_{t=0}\stackrel{\text{chain rule}}=$$
$$=Df(g(0))\cdot Dg(0)=\left(f'_{x_1}(g(0))\,,\,f'_{x_2}(g(0))\,,\ldots,f'_{x_n}(g(0))\right)\cdot(g'_1(0)\,,\,\ldots,g'_n(0))=$$
$$\left(f'_{x_1}(g(0))\,,\,f'_{x_2}(g(0))\,,\ldots,f'_{x_n}(g(0))\right)\cdot(p_1,...,p_n)=\nabla f(x)\cdot p$$
which is what we wrote above...:)
IMPORTANT: You can use the gradient as above to calculate the directional derivative only if you're sure the partial derivatives exist in some neighborhood of $\,x\,$ and are continuous there.
Best Answer
Your limit $(1)$ is valid for any path of approach towards the origin. So choose the path $y = \frac{d}{c}x$, and plug that in:
$$ 0 = \lim_{x \to 0} \frac{f(x,\frac{d}{c}x) - \left(a+ \frac{bd}{c}\right)x - f(0,0)}{\frac{x}{c} \sqrt{c^2+d^2}} = \lim_{x \to 0} \frac{f(x,\frac{d}{c}x)-f(0,0)}{\frac{x}{c}\sqrt{c^2+d^2}} - \frac{ ac + bd }{\sqrt{c^2+d^2}} $$
Move the negative term to the other side, and cancel the $\sqrt{c^2+d^2}$ factor:
$$ ac+bd = \lim_{x \to 0} \frac{f(x,\frac{d}{c}x) - f(0,0)}{x/c} $$
Just re-parameterize $t = \frac{x}{c}$, and this becomes
$$ ac + bd = \lim_{t \to 0} \frac{f(ct,dt)-f(0,0)}{t} $$