Branch data Line data Source code
1 : : /*
2 : : * This file is part of the MicroPython project, http://micropython.org/
3 : : *
4 : : * The MIT License (MIT)
5 : : *
6 : : * Copyright (c) 2013, 2014 Damien P. George
7 : : *
8 : : * Permission is hereby granted, free of charge, to any person obtaining a copy
9 : : * of this software and associated documentation files (the "Software"), to deal
10 : : * in the Software without restriction, including without limitation the rights
11 : : * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 : : * copies of the Software, and to permit persons to whom the Software is
13 : : * furnished to do so, subject to the following conditions:
14 : : *
15 : : * The above copyright notice and this permission notice shall be included in
16 : : * all copies or substantial portions of the Software.
17 : : *
18 : : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 : : * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 : : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 : : * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 : : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 : : * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24 : : * THE SOFTWARE.
25 : : */
26 : :
27 : : #include <stdlib.h>
28 : : #include <stdio.h>
29 : : #include <assert.h>
30 : :
31 : : #include "py/parsenum.h"
32 : : #include "py/runtime.h"
33 : :
34 : : #if MICROPY_PY_BUILTINS_COMPLEX
35 : :
36 : : #include <math.h>
37 : : #include "py/formatfloat.h"
38 : :
39 : : typedef struct _mp_obj_complex_t {
40 : : mp_obj_base_t base;
41 : : mp_float_t real;
42 : : mp_float_t imag;
43 : : } mp_obj_complex_t;
44 : :
45 : 344 : static void complex_print(const mp_print_t *print, mp_obj_t o_in, mp_print_kind_t kind) {
46 : 344 : (void)kind;
47 : 344 : mp_obj_complex_t *o = MP_OBJ_TO_PTR(o_in);
48 : : #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
49 : : char buf[16];
50 : : #if MICROPY_OBJ_REPR == MICROPY_OBJ_REPR_C
51 : : const int precision = 6;
52 : : #else
53 : : const int precision = 7;
54 : : #endif
55 : : #else
56 : 344 : char buf[32];
57 : 344 : const int precision = 16;
58 : : #endif
59 [ + + ]: 344 : if (o->real == 0) {
60 : 72 : mp_format_float(o->imag, buf, sizeof(buf), 'g', precision, '\0');
61 : 72 : mp_printf(print, "%sj", buf);
62 : : } else {
63 : 272 : mp_format_float(o->real, buf, sizeof(buf), 'g', precision, '\0');
64 : 272 : mp_printf(print, "(%s", buf);
65 [ + + + + ]: 272 : if (o->imag >= 0 || isnan(o->imag)) {
66 : 172 : mp_print_str(print, "+");
67 : : }
68 : 272 : mp_format_float(o->imag, buf, sizeof(buf), 'g', precision, '\0');
69 : 272 : mp_printf(print, "%sj)", buf);
70 : : }
71 : 344 : }
72 : :
73 : 508 : static mp_obj_t complex_make_new(const mp_obj_type_t *type_in, size_t n_args, size_t n_kw, const mp_obj_t *args) {
74 : 508 : (void)type_in;
75 : 508 : mp_arg_check_num(n_args, n_kw, 0, 2, false);
76 : :
77 [ + + + ]: 508 : switch (n_args) {
78 : 4 : case 0:
79 : 4 : return mp_obj_new_complex(0, 0);
80 : :
81 : 104 : case 1:
82 [ + + + + : 104 : if (mp_obj_is_str(args[0])) {
+ + ]
83 : : // a string, parse it
84 : 72 : size_t l;
85 : 72 : const char *s = mp_obj_str_get_data(args[0], &l);
86 : 72 : return mp_parse_num_complex(s, l, NULL);
87 [ + + + + ]: 32 : } else if (mp_obj_is_type(args[0], &mp_type_complex)) {
88 : : // a complex, just return it
89 : : return args[0];
90 : : } else {
91 : 28 : mp_float_t real, imag;
92 : 28 : mp_obj_get_complex(args[0], &real, &imag);
93 : 16 : return mp_obj_new_complex(real, imag);
94 : : }
95 : :
96 : : case 2:
97 : : default: {
98 : 400 : mp_float_t real, imag;
99 [ + + + + ]: 400 : if (mp_obj_is_type(args[0], &mp_type_complex)) {
100 : 4 : mp_obj_complex_get(args[0], &real, &imag);
101 : : } else {
102 : 396 : real = mp_obj_get_float(args[0]);
103 : 396 : imag = 0;
104 : : }
105 [ + + + + ]: 400 : if (mp_obj_is_type(args[1], &mp_type_complex)) {
106 : 4 : mp_float_t real2, imag2;
107 : 4 : mp_obj_complex_get(args[1], &real2, &imag2);
108 : 4 : real -= imag2;
109 : 4 : imag += real2;
110 : : } else {
111 : 396 : imag += mp_obj_get_float(args[1]);
112 : : }
113 : 400 : return mp_obj_new_complex(real, imag);
114 : : }
115 : : }
116 : : }
117 : :
118 : 2884 : static mp_obj_t complex_unary_op(mp_unary_op_t op, mp_obj_t o_in) {
119 : 2884 : mp_obj_complex_t *o = MP_OBJ_TO_PTR(o_in);
120 [ + + + + : 2884 : switch (op) {
+ + ]
121 : 8 : case MP_UNARY_OP_BOOL:
122 [ + - + + ]: 8 : return mp_obj_new_bool(o->real != 0 || o->imag != 0);
123 : 8 : case MP_UNARY_OP_HASH:
124 : 8 : return MP_OBJ_NEW_SMALL_INT(mp_float_hash(o->real) ^ mp_float_hash(o->imag));
125 : : case MP_UNARY_OP_POSITIVE:
126 : : return o_in;
127 : 4 : case MP_UNARY_OP_NEGATIVE:
128 : 4 : return mp_obj_new_complex(-o->real, -o->imag);
129 : 8 : case MP_UNARY_OP_ABS:
130 : 8 : return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(o->real * o->real + o->imag * o->imag));
131 : 2852 : default:
132 : 2852 : return MP_OBJ_NULL; // op not supported
133 : : }
134 : : }
135 : :
136 : 1304 : static mp_obj_t complex_binary_op(mp_binary_op_t op, mp_obj_t lhs_in, mp_obj_t rhs_in) {
137 : 1304 : mp_obj_complex_t *lhs = MP_OBJ_TO_PTR(lhs_in);
138 : 1304 : return mp_obj_complex_binary_op(op, lhs->real, lhs->imag, rhs_in);
139 : : }
140 : :
141 : 2484 : static void complex_attr(mp_obj_t self_in, qstr attr, mp_obj_t *dest) {
142 [ + + ]: 2484 : if (dest[0] != MP_OBJ_NULL) {
143 : : // not load attribute
144 : : return;
145 : : }
146 : 2480 : mp_obj_complex_t *self = MP_OBJ_TO_PTR(self_in);
147 [ + + ]: 2480 : if (attr == MP_QSTR_real) {
148 : 1240 : dest[0] = mp_obj_new_float(self->real);
149 [ + - ]: 1240 : } else if (attr == MP_QSTR_imag) {
150 : 1240 : dest[0] = mp_obj_new_float(self->imag);
151 : : }
152 : : }
153 : :
154 : : MP_DEFINE_CONST_OBJ_TYPE(
155 : : mp_type_complex, MP_QSTR_complex, MP_TYPE_FLAG_EQ_NOT_REFLEXIVE | MP_TYPE_FLAG_EQ_CHECKS_OTHER_TYPE,
156 : : make_new, complex_make_new,
157 : : print, complex_print,
158 : : unary_op, complex_unary_op,
159 : : binary_op, complex_binary_op,
160 : : attr, complex_attr
161 : : );
162 : :
163 : 1988 : mp_obj_t mp_obj_new_complex(mp_float_t real, mp_float_t imag) {
164 : 1988 : mp_obj_complex_t *o = mp_obj_malloc(mp_obj_complex_t, &mp_type_complex);
165 : 1988 : o->real = real;
166 : 1988 : o->imag = imag;
167 : 1988 : return MP_OBJ_FROM_PTR(o);
168 : : }
169 : :
170 : 2864 : void mp_obj_complex_get(mp_obj_t self_in, mp_float_t *real, mp_float_t *imag) {
171 [ + - - + ]: 2864 : assert(mp_obj_is_type(self_in, &mp_type_complex));
172 : 2864 : mp_obj_complex_t *self = MP_OBJ_TO_PTR(self_in);
173 : 2864 : *real = self->real;
174 : 2864 : *imag = self->imag;
175 : 2864 : }
176 : :
177 : 1384 : mp_obj_t mp_obj_complex_binary_op(mp_binary_op_t op, mp_float_t lhs_real, mp_float_t lhs_imag, mp_obj_t rhs_in) {
178 : 1384 : mp_float_t rhs_real, rhs_imag;
179 [ + + ]: 1384 : if (!mp_obj_get_complex_maybe(rhs_in, &rhs_real, &rhs_imag)) {
180 : : return MP_OBJ_NULL; // op not supported
181 : : }
182 : :
183 [ + + + + : 1376 : switch (op) {
+ + + + ]
184 : 72 : case MP_BINARY_OP_ADD:
185 : : case MP_BINARY_OP_INPLACE_ADD:
186 : 72 : lhs_real += rhs_real;
187 : 72 : lhs_imag += rhs_imag;
188 : 72 : break;
189 : 8 : case MP_BINARY_OP_SUBTRACT:
190 : : case MP_BINARY_OP_INPLACE_SUBTRACT:
191 : 8 : lhs_real -= rhs_real;
192 : 8 : lhs_imag -= rhs_imag;
193 : 8 : break;
194 : : case MP_BINARY_OP_MULTIPLY:
195 : : case MP_BINARY_OP_INPLACE_MULTIPLY: {
196 : 36 : mp_float_t real;
197 : 36 : multiply:
198 : 36 : real = lhs_real * rhs_real - lhs_imag * rhs_imag;
199 : 36 : lhs_imag = lhs_real * rhs_imag + lhs_imag * rhs_real;
200 : 36 : lhs_real = real;
201 : 36 : break;
202 : : }
203 : : case MP_BINARY_OP_FLOOR_DIVIDE:
204 : : case MP_BINARY_OP_INPLACE_FLOOR_DIVIDE:
205 : 4 : mp_raise_TypeError(MP_ERROR_TEXT("can't truncate-divide a complex number"));
206 : :
207 : 16 : case MP_BINARY_OP_TRUE_DIVIDE:
208 : : case MP_BINARY_OP_INPLACE_TRUE_DIVIDE:
209 [ + + ]: 16 : if (rhs_imag == 0) {
210 [ + + ]: 8 : if (rhs_real == 0) {
211 : 4 : mp_raise_msg(&mp_type_ZeroDivisionError, MP_ERROR_TEXT("complex divide by zero"));
212 : : }
213 : 4 : lhs_real /= rhs_real;
214 : 4 : lhs_imag /= rhs_real;
215 [ + + ]: 8 : } else if (rhs_real == 0) {
216 : 4 : mp_float_t real = lhs_imag / rhs_imag;
217 : 4 : lhs_imag = -lhs_real / rhs_imag;
218 : 4 : lhs_real = real;
219 : : } else {
220 : 4 : mp_float_t rhs_len_sq = rhs_real * rhs_real + rhs_imag * rhs_imag;
221 : 4 : rhs_real /= rhs_len_sq;
222 : 4 : rhs_imag /= -rhs_len_sq;
223 : 4 : goto multiply;
224 : : }
225 : : break;
226 : :
227 : 36 : case MP_BINARY_OP_POWER:
228 : : case MP_BINARY_OP_INPLACE_POWER: {
229 : : // z1**z2 = exp(z2*ln(z1))
230 : : // = exp(z2*(ln(|z1|)+i*arg(z1)))
231 : : // = exp( (x2*ln1 - y2*arg1) + i*(y2*ln1 + x2*arg1) )
232 : : // = exp(x3 + i*y3)
233 : : // = exp(x3)*(cos(y3) + i*sin(y3))
234 : 36 : mp_float_t abs1 = MICROPY_FLOAT_C_FUN(sqrt)(lhs_real * lhs_real + lhs_imag * lhs_imag);
235 [ + + ]: 36 : if (abs1 == 0) {
236 [ + + + + ]: 20 : if (rhs_imag == 0 && rhs_real >= 0) {
237 [ + + ]: 12 : lhs_real = (rhs_real == 0);
238 : : } else {
239 : 8 : mp_raise_msg(&mp_type_ZeroDivisionError, MP_ERROR_TEXT("0.0 to a complex power"));
240 : : }
241 : : } else {
242 : 16 : mp_float_t ln1 = MICROPY_FLOAT_C_FUN(log)(abs1);
243 : 16 : mp_float_t arg1 = MICROPY_FLOAT_C_FUN(atan2)(lhs_imag, lhs_real);
244 : 16 : mp_float_t x3 = rhs_real * ln1 - rhs_imag * arg1;
245 : 16 : mp_float_t y3 = rhs_imag * ln1 + rhs_real * arg1;
246 : 16 : mp_float_t exp_x3 = MICROPY_FLOAT_C_FUN(exp)(x3);
247 : 16 : lhs_real = exp_x3 * MICROPY_FLOAT_C_FUN(cos)(y3);
248 : 16 : lhs_imag = exp_x3 * MICROPY_FLOAT_C_FUN(sin)(y3);
249 : : }
250 : : break;
251 : : }
252 : :
253 : 1200 : case MP_BINARY_OP_EQUAL:
254 [ + + + + ]: 1200 : return mp_obj_new_bool(lhs_real == rhs_real && lhs_imag == rhs_imag);
255 : :
256 : : default:
257 : : return MP_OBJ_NULL; // op not supported
258 : : }
259 : 152 : return mp_obj_new_complex(lhs_real, lhs_imag);
260 : : }
261 : :
262 : : #endif
|
