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 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, write to the Free 18 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 19 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) 34one: .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) 40limit: .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 50ENTRY(__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 611: 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 822: 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) 883: 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) 954: 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) 1095: fldz 110 fdiv %st, %st(0) 111 ret 112END(__ieee754_acoshl) 113