uss_dpfncs.s

VXWORKS源代码

/* Copyright 1991-1992 Wind River Systems, Inc. */
	.data
	.globl	_copyright_wind_river
	.long	_copyright_wind_river

/*
modification history
--------------------
01f,23aug92,jcf  changed bxxx to jxx.
01e,05jun92,kdl  changed external branches to jumps (and bsr's to jsr's).
01d,26may92,rrr  the tree shuffle
01c,30mar92,kdl  added include of "uss_fp.h"; commented-out ".set" directives
		 (SPR #1398).
01b,04oct91,rrr  passed through the ansification filter
		  -changed ASMLANGUAGE to _ASMLANGUAGE
		  -changed copyright notice
01a,28jan91,kdl  original US Software version.
*/


/*
DESCRIPTION

|       ttl     FPAC 68K/DPFNCS: IEEE Double Precision Functions
|DPFNCS idnt    1,0             ; IEEE Double Precision Functions
|                               ; DPFNCS.A68
|
|* * * * * * * * * *
|
|       68000 DPAC  --  Double Precision Functions
|
|       Copyright (C) 1984 - 1990 By
|       United States Software Corporation
|       14215 N.W. Science Park Drive
|       Portland, Oregon 97229
|
|       This software is furnished under a license and may be used
|       and copied only in accordance with the terms of such license
|       and with the inclusion of the above copyright notice.
|       This software or any other copies thereof may not be provided
|       or otherwise made available to any other person.  No title to
|       and ownership of the software is hereby transferred.
|
|       The information in this software is subject to change without
|       notice and should not be construed as a commitment by United
|       States Software Corporation.
|
|       Version:        See VSNLOG.TXT
|       Released:       01 Jan 1984
|
|* * * * * * * * * *
|
|
|  Double Precision Floating Point Function library.
|
|  These functions operate on double precision floating point values,
|  yielding double precision floating point results.  All functions
|  process the value on top of the stack and return a result on the
|  stack.  The XTOI function expects its integer argument in D0.
|
|
|  NAME       ARGUMENT                DESCRIPTION
|  ========== =====================   ====================================
|  EXP        -708.40 < a < 709.78    e to the power specified.
|                                     (if out of range, 0.0 or +INF)
|
|  LN         a > 0.0                 natural log of value specified.
|                                     (invalid opn if neg or zero)
|
|  XTOI       all value pairs         number raised to integer power.
|
|  ATN        any numeric value       arc whose tanget is given value
|                                     (no domain restriction, result in
|                                     the range is -PI/2 to +PI/2)
|
|  SQRT       a >= 0.0                square root of the given value
|                                     (invalid operation if negative)
|
|  COS                                All trigonometric functions
|  SIN                                expect a radian mode argument
|  TAN                                whose magnitude is < 65536 (2^16)
|

NOMANUAL
*/

#define	_ASMLANGUAGE

#include "uss_fp.h"

|       page
|
|       opt     BRS             ; Default forward branches to SHORT
|
        .globl  DPEXP
        .globl  DPLN
        .globl  DPLOG
        .globl  DPXTOI
        .globl  DPSQRT
        .globl  DPATN
        .globl  DPCOS
        .globl  DPSIN
        .globl  DPTAN
|
|
|       xref    DPMUL           ;Floating point operations
|       xref    DPDIV
|       xref    DPRDIV
|       xref    DPADD
|       xref    DAINT           ;Basic floating point functions
|       xref    DFLOAT
|       xref    DINT
|
|       xref    DNANRS
|       xref    DINFRS
|       xref    DUNFRS
|       xref    DZERRS
|
|
|
|
|       .set    DBIAS,1023
|
|
|       xref    FPERR
|
|
|       .set    ERUNF,1         |FPERR UNDERFLOW ERROR CODE
|       .set    EROVF,2         |FPERR OVERFLOW ERROR CODE
|       .set    ERNAN,3         |FPERR INVALID OPERATION ERROR CODE
/*
|
|
|
| ###   SUBTTL  DPEXP: Exponentiation of e Function
|       page
|
|
|  EXP:   Calculate e to the power specified by the stack argument.
|         The result is returned on the stack.
|         placed at DI.
|
|         FPERR is set upon overflow, underflow, or NaN argument.
|
|  Algorithm:
|  The argument is multiplied by 1/LN 2.  The integer portion of the
|  product becomes the two*s exponent.  A Taylor series approximation
|  about zero is used to compute 2**FRAC(ARG * 1/LN 2).
|
*/
        .text
