sysdeps/i386/fpu/e_acoshl.S

   1 /* ix87 specific implementation of arcsinh.
   2    Copyright (C) 1996-2020 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    <https://www.gnu.org/licenses/>.  */
  19
  20 #include <machine/asm.h>
  21 #include <libm-alias-finite.h>
  22
  23         .section .rodata.cst8,"aM",@progbits,8
  24
  25         .p2align 3
  26         /* Please note that we use double value for 1.0.  This number
  27            has an exact representation and so we don't get accuracy
  28            problems.  The advantage is that the code is simpler.  */
  29         .type one,@object
  30 one:    .double 1.0
  31         ASM_SIZE_DIRECTIVE(one)
  32         /* It is not important that this constant is precise.  It is only
  33            a value which is known to be on the safe side for using the
  34            fyl2xp1 instruction.  */
  35         .type limit,@object
  36 limit:  .double 0.29
  37         ASM_SIZE_DIRECTIVE(limit)
  38
  39 #ifdef PIC
  40 #define MO(op) op##@GOTOFF(%edx)
  41 #else
  42 #define MO(op) op
  43 #endif
  44
  45         .text
  46 ENTRY(__ieee754_acoshl)
  47         movl    12(%esp), %ecx
  48         andl    $0xffff, %ecx
  49         cmpl    $0x3fff, %ecx
  50         jl      5f                      // < 1 => invalid
  51         fldln2                          // log(2)
  52         fldt    4(%esp)                 // x : log(2)
  53         cmpl    $0x4020, %ecx
  54         ja      3f                      // x > 2^34
  55 #ifdef  PIC
  56         LOAD_PIC_REG (dx)
  57 #endif
  58         cmpl    $0x4000, %ecx
  59         ja      4f                      // x > 2
  60
  61         // 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
  62         fsubl   MO(one)                 // x-1 : log(2)
  63         fabs                            // acosh(1) is +0 in all rounding modes
  64         fld     %st                     // x-1 : x-1 : log(2)
  65         fmul    %st(1)                  // (x-1)^2 : x-1 : log(2)
  66         fadd    %st(1)                  // x-1+(x-1)^2 : x-1 : log(2)
  67         fadd    %st(1)                  // 2*(x-1)+(x-1)^2 : x-1 : log(2)
  68         fsqrt                           // sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2)
  69         faddp                           // x-1+sqrt(2*(x-1)+(x-1)^2) : log(2)
  70         fcoml   MO(limit)
  71         fnstsw
  72         sahf
  73         ja      2f
  74         fyl2xp1                         // log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
  75         ret
  76
  77 2:      faddl   MO(one)                 // x+sqrt(2*(x-1)+(x-1)^2) : log(2)
  78         fyl2x                           // log(x+sqrt(2*(x-1)+(x-1)^2))
  79         ret
  80
  81         // x > 2^34 => y = log(x) + log(2)
  82         .align ALIGNARG(4)
  83 3:      fyl2x                           // log(x)
  84         fldln2                          // log(2) : log(x)
  85         faddp                           // log(x)+log(2)
  86         ret
  87
  88         // 2^34 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1)))
  89         .align ALIGNARG(4)
  90 4:      fld     %st                     // x : x : log(2)
  91         fadd    %st, %st(1)             // x : 2*x : log(2)
  92         fld     %st                     // x : x : 2*x : log(2)
  93         fmul    %st(1)                  // x^2 : x : 2*x : log(2)
  94         fsubl   MO(one)                 // x^2-1 : x : 2*x : log(2)
  95         fsqrt                           // sqrt(x^2-1) : x : 2*x : log(2)
  96         faddp                           // x+sqrt(x^2-1) : 2*x : log(2)
  97         fdivrl  MO(one)                 // 1/(x+sqrt(x^2-1)) : 2*x : log(2)
  98         fsubrp                          // 2*x+1/(x+sqrt(x^2)-1) : log(2)
  99         fyl2x                           // log(2*x+1/(x+sqrt(x^2-1)))
 100         ret
 101
 102         // x < 1 => NaN
 103         .align ALIGNARG(4)
 104 5:      fldz
 105         fdiv    %st, %st(0)
 106         ret
 107 END(__ieee754_acoshl)
 108 libm_alias_finite (__ieee754_acoshl, __acoshl)