extern "C" {

/*

table[i][...] = cos lo, sin hi, sin lo (24 bits), cos hi (1 bit)

cos hi + cos lo = cos(i*32/pi)

sin hi + sin lo = sin(i*32/pi)

*/

const double table [ 64 ] [ 4 ] = {

{ 0.0 , 0.0 , 0.0 , 1.0 } ,

{ - 0.004815273327803114 , 0.0980171403295606 , - 1.6345824060817233e-18 , 1.0 } ,

{ - 0.019214719596769552 , 0.19509032201612828 , - 7.991079223603255e-18 , 1.0 } ,

{ - 0.043059664267791134 , 0.2902846772544624 , - 1.8927978564718132e-17 , 1.0 } ,

{ - 0.07612046748871325 , 0.3826834323650898 , - 1.0050773010930708e-17 , 1.0 } ,

{ - 0.11807873565164496 , 0.47139673682599764 , 6.516678171793318e-18 , 1.0 } ,

{ - 0.16853038769745476 , 0.5555702330196022 , 4.709410873377549e-17 , 1.0 } ,

{ - 0.22698954663726303 , 0.6343932841636455 , 1.042090159346246e-17 , 1.0 } ,

{ 0.20710678118654752 , 0.7071067811865476 , - 4.833646783938113e-17 , 0.5 } ,

{ 0.1343932841636455 , 0.773010453362737 , - 3.256590550158814e-17 , 0.5 } ,

{ 0.055570233019602226 , 0.8314696123025452 , 1.4073856703154585e-18 , 0.5 } ,

{ - 0.02860326317400235 , 0.881921264348355 , - 1.9843248949424382e-17 , 0.5 } ,

{ - 0.11731656763491023 , 0.9238795325112867 , 1.7645047904427988e-17 , 0.5 } ,

{ 0.040284677254462366 , 0.9569403357322088 , 4.055386956300452e-17 , 0.25 } ,

{ - 0.05490967798387173 , 0.9807852804032304 , 1.8546939469906805e-17 , 0.25 } ,

{ - 0.026982859670439396 , 0.9951847266721969 , - 4.248691455392459e-17 , 0.125 } ,

{ 0.0 , 1.0 , 0.0 , 0.0 } ,

{ 0.026982859670439396 , 0.9951847266721969 , - 4.248691455392459e-17 , - 0.125 } ,

{ 0.05490967798387173 , 0.9807852804032304 , 1.8546939469906805e-17 , - 0.25 } ,

{ - 0.040284677254462366 , 0.9569403357322088 , 4.055386956300452e-17 , - 0.25 } ,

{ 0.11731656763491023 , 0.9238795325112867 , 1.7645047904427988e-17 , - 0.5 } ,

{ 0.02860326317400235 , 0.881921264348355 , - 1.9843248949424382e-17 , - 0.5 } ,

{ - 0.055570233019602226 , 0.8314696123025452 , 1.4073856703154585e-18 , - 0.5 } ,

{ - 0.1343932841636455 , 0.773010453362737 , - 3.256590550158814e-17 , - 0.5 } ,

{ - 0.20710678118654752 , 0.7071067811865476 , - 4.833646783938113e-17 , - 0.5 } ,

{ 0.22698954663726303 , 0.6343932841636455 , 1.042090159346246e-17 , - 1.0 } ,

{ 0.16853038769745476 , 0.5555702330196022 , 4.709410873377549e-17 , - 1.0 } ,

{ 0.11807873565164496 , 0.47139673682599764 , 6.516678171793318e-18 , - 1.0 } ,

{ 0.07612046748871325 , 0.3826834323650898 , - 1.0050773010930708e-17 , - 1.0 } ,

{ 0.043059664267791134 , 0.2902846772544624 , - 1.8927978564718132e-17 , - 1.0 } ,

{ 0.019214719596769552 , 0.19509032201612828 , - 7.991079223603255e-18 , - 1.0 } ,

{ 0.004815273327803114 , 0.0980171403295606 , - 1.6345824060817233e-18 , - 1.0 } ,

{ 0.0 , 0.0 , 0.0 , - 1.0 } ,

{ 0.004815273327803114 , - 0.0980171403295606 , 1.6345824060817233e-18 , - 1.0 } ,

{ 0.019214719596769552 , - 0.19509032201612828 , 7.991079223603255e-18 , - 1.0 } ,

{ 0.043059664267791134 , - 0.2902846772544624 , 1.8927978564718132e-17 , - 1.0 } ,

{ 0.07612046748871325 , - 0.3826834323650898 , 1.0050773010930708e-17 , - 1.0 } ,

{ 0.11807873565164496 , - 0.47139673682599764 , - 6.516678171793318e-18 , - 1.0 } ,

{ 0.16853038769745476 , - 0.5555702330196022 , - 4.709410873377549e-17 , - 1.0 } ,

{ 0.22698954663726303 , - 0.6343932841636455 , - 1.042090159346246e-17 , - 1.0 } ,

{ - 0.20710678118654752 , - 0.7071067811865476 , 4.833646783938113e-17 , - 0.5 } ,

{ - 0.1343932841636455 , - 0.773010453362737 , 3.256590550158814e-17 , - 0.5 } ,

{ - 0.055570233019602226 , - 0.8314696123025452 , - 1.4073856703154585e-18 , - 0.5 } ,

{ 0.02860326317400235 , - 0.881921264348355 , 1.9843248949424382e-17 , - 0.5 } ,

{ 0.11731656763491023 , - 0.9238795325112867 , - 1.7645047904427988e-17 , - 0.5 } ,

{ - 0.040284677254462366 , - 0.9569403357322088 , - 4.055386956300452e-17 , - 0.25 } ,

{ 0.05490967798387173 , - 0.9807852804032304 , - 1.8546939469906805e-17 , - 0.25 } ,

{ 0.026982859670439396 , - 0.9951847266721969 , 4.248691455392459e-17 , - 0.125 } ,

{ 0.0 , - 1.0 , 0.0 , 0.0 } ,

{ - 0.026982859670439396 , - 0.9951847266721969 , 4.248691455392459e-17 , 0.125 } ,

{ - 0.05490967798387173 , - 0.9807852804032304 , - 1.8546939469906805e-17 , 0.25 } ,

{ 0.040284677254462366 , - 0.9569403357322088 , - 4.055386956300452e-17 , 0.25 } ,

{ - 0.11731656763491023 , - 0.9238795325112867 , - 1.7645047904427988e-17 , 0.5 } ,

{ - 0.02860326317400235 , - 0.881921264348355 , 1.9843248949424382e-17 , 0.5 } ,

{ 0.055570233019602226 , - 0.8314696123025452 , - 1.4073856703154585e-18 , 0.5 } ,

{ 0.1343932841636455 , - 0.773010453362737 , 3.256590550158814e-17 , 0.5 } ,

{ 0.20710678118654752 , - 0.7071067811865476 , 4.833646783938113e-17 , 0.5 } ,

{ - 0.22698954663726303 , - 0.6343932841636455 , - 1.042090159346246e-17 , 1.0 } ,

{ - 0.16853038769745476 , - 0.5555702330196022 , - 4.709410873377549e-17 , 1.0 } ,

{ - 0.11807873565164496 , - 0.47139673682599764 , - 6.516678171793318e-18 , 1.0 } ,

{ - 0.07612046748871325 , - 0.3826834323650898 , 1.0050773010930708e-17 , 1.0 } ,

{ - 0.043059664267791134 , - 0.2902846772544624 , 1.8927978564718132e-17 , 1.0 } ,

{ - 0.019214719596769552 , - 0.19509032201612828 , 7.991079223603255e-18 , 1.0 } ,

{ - 0.004815273327803114 , - 0.0980171403295606 , 1.6345824060817233e-18 , 1.0 }

} ;

#define TWO_52_51 6755399441055744.0 // 2^52 + 2^51

#define ROUND(x) (((x) + TWO_52_51) - TWO_52_51)

/*

If |a| > |b| then

xh + xl = a + b (exactly),

xh = a + b (rounded)

*/

#define FAST2SUM(xh, xl, a, b) \

xh = (a) + (b); \

xl = ((a) - (xh)) + (b);

/*

pi/32 0.00011001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101110000000110111000001110011010001

Ch 0.0001100100100001111110110101010001

Cm 0.000000000000000000000000000000000000010000101101000110000100011010011

Cl 0.00000000000000000000000000000000000000000000000000000000000000000000000010011000110011000101000101110000000110111000001110011

*/

#define PI 3.1415926535897931

#define C_h 0.09817477042088285

#define C_m 3.798187816439979e-12

#define C_l 1.2639164054974691e-22

/*

Compute

sin(x) ~ sin(a + y) = sin(a)*cos(y) + cos(a)*sin(y)

where

sin(a), cos(a) from the table above

sin(y), cos(y) computed using a 9th order Taylor series approximation

If |x| >= 90112 then msvcrt uses the x87 FSIN

*/

__declspec ( dllexport ) double __cdecl msvcrt_sin ( double x )

{

// Cody and Waite's range reduction

double k = ROUND ( x * ( 32 / PI ) ) ;

double y_h1 = x - k * C_h ;

double y_h , y_l ;

FAST2SUM ( y_h , y_l , y_h1 , - k * C_m ) ;

y_l = k * C_l - y_l ;

const double * t = table [ ( int ) k & 63 ] ;

double sina_h = t [ 1 ] ;

double sina_l = t [ 2 ] ;

double cosa = t [ 0 ] + t [ 3 ] ;

// sin(a) + yh*cos(a) with low terms

double u , syc_h , syc_l1 , syc_l2 ;

FAST2SUM ( u , syc_l1 , sina_h , y_h * t [ 3 ] ) ;

FAST2SUM ( syc_h , syc_l2 , u , y_h * t [ 0 ] ) ;

double Y2 = y_h * y_h ;

double Y4 = Y2 * Y2 ;

// Taylor series, Horner's rule for evaluation

double ts = y_h * ( y_h1 * ( 1.0 / 362880 ) ) - 1.0 / 5040 ;

double tc = y_h * ( y_h1 * ( 1.0 / 40320 ) ) - 1.0 / 720 ;

ts *= Y4 ;

tc *= Y4 ;

ts += Y2 * ( 1.0 / 120 ) - 1.0 / 6 ;

tc += Y2 * ( 1.0 / 24 ) - 1.0 / 2 ;

ts *= Y2 * ( y_h * cosa ) ;

tc *= Y2 * sina_h ;

/*

At this point,

ts = cos(a)*(yh^8*yh1/9! - yh^7/7! + yh^5/5! - yh^3/3!)

tc = sin(a)*(yh^7*yh1/8! - yh^6/6! + yh^4/4! - yh^2/2!)

The missing terms (sin(a) + yh*cos(a)) will be added later

(why yh^8*yh1 instead of yh^9?)

*/

double s = ts ;

s += ( ( y_l * ( y_h * sina_h - cosa ) + sina_l ) + syc_l1 ) + syc_l2 ; // low terms

/*

sin(a + yh - yl) up to the 2nd order term =

sin(a) + yh*cos(a) - yh^2*sin(a)/2 (high terms, included in syc_h and tc)

- yl*cos(a) + yh*yl*sin(a) (low terms)

- yl^2*sin(a)/2 (ignore)

Thus the term y_l*(y_h*sina_h - cosa)

*/

s += tc ;

s += syc_h ;

return s ;

}