{ power is a function for Virtual Pascal to raise a number to a power,
  like fortran's y**x or C's pow(y, x).  Inspired by Richard Startz's
  '8087 Applications and Programming for the IBM PC and other PCs'

  Your program must include this unit in its 'uses' statement.

  Ex call:      answer:=power(y, x);

  The math coprocessor is assumed present, and the function
  uses the built-in asssember.  You MUST have a coprocessor
  to use this function.

  The number to be raised is not allowed to be negative; if a negative
  number is passed, it's absolute value is used.

  If a zero is passed, the function returns zero,
  so zero to -any- power is zero.  This should be suitable
  for most uses.

  The power may be any valid number.

  There is no other error checking, and no error codes are set.
  Infinity can be generated, and is returned as such; the calling
  code has to handle this as needed.  When INF is returned, it can
  handled by testing against the largest possible floating point number,
  defined in this unit as MAXFLOAT.
  }
unit MATH;
interface
function power(number, power: extended) : extended;
const
   MaxFloat : extended = 1.1E4932;

implementation
uses Use32;
{$N+}

function power(number, power: extended) : extended;
var
   StatWord : smallword;
   ContWord : smallword;
   temp : smallword;
const
   Half : word = $3F000000;

begin
   asm
      fld     power          { put power on top of stack }
      fld     number         { put number on top of stack }
      fabs                   { make sure it's not negative }
      ftst                   { is it's zero .... }
      fstsw   AX
      sahf
      jne     {at}not_zero
      jmp     {at}get_out       { ...yes, so get out with 0.0 on T.O.S }
{at}not_zero:                   { now st(1) = power, st = number }
      fyl2x                  { st = number * log power base 2 }
      fstcw   ContWord       { save the control word to restore later }
      fstcw   temp           { to change rounding control to round-down }
      and     temp, $F33F    { clear out RoundingControl bits }
      or      temp, $0400    { set for Round-down }
      fldcw   temp           { put it in control word }
      fld     st(0)          { push copy of number }
      frndint                { st = z1, st(1) = z }
      fldcw   ContWord       { restore initial control word }
      fsub    st(1), st      { st(1) = z2 = (z-z1) }
      fxch                   { st = z2, st(1) = z1 }
      fld     half           { 1/2 -> top of stack }
      fxch                   { st = z2, st(1) = 1/2 }
      fprem                  { st is z2 or z2=1/2; if z2=1/2, c1=1 }
      fstsw   StatWord
      fstp    st(1)          { now flags are set, so get rid of the half }
      f2xm1                  { st is now  ((2 to the st) - 1)     }
      fld1
      faddp   st(1), st
      test    byte ptr StatWord+1, 00000010B { st has z2 if bit 1 on }
      jz      {at}was_z2                        { else it had z2 - 1/2 }
      fld1
      fadd    st, st(0)
      fsqrt                  { so multiply it by the square root of 2 }
      fmulp   st(1), st
{at}was_z2:
      fscale       { just need to scale by 2**st(1) }
      fstp   st(1)
{at}get_out:
      fstp {at}result
   end;
end;

begin
end.

___
 X KWQ/2 1.2i X I'm not young enough to know everything anymore...
--- Maximus/2 2.02