A comment : I think the best naive way to decide the number of components to retain is to base your estimate on some threshold of sample variation you would like to retain in your reduced dimensionality sample rather than just some arbitrary number eg. 3, 100, 200. As user4959 explained you can check that cumulative variation by checking the relevant field of the list under the $loadingsfield in the list object produced by princomp.

I believe your original question stems from being a bit uncertain about what PCA is doing. Principal Component Analysis allows you to identify the principal mode of variation in your sample. Those modes are emperically calculated as the eigenvectors of your sample's covariance matrix (the "loadings"). Subsequently those vectors serve as the new "coordinate system" of your sample as you project your original sample in the space they define (the "scores"). The proportion of variation associated with the $i$-th eigenvector/ mode of variation/ loading/ principal component equals $\frac{\lambda_i}{\Sigma_{k=1}^{p} \lambda_k}$ where $p$ is your sample's original dimensionality ($p =784$ in your case). [Remember because your covariance matrix is non-negative Definite you'll have no negative eigenvalues $\lambda$.] Now, by definition the eigenvectors are orthogonal to each other. That means that their respective projections are also orthogonal and where originally you had a big possibly correlated sample of variables, now you have a (hopefully significantly) smaller linearly independent sample (the "scores").

Practically with PCA you are using the projections of the PCs (the "scores") as surrogate data for your original sample. You do all your analysis on the scores, and afterwards you reconstruct your original sample back using the PCs to find out out what happened on your original space (that's basically Principal Component Regression). Clearly, if you are able to meaningful interpreter your eigenvectors ("loadings") then you are in an even better position: You can describe what happens to your sample in the mode of variation presented by that loading by doing inference on that loading directly and not care about reconstruction at all. :)

In general what do you "after calculating the PCA" depends on the target of your analysis. PCA just gives you a linearly independent sub-sample of your data that is the optimal under an RSS reconstruction criterion. You might use it for classification, or regression, or both, or as I mentioned you might want to recognise meaningful orthogonal modes of variations in your sample.

I believe your original question stems from being a bit uncertain about what PCA is doing. Principal Component Analysis allows you to identify the principal mode of variation in your sample. Those modes are emperically calculated as the eigenvectors of your sample's covariance matrix (the "loadings"). Subsequently those vectors serve as the new "coordinate system" of your sample as you project your original sample in the space they define (the "scores"). The proportion of variation associated with the $i$-th eigenvector/ mode of variation/ loading/ principal component equals $\frac{\lambda_i}{\Sigma_{k=1}^{p} \lambda_k}$ where $p$ is your sample's original dimensionality ($p =784$ in your case). [Remember because your covariance matrix is non-negative Definite you'll have no negative eigenvalues $\lambda$.] Now, by definition the eigenvectors are orthogonal to each other. That means that their respective projections are also orthogonal and where originally you had a big possibly correlated sample of variables, now you have a (hopefully significantly) smaller linearly independent sample (the "scores").

Practically with PCA you are using the projections of the PCs (the "scores") as surrogate data for your original sample. You do all your analysis on the scores, and afterwards you reconstruct your original sample back using the PCs to find out out what happened on your original space (that's basically Principal Component Regression). Clearly, if you are able to meaningful interpreter your eigenvectors ("loadings") then you are in an even better position: You can describe what happens to your sample in the mode of variation presented by that loading by doing inference on that loading directly and not care about reconstruction at all. :)

In general what do you "after calculating the PCA" depends on the target of your analysis. PCA just gives you a linearly independent sub-sample of your data that is the optimal under an RSS reconstruction criterion. You might use it for classification, or regression, or both, or as I mentioned you might want to recognise meaningful orthogonal modes of variations in your sample.

A comment : I think the best naive way to decide the number of components to retain is to base your estimate on some threshold of sample variation you would like to retain in your reduced dimensionality sample rather than just some arbitrary number eg. 3, 100, 200. As user4959 explained you can check that cumulative variation by checking the relevant field of the list under the $loadingsfield in the list object produced by princomp.