{ Turbo Pascal functions to compute Powers.
  Original written by Paul F. Hultquist, published in PC Tech Journal,
  Vol. 3, No. 5 (May '85), p.213.

  Additions and modifications by Roy Collins, 5/8/85:
    * Error conditions caused the program to display an error message
      and halt.  Changed to display error message, beep, and return
      a value of zero.
    * Since Paul stated that PWRI is faster than PWRR when raising a
      real to an integer power, I added code to PWRR to check for an
      'Y' value and to call PWRI when appropraitte.
    * Added new testing procedures. }

procedure beep;
    { Test use only.
      Draw attention to error conditions. }
begin
  sound(1000); delay(150); nosound;
  delay(100);
  sound(2000); delay(150); nosound;
end; (* proc beep *)

function pwri(x:real; n:integer):real;
    { PWRI performs the tests necessary to eliminate non-computable cases
      of finding X (real) to the power N (integer).  It calls upon function
      RLSCAN to do the actual computation after it has, for example,
      replaced a negative with a positive one (it does a reciprocation
      after return from RLSCAN in that case. Error conditions return a
      zero value. }

  function rlscan(x:real; n:integer):real;
    { This function scans the positive exponent from right to left to
      determine a sequence of multiplications and squarings that produce
      X (real) to the power N (integer) in a near-minimum number of
      multiplications. The algorithm is Algorithm A, p. 442, Vol. 2,
      2nd Ed. of Knuth: "The Art of Computer Programming: Seminumerical
      Algorithms", Addison-Wesley, 1981. }
  var
    y, z : real;
    o    : boolean;
    bign : integer;
  begin
    bign := n;
    y := 1.0;
    z := x;
    while bign > 0 do begin
      o := odd(bign);
      bign := bign div 2;
      if o then begin
        y := y * z;
        rlscan := y;
        end;
      z := z * z;
      end;
  end; (* func rlscan *)

begin
  if n > 0 then
    pwri := rlscan(x,n)
  else
  if (x <> 0.0) and (n < 0) then begin
    n := -n;
    pwri := 1.0 / rlscan(x,n);
    end
  else
  if (n = 0) and (x <> 0) then
    pwri := 1.0
  else
  if (n = 0) and (x = 0) then begin
    writeln('0 to the 0 power.');
    beep;
    pwri := 0.0;
    end
  else
  if (n < 0) and (x = 0) then begin
    writeln('Division by zero.');
    beep;
    pwri := 0.0;
    end;
end; (* func pwri *)

function pwrr(x,y:real):real;
    { PWRR finds X (real) to the power Y (real) using algorithms and
      exponentials. If Y is an integer PWRI is faster.  The function
      eliminates the undefined cases and the case where the result is
      complex. For these error conditions, zero is returned. }
    { Added by Roy Collins, 5/8/85:
         Since PWRI is faster for cases where Y is an integer,
         if Y is an acceptable integer, use PWRI to compute x**y }
begin
  if (y = int(y)) and (abs(y) <= 32767) then
    pwrr := pwri(x,trunc(y))
  else
  if x > 0 then
    pwrr := exp(y*ln(x))
  else
  if x < 0 then begin
    writeln('X < 0.');
    beep;
    pwrr := 0.0;
    end
  else
  if (x=0) and (y=0) then begin
    writeln('0 to the 0 power.');
    beep;
    pwrr := 0.0;
    end
  else
  if (x=0) and (y<0) then begin
    writeln('0 ti a negative power.');
    beep;
    pwrr := 0.0;
    end
  else
    pwrr := 0.0;
end; (* func pwrr *)

var
  ch : char;
  x, y : real;
begin (* test pwri and pwrr functions *)
  repeat
    writeln;
    write('Test R)eal or I)nteger power function (or Q)uit)? (R/N/Q)');
    repeat
      read(kbd,ch);
      ch := upcase(ch);
    until ch in ['R','I','Q'];
    if ch = 'Q' then
      halt;
    writeln(ch);
    if ch = 'R' then
      writeln('Raise X (real) to the Y (real) power')
    else
      writeln('Raise X (real) to the Y (integer) power');
    write('Enter X value: ');
    readln(x);
    write('Enter Y value: ');
    readln(y);
    if ch = 'R' then
      writeln('Real power = ',pwrr(x,y):12:11)
    else
      writeln('Integer power = ',int(pwri(x,trunc(y))):12:11);
  until false;
end.                        