sysdeps/libm-i387/e_acoshl.S

   1 /* ix87 specific implementation of arcsinh.
   2    Copyright (C) 1996, 1997 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 Library General Public License as
   8    published by the Free Software Foundation; either version 2 of the
   9    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    Library General Public License for more details.
  15
  16    You should have received a copy of the GNU Library General Public
  17    License along with the GNU C Library; see the file COPYING.LIB.  If not,
  18    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19    Boston, MA 02111-1307, USA.  */
  20
  21 #include <machine/asm.h>
  22
  23 #ifdef __ELF__
  24         .section .rodata
  25 #else
  26         .text
  27 #endif
  28
  29         .align ALIGNARG(4)
  30         /* Please note that we use double value for 1.0.  This number
  31            has an exact representation and so we don't get accuracy
  32            problems.  The advantage is that the code is simpler.  */
  33         ASM_TYPE_DIRECTIVE(one,@object)
  34 one:    .double 1.0
  35         ASM_SIZE_DIRECTIVE(one)
  36         /* It is not important that this constant is precise.  It is only
  37            a value which is known to be on the safe side for using the
  38            fyl2xp1 instruction.  */
  39         ASM_TYPE_DIRECTIVE(limit,@object)
  40 limit:  .double 0.29
  41         ASM_SIZE_DIRECTIVE(limit)
  42
  43 #ifdef PIC
  44 #define MO(op) op##@GOTOFF(%edx)
  45 #else
  46 #define MO(op) op
  47 #endif
  48
  49         .text
  50 ENTRY(__ieee754_acoshl)
  51         movl    12(%esp), %ecx
  52         andl    $0xffff, %ecx
  53         cmpl    $0x3fff, %ecx
  54         jl      5f                      // < 1 => invalid
  55         fldln2                          // log(2)
  56         fldt    4(%esp)                 // x : log(2)
  57         cmpl    $0x4020, %ecx
  58         ja      3f                      // x > 2^34
  59 #ifdef  PIC
  60         call    1f
  61 1:      popl    %edx
  62         addl    $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
  63 #endif
  64         cmpl    $0x4000, %ecx
  65         ja      4f                      // x > 2
  66
  67         // 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
  68         fsubl   MO(one)                 // x-1 : log(2)
  69         fld     %st                     // x-1 : x-1 : log(2)
  70         fmul    %st(1)                  // (x-1)^2 : x-1 : log(2)
  71         fadd    %st(1)                  // x-1+(x-1)^2 : x-1 : log(2)
  72         fadd    %st(1)                  // 2*(x-1)+(x-1)^2 : x-1 : log(2)
  73         fsqrt                           // sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2)
  74         faddp                           // x-1+sqrt(2*(x-1)+(x-1)^2) : log(2)
  75         fcoml   MO(limit)
  76         fnstsw
  77         sahf
  78         ja      2f
  79         fyl2xp1                         // log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
  80         ret
  81
  82 2:      faddl   MO(one)                 // x+sqrt(2*(x-1)+(x-1)^2) : log(2)
  83         fyl2x                           // log(x+sqrt(2*(x-1)+(x-1)^2))
  84         ret
  85
  86         // x > 2^34 => y = log(x) + log(2)
  87         .align ALIGNARG(4)
  88 3:      fyl2x                           // log(x)
  89         fldln2                          // log(2) : log(x)
  90         faddp                           // log(x)+log(2)
  91         ret
  92
  93         // 2^34 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1)))
  94         .align ALIGNARG(4)
  95 4:      fld     %st                     // x : x : log(2)
  96         fadd    %st, %st(1)             // x : 2*x : log(2)
  97         fld     %st                     // x : x : 2*x : log(2)
  98         fmul    %st(1)                  // x^2 : x : 2*x : log(2)
  99         fsubl   MO(one)                 // x^2-1 : x : 2*x : log(2)
 100         fsqrt                           // sqrt(x^2-1) : x : 2*x : log(2)
 101         faddp                           // x+sqrt(x^2-1) : 2*x : log(2)
 102         fdivrl  MO(one)                 // 1/(x+sqrt(x^2-1)) : 2*x : log(2)
 103         fsubrp                          // 2*x+1/(x+sqrt(x^2)-1) : log(2)
 104         fyl2x                           // log(2*x+1/(x+sqrt(x^2-1)))
 105         ret
 106
 107         // x < 1 => NaN
 108         .align ALIGNARG(4)
 109 5:      fldz
 110         fdiv    %st, %st(0)
 111         ret
 112 END(__ieee754_acoshl)