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.
Here's a real-analytic example: Let $g(t) =1/(1+t^2).$ Define
$$f(x,y) = \begin{cases} xg((y-x^2)/(x^2+y^2)^2), (x,y) \ne (0,0) \\
0, (x,y) =(0,0).\end{cases}$$
Then $f$ is real-analytic on $\mathbb {R}^2 \setminus \{(0,0)\},$ $f$ is continuous at $(0,0),$ and $D_uf(0,0) = 0$ for all unit vectors $u.$ Because $f(x,x^2) = xg(0) = x\ne o(x,x^2),$ $f$ is not differentiable at $(0,0).$
Best Answer
If the directional derivative is zero for $n$ linearly independent vectors, then it is zero in every direction, since the linearly independent vectors span the entire vector space. More explicitly, $D_{\vec{y}} f(\vec{x}) = \langle \nabla f(\vec{x}),\vec{y} \rangle$, and if $\nabla f(\vec{x})$ is orthogonal to $n$ linearly independent vectors then it is orthogonal to $\mathbb{R}^n$ and must be zero. Your mean value theorem now shows that the function is constant.