/* 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/FPFNCS: IEEE Single Precision Functions|FPFNCS idnt    1,0             ; IEEE Single Precision Functions|                               ; FPFNCS.A68||* * * * * * * * * *||       68000 FPAC  --  Single 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||* * * * * * * * * *|||  Single Precision Floating Point Function library.||  These functions operate on single precision floating point values,|  yielding single 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        -87.43 < a < 88.72      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 to forward branches SHORT|        .globl  FPEXP        .globl  FPLN        .globl  FPLOG        .globl  FPXTOI        .globl  FPATN        .globl  FPSQRT        .globl  FPCOS        .globl  FPSIN        .globl  FPTAN||       xref    FPMUL   ;Floating point operations|       xref    FPDIV|       xref    FPRDIV|       xref    FPADD|       xref    FLOAT|       xref    INT||       xref    FNANRS|       xref    FINFRS|       xref    FUNFRS|       xref    FZERRS||       xref    FPERR|||       .set    FBIAS,127||       .set    ERUNF,1 |FPERR UNDERFLOW ERROR CODE|       .set    EROVF,2 |FPERR OVERFLOW ERROR CODE|       .set    ERNAN,3 |FPERR INVALID OPERATION ERROR CODE/*|||| ###   SUBTTL  FPEXP: Exponentiation of e Function|       page|||  EXP:   Calculate e to the power specified by the stack argument.|         The result is returned on the stack.||         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|FEXPCN:        |dsw    0       | EXP CONSTANT LIST (2**X EVALUATION)||         (LN 2)^ 9/ 9!   =   1.0178 08601E-07        .long   0x33DA929E|         (LN 2)^ 8/ 8!   =   1.3215 48679E-06        .long   0x35B16011|         (LN 2)^ 7/ 7!   =   1.5252 73380E-05        .long   0x377FE5FE|         (LN 2)^ 6/ 6!   =   1.5403 53039E-04        .long   0x39218489|         (LN 2)^ 5/ 5!   =   1.3333 55815E-03        .long   0x3AAEC3FF|         (LN 2)^ 4/ 4!   =   9.6181 29108E-03        .long   0x3C1D955B|         (LN 2)^ 3/ 3!   =   5.5504 10866E-02        .long   0x3D635847|         (LN 2)^ 2/ 2!   =   2.4022 65070E-01        .long   0x3E75FDF0|         (LN 2)^ 1/ 1!   =   6.9314 71806E-01FLN2:        .long   0x3F317218|                             1.0000 00000E+00        .long   0x3F800000||       .set    FNEXPC,10               | Number of constants in FEXP poly||         1 / LN 2  =  1.4426 95040 88896 34074 D+00FINLN2:        .long   0x3FB8AA3B|||| ###   PUBLIC  FPEXP|FPEXP:        |dsw    0        bsr     FPFADJ          | Set-up for FP functions|        jvs     FFNANR          | NaN parm -> NaN result        jcc     FEXP01          | J/ arg not INF|        jmi     FFUNFR          | -INF parm -> underflow result        jra     FFPINR          | +INF parm -> positive INF result|FEXP01:        movel   FINLN2,sp@-     | 1/LN 2 on stack        jsr     FPMUL           | Convert to 2^ function|        movew   sp@,d1        rorw    #7,d1           | Exponent in d1.B        cmpib   #FBIAS+6,d1        jhi     FEXP10          | J/ overflow or underflow  >= +/- 2^7        jne     FEXP20          | J/ within range            < +/- 2^6        tstb    sp@        jpl     FEXP20          | J/ result >= 1.0|        cmpiw   #0xF800+0x0100+FBIAS+6,d1       | Check ms 7 man bits (+s+exp)        jcs     FEXP20                  | J/ result >= 2^-126|FEXP10:        tstb    sp@        jpl     FFPINR          | J/ FP function +INF result        jra     FFUNFR          | J/ FP function underflow result|FEXP20:        movel   sp@,sp@-        | Duplicate the argument on the stack|        jsr     INT|        movew   d0,a7@(4)       | Save result of INT|        negl    d0              | Negate integer        jsr     FLOAT        jsr     FPADD           | Special FFRAC|        pea     FEXPCN        movel   sp@+,a6 | Point a6 to constant list        moveq   #FNEXPC,d0      | d0 holds the number of constants        bsr     FXSER           | Polynomial series|        movew   a7@(4),d0       | Fetch power of two scaling        lslw    #7,d0        addw    d0,sp@          | Combine with FXSER result|        movel   sp@+,sp@        | Downshift result        jmp     a4@             | Return/*|| ###   SUBTTL  FPLOG, FPLN: Logarithm Functions|       page||  LOGARITHM (BASE E) FUNCTION.||  The natural logarithm of the single 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 7th|  degree Chebyshey polynomial approximation, to calculate a|  logarithm.||  The center point of each class is adjusted so that its natural|  logarithm expressed in single precision form is accurate to|  more than 40 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]|*/FLNCNS:        |dsw    0               | LN constants list||         c1  =   0.28571 20487        .long   0x3E9248DA|         c2  =   0.40000 00019        .long   0x3ECCCCCD|         c3  =   0.66666 66667        .long   0x3F2AAAAB|         c4  =   2.00000 00000        .long   0x40000000||       .set    FNLNCN,4                | Number of constants in FLN poly||FLNCEN:        |dsw    0               | Center Points|         Center  #1 = 1.0000 00000 00000 00000        .long   0x3F800000|         Center  #2 = 1.1892 07098     (close to 2^(1/4))        .long   0x3F9837F0|         Center  #3 = 1.4142 13521     (close to 2^(1/2))        .long   0x3FB504F3|         Center  #4 = 1.6817 92733     (close to 2^(3/4))        .long   0x3FD744FC|         Center  #5 = 2.0000 00000| ---     shifted to 1.0 center..|FLNLOG:        |dsw    0               | FLNCEN values adj for exact ln values|         LN(Center #2) = 0.17328 67807        .long   0x3E317217|         LN(Center #3) = 0.34657 35614        .long   0x3EB17217|         LN(Center #4) = 0.51986 03272        .long   0x3F051591|||         1 / LN 10  =  0.43429 44819FILN10:        .long   0x3EDE5BD9|||| ###   PUBLIC  FPLN| ###   PUBLIC  FPLOG|FPLOG:        |dsw    0        moveq   #-1,d0          | Signal common log scaling        jra     FLN000|FPLN:        |dsw    0        clrb    d0              | Signal no post calculation scaling|FLN000:        bsr     FPFADJ          | Adjust parameter on stack|        jvs     FFNANR          | J/ NaN arg -> NaN result        jmi     FFNANR          | J/ Neg arg -> NaN result        jeq     FFMINR          | J/ 0.0 arg -> -INF result        jcs     FFPINR          | J/ +INF arg -> +INF result|        moveb   d0,a7@(6)       | Save natural/common flag|        movew   sp@,d1          | Prepare to calc parm sign and exp        lsrw    #7,d1           | Position exponent (right justified)        subiw   #FBIAS,d1        movew   d1,a7@(4)       | /* Save two's exponent value */|        lslw    #7,d1           | Scale parameter        subw    d1,sp@|        clrw    d0              | Set class number to 0        movel   sp@,d1        lsrl    #8,d1           | d1.W = 0xEMMM|        cmpiw   #0x8000+2966,d1        jcs     FLN050          | J/ < 1 +  2966/32768  (1.090515)        addqw   #4,d0        cmpiw   #0x8000+9727,d1        jcs     FLN050          | J/ < 1 +  9727/32768  (1.296844)        addqw   #4,d0        cmpiw   #0x8000+17767,d1        jcs     FLN050          | J/ < 1 + 17767/32768  (1.542206)        addqw   #4,d0        cmpiw   #0x8000+27329,d1        jcs     FLN050          | J/ < 1 + 27329/32768  (1.834015)        addqw   #1,a7@(4)       | /* Class 4: Bump two's exp, then class 0 */        clrw    d0        subib   #0x80,a7@(1)    | Reduce scaled exponent|FLN050:        moveb   d0,a7@(7)       | Save center number|        pea     FLNCEN        movel   sp@+,a2        movel   a2@(0,d0:W),sp@- |Place center on stack|        movel   a7@(4),sp@-     | Copy scaled argument|        movel   a2@(0,d0:W),sp@- |Place center on stack|        jsr     FPADD           | /* Calc  (CENTER + PARM') */|        movel   sp@,d0          | Exch top stack item w/ 3rd        movel   a7@(8),d1        movel   d1,sp@        movel   d0,a7@(8)|        bset    #7,a7@(4)       | Negate CENTER|        jsr     FPADD           | /* Calc  - (CENTER - PARM') */        jsr     FPRDIV          | /* T = - (CENTER-PARM')/(CENTER+PARM) */|        movel   sp@,sp@-        | Duplicate t|        pea     FLNCNS        movel   sp@+,a6 | Poly approx        moveq   #FNLNCN,d0        bsr     FX2SER|        jsr     FPMUL           | Times t        clrw    d0