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
scalbn(double x,int n)19 double scalbn(double x, int n)
20 {
21 union {double f; uint64_t i;} u;
22 long double x_tmp = x;
23 double_t y = x_tmp;
24 int z = n;
25
26 if (z > 1023) { // 1023:byte alignment
27 y *= 0x1p1023;
28 z -= 1023; // 1023:byte alignment
29 if (z > 1023) { // 1023:byte alignment
30 y *= 0x1p1023;
31 z -= 1023; // 1023:byte alignment
32 if (z > 1023) { // 1023:byte alignment
33 z = 1023; // 1023:byte alignment
34 }
35 }
36 } else if (z < -1022) { // -1022:byte alignment
37 /* make sure final z < -53 to avoid double
38 rounding in the subnormal range */
39 y *= 0x1p-1022 * 0x1p53;
40 z += 1022 - 53; // 1022:byte alignment, 53:byte alignment
41 if (z < -1022) { // -1022:byte alignment
42 y *= 0x1p-1022 * 0x1p53;
43 z += 1022 - 53; // 1022:byte alignment, 53:byte alignment
44 if (z < -1022) { // -1022:byte alignment
45 z = -1022; // -1022:byte alignment
46 }
47 }
48 }
49 u.i = (uint64_t)(0x3ff + z) << 52; // 52:byte alignment
50 x_tmp = y * u.f;
51 return x_tmp;
52 }