Skip to content

Commit

Permalink
Fix complex sqrt
Browse files Browse the repository at this point in the history
  • Loading branch information
@paulmasson
paulmasson committed Jun 28, 2021
1 parent d32482b commit 078de16
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 8 deletions.
18 changes: 14 additions & 4 deletions build/math.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -377,7 +377,11 @@ function sqrt( x ) {

if ( isArbitrary(x) ) {

if ( x.im === 0n ) return { re: sqrt(x.re), im: 0n };
if ( x.im === 0n )
if ( x.re < 0n ) return { re: 0n, im: sqrt(-x.re) };
else return { re: sqrt(x.re), im: 0n };

// need evaluation independent of natural logarithm

var c = abs(x);
var sign = x.im < 0n ? -1n : 1n;
Expand All @@ -386,10 +390,16 @@ function sqrt( x ) {

}

var c = Math.hypot( x.re, x.im );
var sign = x.im === 0 ? 1 : Math.sign(x.im);
if ( x.im === 0 )
if ( x.re < 0 ) return { re: 0, im: Math.sqrt(-x.re) };
else return { re: Math.sqrt(x.re), im: 0 };

if ( x.re === 0 ) return { re: 0, im: 0 };

// expression above suffers from rounding errors when not arbitrary,
// especially affecting elliptic integral of third kind

return { re: Math.sqrt( (c + x.re)/2 ), im: sign * Math.sqrt( (c - x.re)/2 ) };
return exp( mul( .5, log(x) ) );

}

Expand Down
18 changes: 14 additions & 4 deletions src/operators.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -344,7 +344,11 @@ function sqrt( x ) {

if ( isArbitrary(x) ) {

if ( x.im === 0n ) return { re: sqrt(x.re), im: 0n };
if ( x.im === 0n )
if ( x.re < 0n ) return { re: 0n, im: sqrt(-x.re) };
else return { re: sqrt(x.re), im: 0n };

// need evaluation independent of natural logarithm

var c = abs(x);
var sign = x.im < 0n ? -1n : 1n;
Expand All @@ -353,10 +357,16 @@ function sqrt( x ) {

}

var c = Math.hypot( x.re, x.im );
var sign = x.im === 0 ? 1 : Math.sign(x.im);
if ( x.im === 0 )
if ( x.re < 0 ) return { re: 0, im: Math.sqrt(-x.re) };
else return { re: Math.sqrt(x.re), im: 0 };

if ( x.re === 0 ) return { re: 0, im: 0 };

// expression above suffers from rounding errors when not arbitrary,
// especially affecting elliptic integral of third kind

return { re: Math.sqrt( (c + x.re)/2 ), im: sign * Math.sqrt( (c - x.re)/2 ) };
return exp( mul( .5, log(x) ) );

}

Expand Down

1 comment on commit 078de16

@paulmasson
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@axkr if you are implementing elliptic integrals then this commit is crucial. The form of the square root for arbitrary precision develops rounding errors at machine precision which particularly affect the integral of the third kind. Took me a few days to figure out this form of the square root is more accurate for machine precision.

Please sign in to comment.