| ... | @@ -19,8 +19,8 @@ pub fn exp(z: var) Complex(@typeOf(z.re)) { | ... | @@ -19,8 +19,8 @@ pub fn exp(z: var) Complex(@typeOf(z.re)) { |
| 19 | fn exp32(z: &const Complex(f32)) Complex(f32) { | 19 | fn exp32(z: &const Complex(f32)) Complex(f32) { |
| 20 | @setFloatMode(this, @import("builtin").FloatMode.Strict); | 20 | @setFloatMode(this, @import("builtin").FloatMode.Strict); |
| 21 | | 21 | |
| 22 | const exp_overflow = 0x42b17218; // max_exp * ln2 ~= 88.72283955 | 22 | const exp_overflow = 0x42b17218; // max_exp * ln2 ~= 88.72283955 |
| 23 | const cexp_overflow = 0x43400074; // (max_exp - min_denom_exp) * ln2 | 23 | const cexp_overflow = 0x43400074; // (max_exp - min_denom_exp) * ln2 |
| 24 | | 24 | |
| 25 | const x = z.re; | 25 | const x = z.re; |
| 26 | const y = z.im; | 26 | const y = z.im; |
| ... | @@ -41,12 +41,10 @@ fn exp32(z: &const Complex(f32)) Complex(f32) { | ... | @@ -41,12 +41,10 @@ fn exp32(z: &const Complex(f32)) Complex(f32) { |
| 41 | // cexp(finite|nan +- i inf|nan) = nan + i nan | 41 | // cexp(finite|nan +- i inf|nan) = nan + i nan |
| 42 | if ((hx & 0x7fffffff) != 0x7f800000) { | 42 | if ((hx & 0x7fffffff) != 0x7f800000) { |
| 43 | return Complex(f32).new(y - y, y - y); | 43 | return Complex(f32).new(y - y, y - y); |
| 44 | } | 44 | } // cexp(-inf +- i inf|nan) = 0 + i0 |
| 45 | // cexp(-inf +- i inf|nan) = 0 + i0 | | |
| 46 | else if (hx & 0x80000000 != 0) { | 45 | else if (hx & 0x80000000 != 0) { |
| 47 | return Complex(f32).new(0, 0); | 46 | return Complex(f32).new(0, 0); |
| 48 | } | 47 | } // cexp(+inf +- i inf|nan) = inf + i nan |
| 49 | // cexp(+inf +- i inf|nan) = inf + i nan | | |
| 50 | else { | 48 | else { |
| 51 | return Complex(f32).new(x, y - y); | 49 | return Complex(f32).new(x, y - y); |
| 52 | } | 50 | } |
| ... | @@ -55,8 +53,7 @@ fn exp32(z: &const Complex(f32)) Complex(f32) { | ... | @@ -55,8 +53,7 @@ fn exp32(z: &const Complex(f32)) Complex(f32) { |
| 55 | // 88.7 <= x <= 192 so must scale | 53 | // 88.7 <= x <= 192 so must scale |
| 56 | if (hx >= exp_overflow and hx <= cexp_overflow) { | 54 | if (hx >= exp_overflow and hx <= cexp_overflow) { |
| 57 | return ldexp_cexp(z, 0); | 55 | return ldexp_cexp(z, 0); |
| 58 | } | 56 | } // - x < exp_overflow => exp(x) won't overflow (common) |
| 59 | // - x < exp_overflow => exp(x) won't overflow (common) | | |
| 60 | // - x > cexp_overflow, so exp(x) * s overflows for s > 0 | 57 | // - x > cexp_overflow, so exp(x) * s overflows for s > 0 |
| 61 | // - x = +-inf | 58 | // - x = +-inf |
| 62 | // - x = nan | 59 | // - x = nan |
| ... | @@ -67,8 +64,8 @@ fn exp32(z: &const Complex(f32)) Complex(f32) { | ... | @@ -67,8 +64,8 @@ fn exp32(z: &const Complex(f32)) Complex(f32) { |
| 67 | } | 64 | } |
| 68 | | 65 | |
| 69 | fn exp64(z: &const Complex(f64)) Complex(f64) { | 66 | fn exp64(z: &const Complex(f64)) Complex(f64) { |
| 70 | const exp_overflow = 0x40862e42; // high bits of max_exp * ln2 ~= 710 | 67 | const exp_overflow = 0x40862e42; // high bits of max_exp * ln2 ~= 710 |
| 71 | const cexp_overflow = 0x4096b8e4; // (max_exp - min_denorm_exp) * ln2 | 68 | const cexp_overflow = 0x4096b8e4; // (max_exp - min_denorm_exp) * ln2 |
| 72 | | 69 | |
| 73 | const x = z.re; | 70 | const x = z.re; |
| 74 | const y = z.im; | 71 | const y = z.im; |
| ... | @@ -95,12 +92,10 @@ fn exp64(z: &const Complex(f64)) Complex(f64) { | ... | @@ -95,12 +92,10 @@ fn exp64(z: &const Complex(f64)) Complex(f64) { |
| 95 | // cexp(finite|nan +- i inf|nan) = nan + i nan | 92 | // cexp(finite|nan +- i inf|nan) = nan + i nan |
| 96 | if (lx != 0 or (hx & 0x7fffffff) != 0x7ff00000) { | 93 | if (lx != 0 or (hx & 0x7fffffff) != 0x7ff00000) { |
| 97 | return Complex(f64).new(y - y, y - y); | 94 | return Complex(f64).new(y - y, y - y); |
| 98 | } | 95 | } // cexp(-inf +- i inf|nan) = 0 + i0 |
| 99 | // cexp(-inf +- i inf|nan) = 0 + i0 | | |
| 100 | else if (hx & 0x80000000 != 0) { | 96 | else if (hx & 0x80000000 != 0) { |
| 101 | return Complex(f64).new(0, 0); | 97 | return Complex(f64).new(0, 0); |
| 102 | } | 98 | } // cexp(+inf +- i inf|nan) = inf + i nan |
| 103 | // cexp(+inf +- i inf|nan) = inf + i nan | | |
| 104 | else { | 99 | else { |
| 105 | return Complex(f64).new(x, y - y); | 100 | return Complex(f64).new(x, y - y); |
| 106 | } | 101 | } |
| ... | @@ -109,9 +104,8 @@ fn exp64(z: &const Complex(f64)) Complex(f64) { | ... | @@ -109,9 +104,8 @@ fn exp64(z: &const Complex(f64)) Complex(f64) { |
| 109 | // 709.7 <= x <= 1454.3 so must scale | 104 | // 709.7 <= x <= 1454.3 so must scale |
| 110 | if (hx >= exp_overflow and hx <= cexp_overflow) { | 105 | if (hx >= exp_overflow and hx <= cexp_overflow) { |
| 111 | const r = ldexp_cexp(z, 0); | 106 | const r = ldexp_cexp(z, 0); |
| 112 | return *r; | 107 | return r.*; |
| 113 | } | 108 | } // - x < exp_overflow => exp(x) won't overflow (common) |
| 114 | // - x < exp_overflow => exp(x) won't overflow (common) | | |
| 115 | // - x > cexp_overflow, so exp(x) * s overflows for s > 0 | 109 | // - x > cexp_overflow, so exp(x) * s overflows for s > 0 |
| 116 | // - x = +-inf | 110 | // - x = +-inf |
| 117 | // - x = nan | 111 | // - x = nan |