aboutsummaryrefslogtreecommitdiff
path: root/src/bigfloat.cpp
blob: c1bfcbbb97d8be4e26b3dedfdbabfe9f335dae6b (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
/*
 * Copyright (c) 2017 Andrew Kelley
 *
 * This file is part of zig, which is MIT licensed.
 * See http://opensource.org/licenses/MIT
 */

#include "bigfloat.hpp"
#include "bigint.hpp"
#include "buffer.hpp"
#include <math.h>
#include <errno.h>

extern "C" {
    __float128 fmodq(__float128 a, __float128 b);
    __float128 ceilq(__float128 a);
    __float128 floorq(__float128 a);
    __float128 strtoflt128 (const char *s, char **sp);
    int quadmath_snprintf (char *s, size_t size, const char *format, ...);
}


void bigfloat_init_float(BigFloat *dest, __float128 x) {
    dest->value = x;
}

void bigfloat_init_bigfloat(BigFloat *dest, const BigFloat *x) {
    dest->value = x->value;
}

void bigfloat_init_bigint(BigFloat *dest, const BigInt *op) {
    dest->value = 0.0;
    if (op->digit_count == 0)
        return;

    __float128 base = (__float128)UINT64_MAX;
    const uint64_t *digits = bigint_ptr(op);

    for (size_t i = op->digit_count - 1;;) {
        uint64_t digit = digits[i];
        dest->value *= base;
        dest->value += (__float128)digit;

        if (i == 0) {
            if (op->is_negative) {
                dest->value = -dest->value;
            }
            return;
        }
        i -= 1;
    }
}

int bigfloat_init_buf_base10(BigFloat *dest, const uint8_t *buf_ptr, size_t buf_len) {
    char *str_begin = (char *)buf_ptr;
    char *str_end;
    errno = 0;
    dest->value = strtoflt128(str_begin, &str_end);
    if (errno) {
        return ErrorOverflow;
    }
    assert(str_end <= ((char*)buf_ptr) + buf_len);
    return 0;
}

void bigfloat_add(BigFloat *dest, const BigFloat *op1, const BigFloat *op2) {
    dest->value = op1->value + op2->value;
}

void bigfloat_negate(BigFloat *dest, const BigFloat *op) {
    dest->value = -op->value;
}

void bigfloat_sub(BigFloat *dest, const BigFloat *op1, const BigFloat *op2) {
    dest->value = op1->value - op2->value;
}

void bigfloat_mul(BigFloat *dest, const BigFloat *op1, const BigFloat *op2) {
    dest->value = op1->value * op2->value;
}

void bigfloat_div(BigFloat *dest, const BigFloat *op1, const BigFloat *op2) {
    dest->value = op1->value / op2->value;
}

void bigfloat_div_trunc(BigFloat *dest, const BigFloat *op1, const BigFloat *op2) {
    dest->value = op1->value / op2->value;
    if (dest->value >= 0.0) {
        dest->value = floorq(dest->value);
    } else {
        dest->value = ceilq(dest->value);
    }
}

void bigfloat_div_floor(BigFloat *dest, const BigFloat *op1, const BigFloat *op2) {
    dest->value = floorq(op1->value / op2->value);
}

void bigfloat_rem(BigFloat *dest, const BigFloat *op1, const BigFloat *op2) {
    dest->value = fmodq(op1->value, op2->value);
}

void bigfloat_mod(BigFloat *dest, const BigFloat *op1, const BigFloat *op2) {
    dest->value = fmodq(fmodq(op1->value, op2->value) + op2->value, op2->value);
}

void bigfloat_write_buf(Buf *buf, const BigFloat *op) {
    buf_resize(buf, 256);
    int len = quadmath_snprintf(buf_ptr(buf), buf_len(buf), "%Qf", op->value);
    assert(len > 0);
    buf_resize(buf, len);
}

Cmp bigfloat_cmp(const BigFloat *op1, const BigFloat *op2) {
    if (op1->value > op2->value) {
        return CmpGT;
    } else if (op1->value < op2->value) {
        return CmpLT;
    } else {
        return CmpEQ;
    }
}

// TODO this is wrong when compiler running on big endian systems. caught by tests
void bigfloat_write_ieee597(const BigFloat *op, uint8_t *buf, size_t bit_count, bool is_big_endian) {
    if (bit_count == 32) {
        float f32 = op->value;
        memcpy(buf, &f32, 4);
    } else if (bit_count == 64) {
        double f64 = op->value;
        memcpy(buf, &f64, 8);
    } else if (bit_count == 128) {
        __float128 f128 = op->value;
        memcpy(buf, &f128, 16);
    } else {
        zig_unreachable();
    }
}

// TODO this is wrong when compiler running on big endian systems. caught by tests
void bigfloat_read_ieee597(BigFloat *dest, const uint8_t *buf, size_t bit_count, bool is_big_endian) {
    if (bit_count == 32) {
        float f32;
        memcpy(&f32, buf, 4);
        dest->value = f32;
    } else if (bit_count == 64) {
        double f64;
        memcpy(&f64, buf, 8);
        dest->value = f64;
    } else if (bit_count == 128) {
        __float128 f128;
        memcpy(&f128, buf, 16);
        dest->value = f128;
    } else {
        zig_unreachable();
    }
}

double bigfloat_to_double(const BigFloat *bigfloat) {
    return bigfloat->value;
}

Cmp bigfloat_cmp_zero(const BigFloat *bigfloat) {
    if (bigfloat->value < 0.0) {
        return CmpLT;
    } else if (bigfloat->value > 0.0) {
        return CmpGT;
    } else {
        return CmpEQ;
    }
}

bool bigfloat_has_fraction(const BigFloat *bigfloat) {
    return floorq(bigfloat->value) != bigfloat->value;
}