Subversion Repositories forwardcom

/*********************************  sin.as  ***********************************
* Author:        Agner Fog
* date created:  2018-03-29
* Last modified: 2021-04-25
* Version:       1.11
* Project:       ForwardCom library math.li
* Description:   sin, cos, and tan functions. Calculate in radians, double precision
*                The argument x can be a scalar or a vector
*                The return value will be a vector with the same length
* C declaration: double sin(double x);
* C declaration: double cos(double x);
* C declaration: double tan(double x);
* C declaration: struct {double s; double c;} sincos(double x);
*
* This code is adapted from C++ vector class library www.github.com/vectorclass
* Copyright 2018-2021 GNU General Public License http://www.gnu.org/licenses
*****************************************************************************/


// define constants
% M_2_PI = 0.636619772367581343076       // 2./pi
% P0sin  = -1.66666666666666307295E-1    // polynomial coefficients for sin
% P1sin  = 8.33333333332211858878E-3
% P2sin  = -1.98412698295895385996E-4
% P3sin  = 2.75573136213857245213E-6
% P4sin  = -2.50507477628578072866E-8
% P5sin  = 1.58962301576546568060E-10

% P0cos  = 4.16666666666665929218E-2     // polynomial coefficients for cos
% P1cos  = -1.38888888888730564116E-3
% P2cos  = 2.48015872888517045348E-5
% P3cos  = -2.75573141792967388112E-7
% P4cos  = 2.08757008419747316778E-9
% P5cos  = -1.13585365213876817300E-11

% DP1    = 7.853981554508209228515625E-1 // modulo pi/2 for extended precision
% DP2    = 7.94662735614792836714E-9     // correction for extended precision modular artithmetic
% DP3    = 3.06161699786838294307E-17    // correction for extended precision modular artithmetic


code section execute align = 4

public _sin:    function, reguse = 0, 0x7BF
public _cos:    function, reguse = 0, 0x7BF
public _sincos: function, reguse = 0, 0x7BF
public _tan:    function, reguse = 0, 0x7BF

// common entry for sin and sincos functions
_sin function
_sincos: 

/* registers:
  v0 = x
  v1 = abs(x)
  v1 = quadrant
  v10 = abs(x) reduced modulo pi/2
  v2 = v10^2
  v3 = v10^4
  v4 = v10^8
  v5 = v10^3
  v6 = unused (vacant flag for calling function)
  v7 = temp
  v8 = sin
  v9 = cos
*/

// Find quadrant:
//      0 -   pi/4 => 0
//   pi/4 - 3*pi/4 => 1
// 3*pi/4 - 5*pi/4 => 2
// 5*pi/4 - 7*pi/4 => 3
// 7*pi/4 - 8*pi/4 => 4

double v1 = clear_bit(v0, 63)                    // abs(x)
double v4 = v1 * M_2_PI
double v4 = round(v4, 0)   // round to integer

// reduce modulo pi/2, with extended precision
// x = ((xa - y * DP1) - y * DP2) - y * DP3;
double v10 = v4 * (-DP1*2) + v1
double v10 = v4 * (-DP2*2) + v10
double v10 = v4 * (-DP3*2) + v10
double v5 = !(v4 > ((1 << 51) + 0.0))              // check for loss of precision and overflow, but not NAN
double v1 = v4 + ((1 << 52) + 0.0)                 // add magic number 2^52 to get integer into lowest bit
double v10 = v5 ? v10 : 0                          // zero if out of range. result will be -1, 0, or 1

// Expansion of sin and cos, valid for -pi/4 <= x <= pi/4
double v2 = v10 * v10                              // x^2
double v3 = v2 * v2                              // x^4
double v4 = v3 * v3                              // x^8

// calculate polynomial P5sin*x2^5 + P4sin*x2^4 + P3sin*x2^3 + P2sin*x2^2 + P1sin*x2 + P0sin
// = (p2+p3*x2)*x4 + ((p4+p5*x2)*x8 + (p0+p1*x2));

double v5 = replace(v10, P0sin)                   // broadcast to same length as x
double v5 = v2 * P1sin + v5
double v7 = replace(v10, P4sin)
double v7 = v2 * P5sin + v7
double v8 = replace(v10, P2sin)
double v8 = v2 * P3sin + v8 
double v7 = v7 * v4 + v5
double v8 = v8 * v3 + v7

// calculate polynomial P5cos*x2^5 + P4cos*x2^4 + P3cos*x2^3 + P2cos*x2^2 + P1cos*x2 + P0cos
// = (p2+p3*x2)*x4 + ((p4+p5*x2)*x8 + (p0+p1*x2));
double v5 = replace(v10, P0cos)
double v5 = v2 * P1cos + v5
double v7 = replace(v10, P4cos)
double v7 = v2 * P5cos + v7
double v9 = replace(v10, P2cos)
double v9 = v2 * P3cos + v9
double v7 = v7 * v4 + v5
double v9 = v9 * v3 + v7

// s = x + (x * x2) * s;
double v5 = v10 * v2
double v8 = v8 * v5 + v10
// c = 1.0 - x2 * 0.5 + (x2 * x2) * c;
double v9 = v9 * v3
double v9 = v2 * (-0.5) + v9
double v9 = 1.0 + v9

// swap sin and cos if odd quadrant
double v3 = v1 ? v9 : v8        // sin
double v4 = v1 ? v8 : v9        // cos

// get sign of sin
int64  v5 = v1 << 62            // get bit 1 into sign bit, x modulo pi/2 = 2 or 3
int64  v5 ^= v0                 // toggle with sign of original x
int64  v5 = and(v5, 1 << 63)    // isolate sign bit
double v0 = v3 ^ v5             // apply sign bit to sin

// get sign of cos
int64  v1 = v1 + 1              // change sign when x modulo pi/2 = 1 or 2
int64  v1 = v1 << 62            // get bit 1 into sign bit
int64  v1 = and(v1, 1 << 63)    // isolate sign bit
double v1 = v4 ^ v1             // apply sign bit to cos

// return sin in v0, cos in v1
return
_sin end

// cosine function
_cos function
call _sincos
double v0 = v1                                   // cos is in v1
return
_cos end

// tangent function
_tan function
call _sincos
double v0 = v0 / v1                              // tan(x) = sin(x)/cos(x)
return
_tan end

code end