sysdeps/i386/fpu/e_acoshl.S

   1 /* ix87 specific implementation of arcsinh.
   2    Copyright (C) 1996-1997, 2005, 2011-2012 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
   5
   6    The GNU C Library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Lesser General Public
   8    License as published by the Free Software Foundation; either
   9    version 2.1 of the License, or (at your option) any later version.
  10
  11    The GNU C Library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with the GNU C Library; if not, see
  18    <http://www.gnu.org/licenses/>.  */
  19
  20 #include <machine/asm.h>
  21
  22         .section .rodata.cst8,"aM",@progbits,8
  23
  24         .p2align 3
  25         /* Please note that we use double value for 1.0.  This number
  26            has an exact representation and so we don't get accuracy
  27            problems.  The advantage is that the code is simpler.  */
  28         ASM_TYPE_DIRECTIVE(one,@object)
  29 one:    .double 1.0
  30         ASM_SIZE_DIRECTIVE(one)
  31         /* It is not important that this constant is precise.  It is only
  32            a value which is known to be on the safe side for using the
  33            fyl2xp1 instruction.  */
  34         ASM_TYPE_DIRECTIVE(limit,@object)
  35 limit:  .double 0.29
  36         ASM_SIZE_DIRECTIVE(limit)
  37
  38 #ifdef PIC
  39 #define MO(op) op##@GOTOFF(%edx)
  40 #else
  41 #define MO(op) op
  42 #endif
  43
  44         .text
  45 ENTRY(__ieee754_acoshl)
  46         movl    12(%esp), %ecx
  47         andl    $0xffff, %ecx
  48         cmpl    $0x3fff, %ecx
  49         jl      5f                      // < 1 => invalid
  50         fldln2                          // log(2)
  51         fldt    4(%esp)                 // x : log(2)
  52         cmpl    $0x4020, %ecx
  53         ja      3f                      // x > 2^34
  54 #ifdef  PIC
  55         LOAD_PIC_REG (dx)
  56 #endif
  57         cmpl    $0x4000, %ecx
  58         ja      4f                      // x > 2
  59
  60         // 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
  61         fsubl   MO(one)                 // x-1 : log(2)
  62         fld     %st                     // x-1 : x-1 : log(2)
  63         fmul    %st(1)                  // (x-1)^2 : x-1 : log(2)
  64         fadd    %st(1)                  // x-1+(x-1)^2 : x-1 : log(2)
  65         fadd    %st(1)                  // 2*(x-1)+(x-1)^2 : x-1 : log(2)
  66         fsqrt                           // sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2)
  67         faddp                           // x-1+sqrt(2*(x-1)+(x-1)^2) : log(2)
  68         fcoml   MO(limit)
  69         fnstsw
  70         sahf
  71         ja      2f
  72         fyl2xp1                         // log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
  73         ret
  74
  75 2:      faddl   MO(one)                 // x+sqrt(2*(x-1)+(x-1)^2) : log(2)
  76         fyl2x                           // log(x+sqrt(2*(x-1)+(x-1)^2))
  77         ret
  78
  79         // x > 2^34 => y = log(x) + log(2)
  80         .align ALIGNARG(4)
  81 3:      fyl2x                           // log(x)
  82         fldln2                          // log(2) : log(x)
  83         faddp                           // log(x)+log(2)
  84         ret
  85
  86         // 2^34 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1)))
  87         .align ALIGNARG(4)
  88 4:      fld     %st                     // x : x : log(2)
  89         fadd    %st, %st(1)             // x : 2*x : log(2)
  90         fld     %st                     // x : x : 2*x : log(2)
  91         fmul    %st(1)                  // x^2 : x : 2*x : log(2)
  92         fsubl   MO(one)                 // x^2-1 : x : 2*x : log(2)
  93         fsqrt                           // sqrt(x^2-1) : x : 2*x : log(2)
  94         faddp                           // x+sqrt(x^2-1) : 2*x : log(2)
  95         fdivrl  MO(one)                 // 1/(x+sqrt(x^2-1)) : 2*x : log(2)
  96         fsubrp                          // 2*x+1/(x+sqrt(x^2)-1) : log(2)
  97         fyl2x                           // log(2*x+1/(x+sqrt(x^2-1)))
  98         ret
  99
 100         // x < 1 => NaN
 101         .align ALIGNARG(4)
 102 5:      fldz
 103         fdiv    %st, %st(0)
 104         ret
 105 END(__ieee754_acoshl)
 106 strong_alias (__ieee754_acoshl, __acoshl_finite)