objcomplex.c 7.27 KB
Newer Older
1
2
3
#include <stdlib.h>
#include <assert.h>

4
#include "mpconfig.h"
5
6
#include "nlr.h"
#include "misc.h"
7
#include "qstr.h"
8
#include "obj.h"
9
#include "parsenum.h"
10
11
12
13
#include "runtime0.h"

#if MICROPY_ENABLE_FLOAT

14
15
#include <math.h>

16
17
18
19
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
#include "formatfloat.h"
#endif

20
21
22
23
24
25
26
27
typedef struct _mp_obj_complex_t {
    mp_obj_base_t base;
    mp_float_t real;
    mp_float_t imag;
} mp_obj_complex_t;

mp_obj_t mp_obj_new_complex(mp_float_t real, mp_float_t imag);

28
STATIC void complex_print(void (*print)(void *env, const char *fmt, ...), void *env, mp_obj_t o_in, mp_print_kind_t kind) {
29
    mp_obj_complex_t *o = o_in;
30
31
32
33
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
    char buf[32];
    if (o->real == 0) {
        format_float(o->imag, buf, sizeof(buf), 'g', 6, '\0');
34
        print(env, "%sj", buf);
35
36
37
    } else {
        format_float(o->real, buf, sizeof(buf), 'g', 6, '\0');
        print(env, "(%s+", buf);
38
        format_float(o->imag, buf, sizeof(buf), 'g', 6, '\0');
39
40
41
        print(env, "%sj)", buf);
    }
#else
42
    if (o->real == 0) {
mux's avatar
mux committed
43
        print(env, "%.8gj",  (double) o->imag);
44
    } else {
mux's avatar
mux committed
45
        print(env, "(%.8g+%.8gj)", (double) o->real, (double) o->imag);
46
    }
47
#endif
48
49
}

50
STATIC mp_obj_t complex_make_new(mp_obj_t type_in, uint n_args, uint n_kw, const mp_obj_t *args) {
51
52
    // TODO check n_kw == 0

53
54
55
56
57
    switch (n_args) {
        case 0:
            return mp_obj_new_complex(0, 0);

        case 1:
58
59
60
61
62
63
64
            if (MP_OBJ_IS_STR(args[0])) {
                // a string, parse it
                uint l;
                const char *s = mp_obj_str_get_data(args[0], &l);
                return mp_parse_num_decimal(s, l, true, true);
            } else if (MP_OBJ_IS_TYPE(args[0], &mp_type_complex)) {
                // a complex, just return it
65
66
                return args[0];
            } else {
67
                // something else, try to cast it to a complex
68
69
70
                return mp_obj_new_complex(mp_obj_get_float(args[0]), 0);
            }

71
        case 2: {
72
            mp_float_t real, imag;
73
            if (MP_OBJ_IS_TYPE(args[0], &mp_type_complex)) {
74
                mp_obj_complex_get(args[0], &real, &imag);
75
            } else {
76
                real = mp_obj_get_float(args[0]);
77
78
                imag = 0;
            }
79
            if (MP_OBJ_IS_TYPE(args[1], &mp_type_complex)) {
80
                mp_float_t real2, imag2;
81
                mp_obj_complex_get(args[1], &real2, &imag2);
82
83
84
                real -= imag2;
                imag += real2;
            } else {
85
                imag += mp_obj_get_float(args[1]);
86
87
88
89
90
            }
            return mp_obj_new_complex(real, imag);
        }

        default:
91
            nlr_raise(mp_obj_new_exception_msg_varg(&mp_type_TypeError, "complex takes at most 2 arguments, %d given", n_args));
92
93
94
    }
}

95
STATIC mp_obj_t complex_unary_op(int op, mp_obj_t o_in) {
96
97
    mp_obj_complex_t *o = o_in;
    switch (op) {
Damien George's avatar
Damien George committed
98
99
100
        case MP_UNARY_OP_BOOL: return MP_BOOL(o->real != 0 || o->imag != 0);
        case MP_UNARY_OP_POSITIVE: return o_in;
        case MP_UNARY_OP_NEGATIVE: return mp_obj_new_complex(-o->real, -o->imag);
101
        default: return MP_OBJ_NOT_SUPPORTED;
102
103
104
    }
}

105
STATIC mp_obj_t complex_binary_op(int op, mp_obj_t lhs_in, mp_obj_t rhs_in) {
106
107
108
109
    mp_obj_complex_t *lhs = lhs_in;
    return mp_obj_complex_binary_op(op, lhs->real, lhs->imag, rhs_in);
}

110
const mp_obj_type_t mp_type_complex = {
111
    { &mp_type_type },
112
    .name = MP_QSTR_complex,
113
114
115
116
    .print = complex_print,
    .make_new = complex_make_new,
    .unary_op = complex_unary_op,
    .binary_op = complex_binary_op,
117
118
119
120
};

mp_obj_t mp_obj_new_complex(mp_float_t real, mp_float_t imag) {
    mp_obj_complex_t *o = m_new_obj(mp_obj_complex_t);
121
    o->base.type = &mp_type_complex;
122
123
124
125
126
127
    o->real = real;
    o->imag = imag;
    return o;
}

