Octave - freqz strange results
What you are seeing are indeed numerical artefacts. The conversion from polynomial roots to polynomial coefficients (and vice-versa) is known to be an ill-conditioned problem. In your case the problem occurs at a relatively low polynomial order because zeros and poles are all clustered around the same point (DC), and the poles and zeros almost cancel each other. By removing one zero, as you tried, the problem's numerical condition is improved, simply because the polynomial order becomes smaller. Another way to avoid the problem is to distribute the poles and zeros more evenly around the circle in the complex plane. Note that the available frequencies extend up to $f_2/2=22050\text{ Hz}.$ That means that just decreasing the sampling frequency would also improve the numerical condition.

Apart from the numerical effects, your approach for converting zeros and poles to numerator and denominator coefficients is correct. However, it can be done in a much simpler way by using the function zp2tf, which transforms given poles and zeros to transfer function coefficients.

Also note that with your approach you generate a double zero at zero frequency because you produce complex conjugate pairs, even for real-valued zeros.

Finally, the filter you computed is not causal and stable. This has nothing to do with the effects you see, but I believe it's still worth mentioning. So you wouldn't be able to implement that filter with the coefficients you computed. The reason for this is that the transfer function is not proper, i.e., its numerator order is greater than the order of its denominator. This means that apart from the given poles, there occur a number poles at infinity, and that number equals the difference between the numerator and denominator orders. This problem can be solved by adding poles at $z=0$, i.e., by multiplying the denominator by $z^K$ where $K$ is the difference between the polynomial orders of numerator and denominator. This will not change the frequency response, it will just make the filter causal and stable.

EDIT:

For maximum precision you shouldn't compute the coefficients of the numerator and denominator polynomials. Instead you should work with second-order sections. Below is a simple example in Octave/Matlab. Note that I didn't optimize the sequence of sections, neither did I optimize the pole-zero pairing per section. But, as you can see in the figure, the numerical artefacts are greatly reduced.


fs = 44100;
rz = 1; rp = 0.9999;
phi_z = 2*pi*[30; 45.8; 66] / fs;
phi_p = 2*pi*[41.7; 60; 95.2 ] / fs;

m = length(phi_p);
L = 2^15;

a = zeros(m,3); b = zeros(m,3);
H = zeros(m,L);

for k = 1:m,
    b(k,:) = [1, -2*rz*cos(phi_z(k)), rz^2];
    a(k,:) = [1, -2*rz*cos(phi_p(k)), rp^2];
    Htmp = freqz( b(k,:), a(k,:), L );
    H(k,:) = Htmp(:).';
end

% add real-valued zero at DC
[H0,w] = freqz( [1, -rz], 1, L );
H = [ H0(:).'; H ];

Htot = prod(H);
semilogx( w/pi+1e-5, abs(Htot) );
axis([1e-5,.01,0,.0005])