function y = sindx(x)
% sind 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, 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 sin(-x) = -sin(x)
  t = y == 0; y(t) = (1-2*signbit(x(t))) * 0;

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