apple / swift-numerics

Advanced mathematical types and functions for Swift
Apache License 2.0
1.67k stars 142 forks source link

Possible improvement in 'pow' Elementary Function when operands are zero #234

Closed albixp closed 1 year ago

albixp commented 2 years ago

Some elementary functions might require to catch when one or both operands are zero to report the correct result. One of those is pow and of course this is applicable for both implementation in Complex and Real Module:

Issue Result of Complex.pow(0, 0) is .infinite instead 1.
Package swift-numerics 1.0.2 Modules Complex Module / Elementary Functions Real Module / Elementary Functions

Example of Complex Module 'pow' possible fix

extension Complex: ElementaryFunctions {
   ...
    // MARK_ - pow-like function catching both inputs equal to zero
    @inlinable
    public static func pow(_ z: Complex, _ w: Complex) -> Complex {
        if z.isZero && w.isZero { return .one} // Proposed improvement
        if z.isZero { return .zero } // Proposed improvement
        return exp(w * log(z))
    }

  @inlinable
  public static func pow(_ z: Complex, _ n: Int) -> Complex {
    if z.isZero && n == 0 { return .one} // Proposed improvement
    if z.isZero { return .zero }
    // TODO: this implementation is not quite correct, because n may be
    // rounded in conversion to RealType. This only effects very extreme
    // cases, so we'll leave it alone for now.
    //
    // Note that this does not have the same problems that a similar
    // implementation for a real type would have, because there's no
    // parity/sign interaction in the complex plane.
    return exp(log(z).multiplied(by: RealType(n)))
  }

}
NevinBR commented 2 years ago

You’ve identified a real issue, and I think we ought to check the sign of the real part of the exponent as well:

if z.isZero {
  if w.isZero { return .one }
  if w.real > 0 { return .zero }
  return .infinity
}
stephentyrone commented 1 year ago

pow(x, y) is defined by exp(y log(x)), which has an essential singularity at (0, 0):

  /// exp(y * log(x)) computed with additional internal precision.
  ///
  /// See also `sqrt()` and `root()`.
  ///
  static func pow(_ x: Self, _ y: Self) -> Self

The suggested change would be a bug, except for the Complex.pow(: Self, : Int) case.

NevinBR commented 1 year ago

pow(x, y) is defined by exp(y log(x))

Does that mean pow(.zero, w) is always infinite for any w: Complex?

And if so, is that actually what we want?

stephentyrone commented 1 year ago

Yes, that's correct. pow(_:Self,_:Self) binds the IEEE 754 powr function for real types, and extends its definition for complex types in the obvious fashion (which means that we should fix this case for real types as well).

pow(_:Self,_:Int) binds the IEEE 754 pown function, with 0, 1, +infinity results.

stephentyrone commented 1 year ago

https://github.com/apple/swift-numerics/pull/235