Skip to content

Commit

Permalink
repairing inexact complex division so it uses Smith's method to avoid…
Browse files Browse the repository at this point in the history
… overflow.
  • Loading branch information
Danny Yoo committed Dec 10, 2010
1 parent 87d64fc commit ce68f8f
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 14 deletions.
49 changes: 38 additions & 11 deletions src/js-numbers.js
Original file line number Original file line Diff line number Diff line change
Expand Up @@ -2165,32 +2165,59 @@ if (typeof(exports) !== 'undefined') {
return Complex.makeInstance(r, i); return Complex.makeInstance(r, i);
}; };






Complex.prototype.divide = function(other){ Complex.prototype.divide = function(other){
var a, b, c, d, r, x, y;
// If the other value is real, just do primitive division // If the other value is real, just do primitive division
if (other.isReal()) { if (other.isReal()) {
return Complex.makeInstance( return Complex.makeInstance(
divide(this.r, other.r), divide(this.r, other.r),
divide(this.i, other.r)); divide(this.i, other.r));
} }


if (this.isInexact() || other.isInexact()) {
// http://portal.acm.org/citation.cfm?id=1039814
// We currently use Smith's method, though we should
// probably switch over to Priest's method.
a = this.r;
b = this.i;
c = other.r;
d = other.i;
if (lessThanOrEqual(abs(d), abs(c))) {
r = divide(d, c);
x = divide(add(a, multiply(b, r)),
add(c, multiply(d, r)));
y = divide(subtract(b, multiply(a, r)),
add(c, multiply(d, r)));
} else {
r = divide(c, d);
x = divide(add(multiply(a, r), b),
add(multiply(c, r), d));
y = divide(subtract(multiply(b, r), a),
add(multiply(c, r), d));
}
return makeComplex(x, y);
} else {
var con = conjugate(other);
var up = multiply(this, con);


var con = conjugate(other); // Down is guaranteed to be real by this point.
var up = multiply(this, con); var down = realPart(multiply(other, con));

// Down is guaranteed to be real by this point.
var down = realPart(multiply(other, con));


var result = Complex.makeInstance( var result = Complex.makeInstance(
divide(realPart(up), down), divide(realPart(up), down),
divide(imaginaryPart(up), down)); divide(imaginaryPart(up), down));
return result; return result;
}
}; };


Complex.prototype.conjugate = function(){ Complex.prototype.conjugate = function(){
var result = Complex.makeInstance( var result = Complex.makeInstance(
this.r, this.r,
subtract(0, subtract(0, this.i));
this.i));


return result; return result;
}; };
Expand Down
12 changes: 9 additions & 3 deletions test/tests.js
Original file line number Original file line Diff line number Diff line change
Expand Up @@ -87,7 +87,9 @@ describe('rational constructions', {


describe('complex construction', { describe('complex construction', {
'polar' : function() { 'polar' : function() {
// FIXME: add tests for polar construction assertTrue(eqv(makeComplexPolar(1, 2),
makeComplex(makeFloat(-0.4161468365471424),
makeFloat(0.9092974268256817))));
}, },
'non-real inputs should raise errors' : function() { 'non-real inputs should raise errors' : function() {
// FIXME: add tests for polar construction // FIXME: add tests for polar construction
Expand Down Expand Up @@ -1714,8 +1716,12 @@ describe('divide', {
makeFloat(1e300)), makeFloat(1e300)),
makeComplex(makeFloat(4e300), makeComplex(makeFloat(4e300),
makeFloat(4e300))), makeFloat(4e300))),
makeFloat(.25))) makeComplex(makeFloat(.25), makeFloat(0.0))));
// FIXME: we're missing this assertTrue(eqv(divide(makeComplex(2, 6),
makeComplex(4, 1)),
makeComplex(makeRational(14, 17),
makeRational(22, 17))));

}, },


'division by zeros' : function() { 'division by zeros' : function() {
Expand Down

0 comments on commit ce68f8f

Please sign in to comment.