Subversion Repositories forwardcom

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


// define constants
% M_2_PI = 0.636619772367581343076       // 2./pi

% P0sinf = -1.6666654611E-1
% P1sinf = 8.3321608736E-3
% P2sinf = -1.9515295891E-4
% P0cosf = 4.166664568298827E-2
% P1cosf = -1.388731625493765E-3
% P2cosf = 2.443315711809948E-5

% DP1F = 0.78515625 * 2.
% DP2F = 2.4187564849853515625E-4 * 2.
% DP3F = 3.77489497744594108E-8 * 2.


code section execute align = 4

public _sinf:    function, reguse = 0, 0x1BF
public _cosf:    function, reguse = 0, 0x1BF
public _sincosf: function, reguse = 0, 0x1BF
public _tanf:    function, reguse = 0, 0x1BF

// common entry for sin and sincos functions
_sinf function
_sincosf: 

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

// 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

// reduce modulo pi/2, with extended precision
//nop
float v1 = clear_bit(v0, 31)                    // abs(x)
float v5 = v1 * M_2_PI
float v5 = round(v5, 0)   // round to integer

// x = ((xa - y * DP1) - y * DP2) - y * DP3;
float v8 = v5 * (-DP1F) + v1 
float v8 = v5 * (-DP2F) + v8
float v8 = v5 * (-DP3F) + v8

float v3 = !(v5 > ((1 << 22) + 0.0))            // check for loss of precision and overflow, but not NAN
float v1 = v5 + ((1 << 23) + 0.0)               // add magic number 2^23 to get integer into lowest bit
float v8 = v3 ? v8 : 0                          // zero if out of range. result will be -1, 0, or 1

// Taylor expansion of sin and cos, valid for -pi/4 <= x <= pi/4
// s = polynomial_2(x^2, P0sinf, P1sinf, P2sinf) * (x*x^2) + x;

float v2 = v8 * v8        // x^2
float v3 = v8 * v2        // x^3
float v4 = v2 * v2        // x^4
float v5 = replace(v8, P0sinf)                   // broadcast to same length as x
float v5 = v2 * P1sinf + v5
float v5 = v4 * P2sinf + v5
float v5 = v5 * v3 + v8                          // sin

// c = polynomial_2(x2, P0cosf, P1cosf, P2cosf) * (x2*x2) + nmul_add(0.5f, x2, 1.0f);
float v7 = replace(v8, P0cosf)                   // broadcast to same length as x
float v7 = v2 * P1cosf + v7
float v7 = v4 * P2cosf + v7
float v3 = replace(v8, 1.0)
float v3 = v2 * (-0.5) + v3                       // 1 - 0.5*x^2
float v7 = v7 * v4 + v3                           // cos

// swap sin and cos if odd quadrant
float v3 = v1 ? v7 : v5        // sin
float v4 = v1 ? v5 : v7        // cos

// get sign of sin
int32  v5 = v1 << 30            // get bit 1 into sign bit, x modulo pi/2 = 2 or 3
int32  v5 ^= v0                 // toggle with sign of original x
int32  v5 = and(v5, 1 << 31)    // isolate sign bit
float  v0 = v3 ^ v5             // apply sign bit to sin

// get sign of cos
int32  v1 = v1 + 1              // change sign when x modulo pi/2 = 1 or 2
int32  v1 = v1 << 30            // get bit 1 into sign bit
int32  v1 = and(v1, 1 << 31)    // isolate sign bit
float  v1 = v4 ^ v1             // apply sign bit to cos

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

// cosine function
_cosf function
call _sincosf
float v0 = v1                                   // cos is in v1
return
_cosf end

// tangent function
_tanf function
call _sincosf
float v0 = v0 / v1                              // tan(x) = sin(x)/cos(x)
return
_tanf end

code end