|
DEXPCN:
        |dsw    0               | EXP CONSTANT LIST (2**X EVALUATION)
|
| ???     (LN 2)^16/16!   =   1.3570 24794 87551 47193 D-16
| ???   DC.L    $3CA38E89,$AE79F8B4
|         (LN 2)^15/15!   =   3.1324 36707 08842 86216 D-15
        .long   0x3CEC36E8,0x43B04022
|         (LN 2)^14/14!   =   6.7787 26354 82254 56335 D-14
        .long   0x3D331496,0x4D5878A9
|         (LN 2)^13/13!   =   1.3691 48885 39041 28881 D-12
        .long   0x3D781619,0x3166D0F9
|         (LN 2)^12/12!   =   2.5678 43599 34882 05142 D-11
        .long   0x3DBC3BD6,0x50FC2986
|         (LN 2)^11/11!   =   4.4455 38271 87081 14976 D-10
        .long   0x3DFE8CAC,0x7351BB25
|         (LN 2)^10/10!   =   7.0549 11620 80112 33299 D-09
        .long   0x3E3E4CF5,0x158B8ECA
|         (LN 2)^ 9/ 9!   =   1.0178 08600 92396 99728 D-07
        .long   0x3E7B5253,0xD395E7C4
|         (LN 2)^ 8/ 8!   =   1.3215 48679 01443 09488 D-06
        .long   0x3EB62C02,0x23A5C824
|         (LN 2)^ 7/ 7!   =   1.5252 73380 40598 40280 D-05
        .long   0x3EEFFCBF,0xC588B0C7
|         (LN 2)^ 6/ 6!   =   1.5403 53039 33816 09954 D-04
        .long   0x3F243091,0x2F86C787
|         (LN 2)^ 5/ 5!   =   1.3333 55814 64284 43423 D-03
        .long   0x3F55D87F,0xE78A6731
|         (LN 2)^ 4/ 4!   =   9.6181 29107 62847 71620 D-03
        .long   0x3F83B2AB,0x6FBA4E77
|         (LN 2)^ 3/ 3!   =   5.5504 10866 48215 79953 D-02
        .long   0x3FAC6B08,0xD704A0C0
|         (LN 2)^ 2/ 2!   =   2.4022 65069 59100 71233 D-01
        .long   0x3FCEBFBD,0xFF82C58F
|         (LN 2)^1/ 1!   =   6.9314 71805 59945 30942 D-01
DLN2:
        .long   0x3FE62E42,0xFEFA39EF
|                             1.0000 00000 00000 00000 D+00
        .long   0x3FF00000,0x00000000
|
|       .set    DNEXPC,16               | Number of constants in DEXP poly
|
|         1 / LN 2  =  1.4426 95040 88896 34074 D+00
DINLN2:
        .long   0x3FF71547,0x652B82FE
|
|
|
| ###   PUBLIC  DPEXP
|
DPEXP:
        |dsw    0
        bsr     DPFADJ          | Set-up for DP functions
|
        jvs     DFNANR          | NaN parm -> NaN result
        jcc     DEXP01          | J/ arg not INF
|
        jmi     DFUNFR          | -INF parm -> underflow result
        jra     DFPINR          | +INF parm -> positive INF result
|
DEXP01:
        movel   DINLN2+4,sp@-   | 1/LN 2 on stack
        movel   DINLN2+0,sp@-
        jsr     DPMUL           | Convert to 2^ function
|
        movew   sp@,d1
        andiw   #0x7FF0,d1
        cmpiw   #16*DBIAS+144,d1
        jhi     DEXP10          | J/ overflow or underflow  >= +/- 2^10
        jne     DEXP20          | J/ within range            < +/- 2^09
        tstb    sp@
        jpl     DEXP20          | J/ result >= 1.0
