overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22/*******************************************************************************
23* *
24* File frexpldexp.c, *
25* *
26* Contains: Functions frexp(x) and ldexp(x), *
27* Implementation of frexp and ldexp functions for the PowerPC. *
28* *
29* Copyright � 1991-2001 Apple Computer, Inc. All rights reserved. *
30* *
31* Written by A. Sazegari, started on January 1991, *
32* Modified and ported by Robert A. Murley (ram) for Mac OS X. *
33* *
34* A MathLib v4 file. *
35* *
36* December 03 1992: first rs6000 implementation. *
37* October 05 1993: added special cases for NaN and � in frexp. *
38* May 27 1997: improved the performance of frexp by eliminating the *
39* switch statement. *
40* June 13 2001: rewrote frexp to eliminate calls to scalb and logb. *
41* replaced DblInHex typedef with hexdouble. *
42* September 06 2001: added #if __ppc__. *
43* added routine descriptions. *
44* September 09 2001: added macros to detect PowerPC and correct compiler. *
45* September 10 2001: added more comments. *
46* September 18 2001: added <Limits.h> and <CoreServices/CoreServices.h>. *
47* October 08 2001: removed <Limits.h> and <CoreServices/CoreServices.h>. *
48* added #define INT16_MAX. *
49* changed compiler errors to warnings. *
50* November 06 2001: commented out warning about Intel architectures. *
51* *
52* W A R N I N G: *
53* These routines require a 64-bit double precision IEEE-754 model. *
54* They are written for PowerPC only and are expecting the compiler *
55* to generate the correct sequence of multiply-add fused instructions. *
56* *
57* These routines are not intended for 32-bit Intel architectures. *
58* *
59* A version of gcc higher than 932 is required. *
60* *
61* GCC compiler options: *
62* optimization level 3 (-O3) *
63* -fschedule-insns -finline-functions -funroll-all-loops *
64* *
65*******************************************************************************/
66
67#include "math.h"
68#include "fp_private.h"
69
70static const double two54 = 0x1.0p+54; // 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
71static const double two25 = 0x1.0p+25; // 33554432.0e0;
72
73/******************************************************************************
74* Function ldexp *
75* Multiplies a floating-point number by an integer power of 2. *
76* *
77* Functions used in this routine: *
78* scalb. *
79******************************************************************************/
80
81double ldexp ( double value, int exp )
82{
83 if (unlikely( exp > INT16_MAX ))
84 exp = INT16_MAX;
85 else if (unlikely( exp < -INT16_MAX ))
86 exp = -INT16_MAX;
87 return scalbn ( value, exp );
88}
89
90float ldexpf ( float value, int exp )
91{
92 if (unlikely( exp > INT16_MAX ))
93 exp = INT16_MAX;
94 else if (unlikely( exp < -INT16_MAX ))
95 exp = -INT16_MAX;
96 return scalbnf ( value, exp );
97}
98
99/******************************************************************************
100* Function frexp *
101* Breaks a floating-point number into a normalized fraction and an *
102* integral power of 2. It stores the integer in the object pointed by *
103* *eptr. *
104******************************************************************************/
105
106double frexp ( double value, int *eptr )
107{
108 hexdouble argument;
109 uint32_t valueHead;
110
111 argument.d = value;
112 __NOOP;
113 __NOOP;
114 __NOOP;
115 valueHead = argument.i.hi & 0x7fffffff; // valueHead <- |x|
116
117 *eptr = 0;
118 if (unlikely( valueHead >= 0x7ff00000 || ( valueHead | argument.i.lo ) == 0 ))
119 return value; // 0, inf, or NaN
120
121 if (unlikely( valueHead < 0x00100000 ))
122 { // denorm
123 argument.d = two54 * value;
124 __NOOP;
125 __NOOP;
126 __NOOP;
127 valueHead = argument.i.hi & 0x7fffffff;
128 *eptr = -54;
129 }
130 *eptr += ( valueHead >> 20 ) - 1022;
131 argument.i.hi = ( argument.i.hi & 0x800fffff ) | 0x3fe00000;
132 __NOOP;
133 __NOOP;
134 __NOOP;
135 return argument.d;
136}
137
138float frexpf ( float value, int *eptr )
139{
140 hexsingle argument;
141 uint32_t valueHead;
142
143 argument.fval = value;
144 __NOOP;
145 __NOOP;
146 __NOOP;
147 valueHead = argument.lval & 0x7fffffff; // valueHead <- |x|
148
149 *eptr = 0;
150 if (unlikely( valueHead >= 0x7f800000 || valueHead == 0 ))
151 return value; // 0, inf, or NaN
152
153 if (unlikely( valueHead < 0x00800000 ))
154 { // denorm
155 argument.fval = (float) two25 * value;
156 __NOOP;
157 __NOOP;
158 __NOOP;
159 valueHead = argument.lval & 0x7fffffff;
160 *eptr = -25;
161 }
162 *eptr += ( valueHead >> 23 ) - 126;
163 argument.lval = ( argument.lval & 0x807fffff ) | 0x3f000000;
164 __NOOP;
165 __NOOP;
166 __NOOP;
167 return argument.fval;
168}