cython

1
# ticket: t305
2

3
from cpython.object cimport Py_EQ, Py_NE
4

5
cimport cython
6

7
DEF C21 = 2-1j
8

9

10
cdef class Complex3j:
11
    """
12
    >>> Complex3j() == 3j
13
    True
14
    >>> Complex3j() == Complex3j()
15
    True
16
    >>> Complex3j() != 3j
17
    False
18
    >>> Complex3j() != 3
19
    True
20
    >>> Complex3j() != Complex3j()
21
    False
22
    """
23
    def __richcmp__(a, b, int op):
24
        if op == Py_EQ or op == Py_NE:
25
            if isinstance(a, Complex3j):
26
                eq = isinstance(b, Complex3j) or b == 3j
27
            else:
28
                eq = isinstance(b, Complex3j) and a == 3j
29
            return eq if op == Py_EQ else not eq
30
        return NotImplemented
31

32

33
def test_object_conversion(o):
34
    """
35
    >>> test_object_conversion(2)
36
    ((2+0j), (2+0j))
37
    >>> test_object_conversion(2j - 0.5)
38
    ((-0.5+2j), (-0.5+2j))
39
    """
40
    cdef float complex a = o
41
    cdef double complex b = o
42
    return (a, b)
43

44

45
def test_arithmetic(double complex z, double complex w):
46
    """
47
    >>> test_arithmetic(2j, 4j)
48
    (2j, -2j, 6j, -2j, (-8+0j), (0.5+0j))
49
    >>> test_arithmetic(6+12j, 3j)
50
    ((6+12j), (-6-12j), (6+15j), (6+9j), (-36+18j), (4-2j))
51
    >>> test_arithmetic(5-10j, 3+4j)
52
    ((5-10j), (-5+10j), (8-6j), (2-14j), (55-10j), (-1-2j))
53
    """
54
    return +z, -z+0, z+w, z-w, z*w, z/w
55

56

57
def test_div(double complex a, double complex b, expected):
58
    """
59
    >>> big = 2.0**1023
60
    >>> test_div(1 + 1j, 1 + big*1j, 1/big - 1j/big)
61
    >>> test_div(1 + 1j, 1/big + 1j/big, big)
62
    """
63
    # Can't count on good c99 complex division :(
64
    if '_c99_' not in __name__:
65
        assert a / b == expected, (a / b, expected)
66

67

68
def test_pow(double complex z, double complex w, tol=None):
69
    """
70
    Various implementations produce slightly different results...
71

72
    >>> a = complex(3, 1)
73
    >>> test_pow(a, 1, 1e-15)
74
    True
75
    >>> test_pow(a, 2, 1e-15)
76
    True
77
    >>> test_pow(a, a, 1e-15)
78
    True
79
    >>> test_pow(complex(0.5, -.25), complex(3, 4), 1e-15)
80
    True
81
    >>> test_pow(-0.5, 1j, tol=1e-15)
82
    True
83
    >>> test_pow(-1, 0.5, tol=1e-15)
84
    True
85
    """
86
    if tol is None:
87
        return z**w
88
    else:
89
        return abs(z**w / <object>z ** <object>w - 1) < tol
90

91

92
def test_int_pow(double complex z, int n, tol=None):
93
    """
94
    >>> [test_int_pow(complex(0, 1), k, 1e-15) for k in range(-4, 5)]
95
    [True, True, True, True, True, True, True, True, True]
96
    >>> [test_int_pow(complex(0, 2), k, 1e-15) for k in range(-4, 5)]
97
    [True, True, True, True, True, True, True, True, True]
98
    >>> [test_int_pow(complex(2, 0.5), k, 1e-14) for k in range(0, 10)]
99
    [True, True, True, True, True, True, True, True, True, True]
100
    >>> test_int_pow(-0.5, 5, tol=1e-15)
101
    True
102
    """
103
    if tol is None:
104
        return z**n + <object>0 # add zero to normalize zero sign
105
    else:
106
        return abs(z**n / <object>z ** <object>n - 1) < tol
107

108

109
@cython.cdivision(False)
110
def test_div_by_zero(double complex z):
111
    """
112
    >>> test_div_by_zero(4j)
113
    -0.25j
114
    >>> test_div_by_zero(0)
115
    Traceback (most recent call last):
116
    ...
117
    ZeroDivisionError: float division
118
    """