|
        movel   sp@,d2          | Examine top 12 bits of mantissa
        rorl    #8,d2
        andiw   #0x0FFF,d2
        cmpiw   #0x0FF8,d2
        jcs     DEXP20          | J/ result >= 2^-1022
|
DEXP10:
        tstb    sp@
        jpl     DFPINR          | J/ DP function +INF result
        jra     DFUNFR          | J/ DP function underflow result
|
DEXP20:
        movel   a7@(4),sp@-
        movel   a7@(4),sp@-     | Double the argument on the stack
|
        jsr     DINT
|
        movew   d1,a7@(8)       | Save result of DINT
|
        negl    d1
        negxl   d0              | Negate integer
        jsr     DFLOAT
        jsr     DPADD           | Special DFRAC
|
        pea     DEXPCN
        movel   sp@+,a6 | Point a6 to constant list
        moveq   #DNEXPC ,d0      | d0 holds the number of constants
        moveq   #DNEXPC, d0      | d0 holds the number of constants
        moveq   #16,d0      | XXX DEBUG DEBUG
        bsr     DXSER           | Polynomial series
|
        movew   a7@(8),d0       | Fetch power of two scaling
        lslw    #4,d0
        swap    d0
        clrw    d0
        addl    sp@+,d0 | Combine with DXSER result
|
        movel   sp@,a7@(4)      | Downshift result
        movel   d0,sp@
        jmp     a4@             | Return

/*
|
| ###   SUBTTL  DPLOG, DPLN: Logarithm Functions
|       page
|
|  LOGARITHM (BASE E) FUNCTION.
|
|  The natural logarithm of the double precision floating point
|  value on the stack is computed using a split domain polynomial
|  approximation.  The common log is computed by scaling the
|  result of the natural log computation.
|
|  TWOS = FLOAT(ARG.EXPONENT) * LOG(2)
|
|  The argument, (after the two's power scaling, 2.0 > arg' >= 1.0),
|  is sorted into one of five classes.  Within that class, a center
|  point with a known log value is used in combination with a 9th
|  degree Chebyshey polynomial approximation, to calculate a
|  logarithm.
|
|  The center point of each class is adjusted so that its natural
|  logarithm expressed in double precision form is accurate to
|  more than 80 bits.
|
|  log(arg) = two's log + center point log + polynomial approximation
|
|  Polynomial approximation:
|
|  log((1-t)/(1+t)) = t*(c1*t^8 + c2*t^6 + c3*t^4 + c4*t^2 + c5)
|         (accurate with 1.0E-17 for t within 0.1 of 1.0)
|
|  t = - (center - arg')/(center + arg')  [to calc log of arg'/center]
|
*/

DLNCNS:
        |dsw    0               | LN constants list
|
|         c1  =   0.22330 022
        .long   0x3FCC951A,0x03000000
|         c2  =   0.28571 20487
        .long   0x3FD2491B,0x30450000
|         c3  =   0.40000 00018 947
        .long   0x3FD99999,0x9BA26900
|         c4  =   0.66666 66666 66105
        .long   0x3FE55555,0x55554192
|         c5  =   2.00000 00000 00000 00000
        .long   0x40000000,0x00000000
|
|       .set    DNLNCN,5                | Number of constants in DLN poly
|
|
DLNCEN:
        |dsw    0               | Center Points
|         Center  #1 = 1.0000 00000 00000 00000
        .long   0x3FF00000,0x00000000
|         Center  #2 = 1.1892 07115 00272 05647  (close to 2^(1/4))
        .long   0x3FF306FE,0x0A31B713
|         Center  #3 = 1.4142 13562 37309 38548  (close to 2^(1/2))
        .long   0x3FF6A09E,0x667F3BC7
|         Center  #4 = 1.6817 92830 50742 54625  (close to 2^(3/4))
        .long   0x3FFAE89F,0x995AD39D
|         Center  #5 = 2.0000 00000 00000 00000
| ---     shifted to 1.0 center..
|
DLNLOG:
        |dsw    0               | DLNCEN values adj for exact ln values
