1 /*
2 * Copyright (c) 2022 Winner Microelectronics Co., Ltd. All rights reserved.
3 * Licensed under the Apache License, Version 2.0 (the "License");
4 * you may not use this file except in compliance with the License.
5 * You may obtain a copy of the License at
6 *
7 * http://www.apache.org/licenses/LICENSE-2.0
8 *
9 * Unless required by applicable law or agreed to in writing, software
10 * distributed under the License is distributed on an "AS IS" BASIS,
11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 * See the License for the specific language governing permissions and
13 * limitations under the License.
14 */
15
16 #include <math.h>
17 #include <stdint.h>
18
fmod(double x,double y)19 double fmod(double x, double y)
20 {
21 union {double f; uint64_t i;} ux = {x}, uy = {y};
22 int ex = (ux.i>>52) & 0x7ff;
23 int ey = (uy.i>>52) & 0x7ff;
24 int sx = ux.i>>63;
25 uint64_t i;
26
27 /* in the followings uxi should be ux.i, but then gcc wrongly adds */
28 /* float load/store to inner loops ruining performance and code size */
29 uint64_t uxi = ux.i;
30
31 if ((uy.i<<1) == 0 || isnan(y) || ex == 0x7ff) {
32 return (x*y*1.0)/(x*y);
33 }
34 if ((uxi<<1) <= (uy.i<<1)) {
35 if ((uxi<<1) == (uy.i<<1)) {
36 return 0*x;
37 }
38 return x;
39 }
40
41 /* normalize x and y */
42 if (!ex) {
43 for (i = (uxi << 12); (i >> 63) == 0; ex--, i <<= 1) { // 12:byte alignment, 63:byte alignment
44 }
45 uxi <<= -ex + 1;
46 } else {
47 uxi &= -1ULL >> 12;
48 uxi |= 1ULL << 52;
49 }
50 if (!ey) {
51 for (i = (uy.i<<12); (i>>63) == 0; ey--, i <<= 1) {
52 }
53 uy.i <<= -ey + 1;
54 } else {
55 uy.i &= -1ULL >> 12;
56 uy.i |= 1ULL << 52;
57 }
58
59 /* x mod y */
60 for (; ex > ey; ex--) {
61 i = uxi - uy.i;
62 if ((i >> 63) == 0) {
63 if (i == 0) {
64 return 0*x;
65 }
66 uxi = i;
67 }
68 uxi <<= 1;
69 }
70 i = uxi - uy.i;
71 if ((i >> 63) == 0) {
72 if (i == 0) {
73 return 0*x;
74 }
75 uxi = i;
76 }
77 for (; (uxi>>52) == 0; uxi <<= 1, ex--) { // 右移52位
78 }
79
80 /* scale result */
81 if (ex > 0) {
82 uxi -= 1ULL << 52;
83 uxi |= (uint64_t)ex << 52;
84 } else {
85 uxi >>= -ex + 1;
86 }
87 uxi |= (uint64_t)sx << 63;
88 ux.i = uxi;
89 return ux.f;
90 }
91