Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

rewrite math.complex #20346

Closed
wants to merge 4 commits into from
Closed

rewrite math.complex #20346

wants to merge 4 commits into from

Conversation

expikr
Copy link
Contributor

@expikr expikr commented Jun 19, 2024

Closes #19207

This decouples the math.complex namespace away from a centralized type definition, instead providing a bring-your-own-type collection of utility functions with no dependency on the specific types' declarations.

The math.Complex type constructor is removed, instead encouraging programmers to pick-and-choose the subset of provided functions to incorporate into their type of choice.

This remains the case even once #16278 is implemented, as the collection already makes no assumptions about the concrete type beyond the presence of a re and im field.

The source files are reorganized such that public-facing dispatch interfaces are all inside complex.zig, the subdirectory instead only containing datatype-specific internal implementations that is only @imported inside the dispatch function body.

The datatype-specific internal implementations are entirely oblivious to any sort of type constructs, receiving inputs as multiple scalar arguments and returning tuples of values.

Future Incremental Improvements:

  • audit the original musl implementation, reimplement from first-principles where appropriate
    • ldexp_cexp (no error)
    • exp (no error)
    • cosh (3 errors)
    • sinh (no error)
    • tanh (no error)
    • sqrt (2 errors)
    • atan (2 errors)
  • fill out fallback branches with naive algebraic expressions
    • atan
    • tanh
    • sqrt
  • create exhaustive numerical precision tests for all datatype-specific impls

@expikr
Copy link
Contributor Author

expikr commented Jul 20, 2024

c.c. @tiehuis since this largely overhauls #949 and also partly progresses #19476

@tiehuis
Copy link
Member

tiehuis commented Jul 22, 2024

Happy to merge the fixes separately as these look good, but will need a bit more time to check over the api changes.

I generally think they look promising, but given the proposal isn't accepted and it changes quite a bit want a closer look.

Comment on lines -35 to +24
const u = @as(f32, @floatFromInt(@as(i32, @intFromFloat(t))));
return ((x - u * DP1) - u * DP2) - t * DP3;
const u = @trunc(t);
return ((x - u * DP1) - u * DP2) - u * DP3;
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Replaced unnecessary foat->int->float cast with a single @trunc.
Correction: fixed mistakenly using t when it should be u.

Comment on lines -52 to -55
if (a == 0.0) {
// overflow
return Complex(f32).init(maxnum, maxnum);
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Overflow branch no longer present in upstream.

Comment on lines -62 to -65
if (a == 0.0) {
// overflow
return Complex(f32).init(maxnum, maxnum);
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Overflow branch no longer present in upstream.

Comment on lines -84 to +55
const u = @as(f64, @floatFromInt(@as(i64, @intFromFloat(t))));
return ((x - u * DP1) - u * DP2) - t * DP3;
const u: f64 = @trunc(t);
return ((x - u * DP1) - u * DP2) - u * DP3;
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Replaced unnecessary foat->int->float cast with a single @trunc.
Correction: fixed mistakenly using t when it should be u.

Copy link
Contributor Author

@expikr expikr Jul 25, 2024

Choose a reason for hiding this comment

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

This file is replaced in its entirety with nothing carried over.

API-wise, the following are changed:

  • various helper additions, and tests are all rewritten.
  • removed: Complex type constructor.
  • changed: unary complex functions that output complex types now inherits the the argument type.
  • changed: binary complex functions that output complex types now receives an explicit type argument for its output type.

Impl-wise, the following are changed:

  • Some functions which used to only support c32 and c64 now have a generic fallback implementations.
  • Inverse complex-trig functions that depended on each other as trivial permutations are now instead rewritten from first principles using the primitive operations (log, etc). This also eliminates some problematic choices in upstream where, instead of permuting argument or output components, they add a pi/2 phase shift to values which is imprecise.

Comment on lines -45 to +35
if (math.signbit(x)) {
return Complex(f32).init(@abs(x - y), math.copysign(x, y));
} else {
return Complex(f32).init(x, math.copysign(y - y, y));
}
return if (math.signbit(x)) .{
@abs(y - y),
math.copysign(x, y),
} else .{
x,
math.copysign(y - y, y),
};
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Correction: the real component of the have-sign-bit case should be @abs(y - y) instead of @abs(x - y).

Comment on lines -96 to +86
if (math.signbit(x)) {
return Complex(f64).init(@abs(x - y), math.copysign(x, y));
} else {
return Complex(f64).init(x, math.copysign(y - y, y));
}
return if (math.signbit(x)) .{
@abs(y - y),
math.copysign(x, y),
} else .{
x,
math.copysign(y - y, y),
};
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Correction: the real component of the have-sign-bit case should be @abs(y - y) instead of @abs(x - y).

Comment on lines -67 to +75
if (iy == 0 and ix >= 0x7f800000) {
if (hx & 0x7fffff == 0) {
return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), x) * y);
}
return Complex(f32).init(x, math.copysign(@as(f32, 0.0), (x + x) * y));
}
if (iy == 0 and ix >= 0x7f800000) return .{
x * x,
if (hx & 0x7fffff == 0)
math.copysign(zero, x) * y
else
math.copysign(zero, (x + x) * y),
};
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Correction: the real component of both cases should be x * x instead of the nonzero-x case being different as x.

Comment on lines -48 to +44
return Complex(f32).init(math.copysign(h, x) * @cos(y), h * @sin(y));
break :ret .{
@cos(y) * h,
@sin(y) * math.copysign(h, x),
};
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Correction: the multipliers on cos/sin of the real/imag parts are swapped.

Comment on lines -35 to +30
if (ix < 0x7f800000 and iy < 0x7f800000) {
if (iy == 0) {
return Complex(f32).init(math.cosh(x), y);
}
if (ix < 0x7f800000 and iy < 0x7f800000) return ret: {
if (iy == 0) break :ret .{
math.cosh(x),
x * y,
};

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Correction: the imag component should be x * y instead of just y.

@tiehuis
Copy link
Member

tiehuis commented Jul 27, 2024

My primary question is what does removing the Complex type actually give us?

A user having to create a structurally matching type seems odd. I can't think of any great use-cases that this approach gives over the dedicated type in 99% of cases. We also lose documentation on a base type we could refer to.

@andrewrk
Copy link
Member

Author has not complied with my request to limit to 1 open PR at a time, is now in timeout for 1 year.

@andrewrk andrewrk closed this Jul 29, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Proposal: nuke the std.math.complex namespace and make math functions generic over Complex
3 participants