lib/libm/arch/i386/e_exp.S

   1 /*
   2  * Copyright (c) 1993,94 Winning Strategies, Inc.
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  * 3. All advertising materials mentioning features or use of this software
  14  *    must display the following acknowledgement:
  15  *      This product includes software developed by Winning Strategies, Inc.
  16  * 4. The name of the author may not be used to endorse or promote products
  17  *    derived from this software without specific prior written permission.
  18  *
  19  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  20  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  21  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  22  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  23  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  24  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  25  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  26  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  27  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  28  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  29  *
  30  * $NetBSD: e_exp.S,v 1.12 2002/02/27 16:32:46 christos Exp $
  31  * $FreeBSD: src/lib/msun/i387/e_exp.S,v 1.8.2.1 2000/07/10 09:16:28 obrien Exp $
  32  * $DragonFly: src/lib/libm/arch/i386/e_exp.S,v 1.1 2005/07/26 21:15:19 joerg Exp $
  33  */
  34
  35 /*
  36  * Written by:
  37  *      J.T. Conklin (jtc@wimsey.com), Winning Strategies, Inc.
  38  */
  39
  40 #include <machine/asm.h>
  41
  42 #include "abi.h"
  43
  44
  45 /* e^x = 2^(x * log2(e)) */
  46 ENTRY(exp)
  47 #ifndef __i386__
  48         /*
  49          * XXX: This code is broken and needs to be merged with the i386 case.
  50          */
  51         fstcw   -12(%rsp)
  52         movw    -12(%rsp),%dx
  53         orw     $0x0180,%dx
  54         movw    %dx,-16(%rsp)
  55         fldcw   -16(%rsp)
  56         movsd   %xmm0,-8(%rsp)
  57         fldl    -8(%rsp)
  58
  59         fldl2e
  60         fmulp                           /* x * log2(e) */
  61         fld     %st(0)
  62         frndint                         /* int(x * log2(e)) */
  63         fxch    %st(1)
  64         fsub    %st(1),%st              /* fract(x * log2(e)) */
  65         f2xm1                           /* 2^(fract(x * log2(e))) - 1 */
  66         fld1
  67         faddp                           /* 2^(fract(x * log2(e))) */
  68         fscale                          /* e^x */
  69         fstp    %st(1)
  70
  71         fstpl   -8(%rsp)
  72         movsd   -8(%rsp),%xmm0
  73         fldcw   -12(%rsp)
  74         ret
  75 #else
  76         /*
  77          * If x is +-Inf, then the subtraction would give Inf-Inf = NaN.
  78          * Avoid this.  Also avoid it if x is NaN for convenience.
  79          */
  80         movl    8(%esp),%eax
  81         andl    $0x7fffffff,%eax
  82         cmpl    $0x7ff00000,%eax
  83         jae     x_Inf_or_NaN
  84
  85         fldl    4(%esp)
  86
  87         /*
  88          * Ensure that the rounding mode is to nearest (to give the smallest
  89          * possible fraction) and that the precision is as high as possible.
  90          * We may as well mask interrupts if we switch the mode.
  91          */
  92         fstcw   4(%esp)
  93         movl    4(%esp),%eax
  94         andl    $0x0300,%eax
  95         cmpl    $0x0300,%eax            /* RC == 0 && PC == 3? */
  96         je      1f                      /* jump if mode is good */
  97         movl    $0x137f,8(%esp)
  98         fldcw   8(%esp)
  99 1:
 100         fldl2e
 101         fmulp                           /* x * log2(e) */
 102         fst     %st(1)
 103         frndint                         /* int(x * log2(e)) */
 104         fst     %st(2)
 105         fsubrp                          /* fract(x * log2(e)) */
 106         f2xm1                           /* 2^(fract(x * log2(e))) - 1 */
 107         fld1
 108         faddp                           /* 2^(fract(x * log2(e))) */
 109         fscale                          /* e^x */
 110         fstp    %st(1)
 111         je      1f
 112         fldcw   4(%esp)
 113 1:
 114         ret
 115
 116 x_Inf_or_NaN:
 117         /*
 118          * Return 0 if x is -Inf.  Otherwise just return x, although the
 119          * C version would return (x + x) (Real Indefinite) if x is a NaN.
 120          */
 121         cmpl    $0xfff00000,8(%esp)
 122         jne     x_not_minus_Inf
 123         cmpl    $0,4(%esp)
 124         jne     x_not_minus_Inf
 125         fldz
 126         ret
 127
 128 x_not_minus_Inf:
 129         fldl    4(%esp)
 130         ret
 131 #endif