Files
TeslaRel410/BORLAND/BC45/SOURCE/RTL/SOURCE/MATHLIB/POWL.CAS
T
CydandClaude Fable 5 63312e07f9 source410: literal 4.10 source reconstruction + BC++ 4.52 fleet toolchain archived
- BORLAND/: Borland C++ 4.52 (chosen over 4.5 by byte-match: CODE/RP/CW32.LIB
  is identical to 4.52's install lib). BCC32/TLINK32/TLIB/MAKE run natively on
  Win11; CODE/BT/OPT.MAK is the shipped BTL4OPT.EXE's exact flag recipe
  (extender = Borland PowerPack DPMI32, not Phar Lap TNT).
- restoration/source410/: the literal 1995-form reconstruction of the missing
  BT game source (never mixed into CODE/). Round 1-3 state:
  * 6 of 10 surviving original TUs COMPILE CLEAN under the period toolchain
    (BTMSSN, BTCNSL, BTSCNRL, BTTEAM, BTL4MODE, BTL4ARND) - first builds
    since 1996.
  * BT_L4/BTL4APP.CPP pilot reconstruction: 12/12 functions, Fail() lands on
    its binary-recorded line 400 exactly.
  * BT/BTCNSL.HPP: console wire IDs recovered from the binary's ctors
    (Killed=9, Damaged=10, ScoreUpdate=13, DeathWithoutHonor=15 [T1];
    TeamScore=12 flagged [T4]).
  * MUNGA/: 8 engine-header backfills back-dated from the BT412 WinTesla tree
    (VDATA numbering decomp-verified; AUDREND's OpenAL-era virtual removed -
    the period compiler is the drift detector).
  * Tooling: backdate.py (WinTesla->1995 header transform), compile410.sh
    (per-TU verification sweep under authentic OPT.MAK flags).
  * README: corrected roadmap - MECH.HPP is the capstone grown with the mech
    TU reconstructions; BTREG.CPP green = the header-family milestone.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
2026-07-19 07:33:26 -05:00

399 lines
13 KiB
Plaintext

/*------------------------------------------------------------------------
* filename - powl.cas
*
* function(s)
* powl - long double power function, x^y
*-----------------------------------------------------------------------*/
/*
* C/C++ Run Time Library - Version 6.5
*
* Copyright (c) 1987, 1994 by Borland International
* All Rights Reserved.
*
*/
#pragma inline
#include <asmrules.h>
#include <_math.h>
#include <math.h>
#include <errno.h>
#include <stddef.h>
#define EXTPROC1(x) asm call near ptr (x)
static unsigned short NANLOGL [5] = {0,0,0,0xC024, 0xFFFF};
/*--------------------------------------------------------------------------*
Name powl - power function, x^y
Usage long double powl(long double x, long double y);
Prototype in math.h
Description Return the value of x to the power of y. If x is zero then
y must be positive, and if x is negative then y must be
integral.
First the special cases Y == 0 or X == 0 or X == infinity
are detected and given standard answers.
Two methods of calculation are used, depending upon whether
Y is an integer of less than 64 bits. If not, then
X^Y = 2^(Log2(X) * Y)
= DoExps ( DoLogs (X, Y), not scaled)
If Y is an integer then it can be represented as a binary
number
Y = B0 + B1.2 + B2.+ .. Bn.2^n
where the B coefficients are 0 or 1, and Bn is always 1,
for some n. The power of X is then calculated as:
Z = X;
while (n-- > 0)
Z *= Z;
if (Bn) Z *= X;
which is the standard trick for fast integral powers. It
works for any X, positive or negative, if Y is not zero. In
practice it will run faster than the DoExps (DoLogs())
method for |Y| < 100, roughly, and slower for larger
powers. Such large powers are very rare in actual usage.
Powers greater than 2^64 in theory may be integers.
However, it is also likely that such large numbers have
lost precision (especially when you consider that the C
data type "long double" is 53 bits precise). These will be
treated as if fractional. If X is positive and very close
to 1.0, then an answer may be possible, but if X is
negative an exception is generated. The rationale for the
exception is that if the least bits of Y have been lost
then it is not possible to be sure whether the result
should be positive or negative, so there is a total loss of
precision.
Return value Return the value of x to the power of y.
When the correct value would overflow, powl returns the
value _LHUGE_VAL.
If the argument x passed to powl is less than or equal to 0
and y is not a whole number, then errno is set to
EDOM Domain error
If x and y are both zero, then the return value is 1.0 and
there is no error. Many C compilers consider this a DOMAIN
error as technically 0^0 is undefined. There are continuous
function f(x) and g(x) such that f(0) = 0 and g(0) = 0, but
with the limit of f(x)/g(x) as x tends to 0 being any real
number. However, there is an elementary theorem that states
that f(x) and g(x) are analytic and nonzero, then the limit
is always 1. Thus in a finite precision computing
environment, it is hard to imagine a situation where a
number other than 1 is desirable.
*---------------------------------------------------------------------------*/
/*
Identical to expl(), but calls __matherrl with different arguments.
*/
static long double near pascal __expl (long double x)
{
asm FLD LONGDOUBLE (x)
asm mov ax, x [8] /* select exponent */
asm and ah, 7Fh /* remove sign bit */
asm cmp ax, 3fffh+13
asm jb exp_OK /* expl (+-2^13) is the limit for long double */
exp_tooBig:
asm mov ax, 0FFFFh /* force extreme */
asm ja exp_excess
asm mov ax, x [6]
exp_excess:
asm test BY0 (x [9]), 80h
asm jnz exp_tooTiny
asm cmp ax, 0B172h
asm jb exp_OK
asm mov si, OVERFLOW
asm jmp short exp_err
exp_tooTiny:
asm cmp ax, 0B16Ch
asm jb exp_OK
asm mov si, UNDERFLOW
exp_err:
asm FSTP ST(0) /* discard ST */
#pragma warn -ret
/* should use args to powl, but have no access */
return __matherrl (_SI, "powl", NULL, NULL,
(UNDERFLOW == _SI) ? 0.0 : _LHUGE_VAL);
#pragma warn .ret
exp_OK:
__expld();
#pragma noretval
return;
}
#pragma warn -powl
#pragma warn -use
long double _FARFUNC powl (long double x, long double y)
{
long double temp; /* also used as a 64-bit integer */
char negate = 0; /* boolean, negate after exp() ? */
volatile unsigned int Control;
volatile unsigned int Status;
asm FLD LONGDOUBLE (x)
asm mov bx, 7FFFh /* mask just the exponent */
asm mov ax, x [8]
asm and ax, bx
asm jz powl_ofZero
asm cmp ax, bx
asm je powl_ofInfinity
asm FLD LONGDOUBLE (y)
asm mov ax, y [8]
asm and ax, bx
asm jnz temp1
asm jmp powl_toZero
temp1:
asm cmp ax, bx
asm je powl_toInfinity
asm jmp powl_normal
/*** Special cases ***/
/*
Raising any number to infinity is treated as a range error.
*/
powl_toInfinity:
asm FSTP LONGDOUBLE (temp) /* propagate Y thru to result */
asm jmp short powl_discard
/*
Powers of infinity are range errors.
*/
powl_ofInfinity:
asm FSTP LONGDOUBLE (temp) /* propagate X thru to result */
asm mov ax, y[8]
asm or ax, ax
asm jge powl_overflow /* jump if exponent nonnegative */
asm mov si, UNDERFLOW
temp = 0.0;
asm jmp short powl_complain
/*
Arrive here if Y is zero. The zero'th power of any number is 1.
*/
powl_toZero:
asm FSTP ST(0) /* discard Y */
asm FSTP ST(0) /* discard X */
asm FLD1
#pragma noretval
return; /* 1.0 */
powl_discard:
asm FSTP ST(0) /* discard X */
powl_overflow:
asm mov si, OVERFLOW
asm jmp short powl_complain
/*
Powers of 0 are (EDOM, 1, 0) as Y ranges over (negative, zero, positive).
*/
powl_ofZero:
asm FSTP ST(0) /* discard X */
asm mov ax, y [8]
asm or ax, ax /* was Y positive ? */
asm jg powl_zero
asm mov si, DOMAIN
asm je powl_zz
temp = _LHUGE_VAL;
asm jmp short powl_complain
powl_zz:
temp = 1.0;
powl_complain:
#pragma warn -ret
return __matherrl (_SI, "powl", &x, &y, temp);
#pragma warn .ret
powl_zero:
asm FLDZ
#pragma noretval
return; /* 0.0 */
/*** End of Special Cases ***/
/*
If arrived here then both x and y seem to be ordinary numbers.
*/
powl_normal:
asm FCLEX
asm FRNDINT
asm FSTSW W0 (Status) /* is Y an integer */
asm FWAIT
asm test BY0 (Status), 20h /* precision error if not */
asm jz powl_integral
asm FSTP ST(0) /* discard Y */
/*
Arrive here if the exponent exceeds integer range or if it contains
a fractional part. Calculate using Log and Exp functions. Just
x is on 87-stack.
*/
powl_fractional:
/* make sure that x > 0 */
asm FTST
asm FSTSW W0 (Status) /* is Y an integer */
asm FWAIT
asm mov ax, Status
asm sahf
asm jae powl_log
asm FSTP ST(0)
temp = *((long double *) NANLOGL);
asm mov si, DOMAIN
asm jmp short powl_complain
powl_log:
/* arg is > 0, so log cannot fail */
#ifdef _Windows
_f87_Log();
#else
asm _FAST_ (_FLOG_)
#endif
/* Disable underflow and overflow exceptions temporarily.
*/
asm fstcw Control /* save old control word */
asm fwait
asm mov ax, Control
asm mov cx, ax /* save copy in cx */
asm or ax, 18h /* mask under/overflow */
asm mov Control, ax
asm fldcw Control
asm FLD LONGDOUBLE (y)
asm FMUL
/* Check for exceptions.
*/
asm fstsw Status /* save the status */
asm fclex /* clear exceptions */
asm mov Control, cx
asm fldcw Control /* reload control word */
asm test word ptr Status, 18h /* did under/overflow occur? */
asm jz powl_noexc /* if so, discard result */
asm jmp powl_discard
powl_noexc:
asm sub sp, 10
asm mov bx, sp
asm FSTP LONGDOUBLE (SS_ [bx])
EXTPROC1 (__expl)
if (negate)
{
asm FCHS
}
#pragma noretval
return;
/*
If arrived here then Y is some integer of up to 64 bits and has
been copied to temp. Y is ST(0), X is ST(1), AX is exponent of Y.
*/
powl_integral:
asm sub ax, 3FFFh /* remove the bias */
asm cmp ax, 63 /* AX = n, the exponent */
asm jb powl_trueIntegral
asm FSTP ST(0) /* discard Y */
powl_fracjmp:
asm jmp short powl_fractional
/*
The shift-and-add method is not accurate for extreme powers since
round off errors are magnified. However, we cannot simply call for
evaluation like fractional powers because X may be negative and
fractional negative powers are treated as exceptions.
*/
powl_trueIntegral:
asm cmp al, 12
asm jb powl_shiftAndAdd
powl_unsafeRange:
asm FISTP qword ptr (temp) /* store an integer copy of Y */
asm test BY0 (x [9]), 80h /* X less than 0 ? */
asm jz powl_fracjmp /* X not signed, so no worry */
asm FCHS /* make X absolute */
asm test BY0 (temp), 01h /* odd or even ? */
asm jz powl_fracjmp /* even powers are positive */
/*
If we arrive here then X was negative and Y was odd. Calculate with
abs(X) and then negate result.
*/
negate = 1;
asm jmp short powl_fracjmp
/*
Arrive here for modest integral powers of any number. We must also
check for overflow, by making a worst-case check on log (X^Y). If
it has a potential to overflow, then we use the exp(log()) method.
*/
powl_shiftAndAdd:
asm mov bx, x [8]
asm shl bx, 1
asm sub bx, 7FFEh /* BX estimates log2 (X) */
asm mov dx, bx
asm xchg cx, ax
asm inc cx /* 2^CL is max possible Y */
asm shl bx, cl /* multiply BX by max Y */
asm sar bx, cl
asm dec cx
asm xchg ax, cx
asm cmp bx, dx /* did BX lose any bits ? */
asm jne powl_unsafeRange
asm FLD ST (1) /* Z = X */
asm mov dx, y [6]
asm shl dx,1
powl_iWhileBit:
asm dec al
asm jl powl_maybeInverse
asm FMUL ST(0), ST(0) /* Z *= Z */
asm shl dx, 1
asm jnc powl_iWhileBit
asm FMUL ST(0), ST(2) /* Z *= X */
asm jmp short powl_iWhileBit
powl_maybeInverse:
asm FSTP ST(1) /* overwrite Y */
asm test BY0 (y [9]), 80h /* was Y a negative power ? */
asm FSTP ST(1) /* overwrite X */
asm jz powl_iDone
asm FLD1
asm FDIVRP ST(1), ST(0) /* if so, invert result. */
powl_iDone:
#pragma noretval
return;
}