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