pow

17 Jan 2002

      (I'm surprised a bit that nobody noticed this before.)
The complex "pow" operator gives completely bogus results with negative
exponents. Example:

[artax:pascal]% cat tst.pas
program Foo (Output);

var z: Complex;

begin
  z := Cmplx (5, 0) pow (-1);
  WriteLn (Re (z) : 5 : 1, Im (z) : 5 : 1)
end.
[artax:pascal]% gpc tst.pas
[artax:pascal]% ./a.out
  1.0  0.0

The "Complex_Pow" function in rts/math.pas computes reciprocals as
1 / (x+iy) := (x-iy) / |x+iy|, which is a nonsense.
Here is a patch (recycling a Maurice's idea to make it work for too big
or too small arguments).

----------------------------------------
function Complex_Pow (z: Complex; y: Integer) = r: Complex;
var a, b: Real;
begin
  if y < 0 then
    begin
      if Abs (Re (z)) > Abs (Im (z)) then begin
	a := Im (z) / Re (z);
	b := Re (z) + Im (z) * a;
	z := Cmplx (1 / b, -a / b)
      end else begin
	if Im (z) = 0 then begin
          SetReturnAddress (ReturnAddress (0));
          RuntimeError (706);  { Executed `x pow y' when complex x is zero and y < 0 }
          RestoreReturnAddress
	end;
	a := Re (z) / Im (z);
	b := Re (z) * a + Im (z);
	z := Cmplx (a / b, -1 / b);
      end;
      y := - y
    end;
  r := 1;
  while y <> 0 do
    begin
      if Odd (y) then r := r * z;
      y := y shr 1;
      if y <> 0 then z := Sqr (z)
    end
end;
-------------------------------------

BTW, are there any serious reasons why the complex "**" operator accepts only
_real_ exponents? As far as the RTS is concerned, it would suffice 
to change "Double" to "Complex" in the declaration of "Complex_Power",
but I expect that it needs some compiler magic as well.

Emil Jerabek