* MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
* M68000 Hi-Performance Microprocessor Division
* M68040 Software Package
*
* M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
* All rights reserved.
*
* THE SOFTWARE is provided on an "AS IS" basis and without warranty.
* To the maximum extent permitted by applicable law,
* MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
* INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
* PARTICULAR PURPOSE and any warranty against infringement with
* regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
* and any accompanying written materials.
*
* To the maximum extent permitted by applicable law,
* IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
* (INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
* PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
* OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
* SOFTWARE. Motorola assumes no responsibility for the maintenance
* and support of the SOFTWARE.
*
* You are hereby granted a copyright license to use, modify, and
* distribute the SOFTWARE so long as this entire notice is retained
* without alteration in any modified and/or redistributed versions,
* and that such modified versions are clearly identified as such.
* No licenses are granted by implication, estoppel or otherwise
* under any patents or trademarks of Motorola, Inc.
*
* stwotox.sa 3.1 12/10/90
*
* stwotox --- 2**X
* stwotoxd --- 2**X for denormalized X
* stentox --- 10**X
* stentoxd --- 10**X for denormalized X
*
* Input: Double-extended number X in location pointed to
* by address register a0.
*
* Output: The function values are returned in Fp0.
*
* Accuracy and Monotonicity: The returned result is within 2 ulps in
* 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
* result is subsequently rounded to double precision. The
* result is provably monotonic in double precision.
*
* Speed: The program stwotox takes approximately 190 cycles and the
* program stentox takes approximately 200 cycles.
*
* Algorithm:
*
* twotox
* 1. If |X| > 16480, go to ExpBig.
*
* 2. If |X| < 2**(-70), go to ExpSm.
*
* 3. Decompose X as X = N/64 + r where |r| <= 1/128. Furthermore
* decompose N as
* N = 64(M + M') + j, j = 0,1,2,...,63.
*
* 4. Overwrite r := r * log2. Then
* 2**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
* Go to expr to compute that expression.
*
* tentox
* 1. If |X| > 16480*log_10(2) (base 10 log of 2), go to ExpBig.
*
* 2. If |X| < 2**(-70), go to ExpSm.
*
* 3. Set y := X*log_2(10)*64 (base 2 log of 10). Set
* N := round-to-int(y). Decompose N as
* N = 64(M + M') + j, j = 0,1,2,...,63.
*
* 4. Define r as
* r := ((X - N*L1)-N*L2) * L10
* where L1, L2 are the leading and trailing parts of log_10(2)/64
* and L10 is the natural log of 10. Then
* 10**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
* Go to expr to compute that expression.
*
* expr
* 1. Fetch 2**(j/64) from table as Fact1 and Fact2.
*
* 2. Overwrite Fact1 and Fact2 by
* Fact1 := 2**(M) * Fact1
* Fact2 := 2**(M) * Fact2
* Thus Fact1 + Fact2 = 2**(M) * 2**(j/64).
*
* 3. Calculate P where 1 + P approximates exp(r):
* P = r + r*r*(A1+r*(A2+...+r*A5)).
*
* 4. Let AdjFact := 2**(M'). Return
* AdjFact * ( Fact1 + ((Fact1*P) + Fact2) ).
* Exit.
*
* ExpBig
* 1. Generate overflow by Huge * Huge if X > 0; otherwise, generate
* underflow by Tiny * Tiny.
*
* ExpSm
* 1. Return 1 + X.
*
STWOTOX IDNT 2,1 Motorola 040 Floating Point Software Package
FMOVE.L FP1,N(a6) ...N = ROUND-TO-INT(64 X)
MOVE.L d2,-(sp)
LEA EXPTBL,a1 ...LOAD ADDRESS OF TABLE OF 2^(J/64)
FMOVE.L N(a6),FP1 ...N --> FLOATING FMT
MOVE.L N(a6),D0
MOVE.L D0,d2
ANDI.L #$3F,D0 ...D0 IS J
ASL.L #4,D0 ...DISPLACEMENT FOR 2^(J/64)
ADDA.L D0,a1 ...ADDRESS FOR 2^(J/64)
ASR.L #6,d2 ...d2 IS L, N = 64L + J
MOVE.L d2,D0
ASR.L #1,D0 ...D0 IS M
SUB.L D0,d2 ...d2 IS M', N = 64(M+M') + J
ADDI.L #$3FFF,d2
MOVE.W d2,ADJFACT(a6) ...ADJFACT IS 2^(M')
MOVE.L (sp)+,d2
*--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
*--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
*--ADJFACT = 2^(M').
*--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
EXPBIG:
*--|X| IS LARGE, GENERATE OVERFLOW IF X > 0; ELSE GENERATE UNDERFLOW
*--REGISTERS SAVE SO FAR ARE FPCR AND D0
MOVE.L X(a6),D0
TST.L D0
BLT.B EXPNEG
bclr.b #7,(a0) ;t_ovfl expects positive value
bra t_ovfl
EXPNEG:
bclr.b #7,(a0) ;t_unfl expects positive value
bra t_unfl
xdef stentoxd
stentoxd:
*--ENTRY POINT FOR 10**(X) FOR DENORMALIZED ARGUMENT
fmove.l d1,fpcr ...set user's rounding mode/precision
Fmove.S #:3F800000,FP0 ...RETURN 1 + X
move.l (a0),d0
or.l #$00800001,d0
fadd.s d0,fp0
bra t_frcinx
xdef stentox
stentox:
*--ENTRY POINT FOR 10**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
FMOVEM.X (a0),FP0 ...LOAD INPUT, do not set cc's
FMOVE.L FP1,N(a6) ...N=INT(X*64*LOG10/LOG2)
MOVE.L d2,-(sp)
LEA EXPTBL,a1 ...LOAD ADDRESS OF TABLE OF 2^(J/64)
FMOVE.L N(a6),FP1 ...N --> FLOATING FMT
MOVE.L N(a6),D0
MOVE.L D0,d2
ANDI.L #$3F,D0 ...D0 IS J
ASL.L #4,D0 ...DISPLACEMENT FOR 2^(J/64)
ADDA.L D0,a1 ...ADDRESS FOR 2^(J/64)
ASR.L #6,d2 ...d2 IS L, N = 64L + J
MOVE.L d2,D0
ASR.L #1,D0 ...D0 IS M
SUB.L D0,d2 ...d2 IS M', N = 64(M+M') + J
ADDI.L #$3FFF,d2
MOVE.W d2,ADJFACT(a6) ...ADJFACT IS 2^(M')
MOVE.L (sp)+,d2
*--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
*--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
*--ADJFACT = 2^(M').
*--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
expr:
*--FPCR, FP2, FP3 ARE SAVED IN ORDER AS SHOWN.
*--ADJFACT CONTAINS 2**(M'), FACT1 + FACT2 = 2**(M) * 2**(J/64).
*--FP0 IS R. THE FOLLOWING CODE COMPUTES
*-- 2**(M'+M) * 2**(J/64) * EXP(R)
FMOVE.X FP0,FP1
FMUL.X FP1,FP1 ...FP1 IS S = R*R
FMOVE.D EXPA5,FP2 ...FP2 IS A5
FMOVE.D EXPA4,FP3 ...FP3 IS A4
FMUL.X FP1,FP2 ...FP2 IS S*A5
FMUL.X FP1,FP3 ...FP3 IS S*A4
FADD.D EXPA3,FP2 ...FP2 IS A3+S*A5
FADD.D EXPA2,FP3 ...FP3 IS A2+S*A4
FMUL.X FP1,FP2 ...FP2 IS S*(A3+S*A5)
FMUL.X FP1,FP3 ...FP3 IS S*(A2+S*A4)
FADD.D EXPA1,FP2 ...FP2 IS A1+S*(A3+S*A5)
FMUL.X FP0,FP3 ...FP3 IS R*S*(A2+S*A4)
FMUL.X FP1,FP2 ...FP2 IS S*(A1+S*(A3+S*A5))
FADD.X FP3,FP0 ...FP0 IS R+R*S*(A2+S*A4)
FADD.X FP2,FP0 ...FP0 IS EXP(R) - 1
*--FINAL RECONSTRUCTION PROCESS
*--EXP(X) = 2^M*2^(J/64) + 2^M*2^(J/64)*(EXP(R)-1) - (1 OR 0)