|         LN(Center #2) = 0.17328 67951 39985 90522
        .long   0x3FC62E42,0xFEFA39E0
|         LN(Center #3) = 0.34657 35902 79971 81045
        .long   0x3FD62E42,0xFEFA39E0
|         LN(Center #4) = 0.51986 03854 19956 82749
        .long   0x3FE0A2B2,0x3F3BAB60
|
|
|         1 / LN 10  =  0.43429 44819 03251 82765
DILN10:
        .long   0x3FDBCB7B,0x1526E50E
|
|
|
| ###   PUBLIC  DPLN
| ###   PUBLIC  DPLOG
|
DPLOG:
        |dsw    0
        moveq   #-1,d0          | Signal common log scaling
        jra     DLN000
|
DPLN:
        |dsw    0
        clrb    d0              | Signal no post calculation scaling
|
DLN000:
        bsr     DPFADJ          | Adjust parameter on stack
|
        jvs     DFNANR          | J/ NaN arg -> NaN result
        jmi     DFNANR          | J/ Neg arg -> NaN result
        jeq     DFMINR          | J/ 0.0 arg -> -INF result
        jcs     DFPINR          | J/ +INF arg -> +INF result
|
        moveb   d0,a7@(10)      | Save natural/common flag
|
        movew   sp@,d1          | Prepare to calc parm sign and exp
        lsrw    #4,d1           | Position exponent (right justified)
        subiw   #DBIAS,d1
        movew   d1,a7@(8)       | /* Save two's exponent value */
|
        lslw    #4,d1           | Scale parameter
        subw    d1,sp@
|
        clrw    d0              | Set class number to 0
        movel   sp@,d1
        lsrl    #8,d1           | d1.W = 0xEMMM
|
        cmpiw   #0xF000+371,d1
        jcs     DLN050          | J/ < 1 +  371/4096  (1.090576)
        addqw   #8,d0
        cmpiw   #0xF000+1216,d1
        jcs     DLN050          | J/ < 1 + 1216/4096  (1.296875)
        addqw   #8,d0
        cmpiw   #0xF000+2221,d1
        jcs     DLN050          | J/ < 1 + 2221/4096  (1.542236)
        addqw   #8,d0
        cmpiw   #0xF000+3416,d1
        jcs     DLN050          | J/ < 1 + 3416/4096  (1.833984)
        addqw   #1,a7@(8)       | /* Class 4: Bump two's exp, then class 0 */
        clrw    d0
        subib   #0x10,a7@(1)    | Reduce scaled exponent
|
DLN050:
        moveb   d0,a7@(11)      | Save center number
|
        pea     DLNCEN
        movel   sp@+,a2
        movel   a2@(4,d0:W),sp@- |Place center on stack
        movel   a2@(0,d0:W),sp@-
|
        movel   a7@(12),sp@-    | Copy scaled argument
        movel   a7@(12),sp@-
|
        movel   a2@(4,d0:W),sp@- |Place center on stack
        movel   a2@(0,d0:W),sp@-
|
        jsr     DPADD           | /* Calc  (CENTER + PARM') */
|
        movel   sp@,d0          | Exch top stack item w/ 3rd
        movel   a7@(16),d1
        movel   d1,sp@
        movel   d0,a7@(16)
        movel   a7@(4),d0
        movel   a7@(20),d1
        movel   d1,a7@(4)
        movel   d0,a7@(20)
|
        bset    #7,a7@(8)       | Negate CENTER
|
        jsr     DPADD           | /* Calc  - (CENTER - PARM') */
        jsr     DPRDIV          | /* T = - (CENTER-PARM')/(CENTER+PARM) */
|
        movel   a7@(4),sp@-     | Duplicate t
        movel   a7@(4),sp@-
|
        pea     DLNCNS
        movel   sp@+,a6 | Poly approx
        moveq   #DNLNCN,d0
        bsr     DX2SER
|
        jsr     DPMUL           | Times t
        clrw    d0
        moveb   a7@(11),d0
        jeq     DLN060          | J/ center = 0 -> no additive log val
|
        pea     DLNLOG-8
        movel   sp@+,a6