void mp_obj_complex_get(mp_obj_t self_in, mp_float_t *real, mp_float_t *imag) {
128
    assert(MP_OBJ_IS_TYPE(self_in, &mp_type_complex));
129
130
131
132
133
134
135
136
    mp_obj_complex_t *self = self_in;
    *real = self->real;
    *imag = self->imag;
}

mp_obj_t mp_obj_complex_binary_op(int op, mp_float_t lhs_real, mp_float_t lhs_imag, mp_obj_t rhs_in) {
    mp_float_t rhs_real, rhs_imag;
    mp_obj_get_complex(rhs_in, &rhs_real, &rhs_imag); // can be any type, this function will convert to float (if possible)
137
    switch (op) {
Damien George's avatar
Damien George committed
138
139
        case MP_BINARY_OP_ADD:
        case MP_BINARY_OP_INPLACE_ADD:
140
141
142
            lhs_real += rhs_real;
            lhs_imag += rhs_imag;
            break;
Damien George's avatar
Damien George committed
143
144
        case MP_BINARY_OP_SUBTRACT:
        case MP_BINARY_OP_INPLACE_SUBTRACT:
145
146
147
            lhs_real -= rhs_real;
            lhs_imag -= rhs_imag;
            break;
Damien George's avatar
Damien George committed
148
        case MP_BINARY_OP_MULTIPLY:
149
150
151
152
        case MP_BINARY_OP_INPLACE_MULTIPLY: {
            mp_float_t real;
            multiply:
            real = lhs_real * rhs_real - lhs_imag * rhs_imag;
153
154
155
156
            lhs_imag = lhs_real * rhs_imag + lhs_imag * rhs_real;
            lhs_real = real;
            break;
        }
Damien George's avatar
Damien George committed
157
        case MP_BINARY_OP_FLOOR_DIVIDE:
158
159
160
        case MP_BINARY_OP_INPLACE_FLOOR_DIVIDE:
            nlr_raise(mp_obj_new_exception_msg_varg(&mp_type_TypeError, "can't do truncated division of a complex number"));

Damien George's avatar
Damien George committed
161
        case MP_BINARY_OP_TRUE_DIVIDE:
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
        case MP_BINARY_OP_INPLACE_TRUE_DIVIDE:
            if (rhs_imag == 0) {
                if (rhs_real == 0) {
                    nlr_raise(mp_obj_new_exception_msg(&mp_type_ZeroDivisionError, "complex division by zero"));
                }
                lhs_real /= rhs_real;
                lhs_imag /= rhs_real;
            } else if (rhs_real == 0) {
                mp_float_t real = lhs_imag / rhs_imag;
                lhs_imag = -lhs_real / rhs_imag;
                lhs_real = real;
            } else {
                mp_float_t rhs_len_sq = rhs_real*rhs_real + rhs_imag*rhs_imag;
                rhs_real /= rhs_len_sq;
                rhs_imag /= -rhs_len_sq;
                goto multiply;
            }
            break;

181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
        case MP_BINARY_OP_POWER:
        case MP_BINARY_OP_INPLACE_POWER: {
            // z1**z2 = exp(z2*ln(z1))
            //        = exp(z2*(ln(|z1|)+i*arg(z1)))
            //        = exp( (x2*ln1 - y2*arg1) + i*(y2*ln1 + x2*arg1) )
            //        = exp(x3 + i*y3)
            //        = exp(x3)*(cos(y3) + i*sin(y3))
            mp_float_t abs1 = MICROPY_FLOAT_C_FUN(sqrt)(lhs_real*lhs_real + lhs_imag*lhs_imag);
            if (abs1 == 0) {
                if (rhs_imag == 0) {
                    lhs_real = 1;
                    rhs_real = 0;
                } else {
                    nlr_raise(mp_obj_new_exception_msg(&mp_type_ZeroDivisionError, "0.0 to a complex power"));
                }
            } else {
                mp_float_t ln1 = MICROPY_FLOAT_C_FUN(log)(abs1);
                mp_float_t arg1 = MICROPY_FLOAT_C_FUN(atan2)(lhs_imag, lhs_real);
                mp_float_t x3 = rhs_real * ln1 - rhs_imag * arg1;
                mp_float_t y3 = rhs_imag * ln1 + rhs_real * arg1;
                mp_float_t exp_x3 = MICROPY_FLOAT_C_FUN(exp)(x3);
                lhs_real = exp_x3 * MICROPY_FLOAT_C_FUN(cos)(y3);
                lhs_imag = exp_x3 * MICROPY_FLOAT_C_FUN(sin)(y3);
            }
            break;
        }

208
209
        case MP_BINARY_OP_EQUAL: return MP_BOOL(lhs_real == rhs_real && lhs_imag == rhs_imag);

210
        default:
211
            return MP_OBJ_NOT_SUPPORTED;
212
213
214
215
216
    }
    return mp_obj_new_complex(lhs_real, lhs_imag);
}

#endif