This code provides math support for the power() function for IEEE 80-bit floating point extended precision. The assembly code provides the basic functionality which is augmented with code in C++ that tests common cases and verifies that the arguments are in the function's domain. There are additional supporting functions which are not shown in this sample. Other than minor formatting for HTML presentation, this code is exactly as found in the original source code.
title Copyright (c) 1993..1996 by Michael Lee Finney. All rights reserved.
;
; Extended basicPower(Extended x, Extended y)
;
; This function returns x ** y where x > 0 and y > 0. The
; algorithm used is...
;
; a <- y * log2(x)
; b <- round(a)
; c <- a - b
; d <- 2**c-1
; e <- d + 1
; f <- e * 2**b
; return f
;
align 16
basicPower__FrT1 equ $
?basicPower__FrT1 equ $
; st(0) = x
; st(1) = y
fyl2x
; st(0) = y * log2(x)
fld st(0)
; st(0) = y * log2(x)
; st(1) = y * log2(x)
frndint
; st(0) = round(y * log2(x))
; st(1) = y * log2(x)
fxch
; st(0) = y * log2(x)
; st(1) = round(y * log2(x))
fsub st,st(1)
; st(0) = y * log2(x) - round(y * log2(x))
; st(1) = round(y * log2(x))
f2xm1
; st(0) = 2 ** [y * log2(x) - round(y * log2(x))] - 1
; st(1) = round(y * log2(x))
fld1
; st(0) = 1.0
; st(1) = 2 ** [y * log2(x) - round(y * log2(x))] - 1
; st(2) = round(y * log2(x))
fadd
; st(0) = 2 ** [y * log2(x) - round(y * log2(x))]
; st(1) = round(y * log2(x))
fscale
; st(0) = 2 ** [y * log2(x)] = x ** y
; st(1) = round(y * log2(x))
fxch
; st(0) = round(y * log2(x))
; st(1) = x ** y
fstp st(0)
; st(0) = x ** y
ret
/* basicPower(extended x, extended y)
*
* This function returns x ** y where x > 0 and y > 0.
*
* Note: This function is implemented in MathSupport.asm
* because it cannot be implemented directly in C++.
*/
// extended basicPower(
// extended x,
// extended y) is
// code
// nil;
// end (0.0)
/* power(extended x, extended y)
*
* This module returns x ** y. If x < 0 then y must be an integer and
* if x == 0 then y must be a positive value.
*/
extended power(
extended x,
extended y) is
local
extended a;
extended r;
code
if (x < 0.0)
if (isIntegral(y))
if (x == -1.0)
r = isOdd(y) ? -1.0 : 1.0;
elif (y < 0.0)
if (y == -1.0)
r = 1.0 / x;
elif (y == -2.0)
r = 1.0 / (x * x);
else
{
a = 1.0 / basicPower(-x, -y);
r = isOdd(y) ? -a : a;
}
elif (y == 0.0)
r = 1.0;
elif (y == 1.0)
r = x;
elif (y == 2.0)
r = x * x;
else
{
a = basicPower(-x, y);
r = isOdd(y) ? -a : a;
}
else
domainError ("power(x, y) -- x negative and y non-integer");
elif (x == 0.0)
if (y > 0.0)
r = 0.0;
else
domainError ("power(x, y) -- x zero and y not positive");
elif (x == 1.0)
r = 1.0;
elif (y < 0.0)
if (y == -0.5)
r = 1.0 / squareRoot(x);
elif (y == -1.0)
r = 1.0 / x;
elif (y == -2.0)
r = 1.0 / (x * x);
else
r = 1.0 / basicPower(x, -y);
elif (y == 0.0)
r = 1.0;
elif (y == 0.5)
r = squareRoot(x);
elif (y == 1.0)
r = x;
elif (y == 2.0)
r = x * x;
else
r = basicPower(x, y);
end (r)