Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/rings/complex_number.pyx
4057 views
1
"""
2
Arbitrary Precision Complex Numbers
3
4
AUTHORS:
5
6
- William Stein (2006-01-26): complete rewrite
7
8
- Joel B. Mohler (2006-12-16): naive rewrite into pyrex
9
10
- William Stein(2007-01): rewrite of Mohler's rewrite
11
12
- Vincent Delecroix (2010-01): plot function
13
"""
14
15
#*****************************************************************************
16
# Copyright (C) 2006 William Stein <[email protected]>
17
#
18
# Distributed under the terms of the GNU General Public License (GPL)
19
# as published by the Free Software Foundation; either version 2 of
20
# the License, or (at your option) any later version.
21
# http://www.gnu.org/licenses/
22
#*****************************************************************************
23
24
25
import math
26
import operator
27
28
from sage.structure.element cimport FieldElement, RingElement, Element, ModuleElement
29
from sage.categories.map cimport Map
30
31
from complex_double cimport ComplexDoubleElement
32
from real_mpfr cimport RealNumber
33
34
import complex_field
35
import sage.misc.misc
36
import integer
37
import infinity
38
39
from sage.libs.mpmath.utils cimport mpfr_to_mpfval
40
41
include "../ext/stdsage.pxi"
42
43
cdef object numpy_complex_interface = {'typestr': '=c16'}
44
cdef object numpy_object_interface = {'typestr': '|O'}
45
46
cdef mp_rnd_t rnd
47
rnd = GMP_RNDN
48
49
cdef double LOG_TEN_TWO_PLUS_EPSILON = 3.321928094887363 # a small overestimate of log(10,2)
50
51
def set_global_complex_round_mode(n):
52
global rnd
53
rnd = n
54
55
#from sage.databases.odlyzko import zeta_zeroes
56
57
def is_ComplexNumber(x):
58
r"""
59
Returns True if x is a complex number. In particular, if x is of
60
the ``ComplexNumber`` type.
61
62
EXAMPLES::
63
64
sage: from sage.rings.complex_number import is_ComplexNumber
65
sage: a = ComplexNumber(1,2); a
66
1.00000000000000 + 2.00000000000000*I
67
sage: is_ComplexNumber(a)
68
True
69
sage: b = ComplexNumber(1); b
70
1.00000000000000
71
sage: is_ComplexNumber(b)
72
True
73
74
Note that the global element I is of type
75
``SymbolicConstant``. However, elements of the class
76
``ComplexField_class`` are of type
77
``ComplexNumber``::
78
79
sage: c = 1 + 2*I
80
sage: is_ComplexNumber(c)
81
False
82
sage: d = CC(1 + 2*I)
83
sage: is_ComplexNumber(d)
84
True
85
"""
86
return isinstance(x, ComplexNumber)
87
88
cdef class ComplexNumber(sage.structure.element.FieldElement):
89
"""
90
A floating point approximation to a complex number using any
91
specified precision. Answers derived from calculations with such
92
approximations may differ from what they would be if those
93
calculations were performed with true complex numbers. This is due
94
to the rounding errors inherent to finite precision calculations.
95
96
EXAMPLES::
97
98
sage: I = CC.0
99
sage: b = 1.5 + 2.5*I
100
sage: loads(b.dumps()) == b
101
True
102
"""
103
104
cdef ComplexNumber _new(self):
105
"""
106
Quickly creates a new initialized complex number with the same
107
parent as self.
108
"""
109
cdef ComplexNumber x
110
x = PY_NEW(ComplexNumber)
111
x._parent = self._parent
112
x._prec = self._prec
113
x._multiplicative_order = None
114
mpfr_init2(x.__re, self._prec)
115
mpfr_init2(x.__im, self._prec)
116
return x
117
118
def __cinit__(self, parent=None, real=None, imag=None):
119
self._prec = -1
120
121
def __init__(self, parent, real, imag=None):
122
r"""
123
Initialize ``ComplexNumber`` instance.
124
125
EXAMPLES::
126
127
sage: a = ComplexNumber(2,1)
128
sage: a.__init__(CC,2,1)
129
sage: a
130
2.00000000000000 + 1.00000000000000*I
131
sage: parent(a)
132
Complex Field with 53 bits of precision
133
sage: real(a)
134
2.00000000000000
135
sage: imag(a)
136
1.00000000000000
137
"""
138
cdef real_mpfr.RealNumber rr, ii
139
self._parent = parent
140
self._prec = self._parent._prec
141
self._multiplicative_order = None
142
143
mpfr_init2(self.__re, self._prec)
144
mpfr_init2(self.__im, self._prec)
145
146
if imag is None:
147
if real is None: return
148
149
if PY_TYPE_CHECK(real, ComplexNumber):
150
real, imag = (<ComplexNumber>real).real(), (<ComplexNumber>real).imag()
151
elif isinstance(real, sage.libs.pari.all.pari_gen):
152
real, imag = real.real(), real.imag()
153
elif isinstance(real, list) or isinstance(real, tuple):
154
re, imag = real
155
real = re
156
elif isinstance(real, complex):
157
real, imag = real.real, real.imag
158
else:
159
imag = 0
160
try:
161
R = parent._real_field()
162
rr = R(real)
163
ii = R(imag)
164
mpfr_set(self.__re, rr.value, rnd)
165
mpfr_set(self.__im, ii.value, rnd)
166
except TypeError:
167
raise TypeError, "unable to coerce to a ComplexNumber: %s" % type(real)
168
169
170
def __dealloc__(self):
171
"""
172
TESTS:
173
174
Check that #12038 is resolved::
175
176
sage: from sage.rings.complex_number import ComplexNumber as CN
177
sage: coerce(CN, 1+I)
178
Traceback (most recent call last):
179
...
180
TypeError: __init__() takes at least 2 positional arguments (1 given)
181
"""
182
if self._prec != -1:
183
mpfr_clear(self.__re)
184
mpfr_clear(self.__im)
185
186
def _interface_init_(self, I=None):
187
"""
188
Returns self formatted as a string, suitable as input to another
189
computer algebra system. (This is the default function used for
190
exporting to other computer algebra systems.)
191
192
EXAMPLES::
193
194
sage: s1 = CC(exp(I)); s1
195
0.540302305868140 + 0.841470984807897*I
196
sage: s1._interface_init_()
197
'0.54030230586813977 + 0.84147098480789650*I'
198
sage: s1 == CC(gp(s1))
199
True
200
"""
201
return self.str(truncate=False)
202
203
property __array_interface__:
204
def __get__(self):
205
"""
206
Used for NumPy conversion.
207
208
EXAMPLES::
209
210
sage: import numpy
211
sage: numpy.array([1.0, 2.5j]).dtype
212
dtype('complex128')
213
sage: numpy.array([1.000000000000000000000000000000000000j]).dtype
214
dtype('object')
215
"""
216
if self._prec <= 57: # max size of repr(float)
217
return numpy_complex_interface
218
else:
219
return numpy_object_interface
220
221
def _sage_input_(self, sib, coerced):
222
r"""
223
Produce an expression which will reproduce this value when evaluated.
224
225
EXAMPLES:
226
sage: for prec in (2, 53, 200):
227
... fld = ComplexField(prec)
228
... var = polygen(fld)
229
... ins = [-20, 0, 1, -2^4000, 2^-4000] + [fld._real_field().random_element() for _ in range(3)]
230
... for v1 in ins:
231
... for v2 in ins:
232
... v = fld(v1, v2)
233
... _ = sage_input(fld(v), verify=True)
234
... _ = sage_input(fld(v) * var, verify=True)
235
sage: x = polygen(CC)
236
sage: for v1 in [-2, 0, 2]:
237
... for v2 in [-2, -1, 0, 1, 2]:
238
... print str(sage_input(x + CC(v1, v2))).splitlines()[1]
239
x + CC(-2 - RR(2)*I)
240
x + CC(-2 - RR(1)*I)
241
x - 2
242
x + CC(-2 + RR(1)*I)
243
x + CC(-2 + RR(2)*I)
244
x - CC(RR(2)*I)
245
x - CC(RR(1)*I)
246
x
247
x + CC(RR(1)*I)
248
x + CC(RR(2)*I)
249
x + CC(2 - RR(2)*I)
250
x + CC(2 - RR(1)*I)
251
x + 2
252
x + CC(2 + RR(1)*I)
253
x + CC(2 + RR(2)*I)
254
sage: from sage.misc.sage_input import SageInputBuilder
255
sage: sib = SageInputBuilder()
256
sage: sib_np = SageInputBuilder(preparse=False)
257
sage: CC(-infinity)._sage_input_(sib, True)
258
{unop:- {call: {atomic:RR}({atomic:Infinity})}}
259
sage: CC(0, infinity)._sage_input_(sib, True)
260
{call: {atomic:CC}({call: {atomic:RR}({atomic:0})}, {call: {atomic:RR}({atomic:Infinity})})}
261
sage: CC(NaN, 5)._sage_input_(sib, True)
262
{call: {atomic:CC}({call: {atomic:RR}({atomic:NaN})}, {call: {atomic:RR}({atomic:5})})}
263
sage: CC(5, NaN)._sage_input_(sib, True)
264
{call: {atomic:CC}({call: {atomic:RR}({atomic:5})}, {call: {atomic:RR}({atomic:NaN})})}
265
sage: CC(12345)._sage_input_(sib, True)
266
{atomic:12345}
267
sage: CC(-12345)._sage_input_(sib, False)
268
{call: {atomic:CC}({binop:+ {unop:- {atomic:12345}} {binop:* {call: {atomic:RR}({atomic:0})} {atomic:I}}})}
269
sage: CC(0, 12345)._sage_input_(sib, True)
270
{call: {atomic:CC}({binop:* {call: {atomic:RR}({atomic:12345})} {atomic:I}})}
271
sage: CC(0, -12345)._sage_input_(sib, False)
272
{unop:- {call: {atomic:CC}({binop:* {call: {atomic:RR}({atomic:12345})} {atomic:I}})}}
273
sage: CC(1.579)._sage_input_(sib, True)
274
{atomic:1.579}
275
sage: CC(1.579)._sage_input_(sib_np, True)
276
{atomic:1.579}
277
sage: ComplexField(150).zeta(37)._sage_input_(sib, True)
278
{call: {call: {atomic:ComplexField}({atomic:150})}({binop:+ {atomic:0.98561591034770846226477029397621845736859851519} {binop:* {call: {call: {atomic:RealField}({atomic:150})}({atomic:0.16900082032184907409303555538443060626072476297})} {atomic:I}}})}
279
sage: ComplexField(150).zeta(37)._sage_input_(sib_np, True)
280
{call: {call: {atomic:ComplexField}({atomic:150})}({binop:+ {call: {call: {atomic:RealField}({atomic:150})}({atomic:'0.98561591034770846226477029397621845736859851519'})} {binop:* {call: {call: {atomic:RealField}({atomic:150})}({atomic:'0.16900082032184907409303555538443060626072476297'})} {atomic:I}}})}
281
"""
282
if coerced and self.imag() == 0:
283
return sib(self.real(), True)
284
285
# The body will be coerced first to symbolics and then to CC.
286
# This works fine if we produce integer or float literals, but
287
# not for infinity or NaN.
288
if not (mpfr_number_p(self.__re) and mpfr_number_p(self.__im)):
289
return sib(self.parent())(self.real(), self.imag())
290
291
# The following uses of .sum() and .prod() will simplify
292
# 3+0*I to 3, 0+1*I to I, etc.
293
real_part = sib(self.real(), 2)
294
imag_part = sib.prod([sib(self.imag()), sib.name('I')],
295
simplify=True)
296
sum = sib.sum([real_part, imag_part], simplify=True)
297
if sum._sie_is_negation():
298
return -sib(self.parent())(sum._sie_operand)
299
else:
300
return sib(self.parent())(sum)
301
302
# The following is an (untested) implementation that produces
303
# CC_I = CC.gen()
304
# 2 + 3*CC_I
305
# instead of CC(2 + 3*I)
306
# cdef int prec
307
308
# if self.real().is_zero() and self.imag() == 1:
309
# v = sib(self.parent()).gen()
310
# prec = self.prec()
311
# if prec == 53:
312
# gen_name = 'CC_I'
313
# else:
314
# gen_name = 'CC%d_I' % prec
315
# sib.cache(self, v, gen_name)
316
317
# real_part = sib(self.real())
318
# imag_part = sib.prod([self.imag(), self.parent().gen()], simplify=True)
319
# return sib.sum([real_part, imag_part], simplify=True)
320
321
def _repr_(self):
322
r"""
323
Returns self formatted as a string.
324
325
EXAMPLES::
326
327
sage: a = ComplexNumber(2,1); a
328
2.00000000000000 + 1.00000000000000*I
329
sage: a._repr_()
330
'2.00000000000000 + 1.00000000000000*I'
331
"""
332
return self.str(10)
333
334
def __hash__(self):
335
"""
336
Returns the hash of self, which coincides with the python complex
337
and float (and often int) types.
338
339
This has the drawback that two very close high precision numbers
340
will have the same hash, but allows them to play nicely with other
341
real types.
342
343
EXAMPLE::
344
345
sage: hash(CC(1.2, 33)) == hash(complex(1.2, 33))
346
True
347
"""
348
return hash(complex(self))
349
350
def __getitem__(self, i):
351
r"""
352
Returns either the real or imaginary component of self depending on
353
the choice of i: real (i=0), imaginary (i=1)
354
355
INPUTS:
356
357
- ``i`` - 0 or 1
358
359
- ``0`` - will return the real component of self
360
361
- ``1`` - will return the imaginary component of self
362
363
EXAMPLES::
364
365
sage: a = ComplexNumber(2,1)
366
sage: a.__getitem__(0)
367
2.00000000000000
368
sage: a.__getitem__(1)
369
1.00000000000000
370
371
::
372
373
sage: b = CC(42,0)
374
sage: b
375
42.0000000000000
376
sage: b.__getitem__(1)
377
0.000000000000000
378
"""
379
if i == 0:
380
return self.real()
381
elif i == 1:
382
return self.imag()
383
raise IndexError, "i must be between 0 and 1."
384
385
def __reduce__( self ):
386
"""
387
Pickling support
388
389
EXAMPLES::
390
391
sage: a = CC(1 + I)
392
sage: loads(dumps(a)) == a
393
True
394
"""
395
# TODO: This is potentially slow -- make a 1 version that
396
# is native and much faster -- doesn't use .real()/.imag()
397
return (make_ComplexNumber0, (self._parent, self._multiplicative_order, self.real(), self.imag()))
398
399
def _set_multiplicative_order(self, n):
400
r"""
401
Function for setting the ``ComplexNumber`` class
402
attribute ``multiplicative_order`` of self.
403
404
INPUTS: n - an integer which will define the multiplicative order
405
of self
406
407
EXAMPLES: Note that it is not advisable to explicitly call
408
``_set_multiplicative_order`` for explicitly declared
409
complex numbers.
410
411
::
412
413
sage: a = ComplexNumber(2,1)
414
sage: a.multiplicative_order()
415
+Infinity
416
sage: a._set_multiplicative_order(5)
417
sage: a.multiplicative_order()
418
5
419
sage: a^5
420
-38.0000000000000 + 41.0000000000000*I
421
"""
422
self._multiplicative_order = integer.Integer(n)
423
424
def str(self, base=10, truncate=True):
425
r"""
426
Return a string representation of this number.
427
428
INPUTS:
429
430
- ``base`` - The base to use for printing (default 10)
431
432
- ``truncate`` - (default: ``True``) Whether to print fewer
433
digits than are available, to mask errors in the last bits.
434
435
EXAMPLES::
436
437
sage: a = CC(pi + I*e)
438
sage: a.str()
439
'3.14159265358979 + 2.71828182845905*I'
440
sage: a.str(truncate=False)
441
'3.1415926535897931 + 2.7182818284590451*I'
442
sage: a.str(base=2)
443
'11.001001000011111101101010100010001000010110100011000 + 10.101101111110000101010001011000101000101011101101001*I'
444
sage: CC(0.5 + 0.625*I).str(base=2)
445
'0.10000000000000000000000000000000000000000000000000000 + 0.10100000000000000000000000000000000000000000000000000*I'
446
sage: a.str(base=16)
447
'3.243f6a8885a30 + 2.b7e151628aed2*I'
448
sage: a.str(base=36)
449
'3.53i5ab8p5fc + 2.puw5nggjf8f*I'
450
sage: CC(0)
451
0.000000000000000
452
"""
453
454
s = ""
455
if self.real() != 0:
456
s = self.real().str(base, truncate=truncate)
457
if self.imag() != 0:
458
y = self.imag()
459
if s!="":
460
if y < 0:
461
s = s+" - "
462
y = -y
463
else:
464
s = s+" + "
465
s = s+"%s*I"%y.str(base, truncate=truncate)
466
if len(s) == 0:
467
s = self.real().str(base, truncate=truncate)
468
return s
469
470
def _latex_(self):
471
r"""
472
Method for converting self to a string with latex formatting.
473
Called by the global function ``latex``.
474
475
EXAMPLES::
476
477
sage: a = ComplexNumber(2,1)
478
sage: a
479
2.00000000000000 + 1.00000000000000*I
480
sage: latex(a)
481
2.00000000000000 + 1.00000000000000i
482
sage: a._latex_()
483
'2.00000000000000 + 1.00000000000000i'
484
485
::
486
487
sage: b = ComplexNumber(7,4,min_prec=16)
488
sage: b
489
7.000 + 4.000*I
490
sage: latex(b)
491
7.000 + 4.000i
492
sage: b._latex_()
493
'7.000 + 4.000i'
494
495
::
496
497
sage: ComplexNumber(0).log()._latex_()
498
'-\\infty'
499
"""
500
import re
501
s = self.str().replace('*I', 'i').replace('infinity','\\infty')
502
return re.sub(r"e(-?\d+)", r" \\times 10^{\1}", s)
503
504
def _pari_(self):
505
r"""
506
Coerces self into a Pari ``complex`` object.
507
508
EXAMPLES: Coerce the object using the ``pari``
509
function::
510
511
sage: a = ComplexNumber(2,1)
512
sage: pari(a)
513
2.00000000000000 + 1.00000000000000*I
514
sage: type(pari(a))
515
<type 'sage.libs.pari.gen.gen'>
516
sage: a._pari_()
517
2.00000000000000 + 1.00000000000000*I
518
sage: type(a._pari_())
519
<type 'sage.libs.pari.gen.gen'>
520
"""
521
return sage.libs.pari.all.pari.complex(self.real()._pari_(), self.imag()._pari_())
522
523
def _mpmath_(self, prec=None, rounding=None):
524
"""
525
Returns an mpmath version of this ComplexNumber.
526
527
.. note::
528
529
Currently, the rounding mode is ignored.
530
531
EXAMPLES::
532
533
sage: CC(1,2)._mpmath_()
534
mpc(real='1.0', imag='2.0')
535
"""
536
if prec is not None:
537
return self.n(prec=prec)._mpmath_()
538
from sage.libs.mpmath.all import make_mpc
539
re = mpfr_to_mpfval(self.__re)
540
im = mpfr_to_mpfval(self.__im)
541
return make_mpc((re, im))
542
543
cpdef ModuleElement _add_(self, ModuleElement right):
544
cdef ComplexNumber x
545
x = self._new()
546
mpfr_add(x.__re, self.__re, (<ComplexNumber>right).__re, rnd)
547
mpfr_add(x.__im, self.__im, (<ComplexNumber>right).__im, rnd)
548
return x
549
550
cpdef ModuleElement _sub_(self, ModuleElement right):
551
cdef ComplexNumber x
552
x = self._new()
553
mpfr_sub(x.__re, self.__re, (<ComplexNumber>right).__re, rnd)
554
mpfr_sub(x.__im, self.__im, (<ComplexNumber>right).__im, rnd)
555
return x
556
557
cpdef RingElement _mul_(self, RingElement right):
558
cdef ComplexNumber x
559
x = self._new()
560
cdef mpfr_t t0, t1
561
mpfr_init2(t0, self._prec)
562
mpfr_init2(t1, self._prec)
563
mpfr_mul(t0, self.__re, (<ComplexNumber>right).__re, rnd)
564
mpfr_mul(t1, self.__im, (<ComplexNumber>right).__im, rnd)
565
mpfr_sub(x.__re, t0, t1, rnd)
566
mpfr_mul(t0, self.__re, (<ComplexNumber>right).__im, rnd)
567
mpfr_mul(t1, self.__im, (<ComplexNumber>right).__re, rnd)
568
mpfr_add(x.__im, t0, t1, rnd)
569
mpfr_clear(t0)
570
mpfr_clear(t1)
571
return x
572
573
def norm(self):
574
r"""
575
Returns the norm of this complex number. If `c = a + bi` is a
576
complex number, then the norm of `c` is defined as the product of
577
`c` and its complex conjugate
578
579
.. MATH::
580
581
\text{norm}(c)
582
=
583
\text{norm}(a + bi)
584
=
585
c \cdot \overline{c}
586
=
587
a^2 + b^2.
588
589
The norm of a complex number is different from its absolute value.
590
The absolute value of a complex number is defined to be the square
591
root of its norm. A typical use of the complex norm is in the
592
integral domain `\ZZ[i]` of Gaussian integers, where the norm of
593
each Gaussian integer `c = a + bi` is defined as its complex norm.
594
595
.. SEEALSO::
596
597
- :func:`sage.misc.functional.norm`
598
599
- :meth:`sage.rings.complex_double.ComplexDoubleElement.norm`
600
601
EXAMPLES:
602
603
This indeed acts as the square function when the
604
imaginary component of self is equal to zero::
605
606
sage: a = ComplexNumber(2,1)
607
sage: a.norm()
608
5.00000000000000
609
sage: b = ComplexNumber(4.2,0)
610
sage: b.norm()
611
17.6400000000000
612
sage: b^2
613
17.6400000000000
614
"""
615
return self.norm_c()
616
617
cdef real_mpfr.RealNumber norm_c(ComplexNumber self):
618
cdef real_mpfr.RealNumber x
619
x = real_mpfr.RealNumber(self._parent._real_field(), None)
620
621
cdef mpfr_t t0, t1
622
mpfr_init2(t0, self._prec)
623
mpfr_init2(t1, self._prec)
624
625
mpfr_mul(t0, self.__re, self.__re, rnd)
626
mpfr_mul(t1, self.__im, self.__im, rnd)
627
628
mpfr_add(x.value, t0, t1, rnd)
629
630
mpfr_clear(t0)
631
mpfr_clear(t1)
632
return x
633
634
cdef real_mpfr.RealNumber abs_c(ComplexNumber self):
635
cdef real_mpfr.RealNumber x
636
x = real_mpfr.RealNumber(self._parent._real_field(), None)
637
638
cdef mpfr_t t0, t1
639
mpfr_init2(t0, self._prec)
640
mpfr_init2(t1, self._prec)
641
642
mpfr_mul(t0, self.__re, self.__re, rnd)
643
mpfr_mul(t1, self.__im, self.__im, rnd)
644
645
mpfr_add(x.value, t0, t1, rnd)
646
mpfr_sqrt(x.value, x.value, rnd)
647
648
mpfr_clear(t0)
649
mpfr_clear(t1)
650
return x
651
652
cpdef RingElement _div_(self, RingElement right):
653
cdef ComplexNumber x
654
x = self._new()
655
cdef mpfr_t a, b, t0, t1, right_nm
656
mpfr_init2(t0, self._prec)
657
mpfr_init2(t1, self._prec)
658
mpfr_init2(a, self._prec)
659
mpfr_init2(b, self._prec)
660
mpfr_init2(right_nm, self._prec)
661
662
mpfr_mul(t0, (<ComplexNumber>right).__re, (<ComplexNumber>right).__re, rnd)
663
mpfr_mul(t1, (<ComplexNumber>right).__im, (<ComplexNumber>right).__im, rnd)
664
mpfr_add(right_nm, t0, t1, rnd)
665
666
mpfr_div(a, (<ComplexNumber>right).__re, right_nm, rnd)
667
mpfr_div(b, (<ComplexNumber>right).__im, right_nm, rnd)
668
669
## Do this: x.__re = a * self.__re + b * self.__im
670
mpfr_mul(t0, a, self.__re, rnd)
671
mpfr_mul(t1, b, self.__im, rnd)
672
mpfr_add(x.__re, t0, t1, rnd)
673
674
## Do this: x.__im = a * self.__im - b * self.__re
675
mpfr_mul(t0, a, self.__im, rnd)
676
mpfr_mul(t1, b, self.__re, rnd)
677
mpfr_sub(x.__im, t0, t1, rnd)
678
mpfr_clear(t0)
679
mpfr_clear(t1)
680
mpfr_clear(a)
681
mpfr_clear(b)
682
mpfr_clear(right_nm)
683
return x
684
685
def __rdiv__(self, left):
686
r"""
687
Returns the quotient of left with self, that is:
688
689
left/self
690
691
as a complex number.
692
693
INPUTS:
694
695
- ``left`` - a complex number to divide by self
696
697
EXAMPLES::
698
699
sage: a = ComplexNumber(2,0)
700
sage: a.__rdiv__(CC(1))
701
0.500000000000000
702
sage: CC(1)/a
703
0.500000000000000
704
"""
705
return ComplexNumber(self._parent, left)/self
706
707
def __pow__(self, right, modulus):
708
"""
709
EXAMPLES::
710
711
sage: C.<i> = ComplexField(20)
712
sage: a = i^2; a
713
-1.0000
714
sage: a.parent()
715
Complex Field with 20 bits of precision
716
sage: a = (1+i)^i; a
717
0.42883 + 0.15487*I
718
sage: (1+i)^(1+i)
719
0.27396 + 0.58370*I
720
sage: a.parent()
721
Complex Field with 20 bits of precision
722
sage: i^i
723
0.20788
724
sage: (2+i)^(0.5)
725
1.4553 + 0.34356*I
726
"""
727
if isinstance(right, (int, long, integer.Integer)):
728
return sage.rings.ring_element.RingElement.__pow__(self, right)
729
730
try:
731
return (self.log()*right).exp()
732
except TypeError:
733
pass
734
735
try:
736
self = right.parent()(self)
737
return self**right
738
except AttributeError:
739
raise TypeError
740
741
def _magma_init_(self, magma):
742
r"""
743
EXAMPLES::
744
745
sage: magma(CC([1, 2])) # optional - magma
746
1.00000000000000 + 2.00000000000000*$.1
747
sage: v = magma(CC([1, 2])).sage(); v # optional - magma
748
1.00000000000000 + 2.00000000000000*I
749
sage: v.parent() # optional - magma
750
Complex Field with 53 bits of precision
751
"""
752
return "%s![%s, %s]" % (self.parent()._magma_init_(magma), self.real(), self.imag())
753
754
def __nonzero__(self):
755
"""
756
Return True if self is not zero. This is an internal function; use
757
self.is_zero() instead.
758
759
EXAMPLES::
760
761
sage: z = 1 + CC(I)
762
sage: z.is_zero()
763
False
764
"""
765
return not (mpfr_zero_p(self.__re) and mpfr_zero_p(self.__im))
766
767
def prec(self):
768
"""
769
Return precision of this complex number.
770
771
EXAMPLES::
772
773
sage: i = ComplexField(2000).0
774
sage: i.prec()
775
2000
776
"""
777
return self._parent.prec()
778
779
def real(self):
780
"""
781
Return real part of self.
782
783
EXAMPLES::
784
785
sage: i = ComplexField(100).0
786
sage: z = 2 + 3*i
787
sage: x = z.real(); x
788
2.0000000000000000000000000000
789
sage: x.parent()
790
Real Field with 100 bits of precision
791
sage: z.real_part()
792
2.0000000000000000000000000000
793
"""
794
cdef real_mpfr.RealNumber x
795
x = real_mpfr.RealNumber(self._parent._real_field(), None)
796
mpfr_set(x.value, self.__re, rnd)
797
return x
798
799
real_part = real
800
801
def imag(self):
802
"""
803
Return imaginary part of self.
804
805
EXAMPLES::
806
807
sage: i = ComplexField(100).0
808
sage: z = 2 + 3*i
809
sage: x = z.imag(); x
810
3.0000000000000000000000000000
811
sage: x.parent()
812
Real Field with 100 bits of precision
813
sage: z.imag_part()
814
3.0000000000000000000000000000
815
"""
816
cdef real_mpfr.RealNumber x
817
x = real_mpfr.RealNumber(self._parent._real_field(), None)
818
mpfr_set(x.value, self.__im, rnd)
819
return x
820
821
imag_part = imag
822
823
def __neg__(self):
824
r"""
825
Method for computing the negative of self.
826
827
-(a + bi) = -a - bi
828
829
EXAMPLES::
830
831
sage: a = ComplexNumber(2,1)
832
sage: -a
833
-2.00000000000000 - 1.00000000000000*I
834
sage: a.__neg__()
835
-2.00000000000000 - 1.00000000000000*I
836
"""
837
cdef ComplexNumber x
838
x = self._new()
839
mpfr_neg(x.__re, self.__re, rnd)
840
mpfr_neg(x.__im, self.__im, rnd)
841
return x
842
843
def __pos__(self):
844
r"""
845
Method for computing the "positive" of self.
846
847
EXAMPLES::
848
849
sage: a = ComplexNumber(2,1)
850
sage: +a
851
2.00000000000000 + 1.00000000000000*I
852
sage: a.__pos__()
853
2.00000000000000 + 1.00000000000000*I
854
"""
855
return self
856
857
def __abs__(self):
858
r"""
859
Method for computing the absolute value or modulus of self
860
861
`|a + bi| = sqrt(a^2 + b^2)`
862
863
EXAMPLES: Note that the absolute value of a complex number with
864
imaginary component equal to zero is the absolute value of the real
865
component.
866
867
::
868
869
sage: a = ComplexNumber(2,1)
870
sage: abs(a)
871
2.23606797749979
872
sage: a.__abs__()
873
2.23606797749979
874
sage: float(sqrt(2^2 + 1^1))
875
2.23606797749979
876
877
::
878
879
sage: b = ComplexNumber(42,0)
880
sage: abs(b)
881
42.0000000000000
882
sage: b.__abs__()
883
42.0000000000000
884
sage: b
885
42.0000000000000
886
"""
887
return self.abs_c()
888
889
def __invert__(self):
890
"""
891
Return the multiplicative inverse.
892
893
EXAMPLES::
894
895
sage: I = CC.0
896
sage: a = ~(5+I)
897
sage: a * (5+I)
898
1.00000000000000
899
"""
900
cdef ComplexNumber x
901
x = self._new()
902
903
cdef mpfr_t t0, t1
904
mpfr_init2(t0, self._prec)
905
mpfr_init2(t1, self._prec)
906
907
mpfr_mul(t0, self.__re, self.__re, rnd)
908
mpfr_mul(t1, self.__im, self.__im, rnd)
909
910
mpfr_add(t0, t0, t1, rnd) # now t0 is the norm
911
mpfr_div(x.__re, self.__re, t0, rnd) # x.__re = self.__re/norm
912
913
mpfr_neg(t1, self.__im, rnd)
914
mpfr_div(x.__im, t1, t0, rnd) # x.__im = -self.__im/norm
915
916
mpfr_clear(t0)
917
mpfr_clear(t1)
918
919
return x
920
921
def __int__(self):
922
r"""
923
Method for converting self to type int. Called by the
924
``int`` function. Note that calling this method returns
925
an error since, in general, complex numbers cannot be coerced into
926
integers.
927
928
EXAMPLES::
929
930
sage: a = ComplexNumber(2,1)
931
sage: int(a)
932
Traceback (most recent call last):
933
...
934
TypeError: can't convert complex to int; use int(abs(z))
935
sage: a.__int__()
936
Traceback (most recent call last):
937
...
938
TypeError: can't convert complex to int; use int(abs(z))
939
"""
940
raise TypeError, "can't convert complex to int; use int(abs(z))"
941
942
def __long__(self):
943
r"""
944
Method for converting self to type long. Called by the
945
``long`` function. Note that calling this method
946
returns an error since, in general, complex numbers cannot be
947
coerced into integers.
948
949
EXAMPLES::
950
951
sage: a = ComplexNumber(2,1)
952
sage: long(a)
953
Traceback (most recent call last):
954
...
955
TypeError: can't convert complex to long; use long(abs(z))
956
sage: a.__long__()
957
Traceback (most recent call last):
958
...
959
TypeError: can't convert complex to long; use long(abs(z))
960
"""
961
raise TypeError, "can't convert complex to long; use long(abs(z))"
962
963
def __float__(self):
964
r"""
965
Method for converting self to type float. Called by the
966
``float`` function. This conversion will throw an error if
967
the number has a nonzero imaginary part.
968
969
EXAMPLES::
970
971
sage: a = ComplexNumber(1, 0)
972
sage: float(a)
973
1.0
974
sage: a = ComplexNumber(2,1)
975
sage: float(a)
976
Traceback (most recent call last):
977
...
978
TypeError: Unable to convert 2.00000000000000 + 1.00000000000000*I to float; use abs() or real_part() as desired
979
sage: a.__float__()
980
Traceback (most recent call last):
981
...
982
TypeError: Unable to convert 2.00000000000000 + 1.00000000000000*I to float; use abs() or real_part() as desired
983
sage: float(abs(ComplexNumber(1,1)))
984
1.4142135623730951
985
"""
986
if mpfr_zero_p(self.__im):
987
return mpfr_get_d(self.__re, rnd)
988
else:
989
raise TypeError, "Unable to convert %s to float; use abs() or real_part() as desired"%self
990
991
def __complex__(self):
992
r"""
993
Method for converting self to type complex. Called by the
994
``complex`` function.
995
996
EXAMPLES::
997
998
sage: a = ComplexNumber(2,1)
999
sage: complex(a)
1000
(2+1j)
1001
sage: type(complex(a))
1002
<type 'complex'>
1003
sage: a.__complex__()
1004
(2+1j)
1005
"""
1006
return complex(mpfr_get_d(self.__re, rnd),
1007
mpfr_get_d(self.__im, rnd))
1008
# return complex(float(self.__re), float(self.__im))
1009
1010
def __richcmp__(left, right, int op):
1011
return (<Element>left)._richcmp(right, op)
1012
1013
cdef int _cmp_c_impl(left, sage.structure.element.Element right) except -2:
1014
cdef int a, b
1015
a = mpfr_nan_p(left.__re)
1016
b = mpfr_nan_p((<ComplexNumber>right).__re)
1017
if a != b:
1018
return -1
1019
1020
cdef int i
1021
i = mpfr_cmp(left.__re, (<ComplexNumber>right).__re)
1022
if i < 0:
1023
return -1
1024
elif i > 0:
1025
return 1
1026
i = mpfr_cmp(left.__im, (<ComplexNumber>right).__im)
1027
if i < 0:
1028
return -1
1029
elif i > 0:
1030
return 1
1031
return 0
1032
1033
def multiplicative_order(self):
1034
"""
1035
Return the multiplicative order of this complex number, if known,
1036
or raise a NotImplementedError.
1037
1038
EXAMPLES::
1039
1040
sage: C.<i> = ComplexField()
1041
sage: i.multiplicative_order()
1042
4
1043
sage: C(1).multiplicative_order()
1044
1
1045
sage: C(-1).multiplicative_order()
1046
2
1047
sage: C(i^2).multiplicative_order()
1048
2
1049
sage: C(-i).multiplicative_order()
1050
4
1051
sage: C(2).multiplicative_order()
1052
+Infinity
1053
sage: w = (1+sqrt(-3.0))/2; w
1054
0.500000000000000 + 0.866025403784439*I
1055
sage: abs(w)
1056
1.00000000000000
1057
sage: w.multiplicative_order()
1058
Traceback (most recent call last):
1059
...
1060
NotImplementedError: order of element not known
1061
"""
1062
if self == 1:
1063
return integer.Integer(1)
1064
elif self == -1:
1065
return integer.Integer(2)
1066
elif self == self._parent.gen():
1067
return integer.Integer(4)
1068
elif self == -self._parent.gen():
1069
return integer.Integer(4)
1070
elif not self._multiplicative_order is None:
1071
return integer.Integer(self._multiplicative_order)
1072
elif abs(abs(self) - 1) > 0.1: # clearly not a root of unity
1073
return infinity.infinity
1074
raise NotImplementedError, "order of element not known"
1075
1076
1077
########################################################################
1078
# Plotting
1079
########################################################################
1080
1081
def plot(self, **kargs):
1082
"""
1083
Plots this complex number as a point in the plane
1084
1085
The accepted options are the ones of :meth:`~sage.plot.point.point2d`.
1086
Type ``point2d.options`` to see all options.
1087
1088
.. note::
1089
1090
Just wraps the sage.plot.point.point2d method
1091
1092
EXAMPLES:
1093
1094
You can either use the indirect::
1095
1096
sage: z = CC(0,1)
1097
sage: plot(z)
1098
1099
or the more direct::
1100
1101
sage: z = CC(0,1)
1102
sage: z.plot()
1103
"""
1104
return sage.plot.point.point2d((self.real(), self.imag()), **kargs)
1105
1106
########################################################################
1107
# Transcendental (and other) functions
1108
########################################################################
1109
1110
1111
1112
1113
# Trig functions
1114
def arccos(self):
1115
"""
1116
EXAMPLES::
1117
1118
sage: (1+CC(I)).arccos()
1119
0.904556894302381 - 1.06127506190504*I
1120
"""
1121
return self._parent(self._pari_().acos())
1122
1123
def arccosh(self):
1124
"""
1125
EXAMPLES::
1126
1127
sage: (1+CC(I)).arccosh()
1128
1.06127506190504 + 0.904556894302381*I
1129
"""
1130
return self._parent(self._pari_().acosh())
1131
1132
def arcsin(self):
1133
"""
1134
EXAMPLES::
1135
1136
sage: (1+CC(I)).arcsin()
1137
0.666239432492515 + 1.06127506190504*I
1138
"""
1139
return self._parent(self._pari_().asin())
1140
1141
def arcsinh(self):
1142
"""
1143
EXAMPLES::
1144
1145
sage: (1+CC(I)).arcsinh()
1146
1.06127506190504 + 0.666239432492515*I
1147
"""
1148
return self._parent(self._pari_().asinh())
1149
1150
def arctan(self):
1151
"""
1152
EXAMPLES::
1153
1154
sage: (1+CC(I)).arctan()
1155
1.01722196789785 + 0.402359478108525*I
1156
"""
1157
return self._parent(self._pari_().atan())
1158
1159
def arctanh(self):
1160
"""
1161
EXAMPLES::
1162
1163
sage: (1+CC(I)).arctanh()
1164
0.402359478108525 + 1.01722196789785*I
1165
"""
1166
return self._parent(self._pari_().atanh())
1167
1168
def coth(self):
1169
"""
1170
EXAMPLES::
1171
1172
sage: ComplexField(100)(1,1).coth()
1173
0.86801414289592494863584920892 - 0.21762156185440268136513424361*I
1174
"""
1175
return ~(self.tanh())
1176
1177
def arccoth(self):
1178
"""
1179
EXAMPLES::
1180
1181
sage: ComplexField(100)(1,1).arccoth()
1182
0.40235947810852509365018983331 - 0.55357435889704525150853273009*I
1183
"""
1184
return (~self).arctanh()
1185
1186
def csc(self):
1187
"""
1188
EXAMPLES::
1189
1190
sage: ComplexField(100)(1,1).csc()
1191
0.62151801717042842123490780586 - 0.30393100162842645033448560451*I
1192
"""
1193
return ~(self.sin())
1194
1195
def csch(self):
1196
"""
1197
EXAMPLES::
1198
1199
sage: ComplexField(100)(1,1).csch()
1200
0.30393100162842645033448560451 - 0.62151801717042842123490780586*I
1201
"""
1202
return ~(self.sinh())
1203
1204
def arccsch(self):
1205
"""
1206
EXAMPLES::
1207
1208
sage: ComplexField(100)(1,1).arccsch()
1209
0.53063753095251782601650945811 - 0.45227844715119068206365839783*I
1210
"""
1211
return (~self).arcsinh()
1212
1213
def sec(self):
1214
"""
1215
EXAMPLES::
1216
1217
sage: ComplexField(100)(1,1).sec()
1218
0.49833703055518678521380589177 + 0.59108384172104504805039169297*I
1219
"""
1220
return ~(self.cos())
1221
1222
def sech(self):
1223
"""
1224
EXAMPLES::
1225
1226
sage: ComplexField(100)(1,1).sech()
1227
0.49833703055518678521380589177 - 0.59108384172104504805039169297*I
1228
"""
1229
return ~(self.cosh())
1230
1231
def arcsech(self):
1232
"""
1233
EXAMPLES::
1234
1235
sage: ComplexField(100)(1,1).arcsech()
1236
0.53063753095251782601650945811 - 1.1185178796437059371676632938*I
1237
"""
1238
return (~self).arccosh()
1239
1240
def cotan(self):
1241
"""
1242
EXAMPLES::
1243
1244
sage: (1+CC(I)).cotan()
1245
0.217621561854403 - 0.868014142895925*I
1246
sage: i = ComplexField(200).0
1247
sage: (1+i).cotan()
1248
0.21762156185440268136513424360523807352075436916785404091068 - 0.86801414289592494863584920891627388827343874994609327121115*I
1249
sage: i = ComplexField(220).0
1250
sage: (1+i).cotan()
1251
0.21762156185440268136513424360523807352075436916785404091068124239 - 0.86801414289592494863584920891627388827343874994609327121115071646*I
1252
"""
1253
return ~(self.tan())
1254
1255
def cos(self):
1256
"""
1257
EXAMPLES::
1258
1259
sage: (1+CC(I)).cos()
1260
0.833730025131149 - 0.988897705762865*I
1261
"""
1262
# write self = a + i*b, then
1263
# cos(self) = cosh(b)*cos(a) - i*sinh(b)*sin(a)
1264
cdef ComplexNumber z
1265
z = self._new()
1266
cdef mpfr_t ch, sh
1267
mpfr_init2(sh, self._prec)
1268
mpfr_sinh(sh, self.__im, rnd)
1269
mpfr_init2(ch, self._prec)
1270
mpfr_sqr(ch, sh, rnd)
1271
mpfr_add_ui(ch, ch, 1, rnd)
1272
mpfr_sqrt(ch, ch, rnd)
1273
mpfr_neg(sh, sh, rnd)
1274
mpfr_sin_cos(z.__im, z.__re, self.__re, rnd)
1275
mpfr_mul(z.__re, z.__re, ch, rnd)
1276
mpfr_mul(z.__im, z.__im, sh, rnd)
1277
mpfr_clear(sh)
1278
mpfr_clear(ch)
1279
return z
1280
1281
def cosh(self):
1282
"""
1283
EXAMPLES::
1284
1285
sage: (1+CC(I)).cosh()
1286
0.833730025131149 + 0.988897705762865*I
1287
"""
1288
# write self = a + i*b, then
1289
# cosh(self) = cosh(a)*cos(b) + i*sinh(a)*sin(b)
1290
cdef ComplexNumber z
1291
z = self._new()
1292
cdef mpfr_t ch, sh
1293
mpfr_init2(sh, self._prec)
1294
mpfr_sinh(sh, self.__re, rnd)
1295
mpfr_init2(ch, self._prec)
1296
mpfr_sqr(ch, sh, rnd)
1297
mpfr_add_ui(ch, ch, 1, rnd)
1298
mpfr_sqrt(ch, ch, rnd)
1299
mpfr_sin_cos(z.__im, z.__re, self.__im, rnd)
1300
mpfr_mul(z.__re, z.__re, ch, rnd)
1301
mpfr_mul(z.__im, z.__im, sh, rnd)
1302
mpfr_clear(sh)
1303
mpfr_clear(ch)
1304
return z
1305
1306
1307
1308
def eta(self, omit_frac=False):
1309
r"""
1310
Return the value of the Dedekind `\eta` function on self,
1311
intelligently computed using `\mathbb{SL}(2,\ZZ)`
1312
transformations.
1313
1314
INPUT:
1315
1316
- ``self`` - element of the upper half plane (if not,
1317
raises a ValueError).
1318
1319
- ``omit_frac`` - (bool, default: False), if True,
1320
omit the `e^{\pi i z / 12}` factor.
1321
1322
1323
OUTPUT: a complex number
1324
1325
The `\eta` function is
1326
1327
.. math::
1328
1329
\eta(z) = e^{\pi i z / 12} \prod_{n=1}^{\infty}(1-e^{2\pi inz})
1330
1331
1332
1333
ALGORITHM: Uses the PARI C library.
1334
1335
EXAMPLES:
1336
1337
First we compute `\eta(1+i)`::
1338
1339
sage: i = CC.0
1340
sage: z = 1+i; z.eta()
1341
0.742048775836565 + 0.198831370229911*I
1342
1343
We compute eta to low precision directly from the definition.
1344
1345
::
1346
1347
sage: z = 1 + i; z.eta()
1348
0.742048775836565 + 0.198831370229911*I
1349
sage: pi = CC(pi) # otherwise we will get a symbolic result.
1350
sage: exp(pi * i * z / 12) * prod([1-exp(2*pi*i*n*z) for n in range(1,10)])
1351
0.742048775836565 + 0.198831370229911*I
1352
1353
The optional argument allows us to omit the fractional part::
1354
1355
sage: z = 1 + i
1356
sage: z.eta(omit_frac=True)
1357
0.998129069925959 - 8.12769318...e-22*I
1358
sage: prod([1-exp(2*pi*i*n*z) for n in range(1,10)])
1359
0.998129069925958 + 4.59099857829247e-19*I
1360
1361
We illustrate what happens when `z` is not in the upper
1362
half plane.
1363
1364
::
1365
1366
sage: z = CC(1)
1367
sage: z.eta()
1368
Traceback (most recent call last):
1369
...
1370
ValueError: value must be in the upper half plane
1371
1372
You can also use functional notation.
1373
1374
::
1375
1376
sage: eta(1+CC(I))
1377
0.742048775836565 + 0.198831370229911*I
1378
"""
1379
try:
1380
return self._parent(self._pari_().eta(not omit_frac))
1381
except sage.libs.pari.all.PariError:
1382
raise ValueError, "value must be in the upper half plane"
1383
1384
1385
def sin(self):
1386
"""
1387
EXAMPLES::
1388
1389
sage: (1+CC(I)).sin()
1390
1.29845758141598 + 0.634963914784736*I
1391
"""
1392
# write self = a + i*b, then
1393
# sin(self) = cosh(b)*sin(a) + i*sinh(b)*cos(a)
1394
cdef ComplexNumber z
1395
z = self._new()
1396
cdef mpfr_t ch, sh
1397
mpfr_init2(sh, self._prec)
1398
mpfr_sinh(sh, self.__im, rnd)
1399
mpfr_init2(ch, self._prec)
1400
mpfr_sqr(ch, sh, rnd)
1401
mpfr_add_ui(ch, ch, 1, rnd)
1402
mpfr_sqrt(ch, ch, rnd)
1403
mpfr_sin_cos(z.__re, z.__im, self.__re, rnd)
1404
mpfr_mul(z.__re, z.__re, ch, rnd)
1405
mpfr_mul(z.__im, z.__im, sh, rnd)
1406
mpfr_clear(sh)
1407
mpfr_clear(ch)
1408
return z
1409
1410
def sinh(self):
1411
"""
1412
EXAMPLES::
1413
1414
sage: (1+CC(I)).sinh()
1415
0.634963914784736 + 1.29845758141598*I
1416
"""
1417
# write self = a + i*b, then
1418
# sinh(self) = sinh(a)*cos(b) + i*cosh(a)*sin(b)
1419
cdef ComplexNumber z
1420
z = self._new()
1421
cdef mpfr_t ch, sh
1422
mpfr_init2(sh, self._prec)
1423
mpfr_sinh(sh, self.__re, rnd)
1424
mpfr_init2(ch, self._prec)
1425
mpfr_sqr(ch, sh, rnd)
1426
mpfr_add_ui(ch, ch, 1, rnd)
1427
mpfr_sqrt(ch, ch, rnd)
1428
mpfr_sin_cos(z.__im, z.__re, self.__im, rnd)
1429
mpfr_mul(z.__re, z.__re, sh, rnd)
1430
mpfr_mul(z.__im, z.__im, ch, rnd)
1431
mpfr_clear(sh)
1432
mpfr_clear(ch)
1433
return z
1434
1435
def tan(self):
1436
"""
1437
EXAMPLES::
1438
1439
sage: (1+CC(I)).tan()
1440
0.271752585319512 + 1.08392332733869*I
1441
"""
1442
# write self = a + i*b, then
1443
# tan(self) = [cos(a)*sin(a) + i*cosh(b)*sinh(b)]/[sinh^2(b)+cos^2(a)]
1444
cdef ComplexNumber z
1445
z = self._new()
1446
cdef mpfr_t ch, sh, c, s, a, b
1447
mpfr_init2(sh, self._prec)
1448
mpfr_sinh(sh, self.__im, rnd)
1449
mpfr_init2(ch, self._prec)
1450
mpfr_init2(a, self._prec)
1451
mpfr_sqr(a, sh, rnd)
1452
mpfr_add_ui(ch, a, 1, rnd)
1453
mpfr_sqrt(ch, ch, rnd)
1454
mpfr_init2(c, self._prec)
1455
mpfr_init2(s, self._prec)
1456
mpfr_sin_cos(s, c, self.__re, rnd)
1457
mpfr_init2(b, self._prec)
1458
mpfr_sqr(b, c, rnd)
1459
mpfr_add(a, a, b, rnd)
1460
mpfr_mul(z.__re, c, s, rnd)
1461
mpfr_div(z.__re, z.__re, a, rnd)
1462
mpfr_mul(z.__im, ch, sh, rnd)
1463
mpfr_div(z.__im, z.__im, a, rnd)
1464
mpfr_clear(sh)
1465
mpfr_clear(ch)
1466
mpfr_clear(c)
1467
mpfr_clear(s)
1468
mpfr_clear(b)
1469
mpfr_clear(a)
1470
return z
1471
1472
1473
def tanh(self):
1474
"""
1475
EXAMPLES::
1476
1477
sage: (1+CC(I)).tanh()
1478
1.08392332733869 + 0.271752585319512*I
1479
"""
1480
# write self = a + i*b, then
1481
# tanh(self) = [cosh(a)*sinh(a) + i*cos(b)*sin(b)]/[sinh^2(a)+cos^2(b)]
1482
cdef ComplexNumber z
1483
z = self._new()
1484
cdef mpfr_t ch, sh, c, s, a, b
1485
mpfr_init2(sh, self._prec)
1486
mpfr_sinh(sh, self.__re, rnd)
1487
mpfr_init2(ch, self._prec)
1488
mpfr_init2(a, self._prec)
1489
mpfr_sqr(a, sh, rnd)
1490
mpfr_add_ui(ch, a, 1, rnd)
1491
mpfr_sqrt(ch, ch, rnd)
1492
mpfr_init2(c, self._prec)
1493
mpfr_init2(s, self._prec)
1494
mpfr_sin_cos(s, c, self.__im, rnd)
1495
mpfr_init2(b, self._prec)
1496
mpfr_sqr(b, c, rnd)
1497
mpfr_add(a, a, b, rnd)
1498
mpfr_mul(z.__im, c, s, rnd)
1499
mpfr_div(z.__im, z.__im, a, rnd)
1500
mpfr_mul(z.__re, ch, sh, rnd)
1501
mpfr_div(z.__re, z.__re, a, rnd)
1502
mpfr_clear(sh)
1503
mpfr_clear(ch)
1504
mpfr_clear(c)
1505
mpfr_clear(s)
1506
mpfr_clear(b)
1507
mpfr_clear(a)
1508
return z
1509
1510
# Other special functions
1511
def agm(self, right, algorithm="optimal"):
1512
"""
1513
Return the Arithmetic-Geometric Mean (AGM) of self and right.
1514
1515
INPUT:
1516
1517
- right (complex)-- another complex number
1518
1519
- algorithm (string, default "optimal")-- the algorithm to use
1520
(see below).
1521
1522
OUTPUT:
1523
1524
(complex) A value of the AGM of self and right. Note that
1525
this is a multi-valued function, and the algorithm used
1526
affects the value returned, as follows:
1527
1528
- "pari": Call the sgm function from the pari library.
1529
1530
- "optimal": Use the AGM sequence such that at each stage
1531
`(a,b)` is replaced by `(a_1,b_1)=((a+b)/2,\pm\sqrt{ab})`
1532
where the sign is chosen so that `|a_1-b_1|\le|a_1+b_1|`, or
1533
equivalently `\Re(b_1/a_1)\ge 0`. The resulting limit is
1534
maximal among all possible values.
1535
1536
- "principal": Use the AGM sequence such that at each stage
1537
`(a,b)` is replaced by `(a_1,b_1)=((a+b)/2,\pm\sqrt{ab})`
1538
where the sign is chosen so that `\Re(b_1)\ge 0` (the
1539
so-called principal branch of the square root).
1540
1541
The values `AGM(a,0)`, `AGM(0,a)`, and `AGM(a,-a)` are all taken to be 0.
1542
1543
EXAMPLES::
1544
1545
sage: a = CC(1,1)
1546
sage: b = CC(2,-1)
1547
sage: a.agm(b)
1548
1.62780548487271 + 0.136827548397369*I
1549
sage: a.agm(b, algorithm="optimal")
1550
1.62780548487271 + 0.136827548397369*I
1551
sage: a.agm(b, algorithm="principal")
1552
1.62780548487271 + 0.136827548397369*I
1553
sage: a.agm(b, algorithm="pari")
1554
1.62780548487271 + 0.136827548397369*I
1555
1556
An example to show that the returned value depends on the algorithm parameter::
1557
1558
sage: a = CC(-0.95,-0.65)
1559
sage: b = CC(0.683,0.747)
1560
sage: a.agm(b, algorithm="optimal")
1561
-0.371591652351761 + 0.319894660206830*I
1562
sage: a.agm(b, algorithm="principal")
1563
0.338175462986180 - 0.0135326969565405*I
1564
sage: a.agm(b, algorithm="pari")
1565
0.0806891850759812 + 0.239036532685557*I
1566
sage: a.agm(b, algorithm="optimal").abs()
1567
0.490319232466314
1568
sage: a.agm(b, algorithm="principal").abs()
1569
0.338446122230459
1570
sage: a.agm(b, algorithm="pari").abs()
1571
0.252287947683910
1572
1573
TESTS:
1574
1575
An example which came up in testing::
1576
1577
sage: I = CC(I)
1578
sage: a = 0.501648970493109 + 1.11877240294744*I
1579
sage: b = 1.05946309435930 + 1.05946309435930*I
1580
sage: a.agm(b)
1581
0.774901870587681 + 1.10254945079875*I
1582
1583
sage: a = CC(-0.32599972608379413, 0.60395514542928641)
1584
sage: b = CC( 0.6062314525690593, 0.1425693337776659)
1585
sage: a.agm(b)
1586
0.199246281325876 + 0.478401702759654*I
1587
sage: a.agm(-a)
1588
0.000000000000000
1589
sage: a.agm(0)
1590
0.000000000000000
1591
sage: CC(0).agm(a)
1592
0.000000000000000
1593
1594
Consistency::
1595
1596
sage: a = 1 + 0.5*I
1597
sage: b = 2 - 0.25*I
1598
sage: a.agm(b) - ComplexField(100)(a).agm(b)
1599
0.000000000000000
1600
sage: ComplexField(200)(a).agm(b) - ComplexField(500)(a).agm(b)
1601
0.00000000000000000000000000000000000000000000000000000000000
1602
sage: ComplexField(500)(a).agm(b) - ComplexField(1000)(a).agm(b)
1603
0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
1604
"""
1605
if algorithm=="pari":
1606
t = self._parent(right)._pari_()
1607
return self._parent(self._pari_().agm(t))
1608
1609
cdef ComplexNumber a, b, a1, b1, d, e, res
1610
cdef mp_exp_t rel_prec
1611
cdef bint optimal = algorithm == "optimal"
1612
1613
if optimal or algorithm == "principal":
1614
1615
if not isinstance(right, ComplexNumber) or (<ComplexNumber>right)._parent is not self._parent:
1616
right = self._parent(right)
1617
1618
res = self._new()
1619
1620
if mpfr_zero_p(self.__re) and mpfr_zero_p(self.__im):
1621
return self
1622
elif mpfr_zero_p((<ComplexNumber>right).__re) and mpfr_zero_p((<ComplexNumber>right).__im):
1623
return right
1624
elif (mpfr_cmpabs(self.__re, (<ComplexNumber>right).__re) == 0 and
1625
mpfr_cmpabs(self.__im, (<ComplexNumber>right).__im) == 0 and
1626
mpfr_cmp(self.__re, (<ComplexNumber>right).__re) != 0 and
1627
mpfr_cmp(self.__im, (<ComplexNumber>right).__im) != 0):
1628
# self = -right
1629
mpfr_set_ui(res.__re, 0, rnd)
1630
mpfr_set_ui(res.__im, 0, rnd)
1631
return res
1632
1633
# Do the computations to a bit higher precicion so rounding error
1634
# won't obscure the termination condition.
1635
a = ComplexNumber(self._parent.to_prec(self._prec+5), None)
1636
b = a._new()
1637
a1 = a._new()
1638
b1 = a._new()
1639
1640
d = a._new()
1641
if optimal:
1642
e = a._new()
1643
1644
# Make copies so we don't mutate self or right.
1645
mpfr_set(a.__re, self.__re, rnd)
1646
mpfr_set(a.__im, self.__im, rnd)
1647
mpfr_set(b.__re, (<ComplexNumber>right).__re, rnd)
1648
mpfr_set(b.__im, (<ComplexNumber>right).__im, rnd)
1649
1650
if optimal:
1651
mpfr_add(e.__re, a.__re, b.__re, rnd)
1652
mpfr_add(e.__im, a.__im, b.__im, rnd)
1653
1654
while True:
1655
1656
# a1 = (a+b)/2
1657
if optimal:
1658
mpfr_swap(a1.__re, e.__re)
1659
mpfr_swap(a1.__im, e.__im)
1660
else:
1661
mpfr_add(a1.__re, a.__re, b.__re, rnd)
1662
mpfr_add(a1.__im, a.__im, b.__im, rnd)
1663
mpfr_mul_2si(a1.__re, a1.__re, -1, rnd)
1664
mpfr_mul_2si(a1.__im, a1.__im, -1, rnd)
1665
1666
# b1 = sqrt(a*b)
1667
mpfr_mul(d.__re, a.__re, b.__re, rnd)
1668
mpfr_mul(d.__im, a.__im, b.__im, rnd)
1669
mpfr_sub(b1.__re, d.__re, d.__im, rnd)
1670
mpfr_mul(d.__re, a.__re, b.__im, rnd)
1671
mpfr_mul(d.__im, a.__im, b.__re, rnd)
1672
mpfr_add(b1.__im, d.__re, d.__im, rnd)
1673
b1 = b1.sqrt() # this would be a *lot* of code duplication
1674
1675
# d = a1 - b1
1676
mpfr_sub(d.__re, a1.__re, b1.__re, rnd)
1677
mpfr_sub(d.__im, a1.__im, b1.__im, rnd)
1678
if mpfr_zero_p(d.__re) and mpfr_zero_p(d.__im):
1679
mpfr_set(res.__re, a1.__re, rnd)
1680
mpfr_set(res.__im, a1.__im, rnd)
1681
return res
1682
1683
if optimal:
1684
# e = a1+b1
1685
mpfr_add(e.__re, a1.__re, b1.__re, rnd)
1686
mpfr_add(e.__im, a1.__im, b1.__im, rnd)
1687
if mpfr_zero_p(e.__re) and mpfr_zero_p(e.__im):
1688
mpfr_set(res.__re, a1.__re, rnd)
1689
mpfr_set(res.__im, a1.__im, rnd)
1690
return res
1691
1692
# |e| < |d|
1693
if cmp_abs(e, d) < 0:
1694
mpfr_swap(d.__re, e.__re)
1695
mpfr_swap(d.__im, e.__im)
1696
mpfr_neg(b1.__re, b1.__re, rnd)
1697
mpfr_neg(b1.__im, b1.__im, rnd)
1698
1699
rel_prec = min_exp_t(max_exp(a1), max_exp(b1)) - max_exp(d)
1700
if rel_prec > self._prec:
1701
mpfr_set(res.__re, a1.__re, rnd)
1702
mpfr_set(res.__im, a1.__im, rnd)
1703
return res
1704
1705
# a, b = a1, b1
1706
mpfr_swap(a.__re, a1.__re)
1707
mpfr_swap(a.__im, a1.__im)
1708
mpfr_swap(b.__re, b1.__re)
1709
mpfr_swap(b.__im, b1.__im)
1710
1711
raise ValueError, "agm algorithm must be one of 'pari', 'optimal', 'principal'"
1712
1713
def argument(self):
1714
r"""
1715
The argument (angle) of the complex number, normalized so that
1716
`-\pi < \theta \leq \pi`.
1717
1718
EXAMPLES::
1719
1720
sage: i = CC.0
1721
sage: (i^2).argument()
1722
3.14159265358979
1723
sage: (1+i).argument()
1724
0.785398163397448
1725
sage: i.argument()
1726
1.57079632679490
1727
sage: (-i).argument()
1728
-1.57079632679490
1729
sage: (RR('-0.001') - i).argument()
1730
-1.57179632646156
1731
"""
1732
cdef real_mpfr.RealNumber x
1733
x = real_mpfr.RealNumber(self._parent._real_field(), None)
1734
mpfr_atan2(x.value, self.__im, self.__re, rnd)
1735
return x
1736
1737
1738
def arg(self):
1739
"""
1740
Same as argument.
1741
1742
EXAMPLES::
1743
1744
sage: i = CC.0
1745
sage: (i^2).arg()
1746
3.14159265358979
1747
"""
1748
return self.argument()
1749
1750
def conjugate(self):
1751
"""
1752
Return the complex conjugate of this complex number.
1753
1754
EXAMPLES::
1755
1756
sage: i = CC.0
1757
sage: (1+i).conjugate()
1758
1.00000000000000 - 1.00000000000000*I
1759
"""
1760
cdef ComplexNumber x
1761
x = self._new()
1762
1763
cdef mpfr_t i
1764
mpfr_init2(i, self._prec)
1765
mpfr_neg(i, self.__im, rnd)
1766
mpfr_set(x.__re, self.__re, rnd)
1767
mpfr_set(x.__im, i, rnd)
1768
mpfr_clear(i)
1769
return x
1770
1771
def dilog(self):
1772
r"""
1773
Returns the complex dilogarithm of self. The complex dilogarithm,
1774
or Spence's function, is defined by
1775
1776
`Li_2(z) = - \int_0^z \frac{\log|1-\zeta|}{\zeta} d(\zeta)`
1777
1778
`= \sum_{k=1}^\infty \frac{z^k}{k}`
1779
1780
Note that the series definition can only be used for
1781
`|z| < 1`
1782
1783
EXAMPLES::
1784
1785
sage: a = ComplexNumber(1,0)
1786
sage: a.dilog()
1787
1.64493406684823
1788
sage: float(pi^2/6)
1789
1.6449340668482262
1790
1791
::
1792
1793
sage: b = ComplexNumber(0,1)
1794
sage: b.dilog()
1795
-0.205616758356028 + 0.915965594177219*I
1796
1797
::
1798
1799
sage: c = ComplexNumber(0,0)
1800
sage: c.dilog()
1801
0.000000000000000
1802
"""
1803
return self._parent(self._pari_().dilog())
1804
1805
def exp(ComplexNumber self):
1806
"""
1807
Compute exp(z).
1808
1809
EXAMPLES::
1810
1811
sage: i = ComplexField(300).0
1812
sage: z = 1 + i
1813
sage: z.exp()
1814
1.46869393991588515713896759732660426132695673662900872279767567631093696585951213872272450 + 2.28735528717884239120817190670050180895558625666835568093865811410364716018934540926734485*I
1815
"""
1816
# write self = a + i*b, then
1817
# exp(self) = exp(a)*(cos(b) + i*sin(b))
1818
cdef ComplexNumber z
1819
z = self._new()
1820
cdef mpfr_t r
1821
mpfr_init2(r, self._prec)
1822
mpfr_exp(r, self.__re, rnd)
1823
mpfr_sin_cos(z.__im, z.__re, self.__im, rnd)
1824
mpfr_mul(z.__re, z.__re, r, rnd)
1825
mpfr_mul(z.__im, z.__im, r, rnd)
1826
mpfr_clear(r)
1827
return z
1828
1829
def gamma(self):
1830
"""
1831
Return the Gamma function evaluated at this complex number.
1832
1833
EXAMPLES::
1834
1835
sage: i = ComplexField(30).0
1836
sage: (1+i).gamma()
1837
0.49801567 - 0.15494983*I
1838
1839
TESTS::
1840
1841
sage: CC(0).gamma()
1842
Infinity
1843
1844
::
1845
1846
sage: CC(-1).gamma()
1847
Infinity
1848
"""
1849
try:
1850
return self._parent(self._pari_().gamma())
1851
except sage.libs.pari.all.PariError:
1852
from sage.rings.infinity import UnsignedInfinityRing
1853
return UnsignedInfinityRing.gen()
1854
1855
def gamma_inc(self, t):
1856
"""
1857
Return the incomplete Gamma function evaluated at this complex
1858
number.
1859
1860
EXAMPLES::
1861
1862
sage: C, i = ComplexField(30).objgen()
1863
sage: (1+i).gamma_inc(2 + 3*i)
1864
0.0020969149 - 0.059981914*I
1865
sage: (1+i).gamma_inc(5)
1866
-0.0013781309 + 0.0065198200*I
1867
sage: C(2).gamma_inc(1 + i)
1868
0.70709210 - 0.42035364*I
1869
sage: CC(2).gamma_inc(5)
1870
0.0404276819945128
1871
"""
1872
return self._parent(self._pari_().incgam(t))
1873
1874
def log(self,base=None):
1875
r"""
1876
Complex logarithm of z with branch chosen as follows: Write
1877
`z = \rho e^{i \theta} with `-\pi < \theta <= pi`. Then
1878
`\mathrm{log}(z) = \mathrm{log}(\rho) + i \theta`.
1879
1880
.. warning::
1881
1882
Currently the real log is computed using floats, so there
1883
is potential precision loss.
1884
1885
EXAMPLES::
1886
1887
sage: a = ComplexNumber(2,1)
1888
sage: a.log()
1889
0.804718956217050 + 0.463647609000806*I
1890
sage: log(a.abs())
1891
0.804718956217050
1892
sage: a.argument()
1893
0.463647609000806
1894
1895
::
1896
1897
sage: b = ComplexNumber(float(exp(42)),0)
1898
sage: b.log()
1899
41.99999999999971
1900
1901
::
1902
1903
sage: c = ComplexNumber(-1,0)
1904
sage: c.log()
1905
3.14159265358979*I
1906
1907
The option of a base is included for compatibility with other logs::
1908
1909
sage: c = ComplexNumber(-1,0)
1910
sage: c.log(2)
1911
4.53236014182719*I
1912
"""
1913
theta = self.argument()
1914
rho = abs(self)
1915
if base is None:
1916
return ComplexNumber(self._parent, rho.log(), theta)
1917
else:
1918
from real_mpfr import RealField
1919
return ComplexNumber(self._parent, rho.log()/RealNumber(RealField(self.prec()),base).log(), theta/RealNumber(RealField(self.prec()),base).log())
1920
1921
def additive_order(self):
1922
"""
1923
EXAMPLES::
1924
1925
sage: CC(0).additive_order()
1926
1
1927
sage: CC.gen().additive_order()
1928
+Infinity
1929
"""
1930
if self == 0:
1931
return 1
1932
else:
1933
return infinity.infinity
1934
1935
def sqrt(self, all=False):
1936
"""
1937
The square root function, taking the branch cut to be the negative
1938
real axis.
1939
1940
INPUT:
1941
1942
1943
- ``all`` - bool (default: False); if True, return a
1944
list of all square roots.
1945
1946
1947
EXAMPLES::
1948
1949
sage: C.<i> = ComplexField(30)
1950
sage: i.sqrt()
1951
0.70710678 + 0.70710678*I
1952
sage: (1+i).sqrt()
1953
1.0986841 + 0.45508986*I
1954
sage: (C(-1)).sqrt()
1955
1.0000000*I
1956
sage: (1 + 1e-100*i).sqrt()^2
1957
1.0000000 + 1.0000000e-100*I
1958
sage: i = ComplexField(200).0
1959
sage: i.sqrt()
1960
0.70710678118654752440084436210484903928483593768847403658834 + 0.70710678118654752440084436210484903928483593768847403658834*I
1961
"""
1962
cdef ComplexNumber z = self._new()
1963
if mpfr_zero_p(self.__im):
1964
if mpfr_sgn(self.__re) >= 0:
1965
mpfr_set_ui(z.__im, 0, rnd)
1966
mpfr_sqrt(z.__re, self.__re, rnd)
1967
else:
1968
mpfr_set_ui(z.__re, 0, rnd)
1969
mpfr_neg(z.__im, self.__re, rnd)
1970
mpfr_sqrt(z.__im, z.__im, rnd)
1971
if all:
1972
return [z, -z] if z else [z]
1973
else:
1974
return z
1975
# self = x + yi = (a+bi)^2
1976
# expand, substitute, solve
1977
# a^2 = (x + sqrt(x^2+y^2))/2
1978
cdef bint avoid_branch = mpfr_sgn(self.__re) < 0 and mpfr_cmpabs(self.__im, self.__re) < 0
1979
cdef mpfr_t a2
1980
mpfr_init2(a2, self._prec)
1981
mpfr_hypot(a2, self.__re, self.__im, rnd)
1982
if avoid_branch:
1983
# x + sqrt(x^2+y^2) numerically unstable for x near negative real axis
1984
# so we compute sqrt of (-z) and shift by i at the end
1985
mpfr_sub(a2, a2, self.__re, rnd)
1986
else:
1987
mpfr_add(a2, a2, self.__re, rnd)
1988
mpfr_mul_2si(a2, a2, -1, rnd)
1989
# a = sqrt(a2)
1990
mpfr_sqrt(z.__re, a2, rnd)
1991
# b = y/(2a)
1992
mpfr_div(z.__im, self.__im, z.__re, rnd)
1993
mpfr_mul_2si(z.__im, z.__im, -1, rnd)
1994
mpfr_clear(a2)
1995
if avoid_branch:
1996
mpfr_swap(z.__re, z.__im)
1997
# Note that y (hence b) was never negated, so we have z=i*sqrt(self).
1998
# if we were below the branch cut, we want the other branch
1999
if mpfr_sgn(self.__im) < 0:
2000
mpfr_neg(z.__re, z.__re, rnd)
2001
mpfr_neg(z.__im, z.__im, rnd)
2002
if all:
2003
return [z, -z]
2004
else:
2005
return z
2006
2007
def nth_root(self, n, all=False):
2008
"""
2009
The n-th root function.
2010
2011
INPUT:
2012
2013
2014
- ``all`` - bool (default: False); if True, return a
2015
list of all n-th roots.
2016
2017
2018
EXAMPLES::
2019
2020
sage: a = CC(27)
2021
sage: a.nth_root(3)
2022
3.00000000000000
2023
sage: a.nth_root(3, all=True)
2024
[3.00000000000000, -1.50000000000000 + 2.59807621135332*I, -1.50000000000000 - 2.59807621135332*I]
2025
sage: a = ComplexField(20)(2,1)
2026
sage: [r^7 for r in a.nth_root(7, all=True)]
2027
[2.0000 + 1.0000*I, 2.0000 + 1.0000*I, 2.0000 + 1.0000*I, 2.0000 + 1.0000*I, 2.0000 + 1.0000*I, 2.0000 + 1.0001*I, 2.0000 + 1.0001*I]
2028
"""
2029
if self.is_zero():
2030
return [self] if all else self
2031
2032
cdef ComplexNumber z
2033
z = self._new()
2034
2035
cdef real_mpfr.RealNumber arg, rho
2036
cdef mpfr_t r
2037
rho = abs(self)
2038
arg = self.argument() / n
2039
mpfr_init2(r, self._prec)
2040
mpfr_root(r, rho.value, n, rnd)
2041
2042
mpfr_sin_cos(z.__im, z.__re, arg.value, rnd)
2043
mpfr_mul(z.__re, z.__re, r, rnd)
2044
mpfr_mul(z.__im, z.__im, r, rnd)
2045
2046
if not all:
2047
mpfr_clear(r)
2048
return z
2049
2050
R = self._parent._real_field()
2051
cdef real_mpfr.RealNumber theta
2052
theta = R.pi()*2/n
2053
zlist = [z]
2054
for k in range(1, n):
2055
z = self._new()
2056
arg += theta
2057
mpfr_sin_cos(z.__im, z.__re, arg.value, rnd)
2058
mpfr_mul(z.__re, z.__re, r, rnd)
2059
mpfr_mul(z.__im, z.__im, r, rnd)
2060
zlist.append(z)
2061
2062
mpfr_clear(r)
2063
return zlist
2064
2065
2066
def is_square(self):
2067
"""
2068
This function always returns true as `\CC` is
2069
algebraically closed.
2070
2071
EXAMPLES::
2072
2073
sage: a = ComplexNumber(2,1)
2074
sage: a.is_square()
2075
True
2076
2077
`\CC` is algebraically closed, hence every element
2078
is a square::
2079
2080
sage: b = ComplexNumber(5)
2081
sage: b.is_square()
2082
True
2083
"""
2084
return True
2085
2086
def is_real(self):
2087
"""
2088
Return True if self is real, i.e. has imaginary part zero.
2089
2090
EXAMPLES::
2091
2092
sage: CC(1.23).is_real()
2093
True
2094
sage: CC(1+i).is_real()
2095
False
2096
"""
2097
return (mpfr_zero_p(self.__im) <> 0)
2098
2099
def is_imaginary(self):
2100
"""
2101
Return True if self is imaginary, i.e. has real part zero.
2102
2103
EXAMPLES::
2104
2105
sage: CC(1.23*i).is_imaginary()
2106
True
2107
sage: CC(1+i).is_imaginary()
2108
False
2109
"""
2110
return (mpfr_zero_p(self.__re) <> 0)
2111
2112
def zeta(self):
2113
"""
2114
Return the Riemann zeta function evaluated at this complex number.
2115
2116
EXAMPLES::
2117
2118
sage: i = ComplexField(30).gen()
2119
sage: z = 1 + i
2120
sage: z.zeta()
2121
0.58215806 - 0.92684856*I
2122
sage: zeta(z)
2123
0.58215806 - 0.92684856*I
2124
2125
sage: CC(1).zeta()
2126
Infinity
2127
"""
2128
if mpfr_zero_p(self.__im) and mpfr_cmp_ui(self.__re, 1) == 0:
2129
import infinity
2130
return infinity.unsigned_infinity
2131
return self._parent(self._pari_().zeta())
2132
2133
def algdep(self, n, **kwds):
2134
"""
2135
Returns a polynomial of degree at most `n` which is
2136
approximately satisfied by this complex number. Note that the
2137
returned polynomial need not be irreducible, and indeed usually
2138
won't be if `z` is a good approximation to an algebraic
2139
number of degree less than `n`.
2140
2141
ALGORITHM: Uses the PARI C-library algdep command.
2142
2143
INPUT: Type algdep? at the top level prompt. All additional
2144
parameters are passed onto the top-level algdep command.
2145
2146
EXAMPLE::
2147
2148
sage: C = ComplexField()
2149
sage: z = (1/2)*(1 + sqrt(3.0) *C.0); z
2150
0.500000000000000 + 0.866025403784439*I
2151
sage: p = z.algdep(5); p
2152
x^3 + 1
2153
sage: p.factor()
2154
(x + 1) * (x^2 - x + 1)
2155
sage: z^2 - z + 1
2156
1.11022302462516e-16
2157
"""
2158
import sage.rings.arith
2159
return sage.rings.arith.algdep(self,n, **kwds)
2160
2161
def algebraic_dependancy( self, n ):
2162
"""
2163
Returns a polynomial of degree at most `n` which is
2164
approximately satisfied by this complex number. Note that the
2165
returned polynomial need not be irreducible, and indeed usually
2166
won't be if `z` is a good approximation to an algebraic
2167
number of degree less than `n`.
2168
2169
ALGORITHM: Uses the PARI C-library algdep command.
2170
2171
INPUT: Type algdep? at the top level prompt. All additional
2172
parameters are passed onto the top-level algdep command.
2173
2174
EXAMPLE::
2175
2176
sage: C = ComplexField()
2177
sage: z = (1/2)*(1 + sqrt(3.0) *C.0); z
2178
0.500000000000000 + 0.866025403784439*I
2179
sage: p = z.algebraic_dependancy(5); p
2180
x^3 + 1
2181
sage: p.factor()
2182
(x + 1) * (x^2 - x + 1)
2183
sage: z^2 - z + 1
2184
1.11022302462516e-16
2185
"""
2186
return self.algdep( n )
2187
2188
def make_ComplexNumber0( fld, mult_order, re, im ):
2189
x = ComplexNumber( fld, re, im )
2190
x._set_multiplicative_order( mult_order )
2191
return x
2192
2193
2194
2195
def create_ComplexNumber(s_real, s_imag=None, int pad=0, min_prec=53):
2196
r"""
2197
Return the complex number defined by the strings s_real and
2198
s_imag as an element of ``ComplexField(prec=n)``,
2199
where n potentially has slightly more (controlled by pad) bits than
2200
given by s.
2201
2202
INPUT:
2203
2204
- ``s_real`` - a string that defines a real number
2205
(or something whose string representation defines a number)
2206
2207
- ``s_imag`` - a string that defines a real number
2208
(or something whose string representation defines a number)
2209
2210
- ``pad`` - an integer >= 0.
2211
2212
- ``min_prec`` - number will have at least this many
2213
bits of precision, no matter what.
2214
2215
2216
EXAMPLES::
2217
2218
sage: ComplexNumber('2.3')
2219
2.30000000000000
2220
sage: ComplexNumber('2.3','1.1')
2221
2.30000000000000 + 1.10000000000000*I
2222
sage: ComplexNumber(10)
2223
10.0000000000000
2224
sage: ComplexNumber(10,10)
2225
10.0000000000000 + 10.0000000000000*I
2226
sage: ComplexNumber(1.000000000000000000000000000,2)
2227
1.00000000000000000000000000 + 2.00000000000000000000000000*I
2228
sage: ComplexNumber(1,2.000000000000000000000)
2229
1.00000000000000000000 + 2.00000000000000000000*I
2230
2231
::
2232
2233
sage: sage.rings.complex_number.create_ComplexNumber(s_real=2,s_imag=1)
2234
2.00000000000000 + 1.00000000000000*I
2235
2236
TESTS:
2237
2238
Make sure we've rounded up log(10,2) enough to guarantee
2239
sufficient precision (trac #10164)::
2240
2241
sage: s = "1." + "0"*10**6 + "1"
2242
sage: sage.rings.complex_number.create_ComplexNumber(s,0).real()-1 == 0
2243
False
2244
sage: sage.rings.complex_number.create_ComplexNumber(0,s).imag()-1 == 0
2245
False
2246
2247
"""
2248
if s_imag is None:
2249
s_imag = 0
2250
2251
if not isinstance(s_real, str):
2252
s_real = str(s_real).strip()
2253
if not isinstance(s_imag, str):
2254
s_imag = str(s_imag).strip()
2255
#if base == 10:
2256
bits = max(int(LOG_TEN_TWO_PLUS_EPSILON*len(s_real)),
2257
int(LOG_TEN_TWO_PLUS_EPSILON*len(s_imag)))
2258
#else:
2259
# bits = max(int(math.log(base,2)*len(s_imag)),int(math.log(base,2)*len(s_imag)))
2260
2261
C = complex_field.ComplexField(prec=max(bits+pad, min_prec))
2262
2263
return ComplexNumber(C, s_real, s_imag)
2264
2265
2266
cdef class RRtoCC(Map):
2267
2268
cdef ComplexNumber _zero
2269
2270
def __init__(self, RR, CC):
2271
"""
2272
EXAMPLES::
2273
2274
sage: from sage.rings.complex_number import RRtoCC
2275
sage: RRtoCC(RR, CC)
2276
Natural map:
2277
From: Real Field with 53 bits of precision
2278
To: Complex Field with 53 bits of precision
2279
"""
2280
Map.__init__(self, RR, CC)
2281
self._zero = ComplexNumber(CC, 0)
2282
self._repr_type_str = "Natural"
2283
2284
cpdef Element _call_(self, x):
2285
"""
2286
EXAMPLES::
2287
2288
sage: from sage.rings.complex_number import RRtoCC
2289
sage: f = RRtoCC(RealField(100), ComplexField(10))
2290
sage: f(1/3)
2291
0.33
2292
"""
2293
cdef ComplexNumber z = self._zero._new()
2294
mpfr_set(z.__re, (<RealNumber>x).value, rnd)
2295
mpfr_set_ui(z.__im, 0, rnd)
2296
return z
2297
2298
2299
cdef class CCtoCDF(Map):
2300
2301
cpdef Element _call_(self, x):
2302
"""
2303
EXAMPLES::
2304
2305
sage: from sage.rings.complex_number import CCtoCDF
2306
sage: f = CCtoCDF(CC, CDF)
2307
sage: f(CC.0)
2308
1.0*I
2309
sage: f(exp(pi*CC.0/4))
2310
0.707106781187 + 0.707106781187*I
2311
"""
2312
cdef ComplexDoubleElement z = <ComplexDoubleElement>PY_NEW(ComplexDoubleElement)
2313
z._complex.dat[0] = mpfr_get_d((<ComplexNumber>x).__re, GMP_RNDN)
2314
z._complex.dat[1] = mpfr_get_d((<ComplexNumber>x).__im, GMP_RNDN)
2315
return z
2316
2317
2318
cdef inline mp_exp_t min_exp_t(mp_exp_t a, mp_exp_t b):
2319
return a if a < b else b
2320
2321
cdef inline mp_exp_t max_exp_t(mp_exp_t a, mp_exp_t b):
2322
return a if a > b else b
2323
2324
cdef inline mp_exp_t max_exp(ComplexNumber z):
2325
"""
2326
Quickly return the maximum exponent of the real and complex parts of z,
2327
which is useful for estimating its magnitude.
2328
"""
2329
if mpfr_zero_p(z.__im):
2330
return mpfr_get_exp(z.__re)
2331
elif mpfr_zero_p(z.__re):
2332
return mpfr_get_exp(z.__im)
2333
return max_exp_t(mpfr_get_exp(z.__re), mpfr_get_exp(z.__im))
2334
2335
cpdef int cmp_abs(ComplexNumber a, ComplexNumber b):
2336
"""
2337
Returns -1, 0, or 1 according to whether `|a|` is less than, equal to, or greater than `|b|`.
2338
2339
Optimized for non-close numbers, where the ordering can be determined by examining exponents.
2340
2341
EXAMPLES::
2342
2343
sage: from sage.rings.complex_number import cmp_abs
2344
sage: cmp_abs(CC(5), CC(1))
2345
1
2346
sage: cmp_abs(CC(5), CC(4))
2347
1
2348
sage: cmp_abs(CC(5), CC(5))
2349
0
2350
sage: cmp_abs(CC(5), CC(6))
2351
-1
2352
sage: cmp_abs(CC(5), CC(100))
2353
-1
2354
sage: cmp_abs(CC(-100), CC(1))
2355
1
2356
sage: cmp_abs(CC(-100), CC(100))
2357
0
2358
sage: cmp_abs(CC(-100), CC(1000))
2359
-1
2360
sage: cmp_abs(CC(1,1), CC(1))
2361
1
2362
sage: cmp_abs(CC(1,1), CC(2))
2363
-1
2364
sage: cmp_abs(CC(1,1), CC(1,0.99999))
2365
1
2366
sage: cmp_abs(CC(1,1), CC(1,-1))
2367
0
2368
sage: cmp_abs(CC(0), CC(1))
2369
-1
2370
sage: cmp_abs(CC(1), CC(0))
2371
1
2372
sage: cmp_abs(CC(0), CC(0))
2373
0
2374
sage: cmp_abs(CC(2,1), CC(1,2))
2375
0
2376
"""
2377
if mpfr_zero_p(b.__re) and mpfr_zero_p(b.__im):
2378
return not ((mpfr_zero_p(a.__re) and mpfr_zero_p(a.__im)))
2379
elif (mpfr_zero_p(a.__re) and mpfr_zero_p(a.__im)):
2380
return -1
2381
cdef mp_exp_t exp_diff = max_exp(a) - max_exp(b)
2382
if exp_diff <= -2:
2383
return -1
2384
elif exp_diff >= 2:
2385
return 1
2386
2387
cdef int res
2388
cdef mpfr_t abs_a, abs_b, tmp
2389
mpfr_init2(abs_a, mpfr_get_prec(a.__re))
2390
mpfr_init2(abs_b, mpfr_get_prec(b.__re))
2391
mpfr_init2(tmp, mpfr_get_prec(a.__re))
2392
2393
mpfr_sqr(abs_a, a.__re, rnd)
2394
mpfr_sqr(tmp, a.__im, rnd)
2395
mpfr_add(abs_a, abs_a, tmp, rnd)
2396
2397
mpfr_sqr(abs_b, b.__re, rnd)
2398
mpfr_sqr(tmp, b.__im, rnd)
2399
mpfr_add(abs_b, abs_b, tmp, rnd)
2400
2401
res = mpfr_cmpabs(abs_a, abs_b)
2402
2403
mpfr_clear(abs_a)
2404
mpfr_clear(abs_b)
2405
mpfr_clear(tmp)
2406
2407
return res
2408
2409