Normalizing posterior distribution

This is from task 5.3. in Bayesian Data Analysis 3rd edition by Gelman et al. It deals with a hierarchical model where I am supposed to simulate the posterior distribution for $tau$, which is the standard deviation for $theta$. The details on the model is not really important here. My question deals with the normalization of the simulated posterior distribution:

# 5.3
# data
y = c(28,8,-3,7,-1,1,18,12)
s = c(15,10,16,11,9,11,10,18)
J = length(y)

# get the likelihood for each value of tau
p.tau = function(tau,y,s) {
  num = sum(1/(s^2+tau^2))
  mu.hat = sum(1/(s^2+tau^2)*y / num)
  V.mu = 1 / num
  p = sqrt(V.mu)*prod(1/sqrt(s^2+tau^2))*exp(-sum((y-mu.hat)^2 / (2*(s^2+tau^2))))
}

# calculate marginal distribution for tau
tau = seq(0.001,30,length=100)
p = rep(NA, length(tau))
for (i in 1:length(tau)) {
  p[i] = p.tau(tau[i],y,s) 
}
delta = diff(tau[1:2])
p = p/sum(p*delta)
plot(tau,p,type="l")

The marginal distribution needs to be normalized, and I know a density is supposed to sum to 1. This would normally be done by

p = p/sum(p)

However, here we were asked to do

p = p/sum(p*delta)

with

delta = diff(tau[1:2])

This does not make the density sum to 1. So why are we supposed to do it this way?