source: trunk/libs/libmath/fmod.c @ 625

Last change on this file since 625 was 469, checked in by alain, 6 years ago

1) Introduce the libsemaphore library.
2) Introduce a small libmath library, required by the "fft" application..
3) Introduce the multithreaded "fft" application.
4) Fix a bad synchronisation bug in the Copy-On-Write mechanism.

File size: 3.0 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/*
17 * Modified for ALMOS-MKH OS at UPMC, France, August 2018. (Alain Greiner)
18 */
19
20
21#include "math.h"
22#include "math_private.h"
23 
24double fmod(double x, double y)
25{
26    union {
27        double d;
28        uint64_t u;
29    } ux = {x};
30    union {
31        double d;
32        uint64_t u;
33    } uy = {y};
34   
35    uint64_t absx = ux.u & ~0x8000000000000000ULL;
36    uint64_t absy = uy.u & ~0x8000000000000000ULL;
37
38    if (absx - 1ULL >= 0x7fefffffffffffffULL || absy - 1ULL >= 0x7fefffffffffffffULL) {
39        double fabsx = fabs(x);
40        double fabsy = fabs(y);
41       
42        // deal with NaN
43        if (x != x || y != y) {
44            return x + y;
45        }
46
47        // x = Inf or y = 0, return Invalid per IEEE-754
48        if (!isfinite(fabsx) || 0.0 == y) {
49            return 0.0 / 0.0; // NaN
50        }
51
52        // handle trivial case
53        if (!isfinite(fabsy) || 0.0 == x) {
54            return x;
55        }
56    }
57 
58    if (absy >= absx) {
59        if (absy == absx) {
60            ux.u ^= absx;
61            return ux.d;
62        }
63       
64        return x;
65    }
66 
67    int32_t expx = absx >> 52;
68    int32_t expy = absy >> 52;
69    int64_t sx = absx & 0x000fffffffffffffULL;
70    int64_t sy = absy & 0x000fffffffffffffULL;
71
72    if (0 == expx) {
73        uint32_t shift = __builtin_clzll(absx) - (64 - 53);
74        sx <<= shift;
75        expx = 1 - shift;
76    }
77
78    if (0 == expy)
79    {
80        uint32_t shift = __builtin_clzll(absy) - (64 - 53);
81        sy <<= shift;
82        expy = 1 - shift;
83    }
84    sx |= 0x0010000000000000ULL;
85    sy |= 0x0010000000000000ULL;
86
87
88    int32_t idiff = expx - expy;
89    int32_t shift = 0;
90    int64_t mask;
91   
92    do {
93        sx <<= shift;
94        idiff += ~shift;
95        sx -= sy;
96        mask = sx >> 63;
97        sx += sx;
98        sx += sy & mask;
99        shift = __builtin_clzll(sx) - (64 - 53);
100    }
101    while (idiff >= shift && sx != 0LL);
102
103    if (idiff < 0) {
104        sx += sy & mask;
105        sx >>= 1;
106        idiff = 0;
107    }
108   
109    sx <<= idiff;
110   
111    if (0 == sx) {
112        ux.u &= 0x8000000000000000ULL;
113        return ux.d;
114    }
115   
116    shift = __builtin_clzll(sx) - (64 - 53);
117    sx <<= shift;
118    expy -= shift;
119    sx &= 0x000fffffffffffffULL;
120    sx |= ux.u & 0x8000000000000000ULL;
121    if (expy > 0) {
122        ux.u = sx | ((int64_t) expy << 52);
123        return ux.d;
124    }
125   
126    expy += 1022;
127    ux.u = sx | ((int64_t) expy << 52);
128    return ux.d * 0x1.0p-1022;
129
130}
131
132
Note: See TracBrowser for help on using the repository browser.