source: soft/giet_vm/giet_libs/math/fmod.c @ 822

Last change on this file since 822 was 809, checked in by alain, 9 years ago

Avoid GCC warnings.

File size: 2.9 KB
Line 
1/*
2 *  fmod.c
3 *  cLibm
4 *
5 *     Written by Jon Okada, started on December 7th, 1992.                     
6 *     Modified by Paul Finlayson (PAF) for MathLib v2.                         
7 *     Modified by A. Sazegari (ali) for MathLib v3.                             
8 *     Modified and ported by Robert A. Murley (ram) for Mac OS X.               
9 *     Modified for cLibm, fixed a few edge cases, rewrote local_ funcs by Ian Ollmann.
10 *     Modified for armLibm, removed code to set flags as this is a soft-float fallback only.
11 *
12 *  Copyright 2007 Apple Inc. All rights reserved.
13 *
14 */
15
16#include "../math.h"
17#include "math_private.h"
18 
19double fmod(double x, double y)
20{
21    union {
22        double d;
23        uint64_t u;
24    } ux = {x};
25    union {
26        double d;
27        uint64_t u;
28    } uy = {y};
29   
30    uint64_t absx = ux.u & ~0x8000000000000000ULL;
31    uint64_t absy = uy.u & ~0x8000000000000000ULL;
32
33    if (absx - 1ULL >= 0x7fefffffffffffffULL || absy - 1ULL >= 0x7fefffffffffffffULL) {
34        double fabsx = fabs(x);
35        double fabsy = fabs(y);
36       
37        // deal with NaN
38        if (x != x || y != y) {
39            return x + y;
40        }
41
42        // x = Inf or y = 0, return Invalid per IEEE-754
43        if (!isfinite(fabsx) || 0.0 == y) {
44            return 0.0 / 0.0; // NaN
45        }
46
47        // handle trivial case
48        if (!isfinite(fabsy) || 0.0 == x) {
49            return x;
50        }
51    }
52 
53    if (absy >= absx) {
54        if (absy == absx) {
55            ux.u ^= absx;
56            return ux.d;
57        }
58       
59        return x;
60    }
61 
62    int32_t expx = absx >> 52;
63    int32_t expy = absy >> 52;
64    int64_t sx = absx & 0x000fffffffffffffULL;
65    int64_t sy = absy & 0x000fffffffffffffULL;
66
67    if (0 == expx) {
68        uint32_t shift = __builtin_clzll(absx) - (64 - 53);
69        sx <<= shift;
70        expx = 1 - shift;
71    }
72
73    if (0 == expy)
74    {
75        uint32_t shift = __builtin_clzll(absy) - (64 - 53);
76        sy <<= shift;
77        expy = 1 - shift;
78    }
79    sx |= 0x0010000000000000ULL;
80    sy |= 0x0010000000000000ULL;
81
82
83    int32_t idiff = expx - expy;
84    int32_t shift = 0;
85    int64_t mask;
86   
87    do {
88        sx <<= shift;
89        idiff += ~shift;
90        sx -= sy;
91        mask = sx >> 63;
92        sx += sx;
93        sx += sy & mask;
94        shift = __builtin_clzll(sx) - (64 - 53);
95    }
96    while (idiff >= shift && sx != 0LL);
97
98    if (idiff < 0) {
99        sx += sy & mask;
100        sx >>= 1;
101        idiff = 0;
102    }
103   
104    sx <<= idiff;
105   
106    if (0 == sx) {
107        ux.u &= 0x8000000000000000ULL;
108        return ux.d;
109    }
110   
111    shift = __builtin_clzll(sx) - (64 - 53);
112    sx <<= shift;
113    expy -= shift;
114    sx &= 0x000fffffffffffffULL;
115    sx |= ux.u & 0x8000000000000000ULL;
116    if (expy > 0) {
117        ux.u = sx | ((int64_t) expy << 52);
118        return ux.d;
119    }
120   
121    expy += 1022;
122    ux.u = sx | ((int64_t) expy << 52);
123    return ux.d * 0x1.0p-1022;
124
125}
126
127
Note: See TracBrowser for help on using the repository browser.