1/* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21/* 22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24/* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 .file "expm1.s" 30 31#include "libm.h" 32LIBM_ANSI_PRAGMA_WEAK(expm1,function) 33 34 .data 35 .align 4 36.mhundred: .float -100.0 37 38 ENTRY(expm1) 39 movl 8(%esp),%ecx / ecx <-- hi_32(x) 40 andl $0x7fffffff,%ecx / ecx <-- hi_32(|x|) 41 cmpl $0x3fe62e42,%ecx / Is |x| < ln(2)? 42 jb .shortcut / If so, take a shortcut. 43 je .check_tail / |x| may be only slightly < ln(2) 44 cmpl $0x7ff00000,%ecx / hi_32(|x|) >= hi_32(INF)? 45 jae .not_finite / if so, x is not finite 46.finite_non_special: / Here, ln(2) < |x| < INF 47 fldl 4(%esp) / push x 48 49 subl $8,%esp / save RP and set round-to-64-bits 50 fstcw (%esp) 51 movw (%esp),%ax 52 movw %ax,4(%esp) 53 orw $0x0300,%ax 54 movw %ax,(%esp) 55 fldcw (%esp) 56 57 fldl2e / push log2e }not for xtndd_dbl 58 fmulp %st,%st(1) / z = x*log2e }not for xtndd_dbl 59 fld %st(0) / duplicate stack top 60 frndint / [z],z 61 / [z] != 0, compute exp(x) and then subtract one to get expm1(x) 62 fxch / z,[z] 63 fsub %st(1),%st / z-[z],[z] 64 f2xm1 / 2**(z-[z])-1,[z] 65 / avoid spurious underflow when scaling to compute exp(x) 66 PIC_SETUP(1) 67 flds PIC_L(.mhundred) 68 PIC_WRAPUP 69 fucom %st(2) / if -100 !< [z], then use -100 70 fstsw %ax 71 sahf 72 jb .got_int_part 73 fxch %st(2) 74.got_int_part: 75 fstp %st(0) / 2**(z-[z])-1,max([z],-100) 76 fld1 / 1,2**(z-[z])-1,max([z],-100) 77 faddp %st,%st(1) / 2**(z-[z]) ,max([z],-100) 78 fscale / exp(x) ,max([z],-100) 79 fld1 / 1,exp(x) ,max([z],-100) 80 fxch / exp(x),1 ,max([z],-100) 81 fsubp %st,%st(1) / exp(x)-1 ,max([z],-100) 82 fstp %st(1) 83 84 fstcw (%esp) / restore old RP 85 movw (%esp),%dx 86 andw $0xfcff,%dx 87 movw 4(%esp),%cx 88 andw $0x0300,%cx 89 orw %dx,%cx 90 movw %cx,(%esp) 91 fldcw (%esp) 92 add $8,%esp 93 94 ret 95 96.check_tail: 97 movl 4(%esp),%edx / edx <-- lo_32(x) 98 cmpl $0xfefa39ef,%edx / Is |x| slightly < ln(2)? 99 ja .finite_non_special / branch if |x| slightly > ln(2) 100.shortcut: 101 / Here, |x| < ln(2), so |z| = |x*log2(e)| < 1, 102 / whence z is in f2xm1's domain. 103 fldl 4(%esp) / push x 104 fldl2e / push log2e }not for xtndd_dbl 105 fmulp %st,%st(1) / z = x*log2e }not for xtndd_dbl 106 f2xm1 / 2**(x*log2(e))-1 = e**x - 1 107 ret 108 109.not_finite: 110 / Here, flags still have settings from execution of 111 / cmpl $0x7ff00000,%ecx / hi_32(|x|) > hi_32(INF)? 112 ja .NaN_or_pinf / if not, x may be +/- INF 113 movl 4(%esp),%edx / edx <-- lo_32(x) 114 cmpl $0,%edx / lo_32(x) = 0? 115 jne .NaN_or_pinf / if not, x is NaN 116 movl 8(%esp),%eax / eax <-- hi_32(x) 117 andl $0x80000000,%eax / here, x is infinite, but +/-? 118 jz .NaN_or_pinf / branch if x = +INF 119 fld1 / Here, x = -inf, so return -1 120 fchs 121 ret 122 123.NaN_or_pinf: 124 / Here, x = NaN or +inf, so load x and return immediately. 125 fldl 4(%esp) 126 fwait 127 ret 128 .align 4 129 SET_SIZE(expm1) 130