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#include "libm_synonyms.h" 34 35 .data 36 .align 4 37.mhundred: .float -100.0 38 39 ENTRY(expm1) 40 movl 8(%esp),%ecx / ecx <-- hi_32(x) 41 andl $0x7fffffff,%ecx / ecx <-- hi_32(|x|) 42 cmpl $0x3fe62e42,%ecx / Is |x| < ln(2)? 43 jb .shortcut / If so, take a shortcut. 44 je .check_tail / |x| may be only slightly < ln(2) 45 cmpl $0x7ff00000,%ecx / hi_32(|x|) >= hi_32(INF)? 46 jae .not_finite / if so, x is not finite 47.finite_non_special: / Here, ln(2) < |x| < INF 48 fldl 4(%esp) / push x 49 50 subl $8,%esp / save RP and set round-to-64-bits 51 fstcw (%esp) 52 movw (%esp),%ax 53 movw %ax,4(%esp) 54 orw $0x0300,%ax 55 movw %ax,(%esp) 56 fldcw (%esp) 57 58 fldl2e / push log2e }not for xtndd_dbl 59 fmulp %st,%st(1) / z = x*log2e }not for xtndd_dbl 60 fld %st(0) / duplicate stack top 61 frndint / [z],z 62 / [z] != 0, compute exp(x) and then subtract one to get expm1(x) 63 fxch / z,[z] 64 fsub %st(1),%st / z-[z],[z] 65 f2xm1 / 2**(z-[z])-1,[z] 66 / avoid spurious underflow when scaling to compute exp(x) 67 PIC_SETUP(1) 68 flds PIC_L(.mhundred) 69 PIC_WRAPUP 70 fucom %st(2) / if -100 !< [z], then use -100 71 fstsw %ax 72 sahf 73 jb .got_int_part 74 fxch %st(2) 75.got_int_part: 76 fstp %st(0) / 2**(z-[z])-1,max([z],-100) 77 fld1 / 1,2**(z-[z])-1,max([z],-100) 78 faddp %st,%st(1) / 2**(z-[z]) ,max([z],-100) 79 fscale / exp(x) ,max([z],-100) 80 fld1 / 1,exp(x) ,max([z],-100) 81 fxch / exp(x),1 ,max([z],-100) 82 fsubp %st,%st(1) / exp(x)-1 ,max([z],-100) 83 fstp %st(1) 84 85 fstcw (%esp) / restore old RP 86 movw (%esp),%dx 87 andw $0xfcff,%dx 88 movw 4(%esp),%cx 89 andw $0x0300,%cx 90 orw %dx,%cx 91 movw %cx,(%esp) 92 fldcw (%esp) 93 add $8,%esp 94 95 ret 96 97.check_tail: 98 movl 4(%esp),%edx / edx <-- lo_32(x) 99 cmpl $0xfefa39ef,%edx / Is |x| slightly < ln(2)? 100 ja .finite_non_special / branch if |x| slightly > ln(2) 101.shortcut: 102 / Here, |x| < ln(2), so |z| = |x*log2(e)| < 1, 103 / whence z is in f2xm1's domain. 104 fldl 4(%esp) / push x 105 fldl2e / push log2e }not for xtndd_dbl 106 fmulp %st,%st(1) / z = x*log2e }not for xtndd_dbl 107 f2xm1 / 2**(x*log2(e))-1 = e**x - 1 108 ret 109 110.not_finite: 111 / Here, flags still have settings from execution of 112 / cmpl $0x7ff00000,%ecx / hi_32(|x|) > hi_32(INF)? 113 ja .NaN_or_pinf / if not, x may be +/- INF 114 movl 4(%esp),%edx / edx <-- lo_32(x) 115 cmpl $0,%edx / lo_32(x) = 0? 116 jne .NaN_or_pinf / if not, x is NaN 117 movl 8(%esp),%eax / eax <-- hi_32(x) 118 andl $0x80000000,%eax / here, x is infinite, but +/-? 119 jz .NaN_or_pinf / branch if x = +INF 120 fld1 / Here, x = -inf, so return -1 121 fchs 122 ret 123 124.NaN_or_pinf: 125 / Here, x = NaN or +inf, so load x and return immediately. 126 fldl 4(%esp) 127 fwait 128 ret 129 .align 4 130 SET_SIZE(expm1) 131