Skip to content

Commit

Permalink
fix: fix C implementation and add dependencies
Browse files Browse the repository at this point in the history
  • Loading branch information
steff456 committed Jul 11, 2023
1 parent e4d06d5 commit e2ac1ed
Show file tree
Hide file tree
Showing 2 changed files with 38 additions and 29 deletions.
12 changes: 9 additions & 3 deletions lib/node_modules/@stdlib/math/base/ops/cdiv/manifest.json
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,9 @@
"@stdlib/math/base/napi/binary",
"@stdlib/complex/float64",
"@stdlib/complex/reim",
"@stdlib/math/base/special/abs"
"@stdlib/math/base/special/abs",
"@stdlib/constants/float64/max",
"@stdlib/constants/float64/eps"
]
},
{
Expand All @@ -55,7 +57,9 @@
"dependencies": [
"@stdlib/complex/float64",
"@stdlib/complex/reim",
"@stdlib/math/base/special/abs"
"@stdlib/math/base/special/abs",
"@stdlib/constants/float64/max",
"@stdlib/constants/float64/eps"
]
},
{
Expand All @@ -71,7 +75,9 @@
"dependencies": [
"@stdlib/complex/float64",
"@stdlib/complex/reim",
"@stdlib/math/base/special/abs"
"@stdlib/math/base/special/abs",
"@stdlib/constants/float64/max",
"@stdlib/constants/float64/eps"
]
}
]
Expand Down
55 changes: 29 additions & 26 deletions lib/node_modules/@stdlib/math/base/ops/cdiv/src/main.c
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,14 @@
#include "stdlib/math/base/special/abs.h"
#include "stdlib/complex/float64.h"
#include "stdlib/complex/reim.h"
#include "stdlib/constants/float64/max.h"
#include "stdlib/constants/float64/eps.h"
#include <stdint.h>

const double LARGE_THRESHOLD = __DBL_MAX__ * 0.5;
const double SMALL_THRESHOLD = __DBL_MIN__ * ( 2.0 / __DBL_EPSILON__ );
const double RECIP_EPS_SQR = 2.0 / ( __DBL_EPSILON__ * __DBL_EPSILON__ );
const double LARGE_THRESHOLD = STDLIB_CONSTANT_FLOAT64_MAX * 0.5;
// TODO: Change to stdlib min constant once implemented
const double SMALL_THRESHOLD = __DBL_MIN__ * ( 2.0 / STDLIB_CONSTANT_FLOAT64_EPS );
const double RECIP_EPS_SQR = 2.0 / ( STDLIB_CONSTANT_FLOAT64_EPS * STDLIB_CONSTANT_FLOAT64_EPS );

/**
* Computes the real part of the quotient.
Expand All @@ -44,8 +47,8 @@ const double RECIP_EPS_SQR = 2.0 / ( __DBL_EPSILON__ * __DBL_EPSILON__ );
* @returns real part of the quotient
*/
static double internalCompreal( const double re1, const double im1, const double re2, const double im2, const double r, const double t ) {
double br;
if ( r == 0.0 ) {
double br;
if ( r == 0.0 ) {
return ( re1 + (im2 * ( im1 / re2 )) ) * t;
}
br = im1 * r;
Expand All @@ -64,22 +67,22 @@ static double internalCompreal( const double re1, const double im1, const double
*
* [@baudin:2012]: https://arxiv.org/abs/1210.4539
*
* @private
* @param {number} re1 - real component
* @param {number} im1 - imaginary component
* @param {number} re2 - real component
* @param {number} im2 - imaginary component
* @returns {Array<number>} result
* @param {number} re - real result
* @param {number} im - imaginary result
*/
static void robustInternal( const double re1, const double im1, const double re2, const double im2, double* re, double* im ) {
double r;
double t;
double r;
double t;

r = im2 / re2;
r = im2 / re2;
t = 1.0 / ( re2 + ( im2 * r ) );

*re = internalCompreal( re1, im1, re2, im2, r, t );
*im = internalCompreal( im1, -re1, re2, im2, r, t );
*re = internalCompreal( re1, im1, re2, im2, r, t );
*im = internalCompreal( im1, -re1, re2, im2, r, t );
}

/**
Expand Down Expand Up @@ -110,29 +113,29 @@ stdlib_complex128_t stdlib_base_cdiv( const stdlib_complex128_t z1, const stdlib
double re2;
double im1;
double im2;
double ab;
double cd;
double ab;
double cd;
double re;
double im;
double s;
double s;

stdlib_reim( z1, &re1, &im1 );
stdlib_reim( z2, &re2, &im2 );

if ( stdlib_base_abs( re1 ) > stdlib_base_abs( im1 ) ){
ab = re1;
} else {
ab = im1;
}
if ( stdlib_base_abs( re2 ) > stdlib_base_abs( im2 ) ){
cd = re2;
} else {
cd = im2;
}
if ( stdlib_base_abs( re1 ) > stdlib_base_abs( im1 ) ){
ab = re1;
} else {
ab = im1;
}
if ( stdlib_base_abs( re2 ) > stdlib_base_abs( im2 ) ){
cd = re2;
} else {
cd = im2;
}

s = 1.0;

if ( ab >= LARGE_THRESHOLD ) {
if ( ab >= LARGE_THRESHOLD ) {
re1 *= 0.5;
im1 *= 0.5;
s *= 2.0;
Expand Down

0 comments on commit e2ac1ed

Please sign in to comment.