119
    return 1/z
120

121

122
def test_coercion(int a, float b, double c, float complex d, double complex e):
123
    """
124
    >>> test_coercion(1, 1.5, 2.5, 4+1j, 10j)
125
    (1+0j)
126
    (1.5+0j)
127
    (2.5+0j)
128
    (4+1j)
129
    10j
130
    (9+21j)
131
    """
132
    cdef double complex z
133
    z = a; print z
134
    z = b; print z
135
    z = c; print z
136
    z = d; print z
137
    z = e; print z
138
    return z + a + b + c + d + e
139

140

141
def test_compare(double complex a, double complex b):
142
    """
143
    >>> test_compare(3, 3)
144
    (True, False, False, False, False, True, False)
145
    >>> test_compare(3j, 3j)
146
    (True, False, True, True, True, False, False)
147
    >>> test_compare(3j, 4j)
148
    (False, True, True, False, True, True, False)
149
    >>> test_compare(3, 4)
150
    (False, True, False, False, False, True, False)
151
    >>> test_compare(2-1j, 4)
152
    (False, True, False, False, False, True, True)
153
    """
154
    return a == b, a != b, a == 3j, 3j == b, a == Complex3j(), Complex3j() != b, a == C21
155

156

157
def test_compare_coerce(double complex a, int b):
158
    """
159
    >>> test_compare_coerce(3, 4)
160
    (False, True, False, False, False, True)
161
    >>> test_compare_coerce(4+1j, 4)
162
    (False, True, False, True, False, True)
163
    >>> test_compare_coerce(4, 4)
164
    (True, False, False, False, False, True)
165
    >>> test_compare_coerce(3j, 4)
166
    (False, True, True, False, True, False)
167
    """
168
    return a == b, a != b, a == 3j, 4+1j == a, a == Complex3j(), Complex3j() != a
169

170

171
def test_literal():
172
    """
173
    >>> test_literal()
174
    (5j, (1-2.5j), (2-1j))
175
    """
176
    return 5j, 1-2.5j, C21
177

178

179
def test_real_imag(double complex z):
180
    """
181
    >>> test_real_imag(1-3j)
182
    (1.0, -3.0)
183
    >>> test_real_imag(5)
184
    (5.0, 0.0)
185
    >>> test_real_imag(1.5j)
186
    (0.0, 1.5)
187
    """
188
    return z.real, z.imag
189

190
def test_real_imag_assignment(object a, double b):
191
    """
192
    >>> test_real_imag_assignment(1, 2)
193
    (1+2j)
194
    >>> test_real_imag_assignment(1.5, -3.5)
195
    (1.5-3.5j)
196
    """
197
    cdef double complex z
198
    z.real = a
199
    z.imag = b
200
    return z
201

202
def test_conjugate(float complex z):
203
    """
204
    >>> test_conjugate(2+3j)
205
    (2-3j)
206
    """
207
    return z.conjugate()
208

209
def test_conjugate_double(double complex z):
210
    """
211
    >>> test_conjugate_double(2+3j)
212
    (2-3j)
213
    """
214
    return z.conjugate()
215

216
ctypedef double complex cdouble
217
def test_conjugate_typedef(cdouble z):
218
    """
219
    >>> test_conjugate_typedef(2+3j)
220
    (2-3j)
221
    """
222
    return z.conjugate()
223

224
cdef cdouble test_conjugate_nogil(cdouble z) nogil:
225
    # Really just a compile test.
226
    return z.conjugate()
227
test_conjugate_nogil(0) # use it
228

229
## cdef extern from "complex_numbers_T305.h":
230
##     ctypedef double double_really_float "myfloat"
231
##     ctypedef float float_really_double "mydouble"
232
##     ctypedef float real_float "myfloat"
233
##     ctypedef double real_double "mydouble"
234

235
## def test_conjugate_nosizeassumptions(double_really_float x,
236
##                                      float_really_double y,
237
##                                      real_float z, real_double w):
238
##     """
239
##     >>> test_conjugate_nosizeassumptions(1, 1, 1, 1)
240
##     (-1j, -1j, -1j, -1j)
241
##     >>> ["%.2f" % x.imag for x in test_conjugate_nosizeassumptions(2e300, 2e300, 2e300, 2e300)]
242
##     ['-inf', '-2e+300', '-inf', '-2e+300']
243
##     """
244
##     cdef double complex I = 1j
245
##     return ((x*I).conjugate(), (y*I).conjugate(), (z*I).conjugate(), (w*I).conjugate())
246

