I would agree with the OP's concerns that the given approach is not best since it is most sensitive to noise. It is reading the absolute peak to peak value, which includes the underlying "truth" of what the peak to peak signal actually is, in addition to any inevitable noise that may be on the sample.
Ideally a standard deviation calculation would be computed or equivalently a time-aligned correlation to the reference 60 Hz signal. For sub-optimum but efficient and far superior implementations when processing resources are a concern, I would suggest an mean absolute deviation (MAD) computation since this can be done with a simple moving average filter, or if further simplication is desired, a CIC filter structure. Conveniently the MAD is related to the standard deviation we seek by a constant $4/(\pi\sqrt{2})$ (The constant $\sqrt{2/\pi}$ that is given in the referenced link is specific to measuring Gaussian distributed noise, but for a sinusoidal distribution the ratio of MAD to RMS is $4/(\pi\sqrt{2}))$. As long as the harmonics are more than 15 dB down or more (very likely as this is a typical requirement for AC supplied to computer equipment but should be verified if accuracy is important), this should result in a very reasonable estimate.
When processing is not a concern, simply compute the RMS value directly (average of the sum of the squares for each sample), or if ultimately a power measurement is desired, a current measurement is also needed and the power will be the average of the products of the current and voltage samples (assuming a resistive load, read further details below on the importance of a phase measurement between the I and V waveforms for the case or reactive loads).
To compute the MAD, simply average the absolute value of each sample over a period of at least several cycles (the longer the average the more accurate the result, so this can be a trade of time versus accuracy). If further simplification is desired, then a CIC filter (which is even simpler) can be implemented. A CIC filter is just an efficient implementation of a moving average.
To convert the resulting moving average to the rms value, simply multiply it by $\pi\sqrt{2}/4$.
This will provide an efficient measurement of the rms voltage, but to measure power you will also need a measure of the current, and the phase relationship between the two (both of which can also be done very efficiently) as the average power is $V_{rms} I_{rms} \cos(\theta)$. There are plenty of low cost current sensors available that will measure the current using inductive coupling (more details on that would be more appropriate on the electronics stack-exchange site so post there if that solution is still elusive). To get an estimate of the phase between the voltage and current readings, you could multiply and scale the voltage and current readings as the average of that product is the $\cos(\theta)$ term desired above:
$$V_p(\cos(\omega t)I_p(\cos(\omega t+\theta)) = \frac{V_pI_p}{2}\cos(\theta)$$
Plus a higher frequency "2f" term at $\cos(2\omega t +\theta)$ that we simply filter out (with another moving average filter!). In cases when the load is resistive or mostly so, $\theta=0$ and $\cos(\theta)=1$. We already know $V_p$ and $I_p$ from the MAD results, since the peak is related to the MAD by another constant $\pi/2$.
If the phase measurement was necessary and done as above with the sample by sample product, then the power can be computed directly by combining the phase and MAD processing, but for power line monitoring knowing if the load is reactive or not may be desirable to see as derived information from this processing.
To simplify the phase estimate, the product phase detector could be implemented with an XOR gate and RC low pass filter of the detected current and voltage signals (converted to pulses), or similarly in the Arduino using an XOR operation from the sign of those two signals. Both of these approaches will produce a linear result vs phase (when normalized linearly going from 1 to 0 as phase goes from $0$ to $\pi/2$, and from 0 to 1 as phase goes from -$\pi/2$ to $0$.
I provide more details on how a CIC filter is structured and its equivalent to a moving average filter here.