function y = cosdx(x)
% cosd with range reduction of x carried out exactly
  if iscomplex(x)
    xr = real(x); xi = imag(x) * (pi / 180);
    y = complex(cosd(xr) .* cosh(xi), -sind(xr) .* sinh(xi));
    return
  end
  y = rem(x, 360);
  q = floor(y / 90 + 0.5);
  y = y - 90 * q;
  y = y * (pi / 180);                   % y in [-45, 45]
  q = mod(q + 1, 4);
  t = bitand(q, 1) == 0; y(t) = sin(y(t)); y(~t) = cos(y(~t));
  t = bitand(q, 2) ~= 0; y(t) = -y(t);
  % enforce cos(-x) = cos(x)
  y = 0 + y;                            

  % Special values follow the specifications for cosPi in IEEE Std 754-2008
  % a=[-inf -540 -450 -360 -270 -180 -90 -0 0 90 180 270 360 450 540 inf nan];
  % cosdx(a) ->
  % [NaN -1 0 1 0 -1 0 1 1 0 -1 0 1 0 -1 NaN NaN]
end