247
ctypedef double mydouble
248
def test_coerce_typedef_multiply(mydouble x, double complex z):
249
    """
250
    >>> test_coerce_typedef_multiply(3, 1+1j)
251
    (3+3j)
252
    """
253
    return x * z
254

255
ctypedef int myint
256
def test_coerce_typedef_multiply_int(myint x, double complex z):
257
    """
258
    >>> test_coerce_typedef_multiply_int(3, 1+1j)
259
    (3+3j)
260
    """
261
    return x * z
262

263
cpdef double complex complex_retval():
264
    """
265
    >>> complex_retval()
266
    1j
267
    """
268
    return 1j
269

270
def stress_test():
271
    """
272
    Run the main operations on 1000 pseudo-random numbers to
273
    try to spot anything accidentally missed from the test cases
274
    (doesn't cover inf and NaN as inputs though)
275
    >>> stress_test()
276
    """
277
    cdef double complex x
278
    cdef double complex y
279

280
    from random import Random
281
    from math import ldexp
282
    r = Random()
283
    r.seed("I'm a seed")  # try to make the test somewhat reproducible
284

285
    # copied from https://docs.python.org/3/library/random.html#recipes
286
    # gets evenly distributed random numbers
287
    def full_random():
288
        mantissa = 0x10_0000_0000_0000 | r.getrandbits(52)
289
        exponent = -53
290
        x = 0
291
        while not x:
292
            x = r.getrandbits(32)
293
            exponent += x.bit_length() - 32
294
        return ldexp(mantissa, exponent)
295

296
    for n in range(1, 1001):
297
        if n % 50 == 0:
298
            # strategical insert some 0 values
299
            a = 0
300
        else:
301
            a = full_random()
302
        if n % 51 == 0:
303
            b = 0
304
        else:
305
            b = full_random()
306
        if n % 52 == 0:
307
            c = 0
308
        else:
309
            c = full_random()
310
        if n % 53 == 0:
311
            d = 0
312
        else:
313
            d = full_random()
314

315
        x= a+1j*b
316
        y = c+1j*d
317
        py_dict = dict(x=x, y=y)
318

319
        sum_ = x+y
320
        sum_py = eval("x+y", py_dict)
321
        delta_sum = abs(sum_/sum_py - 1)
322
        assert delta_sum < 1e-15, f"{x} {y} {sum_} {sum_py} {delta_sum}"
323

324
        minus = x-y
325
        minus_py = eval("x-y", py_dict)
326
        delta_minus = abs(minus/minus_py - 1)
327
        assert delta_minus < 1e-15, f"{x} {y} {minus} {minus_py} {delta_minus}"
328

329
        times = x*y
330
        times_py = eval("x*y", py_dict)
331
        delta_times = abs(times/times_py - 1)
332
        assert delta_times < 1e-15, f"{x} {y} {times} {times_py} {delta_times}"
333

334
        divide = x/y
335
        divide_py = eval("x/y", py_dict)
336
        delta_divide = abs(divide/divide_py - 1)
337
        assert delta_divide < 1e-15, f"{x} {y} {divide} {divide_py} {delta_divide}"
338

339
        divide2 = y/x
340
        divide2_py = eval("y/x", py_dict)
341
        delta_divide2 = abs(divide2/divide2_py - 1)
342
        assert delta_divide2 < 1e-15, f"{x} {y} {divide2} {divide2_py} {delta_divide2}"
343

344
        pow_ = x**y
345
        pow_py = eval("x**y", py_dict)
346
        delta_pow = abs(pow_/pow_py - 1)
347
        assert delta_pow < 1e-15, f"{x} {y} {pow_} {pow_py} {delta_pow}"
348

349
        pow2 = y**x
350
        pow2_py = eval("y**x", py_dict)
351
        delta_pow2 = abs(pow2/pow2_py - 1)
352
        assert delta_pow2 < 1e-15, f"{x} {y} {pow2} {pow2_py} {delta_pow2}"
353