Skip to content

Commit

Permalink
Add complex sqrt thresholds
Browse files Browse the repository at this point in the history 
Update complex test output
  • Loading branch information
Jeremy Clifford committed Apr 21, 2012
1 parent 6a1fb71 commit b6708e1
Show file tree
Hide file tree
Showing 3 changed files with 76 additions and 71 deletions.
23 changes: 14 additions & 9 deletions lib-clay/math/native/sqrt.clay
View file
Expand Up @@ -166,12 +166,17 @@ overload sqrt(x:Float) {
[I when Imaginary?(I) ]
inline overload sqrt(z:I) = sqrt(Complex(z));


alias thresh(#LongDouble) = 0xd.413cccfe77997ffp+16379fl;
alias overload thresh(#Float) = 0x1.a82797f7a7599p+126f;
alias overload thresh(#Double) = 0x1.a827999fcef32p+1022;

[C when Complex?(C)]
overload sqrt(z:C) --> result:C {
alias T = ComplexBaseType(C);
alias ZERO = T(0);
alias QUARTER = T(1/4);
alias HALF = T(1/2);
alias QUARTER = T(0.25);
alias HALF = T(0.5);
alias TWO = T(2);

var a = real(z);
Expand Down Expand Up @@ -204,23 +209,23 @@ overload sqrt(z:C) --> result:C {
else {
var scale = 0;
if (abs(a) >= thresh(T) or abs(b) >= thresh(T)) {
a *: 0.25;
b *: 0.25;
a *: QUARTER;
b *: QUARTER;
scale = 1;
}

// Algorithm 312, CACM vol 10, Oct 1967.
if (a >= ZERO) {
var t = sqrt((a + hypot(a, b)) * 0.5);
result <-- Complex(t, b / (2.0 * t));
var t = sqrt((a + hypot(a, b)) * HALF);
result <-- Complex(t, b / (TWO * t));
}
else {
var t = sqrt((-a + hypot(a, b)) * 0.5);
result <-- Complex(abs(b) / (2.0 * t), copysign(t, b));
var t = sqrt((-a + hypot(a, b)) * HALF);
result <-- Complex(abs(b) / (TWO * t), copysign(t, b));
}

// Rescale.
if (scale == 1)
result *: 2.0;
result *: TWO;
}
}
24 changes: 12 additions & 12 deletions lib-clay/math/native/sqrt.x86.clay
View file
Expand Up @@ -18,16 +18,16 @@ inline overload sqrt(x:Double) = x86_sqrtsd(Vec[Double,2](x,0.))[0];
inline overload sqrt(z:I) = sqrt(Complex(z));


alias thresh(#Double) = 0x1.a827999fcef32p+1022;
// alias overload thresh(#LongDouble) = 0x1.a827999fcef32p+1022;
// alias overload thresh(#Float) = 0x1.a827999fcef32p+1022;
alias thresh(#LongDouble) = 0xd.413cccfe77997ffp+16379fl;
alias overload thresh(#Float) = 0x1.a82797f7a7599p+126f;
alias overload thresh(#Double) = 0x1.a827999fcef32p+1022;

[C when Complex?(C)]
overload sqrt(z:C) --> result:C {
alias T = ComplexBaseType(C);
alias ZERO = T(0);
alias QUARTER = T(1/4);
alias HALF = T(1/2);
alias QUARTER = T(0.25);
alias HALF = T(0.5);
alias TWO = T(2);

var a = real(z);
Expand Down Expand Up @@ -60,23 +60,23 @@ overload sqrt(z:C) --> result:C {
else {
var scale = 0;
if (abs(a) >= thresh(T) or abs(b) >= thresh(T)) {
a *: 0.25;
b *: 0.25;
a *: QUARTER;
b *: QUARTER;
scale = 1;
}

// Algorithm 312, CACM vol 10, Oct 1967.
if (a >= ZERO) {
var t = sqrt((a + hypot(a, b)) * 0.5);
result <-- Complex(t, b / (2.0 * t));
var t = sqrt((a + hypot(a, b)) * HALF);
result <-- Complex(t, b / (TWO * t));
}
else {
var t = sqrt((-a + hypot(a, b)) * 0.5);
result <-- Complex(abs(b) / (2.0 * t), copysign(t, b));
var t = sqrt((-a + hypot(a, b)) * HALF);
result <-- Complex(abs(b) / (TWO * t), copysign(t, b));
}

// Rescale.
if (scale == 1)
result *: 2.0;
result *: TWO;
}
}
100 changes: 50 additions & 50 deletions test/math/complex/out.txt
View file
Expand Up @@ -18,7 +18,7 @@ sin(nan) = nan

cos(0+0j) = 1-0j
cos(-0+0j) = 1+0j
cos(-1-1j) = 0.8337300251311491-0.9888977057628651j
cos(-1-1j) = 0.8337300251311492-0.9888977057628651j
cos(2+2j) = -1.565625835315744-3.297894836311237j
cos(nan) = nan
cos(nan) = nan
Expand All @@ -44,7 +44,7 @@ tan(1+infj) = 0+1j
tan(nan) = nan
tan(-inf+1j) = nan
tan(inf+1j) = nan
tan(-inf+infj) = 0+1j
tan(-inf+infj) = -0+1j
tan(inf+infj) = 0+1j
tan(nan) = nan
tan(nan) = nan
Expand All @@ -54,16 +54,16 @@ tan(nan) = nan

asin(0+0j) = 0+0j
asin(-0+0j) = -0+0j
asin(-1-1j) = -0.09384622049358324-0.825079184917074j
asin(2+2j) = 1.53764515109662-1.356160455761411j
asin(-1-1j) = -0.6662394324925153-1.061275061905036j
asin(2+2j) = 0.7542491446980463+1.734324521487968j
asin(nan) = nan
asin(nan) = nan
asin(1+infj) = nan
asin(1+infj) = 0+infj
asin(nan) = nan
asin(-inf+1j) = nan
asin(inf+1j) = nan
asin(-inf+infj) = nan
asin(inf+infj) = nan
asin(-inf+1j) = -1.570796326794897+infj
asin(inf+1j) = 1.570796326794897+infj
asin(-inf+infj) = -0.7853981633974483+infj
asin(inf+infj) = 0.7853981633974483+infj
asin(nan) = nan
asin(nan) = nan
asin(nan) = nan
Expand All @@ -72,16 +72,16 @@ asin(nan) = nan

acos(0+0j) = 1.570796326794897-0j
acos(-0+0j) = 1.570796326794897-0j
acos(-1-1j) = 1.66464254728848+0.825079184917074j
acos(2+2j) = 0.03315117569827697+1.356160455761411j
acos(-1-1j) = 2.237035759287412+1.061275061905036j
acos(2+2j) = 0.8165471820968503-1.734324521487968j
acos(nan) = nan
acos(nan) = nan
acos(1+infj) = nan
acos(1+infj) = 1.570796326794897-infj
acos(nan) = nan
acos(-inf+1j) = nan
acos(inf+1j) = nan
acos(-inf+infj) = nan
acos(inf+infj) = nan
acos(-inf+1j) = 3.141592653589793-infj
acos(inf+1j) = 0-infj
acos(-inf+infj) = 2.356194490192345-infj
acos(inf+infj) = 0.7853981633974483-infj
acos(nan) = nan
acos(nan) = nan
acos(nan) = nan
Expand All @@ -94,14 +94,14 @@ atan(-1-1j) = -1.017221967897851-0.4023594781085251j
atan(2+2j) = 1.311223269671635+0.2388778612568591j
atan(nan) = nan
atan(nan) = nan
atan(1+infj) = nan
atan(nan) = nan
atan(-inf+1j) = nan
atan(inf+1j) = nan
atan(-inf+infj) = nan
atan(inf+infj) = nan
atan(nan) = nan
atan(1+infj) = 1.570796326794897+0j
atan(nan) = nan
atan(-inf+1j) = -1.570796326794897+0j
atan(inf+1j) = 1.570796326794897+0j
atan(-inf+infj) = -1.570796326794897+0j
atan(inf+infj) = 1.570796326794897+0j
atan(nan) = 1.570796326794897+0j
atan(nan) = -1.570796326794897+0j
atan(nan) = nan
atan(nan) = nan
atan(nan) = nan
Expand Down Expand Up @@ -161,17 +161,17 @@ tanh(nan) = nan
tanh(nan) = nan

asinh(0+0j) = 0+0j
asinh(-0+0j) = 0+0j
asinh(-0+0j) = -0+0j
asinh(-1-1j) = -1.061275061905036-0.6662394324925153j
asinh(2+2j) = 1.734324521487968+0.7542491446980463j
asinh(2+2j) = 1.734324521487966+0.7542491446980462j
asinh(nan) = nan
asinh(nan) = nan
asinh(1+infj) = nan
asinh(1+infj) = inf+1.570796326794897j
asinh(nan) = nan
asinh(-inf+1j) = -inf+0j
asinh(inf+1j) = -inf-2.356194490192345j
asinh(-inf+infj) = nan
asinh(inf+infj) = nan
asinh(inf+1j) = inf+0j
asinh(-inf+infj) = -inf+0.7853981633974483j
asinh(inf+infj) = inf+0.7853981633974483j
asinh(nan) = nan
asinh(nan) = nan
asinh(nan) = nan
Expand All @@ -180,38 +180,38 @@ asinh(nan) = nan

acosh(0+0j) = 0+1.570796326794897j
acosh(-0+0j) = 0+1.570796326794897j
acosh(-1-1j) = -0.825079184917074+1.66464254728848j
acosh(2+2j) = -1.356160455761411+0.03315117569827697j
acosh(-1-1j) = 1.061275061905036-2.237035759287412j
acosh(2+2j) = 1.734324521487966+0.8165471820968504j
acosh(nan) = nan
acosh(nan) = nan
acosh(1+infj) = nan
acosh(1+infj) = inf+1.570796326794897j
acosh(nan) = nan
acosh(-inf+1j) = nan
acosh(inf+1j) = nan
acosh(-inf+infj) = nan
acosh(inf+infj) = nan
acosh(-inf+1j) = inf+3.141592653589793j
acosh(inf+1j) = inf+0j
acosh(-inf+infj) = inf+2.356194490192345j
acosh(inf+infj) = inf+0.7853981633974483j
acosh(nan) = nan
acosh(nan) = nan
acosh(nan) = nan
acosh(nan) = nan
acosh(nan) = nan

atanh(0+0j) = 0+0j
atanh(-0+0j) = 0+0j
atanh(-0+0j) = -0+0j
atanh(-1-1j) = -0.4023594781085251-1.017221967897851j
atanh(2+2j) = 0.2388778612568591+1.311223269671635j
atanh(nan) = nan
atanh(nan) = nan
atanh(1+infj) = nan
atanh(nan) = nan
atanh(-inf+1j) = nan
atanh(inf+1j) = nan
atanh(-inf+infj) = nan
atanh(inf+infj) = nan
atanh(1+infj) = 0+1.570796326794897j
atanh(nan) = nan
atanh(-inf+1j) = -0+1.570796326794897j
atanh(inf+1j) = 0+1.570796326794897j
atanh(-inf+infj) = -0+1.570796326794897j
atanh(inf+infj) = 0+1.570796326794897j
atanh(nan) = nan
atanh(nan) = nan
atanh(nan) = nan
atanh(nan) = 0+1.570796326794897j
atanh(nan) = nan

exp(0+0j) = 1+0j
Expand Down Expand Up @@ -288,20 +288,20 @@ log(nan) = nan

sqrt(0+0j) = 0+0j
sqrt(-0+0j) = 0+0j
sqrt(-1-1j) = 0.4550898605622273+1.09868411346781j
sqrt(-1-1j) = 0.4550898605622273-1.09868411346781j
sqrt(2+2j) = 1.553773974030037+0.6435942529055826j
sqrt(nan) = nan
sqrt(nan) = nan
sqrt(1+infj) = 1.797693134862316e+308+infj
sqrt(1+infj) = inf+infj
sqrt(nan) = nan
sqrt(-inf+1j) = 0+infj
sqrt(inf+1j) = inf+0j
sqrt(-inf+infj) = 1.797693134862316e+308+infj
sqrt(inf+infj) = 1.797693134862316e+308+infj
sqrt(-inf+infj) = inf+infj
sqrt(inf+infj) = inf+infj
sqrt(nan) = nan
sqrt(nan) = nan
sqrt(nan) = nan
sqrt(nan) = 1.797693134862316e+308+infj
sqrt(nan) = inf+infj
sqrt(nan) = nan

proj(0+0j) = 0+0j
Expand All @@ -312,12 +312,12 @@ proj(nan) = nan
proj(nan) = nan
proj(1+infj) = inf+0j
proj(nan) = nan
proj(-inf+1j) = -inf+1j
proj(-inf+1j) = inf+0j
proj(inf+1j) = inf+0j
proj(-inf+infj) = inf+0j
proj(inf+infj) = inf+0j
proj(nan) = inf+0j
proj(nan) = nan
proj(nan) = inf+0j
proj(nan) = nan
proj(nan) = inf+0j
proj(nan) = nan
Expand Down

0 comments on commit b6708e1

Please sign in to comment.