Question on zp2sos gain scaling?

Hello guys, I have a question on matlab zp2sos function.

The question is how matlab do scaling when convert z/p/k to zeros and poles. Based on my understanding, if I choose 'up' and 'inf' option in 'zp2sos' function. The matlab should scale the inter-section gain to be unity max magnitude, i.e., the n-stage sos's forward branch convoluted with the n+1 stage sos's backward branch, this frequency domain magnitude response should be of max value one.

However, I check this in one of my design filter. That's not true. See below for details. The calculated max1 is around 1.17, not 1. So, sometimes, it will introduce overflow.

Is there something wrong in my procedure or a bug of Matlab? (Note: my matlab version is 8.3.0.532 (R2014a))

Fs = 5e6;  % Sampling Frequency
% highpass filter (cheby)
N     = 4;      % Order
Fpass = 13000;  % Passband Frequency
Apass = 1;      % Passband Ripple (dB)

% Construct an FDESIGN object and call its CHEBY1 method.
h  = fdesign.highpass('N,Fp,Ap', N, Fpass, Apass, Fs);
Hd = design(h, 'cheby1');

% get z-p-k
[z_m, p_m, k_m] = Hd.zpk;
[sos_m, g_m] = zp2sos(z_m, p_m, k_m, 'up', 'inf'); 

% check stage gain
sos_0 = sos_m(1, :);
sos_1 = sos_m(:, :);
sos_0(1:3) = [1 0 0];
sos_1(1,4:6) = [1 0 0];
sos_1(2,1:3) = [1 0 0];

Hd0 = dfilt.df2sos(sos_0);
% L-inf norm equals max of frequency domain magnitude.
max0 = norm(Hd0, 'Linf')*g_m; 
Hd1 = dfilt.df2sos(sos_1);
max1 = norm(Hd1, 'Linf');