function y = cosdx(x)
% cosd with range reduction of x carried out exactly
  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);
end
