blob: 530e07c79f051aab246f39b31bfa1571699a8f57 [file] [log] [blame]
Szabolcs Nagy2eaed462013-09-04 15:52:54 +00001#include <math.h>
2#include <stdint.h>
Rich Felkerb69f6952012-03-13 01:17:53 -04003
Rich Felkerb69f6952012-03-13 01:17:53 -04004double scalbn(double x, int n)
5{
Szabolcs Nagy2eaed462013-09-04 15:52:54 +00006 union {double f; uint64_t i;} u;
7 double_t y = x;
Rich Felkerb69f6952012-03-13 01:17:53 -04008
nsz8051e082012-03-19 10:54:07 +01009 if (n > 1023) {
Szabolcs Nagy2eaed462013-09-04 15:52:54 +000010 y *= 0x1p1023;
nsz8051e082012-03-19 10:54:07 +010011 n -= 1023;
12 if (n > 1023) {
Szabolcs Nagy2eaed462013-09-04 15:52:54 +000013 y *= 0x1p1023;
nsz8051e082012-03-19 10:54:07 +010014 n -= 1023;
Szabolcs Nagy1b77b902013-08-15 10:07:46 +000015 if (n > 1023)
16 n = 1023;
nsz8051e082012-03-19 10:54:07 +010017 }
18 } else if (n < -1022) {
Szabolcs Nagy2eaed462013-09-04 15:52:54 +000019 y *= 0x1p-1022;
nsz8051e082012-03-19 10:54:07 +010020 n += 1022;
21 if (n < -1022) {
Szabolcs Nagy2eaed462013-09-04 15:52:54 +000022 y *= 0x1p-1022;
nsz8051e082012-03-19 10:54:07 +010023 n += 1022;
Szabolcs Nagy1b77b902013-08-15 10:07:46 +000024 if (n < -1022)
25 n = -1022;
nsz8051e082012-03-19 10:54:07 +010026 }
Rich Felkerb69f6952012-03-13 01:17:53 -040027 }
Szabolcs Nagy2eaed462013-09-04 15:52:54 +000028 u.i = (uint64_t)(0x3ff+n)<<52;
29 x = y * u.f;
Szabolcs Nagyc4359e02012-11-13 10:55:35 +010030 return x;
Rich Felkerb69f6952012-03-13 01:17:53 -040031}