authorgravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2019-05-02 19:05:26+12:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2019-05-02 19:05:26+12:00
logf950ec0c16de6dba8a541b1a6453ebe431430782
tree47df2ed19eece69d35a60eeff8c26511cfe996e3
parentc00c18de6a5e436b1c362c06d6e5259c8f731a90
parent3370e60dd96828db6f9d1da36e71064303e97de6
signaturebadge-question-mark Signed by PGP key 4AEE18F83AFDEB23

Merge pull request #2397 from ziglang/std.math

Review std/math and update documentation

69 files changed, 806 insertions(+), 549 deletions(-)

std/math/acos.zig+8-2
...@@ -1,11 +1,17 @@...@@ -1,11 +1,17 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - acos(x) = nan if x < -1 or x > 14// https://git.musl-libc.org/cgit/musl/tree/src/math/acosf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/acos.c
46
5const std = @import("../std.zig");7const std = @import("../std.zig");
6const math = std.math;8const math = std.math;
7const expect = std.testing.expect;9const expect = std.testing.expect;
810
11/// Returns the arc-cosine of x.
12///
13/// Special cases:
14/// - acos(x) = nan if x < -1 or x > 1
9pub fn acos(x: var) @typeOf(x) {15pub fn acos(x: var) @typeOf(x) {
10 const T = @typeOf(x);16 const T = @typeOf(x);
11 return switch (T) {17 return switch (T) {
std/math/acosh.zig+9-3
...@@ -1,13 +1,19 @@...@@ -1,13 +1,19 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - acosh(x) = snan if x < 14// https://git.musl-libc.org/cgit/musl/tree/src/math/acoshf.c
4// - acosh(nan) = nan5// https://git.musl-libc.org/cgit/musl/tree/src/math/acosh.c
56
6const builtin = @import("builtin");7const builtin = @import("builtin");
7const std = @import("../std.zig");8const std = @import("../std.zig");
8const math = std.math;9const math = std.math;
9const expect = std.testing.expect;10const expect = std.testing.expect;
1011
12/// Returns the hyperbolic arc-cosine of x.
13///
14/// Special cases:
15/// - acosh(x) = snan if x < 1
16/// - acosh(nan) = nan
11pub fn acosh(x: var) @typeOf(x) {17pub fn acosh(x: var) @typeOf(x) {
12 const T = @typeOf(x);18 const T = @typeOf(x);
13 return switch (T) {19 return switch (T) {
std/math/asin.zig+9-3
...@@ -1,12 +1,18 @@...@@ -1,12 +1,18 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - asin(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/asinf.c
4// - asin(x) = nan if x < -1 or x > 15// https://git.musl-libc.org/cgit/musl/tree/src/math/asin.c
56
6const std = @import("../std.zig");7const std = @import("../std.zig");
7const math = std.math;8const math = std.math;
8const expect = std.testing.expect;9const expect = std.testing.expect;
910
11/// Returns the arc-sin of x.
12///
13/// Special Cases:
14/// - asin(+-0) = +-0
15/// - asin(x) = nan if x < -1 or x > 1
10pub fn asin(x: var) @typeOf(x) {16pub fn asin(x: var) @typeOf(x) {
11 const T = @typeOf(x);17 const T = @typeOf(x);
12 return switch (T) {18 return switch (T) {
std/math/asinh.zig+10-4
...@@ -1,14 +1,20 @@...@@ -1,14 +1,20 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - asinh(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/asinhf.c
4// - asinh(+-inf) = +-inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/asinh.c
5// - asinh(nan) = nan
66
7const std = @import("../std.zig");7const std = @import("../std.zig");
8const math = std.math;8const math = std.math;
9const expect = std.testing.expect;9const expect = std.testing.expect;
10const maxInt = std.math.maxInt;10const maxInt = std.math.maxInt;
1111
12/// Returns the hyperbolic arc-sin of x.
13///
14/// Special Cases:
15/// - asinh(+-0) = +-0
16/// - asinh(+-inf) = +-inf
17/// - asinh(nan) = nan
12pub fn asinh(x: var) @typeOf(x) {18pub fn asinh(x: var) @typeOf(x) {
13 const T = @typeOf(x);19 const T = @typeOf(x);
14 return switch (T) {20 return switch (T) {
std/math/atan.zig+9-3
...@@ -1,12 +1,18 @@...@@ -1,12 +1,18 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - atan(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/atanf.c
4// - atan(+-inf) = +-pi/25// https://git.musl-libc.org/cgit/musl/tree/src/math/atan.c
56
6const std = @import("../std.zig");7const std = @import("../std.zig");
7const math = std.math;8const math = std.math;
8const expect = std.testing.expect;9const expect = std.testing.expect;
910
11/// Returns the arc-tangent of x.
12///
13/// Special Cases:
14/// - atan(+-0) = +-0
15/// - atan(+-inf) = +-pi/2
10pub fn atan(x: var) @typeOf(x) {16pub fn atan(x: var) @typeOf(x) {
11 const T = @typeOf(x);17 const T = @typeOf(x);
12 return switch (T) {18 return switch (T) {
std/math/atan2.zig+24-18
...@@ -1,27 +1,33 @@...@@ -1,27 +1,33 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// atan2(y, nan) = nan4// https://git.musl-libc.org/cgit/musl/tree/src/math/atan2f.c
4// atan2(nan, x) = nan5// https://git.musl-libc.org/cgit/musl/tree/src/math/atan2.c
5// atan2(+0, x>=0) = +0
6// atan2(-0, x>=0) = -0
7// atan2(+0, x<=-0) = +pi
8// atan2(-0, x<=-0) = -pi
9// atan2(y>0, 0) = +pi/2
10// atan2(y<0, 0) = -pi/2
11// atan2(+inf, +inf) = +pi/4
12// atan2(-inf, +inf) = -pi/4
13// atan2(+inf, -inf) = 3pi/4
14// atan2(-inf, -inf) = -3pi/4
15// atan2(y, +inf) = 0
16// atan2(y>0, -inf) = +pi
17// atan2(y<0, -inf) = -pi
18// atan2(+inf, x) = +pi/2
19// atan2(-inf, x) = -pi/2
206
21const std = @import("../std.zig");7const std = @import("../std.zig");
22const math = std.math;8const math = std.math;
23const expect = std.testing.expect;9const expect = std.testing.expect;
2410
11/// Returns the arc-tangent of y/x.
12///
13/// Special Cases:
14/// - atan2(y, nan) = nan
15/// - atan2(nan, x) = nan
16/// - atan2(+0, x>=0) = +0
17/// - atan2(-0, x>=0) = -0
18/// - atan2(+0, x<=-0) = +pi
19/// - atan2(-0, x<=-0) = -pi
20/// - atan2(y>0, 0) = +pi/2
21/// - atan2(y<0, 0) = -pi/2
22/// - atan2(+inf, +inf) = +pi/4
23/// - atan2(-inf, +inf) = -pi/4
24/// - atan2(+inf, -inf) = 3pi/4
25/// - atan2(-inf, -inf) = -3pi/4
26/// - atan2(y, +inf) = 0
27/// - atan2(y>0, -inf) = +pi
28/// - atan2(y<0, -inf) = -pi
29/// - atan2(+inf, x) = +pi/2
30/// - atan2(-inf, x) = -pi/2
25pub fn atan2(comptime T: type, y: T, x: T) T {31pub fn atan2(comptime T: type, y: T, x: T) T {
26 return switch (T) {32 return switch (T) {
27 f32 => atan2_32(y, x),33 f32 => atan2_32(y, x),
std/math/atanh.zig+10-4
...@@ -1,14 +1,20 @@...@@ -1,14 +1,20 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - atanh(+-1) = +-inf with signal4// https://git.musl-libc.org/cgit/musl/tree/src/math/atanhf.c
4// - atanh(x) = nan if |x| > 1 with signal5// https://git.musl-libc.org/cgit/musl/tree/src/math/atanh.c
5// - atanh(nan) = nan
66
7const std = @import("../std.zig");7const std = @import("../std.zig");
8const math = std.math;8const math = std.math;
9const expect = std.testing.expect;9const expect = std.testing.expect;
10const maxInt = std.math.maxInt;10const maxInt = std.math.maxInt;
1111
12/// Returns the hyperbolic arc-tangent of x.
13///
14/// Special Cases:
15/// - atanh(+-1) = +-inf with signal
16/// - atanh(x) = nan if |x| > 1 with signal
17/// - atanh(nan) = nan
12pub fn atanh(x: var) @typeOf(x) {18pub fn atanh(x: var) @typeOf(x) {
13 const T = @typeOf(x);19 const T = @typeOf(x);
14 return switch (T) {20 return switch (T) {
std/math/big/int.zig+112-45
...@@ -21,32 +21,49 @@ comptime {...@@ -21,32 +21,49 @@ comptime {
21 debug.assert(Limb.is_signed == false);21 debug.assert(Limb.is_signed == false);
22}22}
2323
24/// An arbitrary-precision big integer.
25///
26/// Memory is allocated by an Int as needed to ensure operations never overflow. The range of an
27/// Int is bounded only by available memory.
24pub const Int = struct {28pub const Int = struct {
25 const sign_bit: usize = 1 << (usize.bit_count - 1);29 const sign_bit: usize = 1 << (usize.bit_count - 1);
2630
31 /// Default number of limbs to allocate on creation of an Int.
32 pub const default_capacity = 4;
33
34 /// Allocator used by the Int when requesting memory.
27 allocator: ?*Allocator,35 allocator: ?*Allocator,
28 // - little-endian ordered36
29 // - len >= 1 always37 /// Raw digits. These are:
30 // - zero value -> len == 1 with limbs[0] == 038 ///
39 /// * Little-endian ordered
40 /// * limbs.len >= 1
41 /// * Zero is represent as Int.len() == 1 with limbs[0] == 0.
42 ///
43 /// Accessing limbs directly should be avoided.
31 limbs: []Limb,44 limbs: []Limb,
32 // High bit is the sign bit. 1 is negative, 0 positive.
33 // Remaining bits indicate the number of used limbs.
34 //
35 // If Zig gets smarter about packing data, this can be rewritten as a u1 and usize - 1 field.
36 metadata: usize,
3745
38 const default_capacity = 4;46 /// High bit is the sign bit. If set, Int is negative, else Int is positive.
47 /// The remaining bits represent the number of limbs used by Int.
48 metadata: usize,
3949
50 /// Creates a new Int. default_capacity limbs will be allocated immediately.
51 /// Int will be zeroed.
40 pub fn init(allocator: *Allocator) !Int {52 pub fn init(allocator: *Allocator) !Int {
41 return try Int.initCapacity(allocator, default_capacity);53 return try Int.initCapacity(allocator, default_capacity);
42 }54 }
4355
56 /// Creates a new Int. Int will be set to `value`.
57 ///
58 /// This is identical to an `init`, followed by a `set`.
44 pub fn initSet(allocator: *Allocator, value: var) !Int {59 pub fn initSet(allocator: *Allocator, value: var) !Int {
45 var s = try Int.init(allocator);60 var s = try Int.init(allocator);
46 try s.set(value);61 try s.set(value);
47 return s;62 return s;
48 }63 }
4964
65 /// Creates a new Int with a specific capacity. If capacity < default_capacity then the
66 /// default capacity will be used instead.
50 pub fn initCapacity(allocator: *Allocator, capacity: usize) !Int {67 pub fn initCapacity(allocator: *Allocator, capacity: usize) !Int {
51 return Int{68 return Int{
52 .allocator = allocator,69 .allocator = allocator,
...@@ -59,14 +76,17 @@ pub const Int = struct {...@@ -59,14 +76,17 @@ pub const Int = struct {
59 };76 };
60 }77 }
6178
79 /// Returns the number of limbs currently in use.
62 pub fn len(self: Int) usize {80 pub fn len(self: Int) usize {
63 return self.metadata & ~sign_bit;81 return self.metadata & ~sign_bit;
64 }82 }
6583
84 /// Returns whether an Int is positive.
66 pub fn isPositive(self: Int) bool {85 pub fn isPositive(self: Int) bool {
67 return self.metadata & sign_bit == 0;86 return self.metadata & sign_bit == 0;
68 }87 }
6988
89 /// Sets the sign of an Int.
70 pub fn setSign(self: *Int, positive: bool) void {90 pub fn setSign(self: *Int, positive: bool) void {
71 if (positive) {91 if (positive) {
72 self.metadata &= ~sign_bit;92 self.metadata &= ~sign_bit;
...@@ -75,14 +95,17 @@ pub const Int = struct {...@@ -75,14 +95,17 @@ pub const Int = struct {
75 }95 }
76 }96 }
7797
98 /// Sets the length of an Int.
99 ///
100 /// If setLen is used, then the Int must be normalized to suit.
78 pub fn setLen(self: *Int, new_len: usize) void {101 pub fn setLen(self: *Int, new_len: usize) void {
79 self.metadata &= sign_bit;102 self.metadata &= sign_bit;
80 self.metadata |= new_len;103 self.metadata |= new_len;
81 }104 }
82105
83 // Initialize an Int directly from a fixed set of limb values. This is considered read-only106 /// Returns an Int backed by a fixed set of limb values.
84 // and cannot be used as a receiver argument to any functions. If this tries to allocate107 /// This is read-only and cannot be used as a result argument. If the Int tries to allocate
85 // at any point a panic will occur due to the null allocator.108 /// memory a runtime panic will occur.
86 pub fn initFixed(limbs: []const Limb) Int {109 pub fn initFixed(limbs: []const Limb) Int {
87 var self = Int{110 var self = Int{
88 .allocator = null,111 .allocator = null,
...@@ -95,6 +118,9 @@ pub const Int = struct {...@@ -95,6 +118,9 @@ pub const Int = struct {
95 return self;118 return self;
96 }119 }
97120
121 /// Ensures an Int has enough space allocated for capacity limbs. If the Int does not have
122 /// sufficient capacity, the exact amount will be allocated. This occurs even if the requested
123 /// capacity is only greater than the current capacity by one limb.
98 pub fn ensureCapacity(self: *Int, capacity: usize) !void {124 pub fn ensureCapacity(self: *Int, capacity: usize) !void {
99 self.assertWritable();125 self.assertWritable();
100 if (capacity <= self.limbs.len) {126 if (capacity <= self.limbs.len) {
...@@ -110,12 +136,15 @@ pub const Int = struct {...@@ -110,12 +136,15 @@ pub const Int = struct {
110 }136 }
111 }137 }
112138
139 /// Frees all memory associated with an Int.
113 pub fn deinit(self: *Int) void {140 pub fn deinit(self: *Int) void {
114 self.assertWritable();141 self.assertWritable();
115 self.allocator.?.free(self.limbs);142 self.allocator.?.free(self.limbs);
116 self.* = undefined;143 self.* = undefined;
117 }144 }
118145
146 /// Clones an Int and returns a new Int with the same value. The new Int is a deep copy and
147 /// can be modified separately from the original.
119 pub fn clone(other: Int) !Int {148 pub fn clone(other: Int) !Int {
120 other.assertWritable();149 other.assertWritable();
121 return Int{150 return Int{
...@@ -129,6 +158,8 @@ pub const Int = struct {...@@ -129,6 +158,8 @@ pub const Int = struct {
129 };158 };
130 }159 }
131160
161 /// Copies the value of an Int to an existing Int so that they both have the same value.
162 /// Extra memory will be allocated if the receiver does not have enough capacity.
132 pub fn copy(self: *Int, other: Int) !void {163 pub fn copy(self: *Int, other: Int) !void {
133 self.assertWritable();164 self.assertWritable();
134 if (self.limbs.ptr == other.limbs.ptr) {165 if (self.limbs.ptr == other.limbs.ptr) {
...@@ -140,6 +171,8 @@ pub const Int = struct {...@@ -140,6 +171,8 @@ pub const Int = struct {
140 self.metadata = other.metadata;171 self.metadata = other.metadata;
141 }172 }
142173
174 /// Efficiently swap an Int with another. This swaps the limb pointers and a full copy is not
175 /// performed. The address of the limbs field will not be the same after this function.
143 pub fn swap(self: *Int, other: *Int) void {176 pub fn swap(self: *Int, other: *Int) void {
144 self.assertWritable();177 self.assertWritable();
145 mem.swap(Int, self, other);178 mem.swap(Int, self, other);
...@@ -152,35 +185,39 @@ pub const Int = struct {...@@ -152,35 +185,39 @@ pub const Int = struct {
152 debug.warn("\n");185 debug.warn("\n");
153 }186 }
154187
188 /// Negate the sign of an Int.
155 pub fn negate(self: *Int) void {189 pub fn negate(self: *Int) void {
156 self.metadata ^= sign_bit;190 self.metadata ^= sign_bit;
157 }191 }
158192
193 /// Make an Int positive.
159 pub fn abs(self: *Int) void {194 pub fn abs(self: *Int) void {
160 self.metadata &= ~sign_bit;195 self.metadata &= ~sign_bit;
161 }196 }
162197
198 /// Returns true if an Int is odd.
163 pub fn isOdd(self: Int) bool {199 pub fn isOdd(self: Int) bool {
164 return self.limbs[0] & 1 != 0;200 return self.limbs[0] & 1 != 0;
165 }201 }
166202
203 /// Returns true if an Int is even.
167 pub fn isEven(self: Int) bool {204 pub fn isEven(self: Int) bool {
168 return !self.isOdd();205 return !self.isOdd();
169 }206 }
170207
171 // Returns the number of bits required to represent the absolute value of self.208 /// Returns the number of bits required to represent the absolute value an Int.
172 fn bitCountAbs(self: Int) usize {209 fn bitCountAbs(self: Int) usize {
173 return (self.len() - 1) * Limb.bit_count + (Limb.bit_count - @clz(self.limbs[self.len() - 1]));210 return (self.len() - 1) * Limb.bit_count + (Limb.bit_count - @clz(self.limbs[self.len() - 1]));
174 }211 }
175212
176 // Returns the number of bits required to represent the integer in twos-complement form.213 /// Returns the number of bits required to represent the integer in twos-complement form.
177 //214 ///
178 // If the integer is negative the value returned is the number of bits needed by a signed215 /// If the integer is negative the value returned is the number of bits needed by a signed
179 // integer to represent the value. If positive the value is the number of bits for an216 /// integer to represent the value. If positive the value is the number of bits for an
180 // unsigned integer. Any unsigned integer will fit in the signed integer with bitcount217 /// unsigned integer. Any unsigned integer will fit in the signed integer with bitcount
181 // one greater than the returned value.218 /// one greater than the returned value.
182 //219 ///
183 // e.g. -127 returns 8 as it will fit in an i8. 127 returns 7 since it fits in a u7.220 /// e.g. -127 returns 8 as it will fit in an i8. 127 returns 7 since it fits in a u7.
184 fn bitCountTwosComp(self: Int) usize {221 fn bitCountTwosComp(self: Int) usize {
185 var bits = self.bitCountAbs();222 var bits = self.bitCountAbs();
186223
...@@ -203,7 +240,7 @@ pub const Int = struct {...@@ -203,7 +240,7 @@ pub const Int = struct {
203 return bits;240 return bits;
204 }241 }
205242
206 pub fn fitsInTwosComp(self: Int, is_signed: bool, bit_count: usize) bool {243 fn fitsInTwosComp(self: Int, is_signed: bool, bit_count: usize) bool {
207 if (self.eqZero()) {244 if (self.eqZero()) {
208 return true;245 return true;
209 }246 }
...@@ -215,18 +252,20 @@ pub const Int = struct {...@@ -215,18 +252,20 @@ pub const Int = struct {
215 return bit_count >= req_bits;252 return bit_count >= req_bits;
216 }253 }
217254
255 /// Returns whether self can fit into an integer of the requested type.
218 pub fn fits(self: Int, comptime T: type) bool {256 pub fn fits(self: Int, comptime T: type) bool {
219 return self.fitsInTwosComp(T.is_signed, T.bit_count);257 return self.fitsInTwosComp(T.is_signed, T.bit_count);
220 }258 }
221259
222 // Returns the approximate size of the integer in the given base. Negative values accommodate for260 /// Returns the approximate size of the integer in the given base. Negative values accommodate for
223 // the minus sign. This is used for determining the number of characters needed to print the261 /// the minus sign. This is used for determining the number of characters needed to print the
224 // value. It is inexact and will exceed the given value by 1-2 digits.262 /// value. It is inexact and may exceed the given value by ~1-2 bytes.
225 pub fn sizeInBase(self: Int, base: usize) usize {263 pub fn sizeInBase(self: Int, base: usize) usize {
226 const bit_count = usize(@boolToInt(!self.isPositive())) + self.bitCountAbs();264 const bit_count = usize(@boolToInt(!self.isPositive())) + self.bitCountAbs();
227 return (bit_count / math.log2(base)) + 1;265 return (bit_count / math.log2(base)) + 1;
228 }266 }
229267
268 /// Sets an Int to value. Value must be an primitive integer type.
230 pub fn set(self: *Int, value: var) Allocator.Error!void {269 pub fn set(self: *Int, value: var) Allocator.Error!void {
231 self.assertWritable();270 self.assertWritable();
232 const T = @typeOf(value);271 const T = @typeOf(value);
...@@ -290,6 +329,9 @@ pub const Int = struct {...@@ -290,6 +329,9 @@ pub const Int = struct {
290 TargetTooSmall,329 TargetTooSmall,
291 };330 };
292331
332 /// Convert self to type T.
333 ///
334 /// Returns an error if self cannot be narrowed into the requested type without truncation.
293 pub fn to(self: Int, comptime T: type) ConvertError!T {335 pub fn to(self: Int, comptime T: type) ConvertError!T {
294 switch (@typeId(T)) {336 switch (@typeId(T)) {
295 TypeId.Int => {337 TypeId.Int => {
...@@ -353,6 +395,13 @@ pub const Int = struct {...@@ -353,6 +395,13 @@ pub const Int = struct {
353 };395 };
354 }396 }
355397
398 /// Set self from the string representation `value`.
399 ///
400 /// value must contain only digits <= `base`. Base prefixes are not allowed (e.g. 0x43 should
401 /// simply be 43).
402 ///
403 /// Returns an error if memory could not be allocated or `value` has invalid digits for the
404 /// requested base.
356 pub fn setString(self: *Int, base: u8, value: []const u8) !void {405 pub fn setString(self: *Int, base: u8, value: []const u8) !void {
357 self.assertWritable();406 self.assertWritable();
358 if (base < 2 or base > 16) {407 if (base < 2 or base > 16) {
...@@ -380,6 +429,8 @@ pub const Int = struct {...@@ -380,6 +429,8 @@ pub const Int = struct {
380 self.setSign(positive);429 self.setSign(positive);
381 }430 }
382431
432 /// Converts self to a string in the requested base. Memory is allocated from the provided
433 /// allocator and not the one present in self.
383 /// TODO make this call format instead of the other way around434 /// TODO make this call format instead of the other way around
384 pub fn toString(self: Int, allocator: *Allocator, base: u8) ![]const u8 {435 pub fn toString(self: Int, allocator: *Allocator, base: u8) ![]const u8 {
385 if (base < 2 or base > 16) {436 if (base < 2 or base > 16) {
...@@ -463,7 +514,7 @@ pub const Int = struct {...@@ -463,7 +514,7 @@ pub const Int = struct {
463 return s;514 return s;
464 }515 }
465516
466 /// for the std lib format function517 /// To allow `std.fmt.printf` to work with Int.
467 /// TODO make this non-allocating518 /// TODO make this non-allocating
468 pub fn format(519 pub fn format(
469 self: Int,520 self: Int,
...@@ -480,7 +531,7 @@ pub const Int = struct {...@@ -480,7 +531,7 @@ pub const Int = struct {
480 return output(context, str);531 return output(context, str);
481 }532 }
482533
483 // returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively.534 /// Returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively.
484 pub fn cmpAbs(a: Int, b: Int) i8 {535 pub fn cmpAbs(a: Int, b: Int) i8 {
485 if (a.len() < b.len()) {536 if (a.len() < b.len()) {
486 return -1;537 return -1;
...@@ -505,7 +556,7 @@ pub const Int = struct {...@@ -505,7 +556,7 @@ pub const Int = struct {
505 }556 }
506 }557 }
507558
508 // returns -1, 0, 1 if a < b, a == b or a > b respectively.559 /// Returns -1, 0, 1 if a < b, a == b or a > b respectively.
509 pub fn cmp(a: Int, b: Int) i8 {560 pub fn cmp(a: Int, b: Int) i8 {
510 if (a.isPositive() != b.isPositive()) {561 if (a.isPositive() != b.isPositive()) {
511 return if (a.isPositive()) i8(1) else -1;562 return if (a.isPositive()) i8(1) else -1;
...@@ -515,17 +566,17 @@ pub const Int = struct {...@@ -515,17 +566,17 @@ pub const Int = struct {
515 }566 }
516 }567 }
517568
518 // if a == 0569 /// Returns true if a == 0.
519 pub fn eqZero(a: Int) bool {570 pub fn eqZero(a: Int) bool {
520 return a.len() == 1 and a.limbs[0] == 0;571 return a.len() == 1 and a.limbs[0] == 0;
521 }572 }
522573
523 // if |a| == |b|574 /// Returns true if |a| == |b|.
524 pub fn eqAbs(a: Int, b: Int) bool {575 pub fn eqAbs(a: Int, b: Int) bool {
525 return cmpAbs(a, b) == 0;576 return cmpAbs(a, b) == 0;
526 }577 }
527578
528 // if a == b579 /// Returns true if a == b.
529 pub fn eq(a: Int, b: Int) bool {580 pub fn eq(a: Int, b: Int) bool {
530 return cmp(a, b) == 0;581 return cmp(a, b) == 0;
531 }582 }
...@@ -559,7 +610,11 @@ pub const Int = struct {...@@ -559,7 +610,11 @@ pub const Int = struct {
559 };610 };
560 }611 }
561612
562 // r = a + b613 /// r = a + b
614 ///
615 /// r, a and b may be aliases.
616 ///
617 /// Returns an error if memory could not be allocated.
563 pub fn add(r: *Int, a: Int, b: Int) Allocator.Error!void {618 pub fn add(r: *Int, a: Int, b: Int) Allocator.Error!void {
564 r.assertWritable();619 r.assertWritable();
565 if (a.eqZero()) {620 if (a.eqZero()) {
...@@ -617,7 +672,11 @@ pub const Int = struct {...@@ -617,7 +672,11 @@ pub const Int = struct {
617 r[i] = carry;672 r[i] = carry;
618 }673 }
619674
620 // r = a - b675 /// r = a - b
676 ///
677 /// r, a and b may be aliases.
678 ///
679 /// Returns an error if memory could not be allocated.
621 pub fn sub(r: *Int, a: Int, b: Int) !void {680 pub fn sub(r: *Int, a: Int, b: Int) !void {
622 r.assertWritable();681 r.assertWritable();
623 if (a.isPositive() != b.isPositive()) {682 if (a.isPositive() != b.isPositive()) {
...@@ -684,9 +743,11 @@ pub const Int = struct {...@@ -684,9 +743,11 @@ pub const Int = struct {
684 debug.assert(borrow == 0);743 debug.assert(borrow == 0);
685 }744 }
686745
687 // rma = a * b746 /// rma = a * b
688 //747 ///
689 // For greatest efficiency, ensure rma does not alias a or b.748 /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
749 ///
750 /// Returns an error if memory could not be allocated.
690 pub fn mul(rma: *Int, a: Int, b: Int) !void {751 pub fn mul(rma: *Int, a: Int, b: Int) !void {
691 rma.assertWritable();752 rma.assertWritable();
692753
...@@ -759,6 +820,9 @@ pub const Int = struct {...@@ -759,6 +820,9 @@ pub const Int = struct {
759 }820 }
760 }821 }
761822
823 /// q = a / b (rem r)
824 ///
825 /// a / b are floored (rounded towards 0).
762 pub fn divFloor(q: *Int, r: *Int, a: Int, b: Int) !void {826 pub fn divFloor(q: *Int, r: *Int, a: Int, b: Int) !void {
763 try div(q, r, a, b);827 try div(q, r, a, b);
764828
...@@ -771,6 +835,9 @@ pub const Int = struct {...@@ -771,6 +835,9 @@ pub const Int = struct {
771 r.setSign(b.isPositive());835 r.setSign(b.isPositive());
772 }836 }
773837
838 /// q = a / b (rem r)
839 ///
840 /// a / b are truncated (rounded towards -inf).
774 pub fn divTrunc(q: *Int, r: *Int, a: Int, b: Int) !void {841 pub fn divTrunc(q: *Int, r: *Int, a: Int, b: Int) !void {
775 try div(q, r, a, b);842 try div(q, r, a, b);
776 r.setSign(a.isPositive());843 r.setSign(a.isPositive());
...@@ -969,7 +1036,7 @@ pub const Int = struct {...@@ -969,7 +1036,7 @@ pub const Int = struct {
969 r.normalize(r.len());1036 r.normalize(r.len());
970 }1037 }
9711038
972 // r = a << shift, in other words, r = a * 2^shift1039 /// r = a << shift, in other words, r = a * 2^shift
973 pub fn shiftLeft(r: *Int, a: Int, shift: usize) !void {1040 pub fn shiftLeft(r: *Int, a: Int, shift: usize) !void {
974 r.assertWritable();1041 r.assertWritable();
9751042
...@@ -1002,7 +1069,7 @@ pub const Int = struct {...@@ -1002,7 +1069,7 @@ pub const Int = struct {
1002 mem.set(Limb, r[0 .. limb_shift - 1], 0);1069 mem.set(Limb, r[0 .. limb_shift - 1], 0);
1003 }1070 }
10041071
1005 // r = a >> shift1072 /// r = a >> shift
1006 pub fn shiftRight(r: *Int, a: Int, shift: usize) !void {1073 pub fn shiftRight(r: *Int, a: Int, shift: usize) !void {
1007 r.assertWritable();1074 r.assertWritable();
10081075
...@@ -1038,7 +1105,9 @@ pub const Int = struct {...@@ -1038,7 +1105,9 @@ pub const Int = struct {
1038 }1105 }
1039 }1106 }
10401107
1041 // r = a | b1108 /// r = a | b
1109 ///
1110 /// a and b are zero-extended to the longer of a or b.
1042 pub fn bitOr(r: *Int, a: Int, b: Int) !void {1111 pub fn bitOr(r: *Int, a: Int, b: Int) !void {
1043 r.assertWritable();1112 r.assertWritable();
10441113
...@@ -1067,7 +1136,7 @@ pub const Int = struct {...@@ -1067,7 +1136,7 @@ pub const Int = struct {
1067 }1136 }
1068 }1137 }
10691138
1070 // r = a & b1139 /// r = a & b
1071 pub fn bitAnd(r: *Int, a: Int, b: Int) !void {1140 pub fn bitAnd(r: *Int, a: Int, b: Int) !void {
1072 r.assertWritable();1141 r.assertWritable();
10731142
...@@ -1093,7 +1162,7 @@ pub const Int = struct {...@@ -1093,7 +1162,7 @@ pub const Int = struct {
1093 }1162 }
1094 }1163 }
10951164
1096 // r = a ^ b1165 /// r = a ^ b
1097 pub fn bitXor(r: *Int, a: Int, b: Int) !void {1166 pub fn bitXor(r: *Int, a: Int, b: Int) !void {
1098 r.assertWritable();1167 r.assertWritable();
10991168
...@@ -1279,10 +1348,8 @@ test "big.int bitcount/to" {...@@ -1279,10 +1348,8 @@ test "big.int bitcount/to" {
1279 try a.set(0);1348 try a.set(0);
1280 testing.expect(a.bitCountTwosComp() == 0);1349 testing.expect(a.bitCountTwosComp() == 0);
12811350
1282 // TODO: stack smashing1351 testing.expect((try a.to(u0)) == 0);
1283 // testing.expect((try a.to(u0)) == 0);1352 testing.expect((try a.to(i0)) == 0);
1284 // TODO: sigsegv
1285 // testing.expect((try a.to(i0)) == 0);
12861353
1287 try a.set(-1);1354 try a.set(-1);
1288 testing.expect(a.bitCountTwosComp() == 1);1355 testing.expect(a.bitCountTwosComp() == 1);
std/math/big/rational.zig+55-13
...@@ -14,11 +14,22 @@ const Limb = bn.Limb;...@@ -14,11 +14,22 @@ const Limb = bn.Limb;
14const DoubleLimb = bn.DoubleLimb;14const DoubleLimb = bn.DoubleLimb;
15const Int = bn.Int;15const Int = bn.Int;
1616
17/// An arbitrary-precision rational number.
18///
19/// Memory is allocated as needed for operations to ensure full precision is kept. The precision
20/// of a Rational is only bounded by memory.
21///
22/// Rational's are always normalized. That is, for a Rational r = p/q where p and q are integers,
23/// gcd(p, q) = 1 always.
17pub const Rational = struct {24pub const Rational = struct {
18 // Sign of Rational is sign of p. Sign of q is ignored25 /// Numerator. Determines the sign of the Rational.
19 p: Int,26 p: Int,
27
28 /// Denominator. Sign is ignored.
20 q: Int,29 q: Int,
2130
31 /// Create a new Rational. A small amount of memory will be allocated on initialization.
32 /// This will be 2 * Int.default_capacity.
22 pub fn init(a: *Allocator) !Rational {33 pub fn init(a: *Allocator) !Rational {
23 return Rational{34 return Rational{
24 .p = try Int.init(a),35 .p = try Int.init(a),
...@@ -26,18 +37,21 @@ pub const Rational = struct {...@@ -26,18 +37,21 @@ pub const Rational = struct {
26 };37 };
27 }38 }
2839
40 /// Frees all memory associated with a Rational.
29 pub fn deinit(self: *Rational) void {41 pub fn deinit(self: *Rational) void {
30 self.p.deinit();42 self.p.deinit();
31 self.q.deinit();43 self.q.deinit();
32 }44 }
3345
46 /// Set a Rational from a primitive integer type.
34 pub fn setInt(self: *Rational, a: var) !void {47 pub fn setInt(self: *Rational, a: var) !void {
35 try self.p.set(a);48 try self.p.set(a);
36 try self.q.set(1);49 try self.q.set(1);
37 }50 }
3851
39 // TODO: Accept a/b fractions and exponent form52 /// Set a Rational from a string of the form `A/B` where A and B are base-10 integers.
40 pub fn setFloatString(self: *Rational, str: []const u8) !void {53 pub fn setFloatString(self: *Rational, str: []const u8) !void {
54 // TODO: Accept a/b fractions and exponent form
41 if (str.len == 0) {55 if (str.len == 0) {
42 return error.InvalidFloatString;56 return error.InvalidFloatString;
43 }57 }
...@@ -111,8 +125,10 @@ pub const Rational = struct {...@@ -111,8 +125,10 @@ pub const Rational = struct {
111 }125 }
112 }126 }
113127
114 // Translated from golang.go/src/math/big/rat.go.128 /// Set a Rational from a floating-point value. The rational will have enough precision to
129 /// completely represent the provided float.
115 pub fn setFloat(self: *Rational, comptime T: type, f: T) !void {130 pub fn setFloat(self: *Rational, comptime T: type, f: T) !void {
131 // Translated from golang.go/src/math/big/rat.go.
116 debug.assert(@typeId(T) == builtin.TypeId.Float);132 debug.assert(@typeId(T) == builtin.TypeId.Float);
117133
118 const UnsignedIntType = @IntType(false, T.bit_count);134 const UnsignedIntType = @IntType(false, T.bit_count);
...@@ -164,8 +180,13 @@ pub const Rational = struct {...@@ -164,8 +180,13 @@ pub const Rational = struct {
164 try self.reduce();180 try self.reduce();
165 }181 }
166182
167 // Translated from golang.go/src/math/big/rat.go.183 /// Return a floating-point value that is the closest value to a Rational.
184 ///
185 /// The result may not be exact if the Rational is too precise or too large for the
186 /// target type.
168 pub fn toFloat(self: Rational, comptime T: type) !T {187 pub fn toFloat(self: Rational, comptime T: type) !T {
188 // Translated from golang.go/src/math/big/rat.go.
189 // TODO: Indicate whether the result is not exact.
169 debug.assert(@typeId(T) == builtin.TypeId.Float);190 debug.assert(@typeId(T) == builtin.TypeId.Float);
170191
171 const fsize = T.bit_count;192 const fsize = T.bit_count;
...@@ -259,6 +280,7 @@ pub const Rational = struct {...@@ -259,6 +280,7 @@ pub const Rational = struct {
259 return if (self.p.isPositive()) f else -f;280 return if (self.p.isPositive()) f else -f;
260 }281 }
261282
283 /// Set a rational from an integer ratio.
262 pub fn setRatio(self: *Rational, p: var, q: var) !void {284 pub fn setRatio(self: *Rational, p: var, q: var) !void {
263 try self.p.set(p);285 try self.p.set(p);
264 try self.q.set(q);286 try self.q.set(q);
...@@ -273,11 +295,13 @@ pub const Rational = struct {...@@ -273,11 +295,13 @@ pub const Rational = struct {
273 }295 }
274 }296 }
275297
298 /// Set a Rational directly from an Int.
276 pub fn copyInt(self: *Rational, a: Int) !void {299 pub fn copyInt(self: *Rational, a: Int) !void {
277 try self.p.copy(a);300 try self.p.copy(a);
278 try self.q.set(1);301 try self.q.set(1);
279 }302 }
280303
304 /// Set a Rational directly from a ratio of two Int's.
281 pub fn copyRatio(self: *Rational, a: Int, b: Int) !void {305 pub fn copyRatio(self: *Rational, a: Int, b: Int) !void {
282 try self.p.copy(a);306 try self.p.copy(a);
283 try self.q.copy(b);307 try self.q.copy(b);
...@@ -288,23 +312,29 @@ pub const Rational = struct {...@@ -288,23 +312,29 @@ pub const Rational = struct {
288 try self.reduce();312 try self.reduce();
289 }313 }
290314
315 /// Make a Rational positive.
291 pub fn abs(r: *Rational) void {316 pub fn abs(r: *Rational) void {
292 r.p.abs();317 r.p.abs();
293 }318 }
294319
320 /// Negate the sign of a Rational.
295 pub fn negate(r: *Rational) void {321 pub fn negate(r: *Rational) void {
296 r.p.negate();322 r.p.negate();
297 }323 }
298324
325 /// Efficiently swap a Rational with another. This swaps the limb pointers and a full copy is not
326 /// performed. The address of the limbs field will not be the same after this function.
299 pub fn swap(r: *Rational, other: *Rational) void {327 pub fn swap(r: *Rational, other: *Rational) void {
300 r.p.swap(&other.p);328 r.p.swap(&other.p);
301 r.q.swap(&other.q);329 r.q.swap(&other.q);
302 }330 }
303331
332 /// Returns -1, 0, 1 if a < b, a == b or a > b respectively.
304 pub fn cmp(a: Rational, b: Rational) !i8 {333 pub fn cmp(a: Rational, b: Rational) !i8 {
305 return cmpInternal(a, b, true);334 return cmpInternal(a, b, true);
306 }335 }
307336
337 /// Returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively.
308 pub fn cmpAbs(a: Rational, b: Rational) !i8 {338 pub fn cmpAbs(a: Rational, b: Rational) !i8 {
309 return cmpInternal(a, b, false);339 return cmpInternal(a, b, false);
310 }340 }
...@@ -325,9 +355,11 @@ pub const Rational = struct {...@@ -325,9 +355,11 @@ pub const Rational = struct {
325 return if (is_abs) q.cmpAbs(p) else q.cmp(p);355 return if (is_abs) q.cmpAbs(p) else q.cmp(p);
326 }356 }
327357
328 // r/q = ap/aq + bp/bq = (ap*bq + bp*aq) / (aq*bq)358 /// rma = a + b.
329 //359 ///
330 // For best performance, rma should not alias a or b.360 /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
361 ///
362 /// Returns an error if memory could not be allocated.
331 pub fn add(rma: *Rational, a: Rational, b: Rational) !void {363 pub fn add(rma: *Rational, a: Rational, b: Rational) !void {
332 var r = rma;364 var r = rma;
333 var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;365 var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;
...@@ -351,9 +383,11 @@ pub const Rational = struct {...@@ -351,9 +383,11 @@ pub const Rational = struct {
351 try r.reduce();383 try r.reduce();
352 }384 }
353385
354 // r/q = ap/aq - bp/bq = (ap*bq - bp*aq) / (aq*bq)386 /// rma = a - b.
355 //387 ///
356 // For best performance, rma should not alias a or b.388 /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
389 ///
390 /// Returns an error if memory could not be allocated.
357 pub fn sub(rma: *Rational, a: Rational, b: Rational) !void {391 pub fn sub(rma: *Rational, a: Rational, b: Rational) !void {
358 var r = rma;392 var r = rma;
359 var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;393 var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;
...@@ -377,14 +411,22 @@ pub const Rational = struct {...@@ -377,14 +411,22 @@ pub const Rational = struct {
377 try r.reduce();411 try r.reduce();
378 }412 }
379413
380 // r/q = ap/aq * bp/bq = ap*bp / aq*bq414 /// rma = a * b.
415 ///
416 /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
417 ///
418 /// Returns an error if memory could not be allocated.
381 pub fn mul(r: *Rational, a: Rational, b: Rational) !void {419 pub fn mul(r: *Rational, a: Rational, b: Rational) !void {
382 try r.p.mul(a.p, b.p);420 try r.p.mul(a.p, b.p);
383 try r.q.mul(a.q, b.q);421 try r.q.mul(a.q, b.q);
384 try r.reduce();422 try r.reduce();
385 }423 }
386424
387 // r/q = (ap/aq) / (bp/bq) = ap*bq / bp*aq425 /// rma = a / b.
426 ///
427 /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
428 ///
429 /// Returns an error if memory could not be allocated.
388 pub fn div(r: *Rational, a: Rational, b: Rational) !void {430 pub fn div(r: *Rational, a: Rational, b: Rational) !void {
389 if (b.p.eqZero()) {431 if (b.p.eqZero()) {
390 @panic("division by zero");432 @panic("division by zero");
...@@ -395,7 +437,7 @@ pub const Rational = struct {...@@ -395,7 +437,7 @@ pub const Rational = struct {
395 try r.reduce();437 try r.reduce();
396 }438 }
397439
398 // r/q = q/r440 /// Invert the numerator and denominator fields of a Rational. p/q => q/p.
399 pub fn invert(r: *Rational) void {441 pub fn invert(r: *Rational) void {
400 Int.swap(&r.p, &r.q);442 Int.swap(&r.p, &r.q);
401 }443 }
std/math/cbrt.zig+10-4
...@@ -1,13 +1,19 @@...@@ -1,13 +1,19 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - cbrt(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/cbrtf.c
4// - cbrt(+-inf) = +-inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/cbrt.c
5// - cbrt(nan) = nan
66
7const std = @import("../std.zig");7const std = @import("../std.zig");
8const math = std.math;8const math = std.math;
9const expect = std.testing.expect;9const expect = std.testing.expect;
1010
11/// Returns the cube root of x.
12///
13/// Special Cases:
14/// - cbrt(+-0) = +-0
15/// - cbrt(+-inf) = +-inf
16/// - cbrt(nan) = nan
11pub fn cbrt(x: var) @typeOf(x) {17pub fn cbrt(x: var) @typeOf(x) {
12 const T = @typeOf(x);18 const T = @typeOf(x);
13 return switch (T) {19 return switch (T) {
std/math/ceil.zig+10-4
...@@ -1,14 +1,20 @@...@@ -1,14 +1,20 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - ceil(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/ceilf.c
4// - ceil(+-inf) = +-inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/ceil.c
5// - ceil(nan) = nan
66
7const builtin = @import("builtin");7const builtin = @import("builtin");
8const std = @import("../std.zig");8const std = @import("../std.zig");
9const math = std.math;9const math = std.math;
10const expect = std.testing.expect;10const expect = std.testing.expect;
1111
12/// Returns the least integer value greater than of equal to x.
13///
14/// Special Cases:
15/// - ceil(+-0) = +-0
16/// - ceil(+-inf) = +-inf
17/// - ceil(nan) = nan
12pub fn ceil(x: var) @typeOf(x) {18pub fn ceil(x: var) @typeOf(x) {
13 const T = @typeOf(x);19 const T = @typeOf(x);
14 return switch (T) {20 return switch (T) {
std/math/complex.zig+12
...@@ -23,13 +23,18 @@ pub const sqrt = @import("complex/sqrt.zig").sqrt;...@@ -23,13 +23,18 @@ pub const sqrt = @import("complex/sqrt.zig").sqrt;
23pub const tanh = @import("complex/tanh.zig").tanh;23pub const tanh = @import("complex/tanh.zig").tanh;
24pub const tan = @import("complex/tan.zig").tan;24pub const tan = @import("complex/tan.zig").tan;
2525
26/// A complex number consisting of a real an imaginary part. T must be a floating-point value.
26pub fn Complex(comptime T: type) type {27pub fn Complex(comptime T: type) type {
27 return struct {28 return struct {
28 const Self = @This();29 const Self = @This();
2930
31 /// Real part.
30 re: T,32 re: T,
33
34 /// Imaginary part.
31 im: T,35 im: T,
3236
37 /// Create a new Complex number from the given real and imaginary parts.
33 pub fn new(re: T, im: T) Self {38 pub fn new(re: T, im: T) Self {
34 return Self{39 return Self{
35 .re = re,40 .re = re,
...@@ -37,6 +42,7 @@ pub fn Complex(comptime T: type) type {...@@ -37,6 +42,7 @@ pub fn Complex(comptime T: type) type {
37 };42 };
38 }43 }
3944
45 /// Returns the sum of two complex numbers.
40 pub fn add(self: Self, other: Self) Self {46 pub fn add(self: Self, other: Self) Self {
41 return Self{47 return Self{
42 .re = self.re + other.re,48 .re = self.re + other.re,
...@@ -44,6 +50,7 @@ pub fn Complex(comptime T: type) type {...@@ -44,6 +50,7 @@ pub fn Complex(comptime T: type) type {
44 };50 };
45 }51 }
4652
53 /// Returns the subtraction of two complex numbers.
47 pub fn sub(self: Self, other: Self) Self {54 pub fn sub(self: Self, other: Self) Self {
48 return Self{55 return Self{
49 .re = self.re - other.re,56 .re = self.re - other.re,
...@@ -51,6 +58,7 @@ pub fn Complex(comptime T: type) type {...@@ -51,6 +58,7 @@ pub fn Complex(comptime T: type) type {
51 };58 };
52 }59 }
5360
61 /// Returns the product of two complex numbers.
54 pub fn mul(self: Self, other: Self) Self {62 pub fn mul(self: Self, other: Self) Self {
55 return Self{63 return Self{
56 .re = self.re * other.re - self.im * other.im,64 .re = self.re * other.re - self.im * other.im,
...@@ -58,6 +66,7 @@ pub fn Complex(comptime T: type) type {...@@ -58,6 +66,7 @@ pub fn Complex(comptime T: type) type {
58 };66 };
59 }67 }
6068
69 /// Returns the quotient of two complex numbers.
61 pub fn div(self: Self, other: Self) Self {70 pub fn div(self: Self, other: Self) Self {
62 const re_num = self.re * other.re + self.im * other.im;71 const re_num = self.re * other.re + self.im * other.im;
63 const im_num = self.im * other.re - self.re * other.im;72 const im_num = self.im * other.re - self.re * other.im;
...@@ -69,6 +78,7 @@ pub fn Complex(comptime T: type) type {...@@ -69,6 +78,7 @@ pub fn Complex(comptime T: type) type {
69 };78 };
70 }79 }
7180
81 /// Returns the complex conjugate of a number.
72 pub fn conjugate(self: Self) Self {82 pub fn conjugate(self: Self) Self {
73 return Self{83 return Self{
74 .re = self.re,84 .re = self.re,
...@@ -76,6 +86,7 @@ pub fn Complex(comptime T: type) type {...@@ -76,6 +86,7 @@ pub fn Complex(comptime T: type) type {
76 };86 };
77 }87 }
7888
89 /// Returns the reciprocal of a complex number.
79 pub fn reciprocal(self: Self) Self {90 pub fn reciprocal(self: Self) Self {
80 const m = self.re * self.re + self.im * self.im;91 const m = self.re * self.re + self.im * self.im;
81 return Self{92 return Self{
...@@ -84,6 +95,7 @@ pub fn Complex(comptime T: type) type {...@@ -84,6 +95,7 @@ pub fn Complex(comptime T: type) type {
84 };95 };
85 }96 }
8697
98 /// Returns the magnitude of a complex number.
87 pub fn magnitude(self: Self) T {99 pub fn magnitude(self: Self) T {
88 return math.sqrt(self.re * self.re + self.im * self.im);100 return math.sqrt(self.re * self.re + self.im * self.im);
89 }101 }
std/math/complex/abs.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the absolute value (modulus) of z.
7pub fn abs(z: var) @typeOf(z.re) {8pub fn abs(z: var) @typeOf(z.re) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 return math.hypot(T, z.re, z.im);10 return math.hypot(T, z.re, z.im);
std/math/complex/acos.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the arc-cosine of z.
7pub fn acos(z: var) Complex(@typeOf(z.re)) {8pub fn acos(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const q = cmath.asin(z);10 const q = cmath.asin(z);
std/math/complex/acosh.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the hyperbolic arc-cosine of z.
7pub fn acosh(z: var) Complex(@typeOf(z.re)) {8pub fn acosh(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const q = cmath.acos(z);10 const q = cmath.acos(z);
std/math/complex/arg.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the angular component (in radians) of z.
7pub fn arg(z: var) @typeOf(z.re) {8pub fn arg(z: var) @typeOf(z.re) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 return math.atan2(T, z.im, z.re);10 return math.atan2(T, z.im, z.re);
std/math/complex/asin.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7// Returns the arc-sine of z.
7pub fn asin(z: var) Complex(@typeOf(z.re)) {8pub fn asin(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const x = z.re;10 const x = z.re;
std/math/complex/asinh.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the hyperbolic arc-sine of z.
7pub fn asinh(z: var) Complex(@typeOf(z.re)) {8pub fn asinh(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const q = Complex(T).new(-z.im, z.re);10 const q = Complex(T).new(-z.im, z.re);
std/math/complex/atan.zig+7
...@@ -1,9 +1,16 @@...@@ -1,9 +1,16 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/complex/catanf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/catan.c
6
1const std = @import("../../std.zig");7const std = @import("../../std.zig");
2const testing = std.testing;8const testing = std.testing;
3const math = std.math;9const math = std.math;
4const cmath = math.complex;10const cmath = math.complex;
5const Complex = cmath.Complex;11const Complex = cmath.Complex;
612
13/// Returns the arc-tangent of z.
7pub fn atan(z: var) @typeOf(z) {14pub fn atan(z: var) @typeOf(z) {
8 const T = @typeOf(z.re);15 const T = @typeOf(z.re);
9 return switch (T) {16 return switch (T) {
std/math/complex/atanh.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the hyperbolic arc-tangent of z.
7pub fn atanh(z: var) Complex(@typeOf(z.re)) {8pub fn atanh(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const q = Complex(T).new(-z.im, z.re);10 const q = Complex(T).new(-z.im, z.re);
std/math/complex/conj.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the complex conjugate of z.
7pub fn conj(z: var) Complex(@typeOf(z.re)) {8pub fn conj(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 return Complex(T).new(z.re, -z.im);10 return Complex(T).new(z.re, -z.im);
std/math/complex/cos.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the cosine of z.
7pub fn cos(z: var) Complex(@typeOf(z.re)) {8pub fn cos(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const p = Complex(T).new(-z.im, z.re);10 const p = Complex(T).new(-z.im, z.re);
std/math/complex/cosh.zig+7
...@@ -1,3 +1,9 @@...@@ -1,3 +1,9 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/complex/ccoshf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/ccosh.c
6
1const std = @import("../../std.zig");7const std = @import("../../std.zig");
2const testing = std.testing;8const testing = std.testing;
3const math = std.math;9const math = std.math;
...@@ -6,6 +12,7 @@ const Complex = cmath.Complex;...@@ -6,6 +12,7 @@ const Complex = cmath.Complex;
612
7const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;13const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
814
15/// Returns the hyperbolic arc-cosine of z.
9pub fn cosh(z: var) Complex(@typeOf(z.re)) {16pub fn cosh(z: var) Complex(@typeOf(z.re)) {
10 const T = @typeOf(z.re);17 const T = @typeOf(z.re);
11 return switch (T) {18 return switch (T) {
std/math/complex/exp.zig+7
...@@ -1,3 +1,9 @@...@@ -1,3 +1,9 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/complex/cexpf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/cexp.c
6
1const std = @import("../../std.zig");7const std = @import("../../std.zig");
2const testing = std.testing;8const testing = std.testing;
3const math = std.math;9const math = std.math;
...@@ -6,6 +12,7 @@ const Complex = cmath.Complex;...@@ -6,6 +12,7 @@ const Complex = cmath.Complex;
612
7const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;13const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
814
15/// Returns e raised to the power of z (e^z).
9pub fn exp(z: var) @typeOf(z) {16pub fn exp(z: var) @typeOf(z) {
10 const T = @typeOf(z.re);17 const T = @typeOf(z.re);
1118
std/math/complex/ldexp.zig+7
...@@ -1,9 +1,16 @@...@@ -1,9 +1,16 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/complex/__cexpf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/__cexp.c
6
1const std = @import("../../std.zig");7const std = @import("../../std.zig");
2const debug = std.debug;8const debug = std.debug;
3const math = std.math;9const math = std.math;
4const cmath = math.complex;10const cmath = math.complex;
5const Complex = cmath.Complex;11const Complex = cmath.Complex;
612
13/// Returns exp(z) scaled to avoid overflow.
7pub fn ldexp_cexp(z: var, expt: i32) @typeOf(z) {14pub fn ldexp_cexp(z: var, expt: i32) @typeOf(z) {
8 const T = @typeOf(z.re);15 const T = @typeOf(z.re);
916
std/math/complex/log.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the natural logarithm of z.
7pub fn log(z: var) Complex(@typeOf(z.re)) {8pub fn log(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const r = cmath.abs(z);10 const r = cmath.abs(z);
std/math/complex/pow.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns z raised to the complex power of c.
7pub fn pow(comptime T: type, z: T, c: T) T {8pub fn pow(comptime T: type, z: T, c: T) T {
8 const p = cmath.log(z);9 const p = cmath.log(z);
9 const q = c.mul(p);10 const q = c.mul(p);
std/math/complex/proj.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the projection of z onto the riemann sphere.
7pub fn proj(z: var) Complex(@typeOf(z.re)) {8pub fn proj(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
910
std/math/complex/sin.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the sine of z.
7pub fn sin(z: var) Complex(@typeOf(z.re)) {8pub fn sin(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const p = Complex(T).new(-z.im, z.re);10 const p = Complex(T).new(-z.im, z.re);
std/math/complex/sinh.zig+7
...@@ -1,3 +1,9 @@...@@ -1,3 +1,9 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/complex/csinhf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/csinh.c
6
1const std = @import("../../std.zig");7const std = @import("../../std.zig");
2const testing = std.testing;8const testing = std.testing;
3const math = std.math;9const math = std.math;
...@@ -6,6 +12,7 @@ const Complex = cmath.Complex;...@@ -6,6 +12,7 @@ const Complex = cmath.Complex;
612
7const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;13const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
814
15/// Returns the hyperbolic sine of z.
9pub fn sinh(z: var) @typeOf(z) {16pub fn sinh(z: var) @typeOf(z) {
10 const T = @typeOf(z.re);17 const T = @typeOf(z.re);
11 return switch (T) {18 return switch (T) {
std/math/complex/sqrt.zig+8
...@@ -1,9 +1,17 @@...@@ -1,9 +1,17 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/complex/csqrtf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/csqrt.c
6
1const std = @import("../../std.zig");7const std = @import("../../std.zig");
2const testing = std.testing;8const testing = std.testing;
3const math = std.math;9const math = std.math;
4const cmath = math.complex;10const cmath = math.complex;
5const Complex = cmath.Complex;11const Complex = cmath.Complex;
612
13/// Returns the square root of z. The real and imaginary parts of the result have the same sign
14/// as the imaginary part of z.
7pub fn sqrt(z: var) @typeOf(z) {15pub fn sqrt(z: var) @typeOf(z) {
8 const T = @typeOf(z.re);16 const T = @typeOf(z.re);
917
std/math/complex/tan.zig+1
...@@ -4,6 +4,7 @@ const math = std.math;...@@ -4,6 +4,7 @@ const math = std.math;
4const cmath = math.complex;4const cmath = math.complex;
5const Complex = cmath.Complex;5const Complex = cmath.Complex;
66
7/// Returns the tanget of z.
7pub fn tan(z: var) Complex(@typeOf(z.re)) {8pub fn tan(z: var) Complex(@typeOf(z.re)) {
8 const T = @typeOf(z.re);9 const T = @typeOf(z.re);
9 const q = Complex(T).new(-z.im, z.re);10 const q = Complex(T).new(-z.im, z.re);
std/math/complex/tanh.zig+7
...@@ -1,9 +1,16 @@...@@ -1,9 +1,16 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/complex/ctanhf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/ctanh.c
6
1const std = @import("../../std.zig");7const std = @import("../../std.zig");
2const testing = std.testing;8const testing = std.testing;
3const math = std.math;9const math = std.math;
4const cmath = math.complex;10const cmath = math.complex;
5const Complex = cmath.Complex;11const Complex = cmath.Complex;
612
13/// Returns the hyperbolic tangent of z.
7pub fn tanh(z: var) @typeOf(z) {14pub fn tanh(z: var) @typeOf(z) {
8 const T = @typeOf(z.re);15 const T = @typeOf(z.re);
9 return switch (T) {16 return switch (T) {
std/math/copysign.zig+7
...@@ -1,8 +1,15 @@...@@ -1,8 +1,15 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/copysignf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/copysign.c
6
1const std = @import("../std.zig");7const std = @import("../std.zig");
2const math = std.math;8const math = std.math;
3const expect = std.testing.expect;9const expect = std.testing.expect;
4const maxInt = std.math.maxInt;10const maxInt = std.math.maxInt;
511
12/// Returns a value with the magnitude of x and the sign of y.
6pub fn copysign(comptime T: type, x: T, y: T) T {13pub fn copysign(comptime T: type, x: T, y: T) T {
7 return switch (T) {14 return switch (T) {
8 f16 => copysign16(x, y),15 f16 => copysign16(x, y),
std/math/cos.zig+46-100
...@@ -1,18 +1,23 @@...@@ -1,18 +1,23 @@
1// Special Cases:1// Ported from go, which is licensed under a BSD-3 license.
2// https://golang.org/LICENSE
2//3//
3// - cos(+-inf) = nan4// https://golang.org/src/math/sin.go
4// - cos(nan) = nan
55
6const builtin = @import("builtin");6const builtin = @import("builtin");
7const std = @import("../std.zig");7const std = @import("../std.zig");
8const math = std.math;8const math = std.math;
9const expect = std.testing.expect;9const expect = std.testing.expect;
1010
11/// Returns the cosine of the radian value x.
12///
13/// Special Cases:
14/// - cos(+-inf) = nan
15/// - cos(nan) = nan
11pub fn cos(x: var) @typeOf(x) {16pub fn cos(x: var) @typeOf(x) {
12 const T = @typeOf(x);17 const T = @typeOf(x);
13 return switch (T) {18 return switch (T) {
14 f32 => cos32(x),19 f32 => cos_(f32, x),
15 f64 => cos64(x),20 f64 => cos_(f64, x),
16 else => @compileError("cos not implemented for " ++ @typeName(T)),21 else => @compileError("cos not implemented for " ++ @typeName(T)),
17 };22 };
18}23}
...@@ -33,78 +38,24 @@ const C3 = 2.48015872888517045348E-5;...@@ -33,78 +38,24 @@ const C3 = 2.48015872888517045348E-5;
33const C4 = -1.38888888888730564116E-3;38const C4 = -1.38888888888730564116E-3;
34const C5 = 4.16666666666665929218E-2;39const C5 = 4.16666666666665929218E-2;
3540
36// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.41const pi4a = 7.85398125648498535156e-1;
37//42const pi4b = 3.77489470793079817668E-8;
38// This may have slight differences on some edge cases and may need to replaced if so.43const pi4c = 2.69515142907905952645E-15;
39fn cos32(x_: f32) f32 {44const m4pi = 1.273239544735162542821171882678754627704620361328125;
40 const pi4a = 7.85398125648498535156e-1;
41 const pi4b = 3.77489470793079817668E-8;
42 const pi4c = 2.69515142907905952645E-15;
43 const m4pi = 1.273239544735162542821171882678754627704620361328125;
44
45 var x = x_;
46 if (math.isNan(x) or math.isInf(x)) {
47 return math.nan(f32);
48 }
49
50 var sign = false;
51 if (x < 0) {
52 x = -x;
53 }
54
55 var y = math.floor(x * m4pi);
56 var j = @floatToInt(i64, y);
5745
58 if (j & 1 == 1) {46fn cos_(comptime T: type, x_: T) T {
59 j += 1;47 const I = @IntType(true, T.bit_count);
60 y += 1;
61 }
62
63 j &= 7;
64 if (j > 3) {
65 j -= 4;
66 sign = !sign;
67 }
68 if (j > 1) {
69 sign = !sign;
70 }
71
72 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
73 const w = z * z;
74
75 const r = r: {
76 if (j == 1 or j == 2) {
77 break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
78 } else {
79 break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
80 }
81 };
82
83 if (sign) {
84 return -r;
85 } else {
86 return r;
87 }
88}
89
90fn cos64(x_: f64) f64 {
91 const pi4a = 7.85398125648498535156e-1;
92 const pi4b = 3.77489470793079817668E-8;
93 const pi4c = 2.69515142907905952645E-15;
94 const m4pi = 1.273239544735162542821171882678754627704620361328125;
9548
96 var x = x_;49 var x = x_;
97 if (math.isNan(x) or math.isInf(x)) {50 if (math.isNan(x) or math.isInf(x)) {
98 return math.nan(f64);51 return math.nan(T);
99 }52 }
10053
101 var sign = false;54 var sign = false;
102 if (x < 0) {55 x = math.fabs(x);
103 x = -x;
104 }
10556
106 var y = math.floor(x * m4pi);57 var y = math.floor(x * m4pi);
107 var j = @floatToInt(i64, y);58 var j = @floatToInt(I, y);
10859
109 if (j & 1 == 1) {60 if (j & 1 == 1) {
110 j += 1;61 j += 1;
...@@ -123,56 +74,51 @@ fn cos64(x_: f64) f64 {...@@ -123,56 +74,51 @@ fn cos64(x_: f64) f64 {
123 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;74 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
124 const w = z * z;75 const w = z * z;
12576
126 const r = r: {77 const r = if (j == 1 or j == 2)
127 if (j == 1 or j == 2) {78 z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))))
128 break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));79 else
129 } else {80 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
130 break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
131 }
132 };
13381
134 if (sign) {82 return if (sign) -r else r;
135 return -r;
136 } else {
137 return r;
138 }
139}83}
14084
141test "math.cos" {85test "math.cos" {
142 expect(cos(f32(0.0)) == cos32(0.0));86 expect(cos(f32(0.0)) == cos_(f32, 0.0));
143 expect(cos(f64(0.0)) == cos64(0.0));87 expect(cos(f64(0.0)) == cos_(f64, 0.0));
144}88}
14589
146test "math.cos32" {90test "math.cos32" {
147 const epsilon = 0.000001;91 const epsilon = 0.000001;
14892
149 expect(math.approxEq(f32, cos32(0.0), 1.0, epsilon));93 expect(math.approxEq(f32, cos_(f32, 0.0), 1.0, epsilon));
150 expect(math.approxEq(f32, cos32(0.2), 0.980067, epsilon));94 expect(math.approxEq(f32, cos_(f32, 0.2), 0.980067, epsilon));
151 expect(math.approxEq(f32, cos32(0.8923), 0.627623, epsilon));95 expect(math.approxEq(f32, cos_(f32, 0.8923), 0.627623, epsilon));
152 expect(math.approxEq(f32, cos32(1.5), 0.070737, epsilon));96 expect(math.approxEq(f32, cos_(f32, 1.5), 0.070737, epsilon));
153 expect(math.approxEq(f32, cos32(37.45), 0.969132, epsilon));97 expect(math.approxEq(f32, cos_(f32, -1.5), 0.070737, epsilon));
154 expect(math.approxEq(f32, cos32(89.123), 0.400798, epsilon));98 expect(math.approxEq(f32, cos_(f32, 37.45), 0.969132, epsilon));
99 expect(math.approxEq(f32, cos_(f32, 89.123), 0.400798, epsilon));
155}100}
156101
157test "math.cos64" {102test "math.cos64" {
158 const epsilon = 0.000001;103 const epsilon = 0.000001;
159104
160 expect(math.approxEq(f64, cos64(0.0), 1.0, epsilon));105 expect(math.approxEq(f64, cos_(f64, 0.0), 1.0, epsilon));
161 expect(math.approxEq(f64, cos64(0.2), 0.980067, epsilon));106 expect(math.approxEq(f64, cos_(f64, 0.2), 0.980067, epsilon));
162 expect(math.approxEq(f64, cos64(0.8923), 0.627623, epsilon));107 expect(math.approxEq(f64, cos_(f64, 0.8923), 0.627623, epsilon));
163 expect(math.approxEq(f64, cos64(1.5), 0.070737, epsilon));108 expect(math.approxEq(f64, cos_(f64, 1.5), 0.070737, epsilon));
164 expect(math.approxEq(f64, cos64(37.45), 0.969132, epsilon));109 expect(math.approxEq(f64, cos_(f64, -1.5), 0.070737, epsilon));
165 expect(math.approxEq(f64, cos64(89.123), 0.40080, epsilon));110 expect(math.approxEq(f64, cos_(f64, 37.45), 0.969132, epsilon));
111 expect(math.approxEq(f64, cos_(f64, 89.123), 0.40080, epsilon));
166}112}
167113
168test "math.cos32.special" {114test "math.cos32.special" {
169 expect(math.isNan(cos32(math.inf(f32))));115 expect(math.isNan(cos_(f32, math.inf(f32))));
170 expect(math.isNan(cos32(-math.inf(f32))));116 expect(math.isNan(cos_(f32, -math.inf(f32))));
171 expect(math.isNan(cos32(math.nan(f32))));117 expect(math.isNan(cos_(f32, math.nan(f32))));
172}118}
173119
174test "math.cos64.special" {120test "math.cos64.special" {
175 expect(math.isNan(cos64(math.inf(f64))));121 expect(math.isNan(cos_(f64, math.inf(f64))));
176 expect(math.isNan(cos64(-math.inf(f64))));122 expect(math.isNan(cos_(f64, -math.inf(f64))));
177 expect(math.isNan(cos64(math.nan(f64))));123 expect(math.isNan(cos_(f64, math.nan(f64))));
178}124}
std/math/cosh.zig+10-4
...@@ -1,8 +1,8 @@...@@ -1,8 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - cosh(+-0) = 14// https://git.musl-libc.org/cgit/musl/tree/src/math/coshf.c
4// - cosh(+-inf) = +inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/cosh.c
5// - cosh(nan) = nan
66
7const builtin = @import("builtin");7const builtin = @import("builtin");
8const std = @import("../std.zig");8const std = @import("../std.zig");
...@@ -11,6 +11,12 @@ const expo2 = @import("expo2.zig").expo2;...@@ -11,6 +11,12 @@ const expo2 = @import("expo2.zig").expo2;
11const expect = std.testing.expect;11const expect = std.testing.expect;
12const maxInt = std.math.maxInt;12const maxInt = std.math.maxInt;
1313
14/// Returns the hyperbolic cosine of x.
15///
16/// Special Cases:
17/// - cosh(+-0) = 1
18/// - cosh(+-inf) = +inf
19/// - cosh(nan) = nan
14pub fn cosh(x: var) @typeOf(x) {20pub fn cosh(x: var) @typeOf(x) {
15 const T = @typeOf(x);21 const T = @typeOf(x);
16 return switch (T) {22 return switch (T) {
std/math/exp.zig+9-3
...@@ -1,13 +1,19 @@...@@ -1,13 +1,19 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - exp(+inf) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/expf.c
4// - exp(nan) = nan5// https://git.musl-libc.org/cgit/musl/tree/src/math/exp.c
56
6const std = @import("../std.zig");7const std = @import("../std.zig");
7const math = std.math;8const math = std.math;
8const assert = std.debug.assert;9const assert = std.debug.assert;
9const builtin = @import("builtin");10const builtin = @import("builtin");
1011
12/// Returns e raised to the power of x (e^x).
13///
14/// Special Cases:
15/// - exp(+inf) = +inf
16/// - exp(nan) = nan
11pub fn exp(x: var) @typeOf(x) {17pub fn exp(x: var) @typeOf(x) {
12 const T = @typeOf(x);18 const T = @typeOf(x);
13 return switch (T) {19 return switch (T) {
std/math/exp2.zig+9-3
...@@ -1,12 +1,18 @@...@@ -1,12 +1,18 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - exp2(+inf) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/exp2f.c
4// - exp2(nan) = nan5// https://git.musl-libc.org/cgit/musl/tree/src/math/exp2.c
56
6const std = @import("../std.zig");7const std = @import("../std.zig");
7const math = std.math;8const math = std.math;
8const expect = std.testing.expect;9const expect = std.testing.expect;
910
11/// Returns 2 raised to the power of x (2^x).
12///
13/// Special Cases:
14/// - exp2(+inf) = +inf
15/// - exp2(nan) = nan
10pub fn exp2(x: var) @typeOf(x) {16pub fn exp2(x: var) @typeOf(x) {
11 const T = @typeOf(x);17 const T = @typeOf(x);
12 return switch (T) {18 return switch (T) {
std/math/expm1.zig+13-4
...@@ -1,14 +1,23 @@...@@ -1,14 +1,23 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - expm1(+inf) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/expmf.c
4// - expm1(-inf) = -15// https://git.musl-libc.org/cgit/musl/tree/src/math/expm.c
5// - expm1(nan) = nan6
7// TODO: Updated recently.
68
7const builtin = @import("builtin");9const builtin = @import("builtin");
8const std = @import("../std.zig");10const std = @import("../std.zig");
9const math = std.math;11const math = std.math;
10const expect = std.testing.expect;12const expect = std.testing.expect;
1113
14/// Returns e raised to the power of x, minus 1 (e^x - 1). This is more accurate than exp(e, x) - 1
15/// when x is near 0.
16///
17/// Special Cases:
18/// - expm1(+inf) = +inf
19/// - expm1(-inf) = -1
20/// - expm1(nan) = nan
12pub fn expm1(x: var) @typeOf(x) {21pub fn expm1(x: var) @typeOf(x) {
13 const T = @typeOf(x);22 const T = @typeOf(x);
14 return switch (T) {23 return switch (T) {
std/math/expo2.zig+7
...@@ -1,5 +1,12 @@...@@ -1,5 +1,12 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/__expo2f.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/__expo2.c
6
1const math = @import("../math.zig");7const math = @import("../math.zig");
28
9/// Returns exp(x) / 2 for x >= log(maxFloat(T)).
3pub fn expo2(x: var) @typeOf(x) {10pub fn expo2(x: var) @typeOf(x) {
4 const T = @typeOf(x);11 const T = @typeOf(x);
5 return switch (T) {12 return switch (T) {
std/math/fabs.zig+9-3
...@@ -1,13 +1,19 @@...@@ -1,13 +1,19 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - fabs(+-inf) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/fabsf.c
4// - fabs(nan) = nan5// https://git.musl-libc.org/cgit/musl/tree/src/math/fabs.c
56
6const std = @import("../std.zig");7const std = @import("../std.zig");
7const math = std.math;8const math = std.math;
8const expect = std.testing.expect;9const expect = std.testing.expect;
9const maxInt = std.math.maxInt;10const maxInt = std.math.maxInt;
1011
12/// Returns the absolute value of x.
13///
14/// Special Cases:
15/// - fabs(+-inf) = +inf
16/// - fabs(nan) = nan
11pub fn fabs(x: var) @typeOf(x) {17pub fn fabs(x: var) @typeOf(x) {
12 const T = @typeOf(x);18 const T = @typeOf(x);
13 return switch (T) {19 return switch (T) {
std/math/floor.zig+10-4
...@@ -1,14 +1,20 @@...@@ -1,14 +1,20 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - floor(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/floorf.c
4// - floor(+-inf) = +-inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/floor.c
5// - floor(nan) = nan
66
7const builtin = @import("builtin");7const builtin = @import("builtin");
8const expect = std.testing.expect;8const expect = std.testing.expect;
9const std = @import("../std.zig");9const std = @import("../std.zig");
10const math = std.math;10const math = std.math;
1111
12/// Returns the greatest integer value less than or equal to x.
13///
14/// Special Cases:
15/// - floor(+-0) = +-0
16/// - floor(+-inf) = +-inf
17/// - floor(nan) = nan
12pub fn floor(x: var) @typeOf(x) {18pub fn floor(x: var) @typeOf(x) {
13 const T = @typeOf(x);19 const T = @typeOf(x);
14 return switch (T) {20 return switch (T) {
std/math/fma.zig+9-1
...@@ -1,7 +1,14 @@...@@ -1,7 +1,14 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c
6
1const std = @import("../std.zig");7const std = @import("../std.zig");
2const math = std.math;8const math = std.math;
3const expect = std.testing.expect;9const expect = std.testing.expect;
410
11/// Returns x * y + z with a single rounding error.
5pub fn fma(comptime T: type, x: T, y: T, z: T) T {12pub fn fma(comptime T: type, x: T, y: T, z: T) T {
6 return switch (T) {13 return switch (T) {
7 f32 => fma32(x, y, z),14 f32 => fma32(x, y, z),
...@@ -16,7 +23,7 @@ fn fma32(x: f32, y: f32, z: f32) f32 {...@@ -16,7 +23,7 @@ fn fma32(x: f32, y: f32, z: f32) f32 {
16 const u = @bitCast(u64, xy_z);23 const u = @bitCast(u64, xy_z);
17 const e = (u >> 52) & 0x7FF;24 const e = (u >> 52) & 0x7FF;
1825
19 if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or xy_z - xy == z) {26 if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or (xy_z - xy == z and xy_z - z == xy)) {
20 return @floatCast(f32, xy_z);27 return @floatCast(f32, xy_z);
21 } else {28 } else {
22 // TODO: Handle inexact case with double-rounding29 // TODO: Handle inexact case with double-rounding
...@@ -24,6 +31,7 @@ fn fma32(x: f32, y: f32, z: f32) f32 {...@@ -24,6 +31,7 @@ fn fma32(x: f32, y: f32, z: f32) f32 {
24 }31 }
25}32}
2633
34// NOTE: Upstream fma.c has been rewritten completely to raise fp exceptions more accurately.
27fn fma64(x: f64, y: f64, z: f64) f64 {35fn fma64(x: f64, y: f64, z: f64) f64 {
28 if (!math.isFinite(x) or !math.isFinite(y)) {36 if (!math.isFinite(x) or !math.isFinite(y)) {
29 return x * y + z;37 return x * y + z;
std/math/frexp.zig+11-4
...@@ -1,8 +1,8 @@...@@ -1,8 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - frexp(+-0) = +-0, 04// https://git.musl-libc.org/cgit/musl/tree/src/math/frexpf.c
4// - frexp(+-inf) = +-inf, 05// https://git.musl-libc.org/cgit/musl/tree/src/math/frexp.c
5// - frexp(nan) = nan, undefined
66
7const std = @import("../std.zig");7const std = @import("../std.zig");
8const math = std.math;8const math = std.math;
...@@ -17,6 +17,13 @@ fn frexp_result(comptime T: type) type {...@@ -17,6 +17,13 @@ fn frexp_result(comptime T: type) type {
17pub const frexp32_result = frexp_result(f32);17pub const frexp32_result = frexp_result(f32);
18pub const frexp64_result = frexp_result(f64);18pub const frexp64_result = frexp_result(f64);
1919
20/// Breaks x into a normalized fraction and an integral power of two.
21/// f == frac * 2^exp, with |frac| in the interval [0.5, 1).
22///
23/// Special Cases:
24/// - frexp(+-0) = +-0, 0
25/// - frexp(+-inf) = +-inf, 0
26/// - frexp(nan) = nan, undefined
20pub fn frexp(x: var) frexp_result(@typeOf(x)) {27pub fn frexp(x: var) frexp_result(@typeOf(x)) {
21 const T = @typeOf(x);28 const T = @typeOf(x);
22 return switch (T) {29 return switch (T) {
std/math/hypot.zig+11-5
...@@ -1,15 +1,21 @@...@@ -1,15 +1,21 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - hypot(+-inf, y) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/hypotf.c
4// - hypot(x, +-inf) = +inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/hypot.c
5// - hypot(nan, y) = nan
6// - hypot(x, nan) = nan
76
8const std = @import("../std.zig");7const std = @import("../std.zig");
9const math = std.math;8const math = std.math;
10const expect = std.testing.expect;9const expect = std.testing.expect;
11const maxInt = std.math.maxInt;10const maxInt = std.math.maxInt;
1211
12/// Returns sqrt(x * x + y * y), avoiding unncessary overflow and underflow.
13///
14/// Special Cases:
15/// - hypot(+-inf, y) = +inf
16/// - hypot(x, +-inf) = +inf
17/// - hypot(nan, y) = nan
18/// - hypot(x, nan) = nan
13pub fn hypot(comptime T: type, x: T, y: T) T {19pub fn hypot(comptime T: type, x: T, y: T) T {
14 return switch (T) {20 return switch (T) {
15 f32 => hypot32(x, y),21 f32 => hypot32(x, y),
std/math/ilogb.zig+10-4
...@@ -1,8 +1,8 @@...@@ -1,8 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - ilogb(+-inf) = maxInt(i32)4// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogbf.c
4// - ilogb(0) = maxInt(i32)5// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogb.c
5// - ilogb(nan) = maxInt(i32)
66
7const std = @import("../std.zig");7const std = @import("../std.zig");
8const math = std.math;8const math = std.math;
...@@ -10,6 +10,12 @@ const expect = std.testing.expect;...@@ -10,6 +10,12 @@ const expect = std.testing.expect;
10const maxInt = std.math.maxInt;10const maxInt = std.math.maxInt;
11const minInt = std.math.minInt;11const minInt = std.math.minInt;
1212
13/// Returns the binary exponent of x as an integer.
14///
15/// Special Cases:
16/// - ilogb(+-inf) = maxInt(i32)
17/// - ilogb(0) = maxInt(i32)
18/// - ilogb(nan) = maxInt(i32)
13pub fn ilogb(x: var) i32 {19pub fn ilogb(x: var) i32 {
14 const T = @typeOf(x);20 const T = @typeOf(x);
15 return switch (T) {21 return switch (T) {
std/math/inf.zig+1
...@@ -1,6 +1,7 @@...@@ -1,6 +1,7 @@
1const std = @import("../std.zig");1const std = @import("../std.zig");
2const math = std.math;2const math = std.math;
33
4/// Returns value inf for the type T.
4pub fn inf(comptime T: type) T {5pub fn inf(comptime T: type) T {
5 return switch (T) {6 return switch (T) {
6 f16 => math.inf_f16,7 f16 => math.inf_f16,
std/math/isfinite.zig+1
...@@ -3,6 +3,7 @@ const math = std.math;...@@ -3,6 +3,7 @@ const math = std.math;
3const expect = std.testing.expect;3const expect = std.testing.expect;
4const maxInt = std.math.maxInt;4const maxInt = std.math.maxInt;
55
6/// Returns whether x is a finite value.
6pub fn isFinite(x: var) bool {7pub fn isFinite(x: var) bool {
7 const T = @typeOf(x);8 const T = @typeOf(x);
8 switch (T) {9 switch (T) {
std/math/isinf.zig+3
...@@ -3,6 +3,7 @@ const math = std.math;...@@ -3,6 +3,7 @@ const math = std.math;
3const expect = std.testing.expect;3const expect = std.testing.expect;
4const maxInt = std.math.maxInt;4const maxInt = std.math.maxInt;
55
6/// Returns whether x is an infinity, ignoring sign.
6pub fn isInf(x: var) bool {7pub fn isInf(x: var) bool {
7 const T = @typeOf(x);8 const T = @typeOf(x);
8 switch (T) {9 switch (T) {
...@@ -28,6 +29,7 @@ pub fn isInf(x: var) bool {...@@ -28,6 +29,7 @@ pub fn isInf(x: var) bool {
28 }29 }
29}30}
3031
32/// Returns whether x is an infinity with a positive sign.
31pub fn isPositiveInf(x: var) bool {33pub fn isPositiveInf(x: var) bool {
32 const T = @typeOf(x);34 const T = @typeOf(x);
33 switch (T) {35 switch (T) {
...@@ -49,6 +51,7 @@ pub fn isPositiveInf(x: var) bool {...@@ -49,6 +51,7 @@ pub fn isPositiveInf(x: var) bool {
49 }51 }
50}52}
5153
54/// Returns whether x is an infinity with a negative sign.
52pub fn isNegativeInf(x: var) bool {55pub fn isNegativeInf(x: var) bool {
53 const T = @typeOf(x);56 const T = @typeOf(x);
54 switch (T) {57 switch (T) {
std/math/isnan.zig+4-2
...@@ -3,13 +3,15 @@ const math = std.math;...@@ -3,13 +3,15 @@ const math = std.math;
3const expect = std.testing.expect;3const expect = std.testing.expect;
4const maxInt = std.math.maxInt;4const maxInt = std.math.maxInt;
55
6/// Returns whether x is a nan.
6pub fn isNan(x: var) bool {7pub fn isNan(x: var) bool {
7 return x != x;8 return x != x;
8}9}
910
10/// Note: A signalling nan is identical to a standard nan right now but may have a different bit11/// Returns whether x is a signalling nan.
11/// representation in the future when required.
12pub fn isSignalNan(x: var) bool {12pub fn isSignalNan(x: var) bool {
13 // Note: A signalling nan is identical to a standard nan right now but may have a different bit
14 // representation in the future when required.
13 return isNan(x);15 return isNan(x);
14}16}
1517
std/math/isnormal.zig+1
...@@ -3,6 +3,7 @@ const math = std.math;...@@ -3,6 +3,7 @@ const math = std.math;
3const expect = std.testing.expect;3const expect = std.testing.expect;
4const maxInt = std.math.maxInt;4const maxInt = std.math.maxInt;
55
6// Returns whether x has a normalized representation (i.e. integer part of mantissa is 1).
6pub fn isNormal(x: var) bool {7pub fn isNormal(x: var) bool {
7 const T = @typeOf(x);8 const T = @typeOf(x);
8 switch (T) {9 switch (T) {
std/math/ln.zig+11-5
...@@ -1,9 +1,8 @@...@@ -1,9 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - ln(+inf) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/lnf.c
4// - ln(0) = -inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/ln.c
5// - ln(x) = nan if x < 0
6// - ln(nan) = nan
76
8const std = @import("../std.zig");7const std = @import("../std.zig");
9const math = std.math;8const math = std.math;
...@@ -11,6 +10,13 @@ const expect = std.testing.expect;...@@ -11,6 +10,13 @@ const expect = std.testing.expect;
11const builtin = @import("builtin");10const builtin = @import("builtin");
12const TypeId = builtin.TypeId;11const TypeId = builtin.TypeId;
1312
13/// Returns the natural logarithm of x.
14///
15/// Special Cases:
16/// - ln(+inf) = +inf
17/// - ln(0) = -inf
18/// - ln(x) = nan if x < 0
19/// - ln(nan) = nan
14pub fn ln(x: var) @typeOf(x) {20pub fn ln(x: var) @typeOf(x) {
15 const T = @typeOf(x);21 const T = @typeOf(x);
16 switch (@typeId(T)) {22 switch (@typeId(T)) {
std/math/log.zig+7
...@@ -1,9 +1,16 @@...@@ -1,9 +1,16 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/logf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/log.c
6
1const std = @import("../std.zig");7const std = @import("../std.zig");
2const math = std.math;8const math = std.math;
3const builtin = @import("builtin");9const builtin = @import("builtin");
4const TypeId = builtin.TypeId;10const TypeId = builtin.TypeId;
5const expect = std.testing.expect;11const expect = std.testing.expect;
612
13/// Returns the logarithm of x for the provided base.
7pub fn log(comptime T: type, base: T, x: T) T {14pub fn log(comptime T: type, base: T, x: T) T {
8 if (base == 2) {15 if (base == 2) {
9 return math.log2(x);16 return math.log2(x);
std/math/log10.zig+11-5
...@@ -1,9 +1,8 @@...@@ -1,9 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - log10(+inf) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/log10f.c
4// - log10(0) = -inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/log10.c
5// - log10(x) = nan if x < 0
6// - log10(nan) = nan
76
8const std = @import("../std.zig");7const std = @import("../std.zig");
9const math = std.math;8const math = std.math;
...@@ -12,6 +11,13 @@ const builtin = @import("builtin");...@@ -12,6 +11,13 @@ const builtin = @import("builtin");
12const TypeId = builtin.TypeId;11const TypeId = builtin.TypeId;
13const maxInt = std.math.maxInt;12const maxInt = std.math.maxInt;
1413
14/// Returns the base-10 logarithm of x.
15///
16/// Special Cases:
17/// - log10(+inf) = +inf
18/// - log10(0) = -inf
19/// - log10(x) = nan if x < 0
20/// - log10(nan) = nan
15pub fn log10(x: var) @typeOf(x) {21pub fn log10(x: var) @typeOf(x) {
16 const T = @typeOf(x);22 const T = @typeOf(x);
17 switch (@typeId(T)) {23 switch (@typeId(T)) {
std/math/log1p.zig+12-6
...@@ -1,16 +1,22 @@...@@ -1,16 +1,22 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - log1p(+inf) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/log1pf.c
4// - log1p(+-0) = +-05// https://git.musl-libc.org/cgit/musl/tree/src/math/log1p.c
5// - log1p(-1) = -inf
6// - log1p(x) = nan if x < -1
7// - log1p(nan) = nan
86
9const builtin = @import("builtin");7const builtin = @import("builtin");
10const std = @import("../std.zig");8const std = @import("../std.zig");
11const math = std.math;9const math = std.math;
12const expect = std.testing.expect;10const expect = std.testing.expect;
1311
12/// Returns the natural logarithm of 1 + x with greater accuracy when x is near zero.
13///
14/// Special Cases:
15/// - log1p(+inf) = +inf
16/// - log1p(+-0) = +-0
17/// - log1p(-1) = -inf
18/// - log1p(x) = nan if x < -1
19/// - log1p(nan) = nan
14pub fn log1p(x: var) @typeOf(x) {20pub fn log1p(x: var) @typeOf(x) {
15 const T = @typeOf(x);21 const T = @typeOf(x);
16 return switch (T) {22 return switch (T) {
std/math/log2.zig+11-5
...@@ -1,9 +1,8 @@...@@ -1,9 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - log2(+inf) = +inf4// https://git.musl-libc.org/cgit/musl/tree/src/math/log2f.c
4// - log2(0) = -inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/log2.c
5// - log2(x) = nan if x < 0
6// - log2(nan) = nan
76
8const std = @import("../std.zig");7const std = @import("../std.zig");
9const math = std.math;8const math = std.math;
...@@ -12,6 +11,13 @@ const builtin = @import("builtin");...@@ -12,6 +11,13 @@ const builtin = @import("builtin");
12const TypeId = builtin.TypeId;11const TypeId = builtin.TypeId;
13const maxInt = std.math.maxInt;12const maxInt = std.math.maxInt;
1413
14/// Returns the base-2 logarithm of x.
15///
16/// Special Cases:
17/// - log2(+inf) = +inf
18/// - log2(0) = -inf
19/// - log2(x) = nan if x < 0
20/// - log2(nan) = nan
15pub fn log2(x: var) @typeOf(x) {21pub fn log2(x: var) @typeOf(x) {
16 const T = @typeOf(x);22 const T = @typeOf(x);
17 switch (@typeId(T)) {23 switch (@typeId(T)) {
std/math/modf.zig+10-3
...@@ -1,7 +1,8 @@...@@ -1,7 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - modf(+-inf) = +-inf, nan4// https://git.musl-libc.org/cgit/musl/tree/src/math/modff.c
4// - modf(nan) = nan, nan5// https://git.musl-libc.org/cgit/musl/tree/src/math/modf.c
56
6const std = @import("../std.zig");7const std = @import("../std.zig");
7const math = std.math;8const math = std.math;
...@@ -17,6 +18,12 @@ fn modf_result(comptime T: type) type {...@@ -17,6 +18,12 @@ fn modf_result(comptime T: type) type {
17pub const modf32_result = modf_result(f32);18pub const modf32_result = modf_result(f32);
18pub const modf64_result = modf_result(f64);19pub const modf64_result = modf_result(f64);
1920
21/// Returns the integer and fractional floating-point numbers that sum to x. The sign of each
22/// result is the same as the sign of x.
23///
24/// Special Cases:
25/// - modf(+-inf) = +-inf, nan
26/// - modf(nan) = nan, nan
20pub fn modf(x: var) modf_result(@typeOf(x)) {27pub fn modf(x: var) modf_result(@typeOf(x)) {
21 const T = @typeOf(x);28 const T = @typeOf(x);
22 return switch (T) {29 return switch (T) {
std/math/nan.zig+4-2
...@@ -1,5 +1,6 @@...@@ -1,5 +1,6 @@
1const math = @import("../math.zig");1const math = @import("../math.zig");
22
3/// Returns the nan representation for type T.
3pub fn nan(comptime T: type) T {4pub fn nan(comptime T: type) T {
4 return switch (T) {5 return switch (T) {
5 f16 => math.nan_f16,6 f16 => math.nan_f16,
...@@ -10,9 +11,10 @@ pub fn nan(comptime T: type) T {...@@ -10,9 +11,10 @@ pub fn nan(comptime T: type) T {
10 };11 };
11}12}
1213
13// Note: A signalling nan is identical to a standard right now by may have a different bit14/// Returns the signalling nan representation for type T.
14// representation in the future when required.
15pub fn snan(comptime T: type) T {15pub fn snan(comptime T: type) T {
16 // Note: A signalling nan is identical to a standard right now by may have a different bit
17 // representation in the future when required.
16 return switch (T) {18 return switch (T) {
17 f16 => @bitCast(f16, math.nan_u16),19 f16 => @bitCast(f16, math.nan_u16),
18 f32 => @bitCast(f32, math.nan_u32),20 f32 => @bitCast(f32, math.nan_u32),
std/math/pow.zig+57-39
...@@ -1,32 +1,36 @@...@@ -1,32 +1,36 @@
1// Special Cases:1// Ported from go, which is licensed under a BSD-3 license.
2// https://golang.org/LICENSE
2//3//
3// pow(x, +-0) = 1 for any x4// https://golang.org/src/math/pow.go
4// pow(1, y) = 1 for any y
5// pow(x, 1) = x for any x
6// pow(nan, y) = nan
7// pow(x, nan) = nan
8// pow(+-0, y) = +-inf for y an odd integer < 0
9// pow(+-0, -inf) = +inf
10// pow(+-0, +inf) = +0
11// pow(+-0, y) = +inf for finite y < 0 and not an odd integer
12// pow(+-0, y) = +-0 for y an odd integer > 0
13// pow(+-0, y) = +0 for finite y > 0 and not an odd integer
14// pow(-1, +-inf) = 1
15// pow(x, +inf) = +inf for |x| > 1
16// pow(x, -inf) = +0 for |x| > 1
17// pow(x, +inf) = +0 for |x| < 1
18// pow(x, -inf) = +inf for |x| < 1
19// pow(+inf, y) = +inf for y > 0
20// pow(+inf, y) = +0 for y < 0
21// pow(-inf, y) = pow(-0, -y)
22// pow(x, y) = nan for finite x < 0 and finite non-integer y
235
24const builtin = @import("builtin");6const builtin = @import("builtin");
25const std = @import("../std.zig");7const std = @import("../std.zig");
26const math = std.math;8const math = std.math;
27const expect = std.testing.expect;9const expect = std.testing.expect;
2810
29// This implementation is taken from the go stlib, musl is a bit more complex.11/// Returns x raised to the power of y (x^y).
12///
13/// Special Cases:
14/// - pow(x, +-0) = 1 for any x
15/// - pow(1, y) = 1 for any y
16/// - pow(x, 1) = x for any x
17/// - pow(nan, y) = nan
18/// - pow(x, nan) = nan
19/// - pow(+-0, y) = +-inf for y an odd integer < 0
20/// - pow(+-0, -inf) = +inf
21/// - pow(+-0, +inf) = +0
22/// - pow(+-0, y) = +inf for finite y < 0 and not an odd integer
23/// - pow(+-0, y) = +-0 for y an odd integer > 0
24/// - pow(+-0, y) = +0 for finite y > 0 and not an odd integer
25/// - pow(-1, +-inf) = 1
26/// - pow(x, +inf) = +inf for |x| > 1
27/// - pow(x, -inf) = +0 for |x| > 1
28/// - pow(x, +inf) = +0 for |x| < 1
29/// - pow(x, -inf) = +inf for |x| < 1
30/// - pow(+inf, y) = +inf for y > 0
31/// - pow(+inf, y) = +0 for y < 0
32/// - pow(-inf, y) = pow(-0, -y)
33/// - pow(x, y) = nan for finite x < 0 and finite non-integer y
30pub fn pow(comptime T: type, x: T, y: T) T {34pub fn pow(comptime T: type, x: T, y: T) T {
31 if (@typeInfo(T) == builtin.TypeId.Int) {35 if (@typeInfo(T) == builtin.TypeId.Int) {
32 return math.powi(T, x, y) catch unreachable;36 return math.powi(T, x, y) catch unreachable;
...@@ -53,15 +57,6 @@ pub fn pow(comptime T: type, x: T, y: T) T {...@@ -53,15 +57,6 @@ pub fn pow(comptime T: type, x: T, y: T) T {
53 return x;57 return x;
54 }58 }
5559
56 // special case sqrt
57 if (y == 0.5) {
58 return math.sqrt(x);
59 }
60
61 if (y == -0.5) {
62 return 1 / math.sqrt(x);
63 }
64
65 if (x == 0) {60 if (x == 0) {
66 if (y < 0) {61 if (y < 0) {
67 // pow(+-0, y) = +- 0 for y an odd integer62 // pow(+-0, y) = +- 0 for y an odd integer
...@@ -112,14 +107,16 @@ pub fn pow(comptime T: type, x: T, y: T) T {...@@ -112,14 +107,16 @@ pub fn pow(comptime T: type, x: T, y: T) T {
112 }107 }
113 }108 }
114109
115 var ay = y;110 // special case sqrt
116 var flip = false;111 if (y == 0.5) {
117 if (ay < 0) {112 return math.sqrt(x);
118 ay = -ay;113 }
119 flip = true;114
115 if (y == -0.5) {
116 return 1 / math.sqrt(x);
120 }117 }
121118
122 const r1 = math.modf(ay);119 const r1 = math.modf(math.fabs(y));
123 var yi = r1.ipart;120 var yi = r1.ipart;
124 var yf = r1.fpart;121 var yf = r1.fpart;
125122
...@@ -148,8 +145,18 @@ pub fn pow(comptime T: type, x: T, y: T) T {...@@ -148,8 +145,18 @@ pub fn pow(comptime T: type, x: T, y: T) T {
148 var xe = r2.exponent;145 var xe = r2.exponent;
149 var x1 = r2.significand;146 var x1 = r2.significand;
150147
151 var i = @floatToInt(i32, yi);148 var i = @floatToInt(@IntType(true, T.bit_count), yi);
152 while (i != 0) : (i >>= 1) {149 while (i != 0) : (i >>= 1) {
150 const overflow_shift = math.floatExponentBits(T) + 1;
151 if (xe < -(1 << overflow_shift) or (1 << overflow_shift) < xe) {
152 // catch xe before it overflows the left shift below
153 // Since i != 0 it has at least one bit still set, so ae will accumulate xe
154 // on at least one more iteration, ae += xe is a lower bound on ae
155 // the lower bound on ae exceeds the size of a float exp
156 // so the final call to Ldexp will produce under/overflow (0/Inf)
157 ae += xe;
158 break;
159 }
153 if (i & 1 == 1) {160 if (i & 1 == 1) {
154 a1 *= x1;161 a1 *= x1;
155 ae += xe;162 ae += xe;
...@@ -163,7 +170,7 @@ pub fn pow(comptime T: type, x: T, y: T) T {...@@ -163,7 +170,7 @@ pub fn pow(comptime T: type, x: T, y: T) T {
163 }170 }
164171
165 // a *= a1 * 2^ae172 // a *= a1 * 2^ae
166 if (flip) {173 if (y < 0) {
167 a1 = 1 / a1;174 a1 = 1 / a1;
168 ae = -ae;175 ae = -ae;
169 }176 }
...@@ -202,6 +209,9 @@ test "math.pow.special" {...@@ -202,6 +209,9 @@ test "math.pow.special" {
202 expect(pow(f32, 45, 1.0) == 45);209 expect(pow(f32, 45, 1.0) == 45);
203 expect(pow(f32, -45, 1.0) == -45);210 expect(pow(f32, -45, 1.0) == -45);
204 expect(math.isNan(pow(f32, math.nan(f32), 5.0)));211 expect(math.isNan(pow(f32, math.nan(f32), 5.0)));
212 expect(math.isPositiveInf(pow(f32, -math.inf(f32), 0.5)));
213 expect(math.isPositiveInf(pow(f32, -0, -0.5)));
214 expect(pow(f32, -0, 0.5) == 0);
205 expect(math.isNan(pow(f32, 5.0, math.nan(f32))));215 expect(math.isNan(pow(f32, 5.0, math.nan(f32))));
206 expect(math.isPositiveInf(pow(f32, 0.0, -1.0)));216 expect(math.isPositiveInf(pow(f32, 0.0, -1.0)));
207 //expect(math.isNegativeInf(pow(f32, -0.0, -3.0))); TODO is this required?217 //expect(math.isNegativeInf(pow(f32, -0.0, -3.0))); TODO is this required?
...@@ -232,3 +242,11 @@ test "math.pow.special" {...@@ -232,3 +242,11 @@ test "math.pow.special" {
232 expect(math.isNan(pow(f32, -1.0, 1.2)));242 expect(math.isNan(pow(f32, -1.0, 1.2)));
233 expect(math.isNan(pow(f32, -12.4, 78.5)));243 expect(math.isNan(pow(f32, -12.4, 78.5)));
234}244}
245
246test "math.pow.overflow" {
247 expect(math.isPositiveInf(pow(f64, 2, 1 << 32)));
248 expect(pow(f64, 2, -(1 << 32)) == 0);
249 expect(math.isNegativeInf(pow(f64, -2, (1 << 32) + 1)));
250 expect(pow(f64, 0.5, 1 << 45) == 0);
251 expect(math.isPositiveInf(pow(f64, 0.5, -(1 << 45))));
252}
std/math/powi.zig+13-9
...@@ -1,12 +1,7 @@...@@ -1,12 +1,7 @@
1// Special Cases:1// Based on Rust, which is licensed under the MIT license.
2// https://github.com/rust-lang/rust/blob/360432f1e8794de58cd94f34c9c17ad65871e5b5/LICENSE-MIT
2//3//
3// powi(x, +-0) = 1 for any x4// https://github.com/rust-lang/rust/blob/360432f1e8794de58cd94f34c9c17ad65871e5b5/src/libcore/num/mod.rs#L3423
4// powi(0, y) = 0 for any y
5// powi(1, y) = 1 for any y
6// powi(-1, y) = -1 for for y an odd integer
7// powi(-1, y) = 1 for for y an even integer
8// powi(x, y) = Overflow for for y >= @sizeOf(x) - 1 y > 0
9// powi(x, y) = Underflow for for y > @sizeOf(x) - 1 y < 0
105
11const builtin = @import("builtin");6const builtin = @import("builtin");
12const std = @import("../std.zig");7const std = @import("../std.zig");
...@@ -14,7 +9,16 @@ const math = std.math;...@@ -14,7 +9,16 @@ const math = std.math;
14const assert = std.debug.assert;9const assert = std.debug.assert;
15const testing = std.testing;10const testing = std.testing;
1611
17// This implementation is based on that from the rust stlib12/// Returns the power of x raised by the integer y (x^y).
13///
14/// Special Cases:
15/// - powi(x, +-0) = 1 for any x
16/// - powi(0, y) = 0 for any y
17/// - powi(1, y) = 1 for any y
18/// - powi(-1, y) = -1 for y an odd integer
19/// - powi(-1, y) = 1 for y an even integer
20/// - powi(x, y) = Overflow for y >= @sizeOf(x) - 1 or y > 0
21/// - powi(x, y) = Underflow for y > @sizeOf(x) - 1 or y < 0
18pub fn powi(comptime T: type, x: T, y: T) (error{22pub fn powi(comptime T: type, x: T, y: T) (error{
19 Overflow,23 Overflow,
20 Underflow,24 Underflow,
std/math/round.zig+10-4
...@@ -1,14 +1,20 @@...@@ -1,14 +1,20 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - round(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/roundf.c
4// - round(+-inf) = +-inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/round.c
5// - round(nan) = nan
66
7const builtin = @import("builtin");7const builtin = @import("builtin");
8const expect = std.testing.expect;8const expect = std.testing.expect;
9const std = @import("../std.zig");9const std = @import("../std.zig");
10const math = std.math;10const math = std.math;
1111
12/// Returns x rounded to the nearest integer, rounding half away from zero.
13///
14/// Special Cases:
15/// - round(+-0) = +-0
16/// - round(+-inf) = +-inf
17/// - round(nan) = nan
12pub fn round(x: var) @typeOf(x) {18pub fn round(x: var) @typeOf(x) {
13 const T = @typeOf(x);19 const T = @typeOf(x);
14 return switch (T) {20 return switch (T) {
std/math/scalbn.zig+7
...@@ -1,7 +1,14 @@...@@ -1,7 +1,14 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbnf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbn.c
6
1const std = @import("../std.zig");7const std = @import("../std.zig");
2const math = std.math;8const math = std.math;
3const expect = std.testing.expect;9const expect = std.testing.expect;
410
11/// Returns x * 2^n.
5pub fn scalbn(x: var, n: i32) @typeOf(x) {12pub fn scalbn(x: var, n: i32) @typeOf(x) {
6 const T = @typeOf(x);13 const T = @typeOf(x);
7 return switch (T) {14 return switch (T) {
std/math/signbit.zig+1
...@@ -2,6 +2,7 @@ const std = @import("../std.zig");...@@ -2,6 +2,7 @@ const std = @import("../std.zig");
2const math = std.math;2const math = std.math;
3const expect = std.testing.expect;3const expect = std.testing.expect;
44
5/// Returns whether x is negative or negative 0.
5pub fn signbit(x: var) bool {6pub fn signbit(x: var) bool {
6 const T = @typeOf(x);7 const T = @typeOf(x);
7 return switch (T) {8 return switch (T) {
std/math/sin.zig+52-108
...@@ -1,19 +1,24 @@...@@ -1,19 +1,24 @@
1// Special Cases:1// Ported from go, which is licensed under a BSD-3 license.
2// https://golang.org/LICENSE
2//3//
3// - sin(+-0) = +-04// https://golang.org/src/math/sin.go
4// - sin(+-inf) = nan
5// - sin(nan) = nan
65
7const builtin = @import("builtin");6const builtin = @import("builtin");
8const std = @import("../std.zig");7const std = @import("../std.zig");
9const math = std.math;8const math = std.math;
10const expect = std.testing.expect;9const expect = std.testing.expect;
1110
11/// Returns the sine of the radian value x.
12///
13/// Special Cases:
14/// - sin(+-0) = +-0
15/// - sin(+-inf) = nan
16/// - sin(nan) = nan
12pub fn sin(x: var) @typeOf(x) {17pub fn sin(x: var) @typeOf(x) {
13 const T = @typeOf(x);18 const T = @typeOf(x);
14 return switch (T) {19 return switch (T) {
15 f32 => sin32(x),20 f32 => sin_(T, x),
16 f64 => sin64(x),21 f64 => sin_(T, x),
17 else => @compileError("sin not implemented for " ++ @typeName(T)),22 else => @compileError("sin not implemented for " ++ @typeName(T)),
18 };23 };
19}24}
...@@ -34,83 +39,27 @@ const C3 = 2.48015872888517045348E-5;...@@ -34,83 +39,27 @@ const C3 = 2.48015872888517045348E-5;
34const C4 = -1.38888888888730564116E-3;39const C4 = -1.38888888888730564116E-3;
35const C5 = 4.16666666666665929218E-2;40const C5 = 4.16666666666665929218E-2;
3641
37// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.42const pi4a = 7.85398125648498535156e-1;
38//43const pi4b = 3.77489470793079817668E-8;
39// This may have slight differences on some edge cases and may need to replaced if so.44const pi4c = 2.69515142907905952645E-15;
40fn sin32(x_: f32) f32 {45const m4pi = 1.273239544735162542821171882678754627704620361328125;
41 const pi4a = 7.85398125648498535156e-1;
42 const pi4b = 3.77489470793079817668E-8;
43 const pi4c = 2.69515142907905952645E-15;
44 const m4pi = 1.273239544735162542821171882678754627704620361328125;
45
46 var x = x_;
47 if (x == 0 or math.isNan(x)) {
48 return x;
49 }
50 if (math.isInf(x)) {
51 return math.nan(f32);
52 }
53
54 var sign = false;
55 if (x < 0) {
56 x = -x;
57 sign = true;
58 }
59
60 var y = math.floor(x * m4pi);
61 var j = @floatToInt(i64, y);
62
63 if (j & 1 == 1) {
64 j += 1;
65 y += 1;
66 }
6746
68 j &= 7;47fn sin_(comptime T: type, x_: T) T {
69 if (j > 3) {48 const I = @IntType(true, T.bit_count);
70 j -= 4;
71 sign = !sign;
72 }
73
74 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
75 const w = z * z;
76
77 const r = r: {
78 if (j == 1 or j == 2) {
79 break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
80 } else {
81 break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
82 }
83 };
84
85 if (sign) {
86 return -r;
87 } else {
88 return r;
89 }
90}
91
92fn sin64(x_: f64) f64 {
93 const pi4a = 7.85398125648498535156e-1;
94 const pi4b = 3.77489470793079817668E-8;
95 const pi4c = 2.69515142907905952645E-15;
96 const m4pi = 1.273239544735162542821171882678754627704620361328125;
9749
98 var x = x_;50 var x = x_;
99 if (x == 0 or math.isNan(x)) {51 if (x == 0 or math.isNan(x)) {
100 return x;52 return x;
101 }53 }
102 if (math.isInf(x)) {54 if (math.isInf(x)) {
103 return math.nan(f64);55 return math.nan(T);
104 }56 }
10557
106 var sign = false;58 var sign = x < 0;
107 if (x < 0) {59 x = math.fabs(x);
108 x = -x;
109 sign = true;
110 }
11160
112 var y = math.floor(x * m4pi);61 var y = math.floor(x * m4pi);
113 var j = @floatToInt(i64, y);62 var j = @floatToInt(I, y);
11463
115 if (j & 1 == 1) {64 if (j & 1 == 1) {
116 j += 1;65 j += 1;
...@@ -126,61 +75,56 @@ fn sin64(x_: f64) f64 {...@@ -126,61 +75,56 @@ fn sin64(x_: f64) f64 {
126 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;75 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
127 const w = z * z;76 const w = z * z;
12877
129 const r = r: {78 const r = if (j == 1 or j == 2)
130 if (j == 1 or j == 2) {79 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))))
131 break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));80 else
132 } else {81 z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
133 break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
134 }
135 };
13682
137 if (sign) {83 return if (sign) -r else r;
138 return -r;
139 } else {
140 return r;
141 }
142}84}
14385
144test "math.sin" {86test "math.sin" {
145 expect(sin(f32(0.0)) == sin32(0.0));87 expect(sin(f32(0.0)) == sin_(f32, 0.0));
146 expect(sin(f64(0.0)) == sin64(0.0));88 expect(sin(f64(0.0)) == sin_(f64, 0.0));
147 expect(comptime (math.sin(f64(2))) == math.sin(f64(2)));89 expect(comptime (math.sin(f64(2))) == math.sin(f64(2)));
148}90}
14991
150test "math.sin32" {92test "math.sin32" {
151 const epsilon = 0.000001;93 const epsilon = 0.000001;
15294
153 expect(math.approxEq(f32, sin32(0.0), 0.0, epsilon));95 expect(math.approxEq(f32, sin_(f32, 0.0), 0.0, epsilon));
154 expect(math.approxEq(f32, sin32(0.2), 0.198669, epsilon));96 expect(math.approxEq(f32, sin_(f32, 0.2), 0.198669, epsilon));
155 expect(math.approxEq(f32, sin32(0.8923), 0.778517, epsilon));97 expect(math.approxEq(f32, sin_(f32, 0.8923), 0.778517, epsilon));
156 expect(math.approxEq(f32, sin32(1.5), 0.997495, epsilon));98 expect(math.approxEq(f32, sin_(f32, 1.5), 0.997495, epsilon));
157 expect(math.approxEq(f32, sin32(37.45), -0.246544, epsilon));99 expect(math.approxEq(f32, sin_(f32, -1.5), -0.997495, epsilon));
158 expect(math.approxEq(f32, sin32(89.123), 0.916166, epsilon));100 expect(math.approxEq(f32, sin_(f32, 37.45), -0.246544, epsilon));
101 expect(math.approxEq(f32, sin_(f32, 89.123), 0.916166, epsilon));
159}102}
160103
161test "math.sin64" {104test "math.sin64" {
162 const epsilon = 0.000001;105 const epsilon = 0.000001;
163106
164 expect(math.approxEq(f64, sin64(0.0), 0.0, epsilon));107 expect(math.approxEq(f64, sin_(f64, 0.0), 0.0, epsilon));
165 expect(math.approxEq(f64, sin64(0.2), 0.198669, epsilon));108 expect(math.approxEq(f64, sin_(f64, 0.2), 0.198669, epsilon));
166 expect(math.approxEq(f64, sin64(0.8923), 0.778517, epsilon));109 expect(math.approxEq(f64, sin_(f64, 0.8923), 0.778517, epsilon));
167 expect(math.approxEq(f64, sin64(1.5), 0.997495, epsilon));110 expect(math.approxEq(f64, sin_(f64, 1.5), 0.997495, epsilon));
168 expect(math.approxEq(f64, sin64(37.45), -0.246543, epsilon));111 expect(math.approxEq(f64, sin_(f64, -1.5), -0.997495, epsilon));
169 expect(math.approxEq(f64, sin64(89.123), 0.916166, epsilon));112 expect(math.approxEq(f64, sin_(f64, 37.45), -0.246543, epsilon));
113 expect(math.approxEq(f64, sin_(f64, 89.123), 0.916166, epsilon));
170}114}
171115
172test "math.sin32.special" {116test "math.sin32.special" {
173 expect(sin32(0.0) == 0.0);117 expect(sin_(f32, 0.0) == 0.0);
174 expect(sin32(-0.0) == -0.0);118 expect(sin_(f32, -0.0) == -0.0);
175 expect(math.isNan(sin32(math.inf(f32))));119 expect(math.isNan(sin_(f32, math.inf(f32))));
176 expect(math.isNan(sin32(-math.inf(f32))));120 expect(math.isNan(sin_(f32, -math.inf(f32))));
177 expect(math.isNan(sin32(math.nan(f32))));121 expect(math.isNan(sin_(f32, math.nan(f32))));
178}122}
179123
180test "math.sin64.special" {124test "math.sin64.special" {
181 expect(sin64(0.0) == 0.0);125 expect(sin_(f64, 0.0) == 0.0);
182 expect(sin64(-0.0) == -0.0);126 expect(sin_(f64, -0.0) == -0.0);
183 expect(math.isNan(sin64(math.inf(f64))));127 expect(math.isNan(sin_(f64, math.inf(f64))));
184 expect(math.isNan(sin64(-math.inf(f64))));128 expect(math.isNan(sin_(f64, -math.inf(f64))));
185 expect(math.isNan(sin64(math.nan(f64))));129 expect(math.isNan(sin_(f64, math.nan(f64))));
186}130}
std/math/sinh.zig+10-4
...@@ -1,8 +1,8 @@...@@ -1,8 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - sinh(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/sinhf.c
4// - sinh(+-inf) = +-inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/sinh.c
5// - sinh(nan) = nan
66
7const builtin = @import("builtin");7const builtin = @import("builtin");
8const std = @import("../std.zig");8const std = @import("../std.zig");
...@@ -11,6 +11,12 @@ const expect = std.testing.expect;...@@ -11,6 +11,12 @@ const expect = std.testing.expect;
11const expo2 = @import("expo2.zig").expo2;11const expo2 = @import("expo2.zig").expo2;
12const maxInt = std.math.maxInt;12const maxInt = std.math.maxInt;
1313
14/// Returns the hyperbolic sine of x.
15///
16/// Special Cases:
17/// - sinh(+-0) = +-0
18/// - sinh(+-inf) = +-inf
19/// - sinh(nan) = nan
14pub fn sinh(x: var) @typeOf(x) {20pub fn sinh(x: var) @typeOf(x) {
15 const T = @typeOf(x);21 const T = @typeOf(x);
16 return switch (T) {22 return switch (T) {
std/math/sqrt.zig+7-7
...@@ -1,10 +1,3 @@...@@ -1,10 +1,3 @@
1// Special Cases:
2//
3// - sqrt(+inf) = +inf
4// - sqrt(+-0) = +-0
5// - sqrt(x) = nan if x < 0
6// - sqrt(nan) = nan
7
8const std = @import("../std.zig");1const std = @import("../std.zig");
9const math = std.math;2const math = std.math;
10const expect = std.testing.expect;3const expect = std.testing.expect;
...@@ -12,6 +5,13 @@ const builtin = @import("builtin");...@@ -12,6 +5,13 @@ const builtin = @import("builtin");
12const TypeId = builtin.TypeId;5const TypeId = builtin.TypeId;
13const maxInt = std.math.maxInt;6const maxInt = std.math.maxInt;
147
8/// Returns the square root of x.
9///
10/// Special Cases:
11/// - sqrt(+inf) = +inf
12/// - sqrt(+-0) = +-0
13/// - sqrt(x) = nan if x < 0
14/// - sqrt(nan) = nan
15pub fn sqrt(x: var) (if (@typeId(@typeOf(x)) == TypeId.Int) @IntType(false, @typeOf(x).bit_count / 2) else @typeOf(x)) {15pub fn sqrt(x: var) (if (@typeId(@typeOf(x)) == TypeId.Int) @IntType(false, @typeOf(x).bit_count / 2) else @typeOf(x)) {
16 const T = @typeOf(x);16 const T = @typeOf(x);
17 switch (@typeId(T)) {17 switch (@typeId(T)) {
std/math/tan.zig+50-104
...@@ -1,19 +1,24 @@...@@ -1,19 +1,24 @@
1// Special Cases:1// Ported from go, which is licensed under a BSD-3 license.
2// https://golang.org/LICENSE
2//3//
3// - tan(+-0) = +-04// https://golang.org/src/math/tan.go
4// - tan(+-inf) = nan
5// - tan(nan) = nan
65
7const builtin = @import("builtin");6const builtin = @import("builtin");
8const std = @import("../std.zig");7const std = @import("../std.zig");
9const math = std.math;8const math = std.math;
10const expect = std.testing.expect;9const expect = std.testing.expect;
1110
11/// Returns the tangent of the radian value x.
12///
13/// Special Cases:
14/// - tan(+-0) = +-0
15/// - tan(+-inf) = nan
16/// - tan(nan) = nan
12pub fn tan(x: var) @typeOf(x) {17pub fn tan(x: var) @typeOf(x) {
13 const T = @typeOf(x);18 const T = @typeOf(x);
14 return switch (T) {19 return switch (T) {
15 f32 => tan32(x),20 f32 => tan_(f32, x),
16 f64 => tan64(x),21 f64 => tan_(f64, x),
17 else => @compileError("tan not implemented for " ++ @typeName(T)),22 else => @compileError("tan not implemented for " ++ @typeName(T)),
18 };23 };
19}24}
...@@ -27,80 +32,27 @@ const Tq2 = -1.32089234440210967447E6;...@@ -27,80 +32,27 @@ const Tq2 = -1.32089234440210967447E6;
27const Tq3 = 2.50083801823357915839E7;32const Tq3 = 2.50083801823357915839E7;
28const Tq4 = -5.38695755929454629881E7;33const Tq4 = -5.38695755929454629881E7;
2934
30// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.35const pi4a = 7.85398125648498535156e-1;
31//36const pi4b = 3.77489470793079817668E-8;
32// This may have slight differences on some edge cases and may need to replaced if so.37const pi4c = 2.69515142907905952645E-15;
33fn tan32(x_: f32) f32 {38const m4pi = 1.273239544735162542821171882678754627704620361328125;
34 const pi4a = 7.85398125648498535156e-1;
35 const pi4b = 3.77489470793079817668E-8;
36 const pi4c = 2.69515142907905952645E-15;
37 const m4pi = 1.273239544735162542821171882678754627704620361328125;
38
39 var x = x_;
40 if (x == 0 or math.isNan(x)) {
41 return x;
42 }
43 if (math.isInf(x)) {
44 return math.nan(f32);
45 }
46
47 var sign = false;
48 if (x < 0) {
49 x = -x;
50 sign = true;
51 }
52
53 var y = math.floor(x * m4pi);
54 var j = @floatToInt(i64, y);
55
56 if (j & 1 == 1) {
57 j += 1;
58 y += 1;
59 }
60
61 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
62 const w = z * z;
63
64 var r = r: {
65 if (w > 1e-14) {
66 break :r z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4));
67 } else {
68 break :r z;
69 }
70 };
71
72 if (j & 2 == 2) {
73 r = -1 / r;
74 }
75 if (sign) {
76 r = -r;
77 }
78
79 return r;
80}
8139
82fn tan64(x_: f64) f64 {40fn tan_(comptime T: type, x_: T) T {
83 const pi4a = 7.85398125648498535156e-1;41 const I = @IntType(true, T.bit_count);
84 const pi4b = 3.77489470793079817668E-8;
85 const pi4c = 2.69515142907905952645E-15;
86 const m4pi = 1.273239544735162542821171882678754627704620361328125;
8742
88 var x = x_;43 var x = x_;
89 if (x == 0 or math.isNan(x)) {44 if (x == 0 or math.isNan(x)) {
90 return x;45 return x;
91 }46 }
92 if (math.isInf(x)) {47 if (math.isInf(x)) {
93 return math.nan(f64);48 return math.nan(T);
94 }49 }
9550
96 var sign = false;51 var sign = x < 0;
97 if (x < 0) {52 x = math.fabs(x);
98 x = -x;
99 sign = true;
100 }
10153
102 var y = math.floor(x * m4pi);54 var y = math.floor(x * m4pi);
103 var j = @floatToInt(i64, y);55 var j = @floatToInt(I, y);
10456
105 if (j & 1 == 1) {57 if (j & 1 == 1) {
106 j += 1;58 j += 1;
...@@ -110,63 +62,57 @@ fn tan64(x_: f64) f64 {...@@ -110,63 +62,57 @@ fn tan64(x_: f64) f64 {
110 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;62 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
111 const w = z * z;63 const w = z * z;
11264
113 var r = r: {65 var r = if (w > 1e-14)
114 if (w > 1e-14) {66 z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4))
115 break :r z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4));67 else
116 } else {68 z;
117 break :r z;
118 }
119 };
12069
121 if (j & 2 == 2) {70 if (j & 2 == 2) {
122 r = -1 / r;71 r = -1 / r;
123 }72 }
124 if (sign) {
125 r = -r;
126 }
12773
128 return r;74 return if (sign) -r else r;
129}75}
13076
131test "math.tan" {77test "math.tan" {
132 expect(tan(f32(0.0)) == tan32(0.0));78 expect(tan(f32(0.0)) == tan_(f32, 0.0));
133 expect(tan(f64(0.0)) == tan64(0.0));79 expect(tan(f64(0.0)) == tan_(f64, 0.0));
134}80}
13581
136test "math.tan32" {82test "math.tan32" {
137 const epsilon = 0.000001;83 const epsilon = 0.000001;
13884
139 expect(math.approxEq(f32, tan32(0.0), 0.0, epsilon));85 expect(math.approxEq(f32, tan_(f32, 0.0), 0.0, epsilon));
140 expect(math.approxEq(f32, tan32(0.2), 0.202710, epsilon));86 expect(math.approxEq(f32, tan_(f32, 0.2), 0.202710, epsilon));
141 expect(math.approxEq(f32, tan32(0.8923), 1.240422, epsilon));87 expect(math.approxEq(f32, tan_(f32, 0.8923), 1.240422, epsilon));
142 expect(math.approxEq(f32, tan32(1.5), 14.101420, epsilon));88 expect(math.approxEq(f32, tan_(f32, 1.5), 14.101420, epsilon));
143 expect(math.approxEq(f32, tan32(37.45), -0.254397, epsilon));89 expect(math.approxEq(f32, tan_(f32, 37.45), -0.254397, epsilon));
144 expect(math.approxEq(f32, tan32(89.123), 2.285852, epsilon));90 expect(math.approxEq(f32, tan_(f32, 89.123), 2.285852, epsilon));
145}91}
14692
147test "math.tan64" {93test "math.tan64" {
148 const epsilon = 0.000001;94 const epsilon = 0.000001;
14995
150 expect(math.approxEq(f64, tan64(0.0), 0.0, epsilon));96 expect(math.approxEq(f64, tan_(f64, 0.0), 0.0, epsilon));
151 expect(math.approxEq(f64, tan64(0.2), 0.202710, epsilon));97 expect(math.approxEq(f64, tan_(f64, 0.2), 0.202710, epsilon));
152 expect(math.approxEq(f64, tan64(0.8923), 1.240422, epsilon));98 expect(math.approxEq(f64, tan_(f64, 0.8923), 1.240422, epsilon));
153 expect(math.approxEq(f64, tan64(1.5), 14.101420, epsilon));99 expect(math.approxEq(f64, tan_(f64, 1.5), 14.101420, epsilon));
154 expect(math.approxEq(f64, tan64(37.45), -0.254397, epsilon));100 expect(math.approxEq(f64, tan_(f64, 37.45), -0.254397, epsilon));
155 expect(math.approxEq(f64, tan64(89.123), 2.2858376, epsilon));101 expect(math.approxEq(f64, tan_(f64, 89.123), 2.2858376, epsilon));
156}102}
157103
158test "math.tan32.special" {104test "math.tan32.special" {
159 expect(tan32(0.0) == 0.0);105 expect(tan_(f32, 0.0) == 0.0);
160 expect(tan32(-0.0) == -0.0);106 expect(tan_(f32, -0.0) == -0.0);
161 expect(math.isNan(tan32(math.inf(f32))));107 expect(math.isNan(tan_(f32, math.inf(f32))));
162 expect(math.isNan(tan32(-math.inf(f32))));108 expect(math.isNan(tan_(f32, -math.inf(f32))));
163 expect(math.isNan(tan32(math.nan(f32))));109 expect(math.isNan(tan_(f32, math.nan(f32))));
164}110}
165111
166test "math.tan64.special" {112test "math.tan64.special" {
167 expect(tan64(0.0) == 0.0);113 expect(tan_(f64, 0.0) == 0.0);
168 expect(tan64(-0.0) == -0.0);114 expect(tan_(f64, -0.0) == -0.0);
169 expect(math.isNan(tan64(math.inf(f64))));115 expect(math.isNan(tan_(f64, math.inf(f64))));
170 expect(math.isNan(tan64(-math.inf(f64))));116 expect(math.isNan(tan_(f64, -math.inf(f64))));
171 expect(math.isNan(tan64(math.nan(f64))));117 expect(math.isNan(tan_(f64, math.nan(f64))));
172}118}
std/math/tanh.zig+10-4
...@@ -1,8 +1,8 @@...@@ -1,8 +1,8 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - sinh(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/tanhf.c
4// - sinh(+-inf) = +-15// https://git.musl-libc.org/cgit/musl/tree/src/math/tanh.c
5// - sinh(nan) = nan
66
7const builtin = @import("builtin");7const builtin = @import("builtin");
8const std = @import("../std.zig");8const std = @import("../std.zig");
...@@ -11,6 +11,12 @@ const expect = std.testing.expect;...@@ -11,6 +11,12 @@ const expect = std.testing.expect;
11const expo2 = @import("expo2.zig").expo2;11const expo2 = @import("expo2.zig").expo2;
12const maxInt = std.math.maxInt;12const maxInt = std.math.maxInt;
1313
14/// Returns the hyperbolic tangent of x.
15///
16/// Special Cases:
17/// - sinh(+-0) = +-0
18/// - sinh(+-inf) = +-1
19/// - sinh(nan) = nan
14pub fn tanh(x: var) @typeOf(x) {20pub fn tanh(x: var) @typeOf(x) {
15 const T = @typeOf(x);21 const T = @typeOf(x);
16 return switch (T) {22 return switch (T) {
std/math/trunc.zig+10-4
...@@ -1,14 +1,20 @@...@@ -1,14 +1,20 @@
1// Special Cases:1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
2//3//
3// - trunc(+-0) = +-04// https://git.musl-libc.org/cgit/musl/tree/src/math/truncf.c
4// - trunc(+-inf) = +-inf5// https://git.musl-libc.org/cgit/musl/tree/src/math/trunc.c
5// - trunc(nan) = nan
66
7const std = @import("../std.zig");7const std = @import("../std.zig");
8const math = std.math;8const math = std.math;
9const expect = std.testing.expect;9const expect = std.testing.expect;
10const maxInt = std.math.maxInt;10const maxInt = std.math.maxInt;
1111
12/// Returns the integer value of x.
13///
14/// Special Cases:
15/// - trunc(+-0) = +-0
16/// - trunc(+-inf) = +-inf
17/// - trunc(nan) = nan
12pub fn trunc(x: var) @typeOf(x) {18pub fn trunc(x: var) @typeOf(x) {
13 const T = @typeOf(x);19 const T = @typeOf(x);
14 return switch (T) {20 return switch (T) {