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 hypot.c, *
25* Function hypot(x,y), *
26* Implementation of sqrt(x^2+y^2) for the PowerPC. *
27* *
28* Copyright � 1991 Apple Computer, Inc. All rights reserved. *
29* *
30* Written by Ali Sazegari, started on November 1991, *
31* *
32* based on math.h, library code for Macintoshes with a 68881/68882 *
33* by Jim Thomas. *
34* *
35* W A R N I N G: This routine expects a 64 bit double model. *
36* *
37* December 03 1992: first rs6000 port. *
38* January 05 1993: added the environmental controls. *
39* July 14 1993: added #pragma fenv_access. changed feholdexcept *
40* to the internal routine _feprocentry. *
41* September22 1993: conforming to nceg-fpce specification for �� values. *
42* September19 1994: revamp of the algorithm for performance. *
43* *
44*******************************************************************************/
45
46#include "math.h"
47#include "fenv.h"
48#include "fp_private.h"
49#include "fenv_private.h"
50#include <System/ppc/cpu_capabilities.h>
51
52inline static uint32_t __attribute__((always_inline))
53_get_cpu_capabilities(void)
54{
55 uint32_t caps;
56 asm("lwz %0, %1(0)" : "=r" (caps) : "I" (_COMM_PAGE_CPU_CAPABILITIES));
57 return caps;
58}
59
60/*******************************************************************************
61* Functions needed for the computation. *
62*******************************************************************************/
63
64/* the following fp.h functions are used: */
65/* fabs and sqrt; */
66/* the following environmental functions are used: */
67/* feclearexcept, _feprocentry and feupdateenv. */
68
69
70#pragma STDC FENV_ACCESS ON
71
72static const hexdouble Huge = HEXDOUBLE(0x7FF00000, 0x00000000);
73static const hexdouble NegHuge = HEXDOUBLE(0xFFF00000, 0x00000000);
74
75double hypot ( double x, double y )
76{
77 register double temp;
78 hexdouble OldEnvironment, CurrentEnvironment;
79
80 register double FPR_env, FPR_z, FPR_one, FPR_inf, FPR_Minf,
81 FPR_absx, FPR_absy, FPR_big, FPR_small;
82
83 FPR_z = 0.0; FPR_one = 1.0;
84 FPR_inf = Huge.d; FPR_Minf = NegHuge.d;
85 FPR_absx = __FABS( x ); FPR_absy = __FABS( y );
86
87/*******************************************************************************
88* If argument is SNaN then a QNaN has to be returned and the invalid *
89* flag signaled. *
90*******************************************************************************/
91
92 if (unlikely( ( x == FPR_inf ) || ( y == FPR_inf ) || ( x == FPR_Minf ) || ( y == FPR_Minf ) ))
93 return FPR_inf;
94
95 if (unlikely( ( x != x ) || ( y != y ) ))
96 {
97 x = __FABS ( x + y );
98 return x;
99 }
100
101 if ( FPR_absx > FPR_absy )
102 {
103 FPR_big = FPR_absx;
104 FPR_small = FPR_absy;
105 }
106 else
107 {
108 FPR_big = FPR_absy;
109 FPR_small = FPR_absx;
110 }
111
112 // Now +0.0 <= FPR_small <= FPR_big < INFINITY
113
114 if (unlikely( FPR_small == FPR_z ))
115 return FPR_big;
116
117 FEGETENVD( FPR_env ); // save environment, set default
118 FESETENVD( FPR_z );
119
120 temp = FPR_small / FPR_big;
121 OldEnvironment.d = FPR_env;
122
123 temp = FPR_one + temp * temp;
124 if (likely((_get_cpu_capabilities() & kHasFsqrt)))
125 temp = __fsqrt ( temp );
126 else
127 temp = sqrt ( temp );
128
129 FEGETENVD_GRP( CurrentEnvironment.d );
130 CurrentEnvironment.i.lo &= ~FE_UNDERFLOW; // Clear any inconsequential underflow
131 FESETENVD_GRP( CurrentEnvironment.d );
132
133 temp = FPR_big * temp; // Might raise UNDERFLOW or OVERFLOW
134
135 FEGETENVD_GRP( CurrentEnvironment.d );
136 OldEnvironment.i.lo |= CurrentEnvironment.i.lo; // Pick up any UF or OF
137 FESETENVD_GRP( OldEnvironment.d ); // restore caller's environment
138
139 return temp;
140}