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 
complex_print(const mp_print_t * print,mp_obj_t o_in,mp_print_kind_t kind)45 STATIC void complex_print(const mp_print_t *print, mp_obj_t o_in, mp_print_kind_t kind) {
46     (void)kind;
47     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     char buf[32];
57     const int precision = 16;
58     #endif
59     if (o->real == 0) {
60         mp_format_float(o->imag, buf, sizeof(buf), 'g', precision, '\0');
61         mp_printf(print, "%sj", buf);
62     } else {
63         mp_format_float(o->real, buf, sizeof(buf), 'g', precision, '\0');
64         mp_printf(print, "(%s", buf);
65         if (o->imag >= 0 || isnan(o->imag)) {
66             mp_print_str(print, "+");
67         }
68         mp_format_float(o->imag, buf, sizeof(buf), 'g', precision, '\0');
69         mp_printf(print, "%sj)", buf);
70     }
71 }
72 
complex_make_new(const mp_obj_type_t * type_in,size_t n_args,size_t n_kw,const mp_obj_t * args)73 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     (void)type_in;
75     mp_arg_check_num(n_args, n_kw, 0, 2, false);
76 
77     switch (n_args) {
78         case 0:
79             return mp_obj_new_complex(0, 0);
80 
81         case 1:
82             if (mp_obj_is_str(args[0])) {
83                 // a string, parse it
84                 size_t l;
85                 const char *s = mp_obj_str_get_data(args[0], &l);
86                 return mp_parse_num_decimal(s, l, true, true, NULL);
87             } else if (mp_obj_is_type(args[0], &mp_type_complex)) {
88                 // a complex, just return it
89                 return args[0];
90             } else {
91                 // something else, try to cast it to a complex
92                 return mp_obj_new_complex(mp_obj_get_float(args[0]), 0);
93             }
94 
95         case 2:
96         default: {
97             mp_float_t real, imag;
98             if (mp_obj_is_type(args[0], &mp_type_complex)) {
99                 mp_obj_complex_get(args[0], &real, &imag);
100             } else {
101                 real = mp_obj_get_float(args[0]);
102                 imag = 0;
103             }
104             if (mp_obj_is_type(args[1], &mp_type_complex)) {
105                 mp_float_t real2, imag2;
106                 mp_obj_complex_get(args[1], &real2, &imag2);
107                 real -= imag2;
108                 imag += real2;
109             } else {
110                 imag += mp_obj_get_float(args[1]);
111             }
112             return mp_obj_new_complex(real, imag);
113         }
114     }
115 }
116 
complex_unary_op(mp_unary_op_t op,mp_obj_t o_in)117 STATIC mp_obj_t complex_unary_op(mp_unary_op_t op, mp_obj_t o_in) {
118     mp_obj_complex_t *o = MP_OBJ_TO_PTR(o_in);
119     switch (op) {
120         case MP_UNARY_OP_BOOL:
121             return mp_obj_new_bool(o->real != 0 || o->imag != 0);
122         case MP_UNARY_OP_HASH:
123             return MP_OBJ_NEW_SMALL_INT(mp_float_hash(o->real) ^ mp_float_hash(o->imag));
124         case MP_UNARY_OP_POSITIVE:
125             return o_in;
126         case MP_UNARY_OP_NEGATIVE:
127             return mp_obj_new_complex(-o->real, -o->imag);
128         case MP_UNARY_OP_ABS:
129             return mp_obj_new_float(MICROPY_FLOAT_C_FUN(sqrt)(o->real * o->real + o->imag * o->imag));
130         default:
131             return MP_OBJ_NULL;      // op not supported
132     }
133 }
134 
complex_binary_op(mp_binary_op_t op,mp_obj_t lhs_in,mp_obj_t rhs_in)135 STATIC mp_obj_t complex_binary_op(mp_binary_op_t op, mp_obj_t lhs_in, mp_obj_t rhs_in) {
136     mp_obj_complex_t *lhs = MP_OBJ_TO_PTR(lhs_in);
137     return mp_obj_complex_binary_op(op, lhs->real, lhs->imag, rhs_in);
138 }
139 
complex_attr(mp_obj_t self_in,qstr attr,mp_obj_t * dest)140 STATIC void complex_attr(mp_obj_t self_in, qstr attr, mp_obj_t *dest) {
141     if (dest[0] != MP_OBJ_NULL) {
142         // not load attribute
143         return;
144     }
145     mp_obj_complex_t *self = MP_OBJ_TO_PTR(self_in);
146     if (attr == MP_QSTR_real) {
147         dest[0] = mp_obj_new_float(self->real);
148     } else if (attr == MP_QSTR_imag) {
149         dest[0] = mp_obj_new_float(self->imag);
150     }
151 }
152 
153 const mp_obj_type_t mp_type_complex = {
154     { &mp_type_type },
155     .flags = MP_TYPE_FLAG_EQ_NOT_REFLEXIVE | MP_TYPE_FLAG_EQ_CHECKS_OTHER_TYPE,
156     .name = MP_QSTR_complex,
157     .print = complex_print,
158     .make_new = complex_make_new,
159     .unary_op = complex_unary_op,
160     .binary_op = complex_binary_op,
161     .attr = complex_attr,
162 };
163 
mp_obj_new_complex(mp_float_t real,mp_float_t imag)164 mp_obj_t mp_obj_new_complex(mp_float_t real, mp_float_t imag) {
165     mp_obj_complex_t *o = m_new_obj(mp_obj_complex_t);
166     o->base.type = &mp_type_complex;
167     o->real = real;
168     o->imag = imag;
169     return MP_OBJ_FROM_PTR(o);
170 }
171 
mp_obj_complex_get(mp_obj_t self_in,mp_float_t * real,mp_float_t * imag)172 void mp_obj_complex_get(mp_obj_t self_in, mp_float_t *real, mp_float_t *imag) {
173     assert(mp_obj_is_type(self_in, &mp_type_complex));
174     mp_obj_complex_t *self = MP_OBJ_TO_PTR(self_in);
175     *real = self->real;
176     *imag = self->imag;
177 }
178 
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)179 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) {
180     mp_float_t rhs_real, rhs_imag;
181     if (!mp_obj_get_complex_maybe(rhs_in, &rhs_real, &rhs_imag)) {
182         return MP_OBJ_NULL; // op not supported
183     }
184 
185     switch (op) {
186         case MP_BINARY_OP_ADD:
187         case MP_BINARY_OP_INPLACE_ADD:
188             lhs_real += rhs_real;
189             lhs_imag += rhs_imag;
190             break;
191         case MP_BINARY_OP_SUBTRACT:
192         case MP_BINARY_OP_INPLACE_SUBTRACT:
193             lhs_real -= rhs_real;
194             lhs_imag -= rhs_imag;
195             break;
196         case MP_BINARY_OP_MULTIPLY:
197         case MP_BINARY_OP_INPLACE_MULTIPLY: {
198             mp_float_t real;
199         multiply:
200             real = lhs_real * rhs_real - lhs_imag * rhs_imag;
201             lhs_imag = lhs_real * rhs_imag + lhs_imag * rhs_real;
202             lhs_real = real;
203             break;
204         }
205         case MP_BINARY_OP_FLOOR_DIVIDE:
206         case MP_BINARY_OP_INPLACE_FLOOR_DIVIDE:
207             mp_raise_TypeError(MP_ERROR_TEXT("can't truncate-divide a complex number"));
208 
209         case MP_BINARY_OP_TRUE_DIVIDE:
210         case MP_BINARY_OP_INPLACE_TRUE_DIVIDE:
211             if (rhs_imag == 0) {
212                 if (rhs_real == 0) {
213                     mp_raise_msg(&mp_type_ZeroDivisionError, MP_ERROR_TEXT("complex divide by zero"));
214                 }
215                 lhs_real /= rhs_real;
216                 lhs_imag /= rhs_real;
217             } else if (rhs_real == 0) {
218                 mp_float_t real = lhs_imag / rhs_imag;
219                 lhs_imag = -lhs_real / rhs_imag;
220                 lhs_real = real;
221             } else {
222                 mp_float_t rhs_len_sq = rhs_real * rhs_real + rhs_imag * rhs_imag;
223                 rhs_real /= rhs_len_sq;
224                 rhs_imag /= -rhs_len_sq;
225                 goto multiply;
226             }
227             break;
228 
229         case MP_BINARY_OP_POWER:
230         case MP_BINARY_OP_INPLACE_POWER: {
231             // z1**z2 = exp(z2*ln(z1))
232             //        = exp(z2*(ln(|z1|)+i*arg(z1)))
233             //        = exp( (x2*ln1 - y2*arg1) + i*(y2*ln1 + x2*arg1) )
234             //        = exp(x3 + i*y3)
235             //        = exp(x3)*(cos(y3) + i*sin(y3))
236             mp_float_t abs1 = MICROPY_FLOAT_C_FUN(sqrt)(lhs_real * lhs_real + lhs_imag * lhs_imag);
237             if (abs1 == 0) {
238                 if (rhs_imag == 0 && rhs_real >= 0) {
239                     lhs_real = (rhs_real == 0);
240                 } else {
241                     mp_raise_msg(&mp_type_ZeroDivisionError, MP_ERROR_TEXT("0.0 to a complex power"));
242                 }
243             } else {
244                 mp_float_t ln1 = MICROPY_FLOAT_C_FUN(log)(abs1);
245                 mp_float_t arg1 = MICROPY_FLOAT_C_FUN(atan2)(lhs_imag, lhs_real);
246                 mp_float_t x3 = rhs_real * ln1 - rhs_imag * arg1;
247                 mp_float_t y3 = rhs_imag * ln1 + rhs_real * arg1;
248                 mp_float_t exp_x3 = MICROPY_FLOAT_C_FUN(exp)(x3);
249                 lhs_real = exp_x3 * MICROPY_FLOAT_C_FUN(cos)(y3);
250                 lhs_imag = exp_x3 * MICROPY_FLOAT_C_FUN(sin)(y3);
251             }
252             break;
253         }
254 
255         case MP_BINARY_OP_EQUAL:
256             return mp_obj_new_bool(lhs_real == rhs_real && lhs_imag == rhs_imag);
257 
258         default:
259             return MP_OBJ_NULL; // op not supported
260     }
261     return mp_obj_new_complex(lhs_real, lhs_imag);
262 }
263 
264 #endif
265