@stdlib/math-base-special-rempio2
Advanced tools
Comparing version 0.2.1 to 0.3.0
"use strict";var m=function(x,r){return function(){return r||x((r={exports:{}}).exports,r),r.exports}};var g=m(function(c0,S){ | ||
var L=require('@stdlib/math-base-special-floor/dist'),T=require('@stdlib/math-base-special-ldexp/dist'),W=require('@stdlib/array-base-zeros/dist'),w=[10680707,7228996,1387004,2578385,16069853,12639074,9804092,4427841,16666979,11263675,12935607,2387514,4345298,14681673,3074569,13734428,16653803,1880361,10960616,8533493,3062596,8710556,7349940,6258241,3772886,3769171,3798172,8675211,12450088,3874808,9961438,366607,15675153,9132554,7151469,3571407,2607881,12013382,4155038,6285869,7677882,13102053,15825725,473591,9065106,15363067,6271263,9264392,5636912,4652155,7056368,13614112,10155062,1944035,9527646,15080200,6658437,6231200,6832269,16767104,5075751,3212806,1398474,7579849,6349435,12618859],Y=[1.570796251296997,7549789415861596e-23,5390302529957765e-30,3282003415807913e-37,1270655753080676e-44,12293330898111133e-52,27337005381646456e-60,21674168387780482e-67],G=16777216,h=5960464477539063e-23,c=W(20),b=W(20),N=W(20),F=W(20);function M(x,r,a,_,u,f,t,E,n){var v,i,o,A,e,D,P,O,I;for(A=f,I=_[a],O=a,e=0;O>0;e++)i=h*I|0,F[e]=I-G*i|0,I=_[O-1]+i,O-=1;if(I=T(I,u),I-=8*L(I*.125),P=I|0,I-=P,o=0,u>0?(e=F[a-1]>>24-u,P+=e,F[a-1]-=e<<24-u,o=F[a-1]>>23-u):u===0?o=F[a-1]>>23:I>=.5&&(o=2),o>0){for(P+=1,v=0,e=0;e<a;e++)O=F[e],v===0?O!==0&&(v=1,F[e]=16777216-O):F[e]=16777215-O;if(u>0)switch(u){case 1:F[a-1]&=8388607;break;case 2:F[a-1]&=4194303;break}o===2&&(I=1-I,v!==0&&(I-=T(1,u)))}if(I===0){for(O=0,e=a-1;e>=f;e--)O|=F[e];if(O===0){for(D=1;F[f-D]===0;D++);for(e=a+1;e<=a+D;e++){for(n[E+e]=w[t+e],i=0,O=0;O<=E;O++)i+=x[O]*n[E+(e-O)];_[e]=i}return a+=D,M(x,r,a,_,u,f,t,E,n)}}if(I===0)for(a-=1,u-=24;F[a]===0;)a-=1,u-=24;else I=T(I,-u),I>=G?(i=h*I|0,F[a]=I-G*i|0,a+=1,u+=24,F[a]=i):F[a]=I|0;for(i=T(1,u),e=a;e>=0;e--)_[e]=i*F[e],i*=h;for(e=a;e>=0;e--){for(i=0,D=0;D<=A&&D<=a-e;D++)i+=Y[D]*_[e+D];N[a-e]=i}for(i=0,e=a;e>=0;e--)i+=N[e];for(o===0?r[0]=i:r[0]=-i,i=N[0]-i,e=1;e<=a;e++)i+=N[e];return o===0?r[1]=i:r[1]=-i,P&7}function Z(x,r,a,_){var u,f,t,E,n,v,i,o,A;for(f=4,E=_-1,t=(a-3)/24|0,t<0&&(t=0),v=a-24*(t+1),o=t-E,A=E+f,i=0;i<=A;i++)o<0?c[i]=0:c[i]=w[o],o+=1;for(i=0;i<=f;i++){for(u=0,o=0;o<=E;o++)u+=x[o]*c[E+(i-o)];b[i]=u}return n=f,M(x,r,n,b,v,f,t,E,c)}S.exports=Z | ||
var L=require('@stdlib/math-base-special-floor/dist'),T=require('@stdlib/math-base-special-ldexp/dist'),W=require('@stdlib/array-base-zeros/dist'),w=[10680707,7228996,1387004,2578385,16069853,12639074,9804092,4427841,16666979,11263675,12935607,2387514,4345298,14681673,3074569,13734428,16653803,1880361,10960616,8533493,3062596,8710556,7349940,6258241,3772886,3769171,3798172,8675211,12450088,3874808,9961438,366607,15675153,9132554,7151469,3571407,2607881,12013382,4155038,6285869,7677882,13102053,15825725,473591,9065106,15363067,6271263,9264392,5636912,4652155,7056368,13614112,10155062,1944035,9527646,15080200,6658437,6231200,6832269,16767104,5075751,3212806,1398474,7579849,6349435,12618859],Y=[1.570796251296997,7549789415861596e-23,5390302529957765e-30,3282003415807913e-37,1270655753080676e-44,12293330898111133e-52,27337005381646456e-60,21674168387780482e-67],G=16777216,h=5960464477539063e-23,c=W(20),b=W(20),N=W(20),F=W(20);function M(x,r,a,_,u,f,t,E,n){var v,i,o,A,e,D,P,O,I;for(A=f,I=_[a],O=a,e=0;O>0;e++)i=h*I|0,F[e]=I-G*i|0,I=_[O-1]+i,O-=1;if(I=T(I,u),I-=8*L(I*.125),P=I|0,I-=P,o=0,u>0?(e=F[a-1]>>24-u,P+=e,F[a-1]-=e<<24-u,o=F[a-1]>>23-u):u===0?o=F[a-1]>>23:I>=.5&&(o=2),o>0){for(P+=1,v=0,e=0;e<a;e++)O=F[e],v===0?O!==0&&(v=1,F[e]=16777216-O):F[e]=16777215-O;if(u>0)switch(u){case 1:F[a-1]&=8388607;break;case 2:F[a-1]&=4194303;break}o===2&&(I=1-I,v!==0&&(I-=T(1,u)))}if(I===0){for(O=0,e=a-1;e>=f;e--)O|=F[e];if(O===0){for(D=1;F[f-D]===0;D++);for(e=a+1;e<=a+D;e++){for(n[E+e]=w[t+e],i=0,O=0;O<=E;O++)i+=x[O]*n[E+(e-O)];_[e]=i}return a+=D,M(x,r,a,_,u,f,t,E,n)}for(a-=1,u-=24;F[a]===0;)a-=1,u-=24}else I=T(I,-u),I>=G?(i=h*I|0,F[a]=I-G*i|0,a+=1,u+=24,F[a]=i):F[a]=I|0;for(i=T(1,u),e=a;e>=0;e--)_[e]=i*F[e],i*=h;for(e=a;e>=0;e--){for(i=0,D=0;D<=A&&D<=a-e;D++)i+=Y[D]*_[e+D];N[a-e]=i}for(i=0,e=a;e>=0;e--)i+=N[e];for(o===0?r[0]=i:r[0]=-i,i=N[0]-i,e=1;e<=a;e++)i+=N[e];return o===0?r[1]=i:r[1]=-i,P&7}function Z(x,r,a,_){var u,f,t,E,n,v,i,o,A;for(f=4,E=_-1,t=(a-3)/24|0,t<0&&(t=0),v=a-24*(t+1),o=t-E,A=E+f,i=0;i<=A;i++)o<0?c[i]=0:c[i]=w[o],o+=1;for(i=0;i<=f;i++){for(u=0,o=0;o<=E;o++)u+=x[o]*c[E+(i-o)];b[i]=u}return n=f,M(x,r,n,b,v,f,t,E,c)}S.exports=Z | ||
});var X=m(function(N0,V){ | ||
@@ -4,0 +4,0 @@ var J=require('@stdlib/math-base-special-round/dist'),K=require('@stdlib/number-float64-base-get-high-word/dist'),$=.6366197723675814,q=1.5707963267341256,z=6077100506506192e-26,j=6077100506303966e-26,y=20222662487959506e-37,r0=20222662487111665e-37,e0=84784276603689e-45,Q=2047;function a0(x,r,a){var _,u,f,t,E,n,v;return u=J(x*$),t=x-u*q,E=u*z,v=r>>20|0,a[0]=t-E,_=K(a[0]),n=v-(_>>20&Q),n>16&&(f=t,E=u*j,t=f-E,E=u*y-(f-t-E),a[0]=t-E,_=K(a[0]),n=v-(_>>20&Q),n>49&&(f=t,E=u*r0,t=f-E,E=u*e0-(f-t-E),a[0]=t-E)),a[1]=t-a[0]-E,u}V.exports=a0 |
@@ -209,5 +209,3 @@ /** | ||
} | ||
} | ||
// Chop off zero terms... | ||
if ( z === 0.0 ) { | ||
// Chop off zero terms... | ||
jz -= 1; | ||
@@ -276,4 +274,102 @@ q0 -= 24; | ||
* | ||
* - The method is to compute the integer (`mod 8`) and fraction parts of `2x/π` without doing the full multiplication. In general, we skip the part of the product that is known to be a huge integer (more accurately, equals `0 mod 8` ). Thus, the number of operations is independent of the exponent of the input. | ||
* - The method is to compute the integer (mod 8) and fraction parts of (2/π) * x without doing the full multiplication. In general, we skip the part of the product that are known to be a huge integer (more accurately, = 0 mod 8 ). Thus the number of operations are independent of the exponent of the input. | ||
* | ||
* - (2/π) is represented by an array of 24-bit integers in `ipio2[]`. | ||
* | ||
* - Input parameters: | ||
* | ||
* - `x[]` The input value (must be positive) is broken into `nx` pieces of 24-bit integers in double precision format. `x[i]` will be the i-th 24 bit of x. The scaled exponent of `x[0]` is given in input parameter `e0` (i.e., `x[0]*2^e0` match x's up to 24 bits). | ||
* | ||
* Example of breaking a double positive `z` into `x[0]+x[1]+x[2]`: | ||
* | ||
* ```tex | ||
* e0 = \mathrm{ilogb}(z) - 23 | ||
* z = \mathrm{scalbn}(z, -e0) | ||
* ``` | ||
* | ||
* for `i = 0,1,2` | ||
* | ||
* ```tex | ||
* x[i] = \lfloor z \rfloor | ||
* z = (z - x[i]) \times 2^{24} | ||
* ``` | ||
* | ||
* - `y[]` output result in an array of double precision numbers. | ||
* | ||
* The dimension of `y[]` is: | ||
* 24-bit precision 1 | ||
* 53-bit precision 2 | ||
* 64-bit precision 2 | ||
* 113-bit precision 3 | ||
* | ||
* The actual value is the sum of them. Thus, for 113-bit precision, one may have to do something like: | ||
* | ||
* ```tex | ||
* \mathrm{long\ double} \: t, w, r_{\text{head}}, r_{\text{tail}}; \\ | ||
* t &= (\mathrm{long\ double}) y[2] + (\mathrm{long\ double}) y[1]; \\ | ||
* w &= (\mathrm{long\ double}) y[0]; \\ | ||
* r_{\text{head}} &= t + w; \\ | ||
* r_{\text{tail}} &= w - (r_{\text{head}} - t); | ||
* ``` | ||
* | ||
* - `e0` The exponent of `x[0]`. Must be <= 16360 or you need to expand the `ipio2` table. | ||
* | ||
* - `nx` dimension of `x[]` | ||
* | ||
* - `prec` an integer indicating the precision: | ||
* 0 24 bits (single) | ||
* 1 53 bits (double) | ||
* 2 64 bits (extended) | ||
* 3 113 bits (quad) | ||
* | ||
* - External function: | ||
* | ||
* - double `scalbn()`, `floor()`; | ||
* | ||
* - Here is the description of some local variables: | ||
* | ||
* - `jk` `jk+1` is the initial number of terms of `ipio2[]` needed in the computation. The minimum and recommended value for `jk` is 3,4,4,6 for single, double, extended, and quad. `jk+1` must be 2 larger than you might expect so that our recomputation test works. (Up to 24 bits in the integer part (the 24 bits of it that we compute) and 23 bits in the fraction part may be lost to cancellation before we recompute.) | ||
* | ||
* - `jz` local integer variable indicating the number of terms of `ipio2[]` used. | ||
* | ||
* - `jx` `nx - 1` | ||
* | ||
* - `jv` index for pointing to the suitable `ipio2[]` for the computation. In general, we want | ||
* | ||
* ```tex | ||
* \frac{{2^{e0} \cdot x[0] \cdot \mathrm{ipio2}[jv-1] \cdot 2^{-24jv}}}{{8}} | ||
* ``` | ||
* | ||
* to be an integer. Thus | ||
* | ||
* ```tex | ||
* e0 - 3 - 24 \cdot jv \geq 0 \quad \text{or} \quad \frac{{e0 - 3}}{{24}} \geq jv | ||
* ``` | ||
* | ||
* Hence | ||
* | ||
* ```tex | ||
* jv = \max(0, \frac{{e0 - 3}}{{24}}) | ||
* ``` | ||
* | ||
* - `jp` `jp+1` is the number of terms in `PIo2[]` needed, `jp = jk`. | ||
* | ||
* - `q[]` double array with integral value, representing the 24-bits chunk of the product of `x` and `2/π`. | ||
* | ||
* - `q0` the corresponding exponent of `q[0]`. Note that the exponent for `q[i]` would be `q0-24*i`. | ||
* | ||
* - `PIo2[]` double precision array, obtained by cutting `π/2` into 24 bits chunks. | ||
* | ||
* - `f[]` `ipso2[]` in floating point | ||
* | ||
* - `iq[]` integer array by breaking up `q[]` in 24-bits chunk. | ||
* | ||
* - `fq[]` final product of `x*(2/π)` in `fq[0],..,fq[jk]` | ||
* | ||
* - `ih` integer. If >0 it indicates `q[]` is >= 0.5, hence it also indicates the _sign_ of the result. | ||
* | ||
* - Constants: | ||
* | ||
* - The hexadecimal values are the intended ones for the following constants. The decimal values may be used, provided that the compiler will convert from decimal to binary accurately enough to produce the hexadecimal values shown. | ||
* | ||
* @private | ||
@@ -280,0 +376,0 @@ * @param {PositiveNumber} x - input value |
{ | ||
"name": "@stdlib/math-base-special-rempio2", | ||
"version": "0.2.1", | ||
"version": "0.3.0", | ||
"description": "Compute `x - nπ/2 = r`.", | ||
@@ -17,5 +17,8 @@ "license": "Apache-2.0", | ||
"main": "./lib", | ||
"gypfile": false, | ||
"directories": { | ||
"doc": "./docs", | ||
"include": "./include", | ||
"lib": "./lib", | ||
"src": "./src", | ||
"dist": "./dist" | ||
@@ -35,11 +38,17 @@ }, | ||
"@stdlib/array-base-zeros": "^0.2.1", | ||
"@stdlib/constants-float64-high-word-abs-mask": "^0.2.1", | ||
"@stdlib/constants-float64-high-word-exponent-mask": "^0.2.1", | ||
"@stdlib/constants-float64-high-word-significand-mask": "^0.2.1", | ||
"@stdlib/math-base-special-floor": "^0.2.1", | ||
"@stdlib/math-base-special-ldexp": "^0.2.1", | ||
"@stdlib/math-base-special-round": "^0.2.1", | ||
"@stdlib/number-float64-base-from-words": "^0.2.1", | ||
"@stdlib/number-float64-base-get-high-word": "^0.2.1", | ||
"@stdlib/number-float64-base-get-low-word": "^0.2.1" | ||
"@stdlib/constants-float64-high-word-abs-mask": "^0.2.2", | ||
"@stdlib/constants-float64-high-word-exponent-mask": "^0.2.2", | ||
"@stdlib/constants-float64-high-word-significand-mask": "^0.2.2", | ||
"@stdlib/math-base-special-floor": "^0.2.3", | ||
"@stdlib/math-base-special-ldexp": "^0.2.3", | ||
"@stdlib/math-base-special-round": "^0.3.0", | ||
"@stdlib/napi-argv": "^0.2.2", | ||
"@stdlib/napi-argv-double": "^0.2.1", | ||
"@stdlib/napi-argv-float64array": "^0.2.2", | ||
"@stdlib/napi-create-double": "^0.0.2", | ||
"@stdlib/napi-export": "^0.2.2", | ||
"@stdlib/number-float64-base-from-words": "^0.2.2", | ||
"@stdlib/number-float64-base-get-high-word": "^0.2.2", | ||
"@stdlib/number-float64-base-get-low-word": "^0.2.2", | ||
"@stdlib/utils-library-manifest": "^0.2.2" | ||
}, | ||
@@ -46,0 +55,0 @@ "devDependencies": {}, |
108
README.md
@@ -58,3 +58,3 @@ <!-- | ||
Computes `x - nπ/2 = r`. The function returns `n` and stores the remainder `r` as two numbers in `y`, such that `y[0]+y[1] = r`. | ||
Computes `x - nπ/2 = r`. | ||
@@ -73,3 +73,3 @@ ```javascript | ||
When `x` is `NaN` or infinite, the function returns zero and sets the elements of `y` to `NaN`. | ||
When `x` is `NaN` or infinite, the function returns `0` and sets the elements of `y` to `NaN`. | ||
@@ -108,2 +108,3 @@ ```javascript | ||
- The function returns `n` and stores the remainder `r` as two numbers in `y`, such that `y[0]+y[1] = r`. | ||
- For input values larger than `2^20*π/2` in magnitude, the function **only** returns the last three binary digits of `n` and not the full result. | ||
@@ -140,2 +141,101 @@ | ||
<!-- C interface documentation. --> | ||
* * * | ||
<section class="c"> | ||
## C APIs | ||
<!-- Section to include introductory text. Make sure to keep an empty line after the intro `section` element and another before the `/section` close. --> | ||
<section class="intro"> | ||
</section> | ||
<!-- /.intro --> | ||
<!-- C usage documentation. --> | ||
<section class="usage"> | ||
### Usage | ||
```c | ||
#include "stdlib/math/base/special/rempio2.h" | ||
``` | ||
#### stdlib_base_rempio2( x, &rem1, &rem2 ) | ||
Computes `x - nπ/2 = r`. | ||
```c | ||
#include <stdint.h> | ||
double rem1; | ||
double rem2; | ||
int32_t n = stdlib_base_rempio2( 4.0, &rem1, &rem2 ); | ||
``` | ||
The function accepts the following arguments: | ||
- **x**: `[in] double` input value. | ||
- **rem1**: `[out] double*` destination for first remainder number. | ||
- **rem2**: `[out] double*` destination for second remainder number. | ||
```c | ||
int32_t stdlib_base_rempio2( const double x, double *rem1, double *rem2 ); | ||
``` | ||
</section> | ||
<!-- /.usage --> | ||
<!-- C API usage notes. Make sure to keep an empty line after the `section` element and another before the `/section` close. --> | ||
<section class="notes"> | ||
### Notes | ||
- The function returns `n` and stores the remainder `r` as two numbers in `rem1` and `rem2`, respectively, such that `rem1+rem2 = r`. | ||
</section> | ||
<!-- /.notes --> | ||
<!-- C API usage examples. --> | ||
<section class="examples"> | ||
### Examples | ||
```c | ||
#include "stdlib/math/base/special/rempio2.h" | ||
#include <stdio.h> | ||
#include <stdint.h> | ||
#include <inttypes.h> | ||
int main( void ) { | ||
const double x[] = { 0.0, 1.0, 4.0, 128.0 }; | ||
double rem1; | ||
double rem2; | ||
int32_t n; | ||
int i; | ||
for ( i = 0; i < 4; i++ ) { | ||
n = stdlib_base_rempio2( x[ i ], &rem1, &rem2 ); | ||
printf( "%lf - %"PRId32"π/2 = %lf + %lf\n", x[ i ], n, rem1, rem2 ); | ||
} | ||
} | ||
``` | ||
</section> | ||
<!-- /.examples --> | ||
</section> | ||
<!-- /.c --> | ||
<!-- Section for related `stdlib` packages. Do not manually edit this section, as it is automatically populated. --> | ||
@@ -183,4 +283,4 @@ | ||
[test-image]: https://github.com/stdlib-js/math-base-special-rempio2/actions/workflows/test.yml/badge.svg?branch=v0.2.1 | ||
[test-url]: https://github.com/stdlib-js/math-base-special-rempio2/actions/workflows/test.yml?query=branch:v0.2.1 | ||
[test-image]: https://github.com/stdlib-js/math-base-special-rempio2/actions/workflows/test.yml/badge.svg?branch=v0.3.0 | ||
[test-url]: https://github.com/stdlib-js/math-base-special-rempio2/actions/workflows/test.yml?query=branch:v0.3.0 | ||
@@ -187,0 +287,0 @@ [coverage-image]: https://img.shields.io/codecov/c/github/stdlib-js/math-base-special-rempio2/main.svg |
Sorry, the diff of this file is not supported yet
License Policy Violation
LicenseThis package is not allowed per your license policy. Review the package's license to ensure compliance.
Found 1 instance in 1 package
License Policy Violation
LicenseThis package is not allowed per your license policy. Review the package's license to ensure compliance.
Found 1 instance in 1 package
110041
18
1004
313
16
+ Added@stdlib/napi-argv@^0.2.2
+ Added@stdlib/napi-export@^0.2.2
+ Added@stdlib/assert-napi-equal-typedarray-types@0.2.2(transitive)
+ Added@stdlib/assert-napi-equal-types@0.2.2(transitive)
+ Added@stdlib/assert-napi-is-type@0.2.2(transitive)
+ Added@stdlib/assert-napi-is-typedarray@0.2.2(transitive)
+ Added@stdlib/assert-napi-status-ok@0.2.2(transitive)
+ Added@stdlib/math-base-assert-is-negative-zero@0.2.2(transitive)
+ Added@stdlib/math-base-special-round@0.3.0(transitive)
+ Added@stdlib/napi-argv@0.2.2(transitive)
+ Added@stdlib/napi-argv-double@0.2.1(transitive)
+ Added@stdlib/napi-argv-float64array@0.2.2(transitive)
+ Added@stdlib/napi-create-double@0.0.2(transitive)
+ Added@stdlib/napi-export@0.2.2(transitive)
- Removed@stdlib/math-base-special-round@0.2.1(transitive)
Updated@stdlib/constants-float64-high-word-exponent-mask@^0.2.2
Updated@stdlib/constants-float64-high-word-significand-mask@^0.2.2