delphidabbler / code-snippets

Collections of reusable code snippets, mainly in Pascal.
MIT License
27 stars 1 forks source link

Add implementations of IEEE `pown` function. #20

Open delphidabbler opened 1 year ago

delphidabbler commented 1 year ago

This function raises real numbers to integral powers.

According to IEEE Std 754-2008 for Floating-Point Arithmetic, page 44:

  • pown(x, 0) is 1 for any x (even a zero, quiet NaN, or infinity)
  • pown(±0, n) is ±∞ and signals the divideByZero exception for odd integral n<0
  • pown(±0, n) is +∞ and signals the divideByZero exception for even integral n<0
  • pown(±0, n) is +0 for even integral n>0
  • pown(±0, n) is ±0 for odd integral n>0.

In addition to actual pown function more restricted versions would be useful:

delphidabbler commented 1 year ago

Here's a version that raises integer X to non-negative integral number power N.

In this function we have n≥0 so the IEEE rules above reduce to:

Since n≥0 and both x and n are integral, pown(x, n) is always integral. So the function can return an integral type.

function PowNZN(X: Integer; N: Cardinal): UInt64;
begin
  if (X = 0) and (N > 0) then
    Exit(0);
  Result := 1;
  for var I := 1 to N do
    Result := Result * X;
end;

Note that PowNZN was chosen as a name rather than PowN since the IEEE pown function has parameters x ∊ [−∞, +∞] and n ∊ ℤ, while PowNZN has parameters x ∊ ℤ and n ∊ ℕ.

delphidabbler commented 1 year ago

Here's a version that raises integer X to integral power N.

The IEEE rules can be simplified as:

Since we can have n<0 and where pown(x, n) = 1 / (pown(x, |n|), the return type must be real.

⚠️ Despite returning an real type, the function only supports results in the range -High(UInt64)..High(UInt64).

function PowNZZ(X: Integer; N: Integer): Extended;
begin
  if N = 0 then
    Exit(1.0);
  if X = 0 then
  begin
    if N < 0 then
      raise EDivideByZero.Create('Negative exponent for X = 0');
    Exit(0.0);
  end;
  ResultInt: UInt64 := 1;
  for var I := 1 to Abs(N) do
    ResultInt := ResultInt * X;
  if N > 0 then
    Result := ResultInt
  else // N < 0
    Result := 1 / ResultInt;
end;

Note that PowNZZ was chosen as a name rather than PowN since the IEEE pown function has parameters x ∊ [−∞, +∞] and n ∊ ℤ, while PowNZZ has parameters x ∊ ℤ and n ∊ ℤ.