/* 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 shuffle01c,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 notice01a,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-01DLN2:        .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+00DINLN2:        .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 82765DILN10:        .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