The local zeta-functions of some cubic plane curves This question is motivated by the work presented in article 358 of Gauss' Disquisitiones Arithmeticae. For the sake of completeness, let me say something about the background and present the question at the end. Any comment or correction is appreciated.
Let $n$ be an odd prime, and suppose further that $n=3m+1$ for some integer $m$. Let $g$ be a generator of the cyclic group $(\mathbb{F}_n)^{\times}$, and let $H$ be its unique subgroup of index $3$. In the above mentioned work, Gauss denotes the cosets $H$, $gH$ and $g^2H$ respectively by $R$, $R'$ and $R''$. He further denotes by $p$ (respectively, $p'$ and $p''$) the Gauss sums
$$\sum_{w}\exp{2\pi iw},\textrm{ for $w$ running over $R$ (respectively, $R'$ and $R''$).}$$
He then studies the polynomial equation
\begin{equation}\label{question}
x^3-Ax^2+Bx-C=0
\end{equation}
of which the roots are $p$, $p'$ and $p''$. 
One easily find $A=-1$ since $A=p+p'+p''$ is the sum of all $n$-th roots of unity except $1$. 
To determine $B$ and $C$, Gauss considers the set
$$RR:=\{w\in (\mathbb{F}_n)^{\times}|w\in R, w+1\in R'\}$$
and similarly defined sets $RR'$, $RR''$, $R'R''$, etc., and denotes by $(RR)$, $(RR')$, etc. their cardinalities. He further determines
\begin{equation}\label{abc}
\begin{cases}
(R'R'')=(R''R')=:a,\\
(R''R'')=(R'R)=(RR')=:b,\\
(R'R')=(R''R)=(RR'')=:c.
\end{cases}
\end{equation}
and
\begin{equation}\label{a+b+c}
a+b+c=m.
\end{equation}
To determine $B$, Gauss uses elementary yet sophisticated tricks to express the products $pp'$, $pp'$ and $p'p''$ in terms of linear combinations of $p$, $p'$ and $p''$ with coefficients $a$, $b$ and $c$. Then using a symmetry argument, Gauss proves
$$B=m(p+p'+p'')=-m.$$
By various computational exertions, Gauss proves
$$C=a^2-bc,$$
and 
\begin{equation}\label{long}
a^2+b^2+c^2-a=ab+bc+ac,
\end{equation}
which yields
\begin{equation}\label{4n}
4n=12a+12b+12c+4=(6a-3b-3c-2)^2+27(b-c)^2.
\end{equation}
Gauss then shows that the number $4n$ (for $n$ prime, or see @Will Sawin's comment) can be written in a unique way as $M^2+27N^2$ for integers $M$, $N$. Hence the values $a$, $b$, $c$, and therefore $C$, are completely determined.
Now notice that $a$, $b$ and $c$ give rise to numbers of solutions of some cubic equations over $\mathbb{F}_n$, i.e., numbers of points of some curves over $\mathbb{F}_n$. For example, $a+3=(RR)+4$ is the number of solutions of the equation 
\begin{equation}\label{final}
x^3-y^3=1.
\end{equation}
The above argument, in principle, determines all the coefficients $N_1$ of the local zeta-function of (the associated projective variety of) the above equation
$$Z(u)=\exp{\sum_{\nu=1}^{\infty}N_{\nu}\frac{u^{\nu}}{\nu}}.$$
As mentioned in this wikipeida page, this is the earliest known non-trivial cases of local zeta-functions. 
My question is: How does one (if possible) start from the above work of Gauss and express $Z(u)$ as a rational function in $u$?
 A: We have $3a+3b+3c+1 = n$ so $M + n = 9a+1$. Probably you mean $9a+3$ is the number of solutions of the equation, and not $a$, as each nonzero number has three cube roots. Except shouldn't it be $9a+6$ as there are also the solutions with $x$ or $y$ zero? The constant term will not affect whether the function is rational, but it's good to get it right - I think the correct formula should end up $M+n-2$, based on what I already know about the rational function.
Given the unique solution $4p = M_1^2 + 27 N_1^2$, we can write $ \alpha=  \frac{-M_1 + \sqrt{-27} N_1 }{ 2} $ and observe that $|\alpha| = p$, and so $|\alpha^r| = p^r$. Then $\alpha$ is an algebraic integer because $M_1$ is congruent to $N$ mod $2$, so $\alpha^r$ can be written as $\frac{ -M_r + \sqrt{-27} N_r } {2}$ and thus $M_r =- ( \alpha^r + \overline{\alpha}^r)$ (if we ignore the non-uniqueness issues, which are discussed below). This gives $M_r + p^r = - (\alpha^r +\overline{\alpha}^r) + p^r$, and plugging $N_r = \alpha^r +\overline{\alpha}^r + p^r+c$ for a constant $c$ into your equation gives you the rational function $$\frac{ (1- \alpha u) (1- \overline{\alpha} u) }{ (1- u)^c (1-pu) }.$$
Presumably Gauss would have done this, not with algebraic numbers, but with writing down explicit polynomials and possibly something to do with quadratic forms.
Now what's not justified here is how to find the right value for $M_r$. We have the equation $M_r^2 + 27 N_r^2 = 4p^r$ and we know that one solution to this equation gives what we want. It seems to me (using unique factorization in $\mathbb Z[\mu_3]$) that there is a unique solution, up to sign, with $M_r$ and $N_r$ prime to $p$. I think we can check that $M_r$ and $N_r$ are prime to $p$ with a separate argument. Possibly the sign can be controlled with congruence conditions mod $2$ or mod $3$, and I guess this is how Gauss would have done it.
The easiest way for us to do this nowadays is probably the Hasse-Davenport relations between Gauss sums over different finite fields. But because neither "Hasse" nor "Davenport" is "Gauss", this may not have been available to Gauss.
