Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/functions/other.py
4056 views
1
"""
2
Other functions
3
"""
4
from sage.symbolic.function import GinacFunction, BuiltinFunction
5
from sage.symbolic.expression import Expression
6
from sage.symbolic.pynac import register_symbol, symbol_table
7
from sage.symbolic.pynac import py_factorial_py
8
from sage.libs.pari.gen import pari
9
from sage.symbolic.all import SR
10
from sage.rings.all import Integer, Rational, RealField, RR, \
11
is_ComplexNumber, ComplexField
12
from sage.misc.latex import latex
13
import math
14
15
import sage.structure.element
16
coercion_model = sage.structure.element.get_coercion_model()
17
18
from sage.structure.coerce import parent
19
from sage.symbolic.constants import pi
20
from sage.symbolic.function import is_inexact
21
from sage.functions.log import exp
22
from sage.functions.trig import arctan2
23
from sage.functions.transcendental import Ei
24
from sage.libs.mpmath import utils as mpmath_utils
25
26
one_half = ~SR(2)
27
28
class Function_erf(BuiltinFunction):
29
r"""
30
The error function, defined for real values as
31
32
`\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt`.
33
34
This function is also defined for complex values, via analytic
35
continuation.
36
37
38
EXAMPLES:
39
40
We can evaluate numerically::
41
42
sage: erf(2)
43
erf(2)
44
sage: erf(2).n()
45
0.995322265018953
46
sage: erf(2).n(100)
47
0.99532226501895273416206925637
48
sage: erf(ComplexField(100)(2+3j))
49
-20.829461427614568389103088452 + 8.6873182714701631444280787545*I
50
51
Basic symbolic properties are handled by Sage and Maxima::
52
53
sage: x = var("x")
54
sage: diff(erf(x),x)
55
2*e^(-x^2)/sqrt(pi)
56
sage: integrate(erf(x),x)
57
x*erf(x) + e^(-x^2)/sqrt(pi)
58
59
ALGORITHM:
60
61
Sage implements numerical evaluation of the error function via the ``erfc()``
62
function in PARI. Symbolics are handled by Sage and Maxima.
63
64
REFERENCES:
65
66
- http://en.wikipedia.org/wiki/Error_function
67
68
TESTS:
69
70
Check limits::
71
72
sage: limit(erf(x),x=0)
73
0
74
sage: limit(erf(x),x=infinity)
75
1
76
77
Check that it's odd::
78
79
sage: erf(1.0)
80
0.842700792949715
81
sage: erf(-1.0)
82
-0.842700792949715
83
84
Check against other implementations and against the definition::
85
86
sage: erf(3).n()
87
0.999977909503001
88
sage: maxima.erf(3).n()
89
0.999977909503001
90
sage: (1-pari(3).erfc())
91
0.999977909503001
92
sage: RR(3).erf()
93
0.999977909503001
94
sage: (integrate(exp(-x**2),(x,0,3))*2/sqrt(pi)).n()
95
0.999977909503001
96
97
:trac:`9044`::
98
99
sage: N(erf(sqrt(2)),200)
100
0.95449973610364158559943472566693312505644755259664313203267
101
102
:trac:`11626`::
103
104
sage: n(erf(2),100)
105
0.99532226501895273416206925637
106
sage: erf(2).n(100)
107
0.99532226501895273416206925637
108
109
Test (indirectly) :trac:`11885`::
110
111
sage: erf(float(0.5))
112
0.5204998778130465
113
sage: erf(complex(0.5))
114
(0.5204998778130465+0j)
115
116
Ensure conversion from maxima elements works::
117
118
sage: merf = maxima(erf(x)).sage().operator()
119
sage: merf == erf
120
True
121
122
Make sure we can dump and load it::
123
124
sage: loads(dumps(erf(2)))
125
erf(2)
126
127
Special-case 0 for immediate evaluation::
128
129
sage: erf(0)
130
0
131
sage: solve(erf(x)==0,x)
132
[x == 0]
133
134
Make sure that we can hold::
135
136
sage: erf(0,hold=True)
137
erf(0)
138
sage: simplify(erf(0,hold=True))
139
0
140
141
Check that high-precision ComplexField inputs work::
142
143
sage: CC(erf(ComplexField(1000)(2+3j)))
144
-20.8294614276146 + 8.68731827147016*I
145
"""
146
147
def __init__(self):
148
r"""
149
See docstring for :meth:`Function_erf`.
150
151
EXAMPLES::
152
153
sage: erf(2)
154
erf(2)
155
"""
156
BuiltinFunction.__init__(self, "erf", latex_name=r"\text{erf}")
157
158
def _eval_(self, x):
159
"""
160
EXAMPLES:
161
162
Input is not an expression but is exact::
163
164
sage: erf(0)
165
0
166
sage: erf(1)
167
erf(1)
168
169
Input is not an expression and is not exact::
170
171
sage: erf(0.0)
172
0.000000000000000
173
174
Input is an expression but not a trivial zero::
175
176
sage: erf(x)
177
erf(x)
178
179
Input is an expression which is trivially zero::
180
181
sage: erf(SR(0))
182
0
183
"""
184
if not isinstance(x, Expression):
185
if is_inexact(x):
186
return self._evalf_(x, parent=parent(x))
187
elif x == Integer(0):
188
return Integer(0)
189
elif x.is_trivial_zero():
190
return x
191
return None
192
193
def _evalf_(self, x, parent):
194
"""
195
EXAMPLES::
196
197
sage: erf(2).n()
198
0.995322265018953
199
sage: erf(2).n(200)
200
0.99532226501895273416206925636725292861089179704006007673835
201
sage: erf(pi - 1/2*I).n(100)
202
1.0000111669099367825726058952 + 1.6332655417638522934072124548e-6*I
203
204
TESTS:
205
206
Check that PARI/GP through the GP interface gives the same answer::
207
208
sage: gp.set_real_precision(59) # random
209
38
210
sage: print gp.eval("1 - erfc(1)"); print erf(1).n(200);
211
0.84270079294971486934122063508260925929606699796630290845994
212
0.84270079294971486934122063508260925929606699796630290845994
213
"""
214
try:
215
prec = parent.prec()
216
except AttributeError: # not a Sage parent
217
prec = 0
218
return parent(1) - parent(pari(x).erfc(prec))
219
220
def _derivative_(self, x, diff_param=None):
221
"""
222
Derivative of erf function
223
224
EXAMPLES::
225
226
sage: erf(x).diff(x)
227
2*e^(-x^2)/sqrt(pi)
228
229
TESTS:
230
231
Check if #8568 is fixed::
232
233
sage: var('c,x')
234
(c, x)
235
sage: derivative(erf(c*x),x)
236
2*c*e^(-c^2*x^2)/sqrt(pi)
237
sage: erf(c*x).diff(x)._maxima_init_()
238
'((%pi)^(-1/2))*(c)*(exp(((c)^(2))*((x)^(2))*(-1)))*(2)'
239
"""
240
return 2*exp(-x**2)/sqrt(pi)
241
242
erf = Function_erf()
243
244
class Function_abs(GinacFunction):
245
def __init__(self):
246
r"""
247
The absolute value function.
248
249
EXAMPLES::
250
251
sage: var('x y')
252
(x, y)
253
sage: abs(x)
254
abs(x)
255
sage: abs(x^2 + y^2)
256
abs(x^2 + y^2)
257
sage: abs(-2)
258
2
259
sage: sqrt(x^2)
260
sqrt(x^2)
261
sage: abs(sqrt(x))
262
abs(sqrt(x))
263
sage: complex(abs(3*I))
264
(3+0j)
265
266
sage: f = sage.functions.other.Function_abs()
267
sage: latex(f)
268
\mathrm{abs}
269
sage: latex(abs(x))
270
{\left| x \right|}
271
272
Test pickling::
273
274
sage: loads(dumps(abs(x)))
275
abs(x)
276
"""
277
GinacFunction.__init__(self, "abs", latex_name=r"\mathrm{abs}")
278
279
abs = abs_symbolic = Function_abs()
280
281
class Function_ceil(BuiltinFunction):
282
def __init__(self):
283
r"""
284
The ceiling function.
285
286
The ceiling of `x` is computed in the following manner.
287
288
289
#. The ``x.ceil()`` method is called and returned if it
290
is there. If it is not, then Sage checks if `x` is one of
291
Python's native numeric data types. If so, then it calls and
292
returns ``Integer(int(math.ceil(x)))``.
293
294
#. Sage tries to convert `x` into a
295
``RealIntervalField`` with 53 bits of precision. Next,
296
the ceilings of the endpoints are computed. If they are the same,
297
then that value is returned. Otherwise, the precision of the
298
``RealIntervalField`` is increased until they do match
299
up or it reaches ``maximum_bits`` of precision.
300
301
#. If none of the above work, Sage returns a
302
``Expression`` object.
303
304
305
EXAMPLES::
306
307
sage: a = ceil(2/5 + x)
308
sage: a
309
ceil(x + 2/5)
310
sage: a(x=4)
311
5
312
sage: a(x=4.0)
313
5
314
sage: ZZ(a(x=3))
315
4
316
sage: a = ceil(x^3 + x + 5/2); a
317
ceil(x^3 + x + 5/2)
318
sage: a.simplify()
319
ceil(x^3 + x + 1/2) + 2
320
sage: a(x=2)
321
13
322
323
::
324
325
sage: ceil(log(8)/log(2))
326
3
327
328
::
329
330
sage: ceil(5.4)
331
6
332
sage: type(ceil(5.4))
333
<type 'sage.rings.integer.Integer'>
334
335
::
336
337
sage: ceil(factorial(50)/exp(1))
338
11188719610782480504630258070757734324011354208865721592720336801
339
sage: ceil(SR(10^50 + 10^(-50)))
340
100000000000000000000000000000000000000000000000001
341
sage: ceil(SR(10^50 - 10^(-50)))
342
100000000000000000000000000000000000000000000000000
343
344
sage: ceil(sec(e))
345
-1
346
347
sage: latex(ceil(x))
348
\left \lceil x \right \rceil
349
350
::
351
352
sage: import numpy
353
sage: a = numpy.linspace(0,2,6)
354
sage: ceil(a)
355
array([ 0., 1., 1., 2., 2., 2.])
356
357
Test pickling::
358
359
sage: loads(dumps(ceil))
360
ceil
361
"""
362
BuiltinFunction.__init__(self, "ceil",
363
conversions=dict(maxima='ceiling'))
364
365
def _print_latex_(self, x):
366
r"""
367
EXAMPLES::
368
369
sage: latex(ceil(x)) # indirect doctest
370
\left \lceil x \right \rceil
371
"""
372
return r"\left \lceil %s \right \rceil"%latex(x)
373
374
#FIXME: this should be moved to _eval_
375
def __call__(self, x, maximum_bits=20000):
376
"""
377
Allows an object of this class to behave like a function. If
378
``ceil`` is an instance of this class, we can do ``ceil(n)`` to get
379
the ceiling of ``n``.
380
381
TESTS::
382
383
sage: ceil(SR(10^50 + 10^(-50)))
384
100000000000000000000000000000000000000000000000001
385
sage: ceil(SR(10^50 - 10^(-50)))
386
100000000000000000000000000000000000000000000000000
387
sage: ceil(int(10^50))
388
100000000000000000000000000000000000000000000000000
389
"""
390
try:
391
return x.ceil()
392
except AttributeError:
393
if isinstance(x, (int, long)):
394
return Integer(x)
395
elif isinstance(x, (float, complex)):
396
return Integer(int(math.ceil(x)))
397
elif type(x).__module__ == 'numpy':
398
import numpy
399
return numpy.ceil(x)
400
401
x_original = x
402
403
from sage.rings.all import RealIntervalField
404
# If x can be coerced into a real interval, then we should
405
# try increasing the number of bits of precision until
406
# we get the ceiling at each of the endpoints is the same.
407
# The precision will continue to be increased up to maximum_bits
408
# of precision at which point it will raise a value error.
409
bits = 53
410
try:
411
x_interval = RealIntervalField(bits)(x)
412
upper_ceil = x_interval.upper().ceil()
413
lower_ceil = x_interval.lower().ceil()
414
415
while upper_ceil != lower_ceil and bits < maximum_bits:
416
bits += 100
417
x_interval = RealIntervalField(bits)(x)
418
upper_ceil = x_interval.upper().ceil()
419
lower_ceil = x_interval.lower().ceil()
420
421
if bits < maximum_bits:
422
return lower_ceil
423
else:
424
try:
425
return ceil(SR(x).full_simplify())
426
except ValueError:
427
pass
428
raise ValueError, "x (= %s) requires more than %s bits of precision to compute its ceiling"%(x, maximum_bits)
429
430
except TypeError:
431
# If x cannot be coerced into a RealField, then
432
# it should be left as a symbolic expression.
433
return BuiltinFunction.__call__(self, SR(x_original))
434
435
def _eval_(self, x):
436
"""
437
EXAMPLES::
438
439
sage: ceil(x).subs(x==7.5)
440
8
441
sage: ceil(x)
442
ceil(x)
443
"""
444
try:
445
return x.ceil()
446
except AttributeError:
447
if isinstance(x, (int, long)):
448
return Integer(x)
449
elif isinstance(x, (float, complex)):
450
return Integer(int(math.ceil(x)))
451
return None
452
453
ceil = Function_ceil()
454
455
456
class Function_floor(BuiltinFunction):
457
def __init__(self):
458
r"""
459
The floor function.
460
461
The floor of `x` is computed in the following manner.
462
463
464
#. The ``x.floor()`` method is called and returned if
465
it is there. If it is not, then Sage checks if `x` is one
466
of Python's native numeric data types. If so, then it calls and
467
returns ``Integer(int(math.floor(x)))``.
468
469
#. Sage tries to convert `x` into a
470
``RealIntervalField`` with 53 bits of precision. Next,
471
the floors of the endpoints are computed. If they are the same,
472
then that value is returned. Otherwise, the precision of the
473
``RealIntervalField`` is increased until they do match
474
up or it reaches ``maximum_bits`` of precision.
475
476
#. If none of the above work, Sage returns a
477
symbolic ``Expression`` object.
478
479
480
EXAMPLES::
481
482
sage: floor(5.4)
483
5
484
sage: type(floor(5.4))
485
<type 'sage.rings.integer.Integer'>
486
sage: var('x')
487
x
488
sage: a = floor(5.4 + x); a
489
floor(x + 5.40000000000000)
490
sage: a.simplify()
491
floor(x + 0.4) + 5
492
sage: a(x=2)
493
7
494
495
::
496
497
sage: floor(log(8)/log(2))
498
3
499
500
::
501
502
sage: floor(factorial(50)/exp(1))
503
11188719610782480504630258070757734324011354208865721592720336800
504
sage: floor(SR(10^50 + 10^(-50)))
505
100000000000000000000000000000000000000000000000000
506
sage: floor(SR(10^50 - 10^(-50)))
507
99999999999999999999999999999999999999999999999999
508
sage: floor(int(10^50))
509
100000000000000000000000000000000000000000000000000
510
511
::
512
513
sage: import numpy
514
sage: a = numpy.linspace(0,2,6)
515
sage: floor(a)
516
array([ 0., 0., 0., 1., 1., 2.])
517
518
Test pickling::
519
520
sage: loads(dumps(floor))
521
floor
522
"""
523
BuiltinFunction.__init__(self, "floor")
524
525
def _print_latex_(self, x):
526
r"""
527
EXAMPLES::
528
529
sage: latex(floor(x))
530
\left \lfloor x \right \rfloor
531
"""
532
return r"\left \lfloor %s \right \rfloor"%latex(x)
533
534
#FIXME: this should be moved to _eval_
535
def __call__(self, x, maximum_bits=20000):
536
"""
537
Allows an object of this class to behave like a function. If
538
``floor`` is an instance of this class, we can do ``floor(n)`` to
539
obtain the floor of ``n``.
540
541
TESTS::
542
543
sage: floor(SR(10^50 + 10^(-50)))
544
100000000000000000000000000000000000000000000000000
545
sage: floor(SR(10^50 - 10^(-50)))
546
99999999999999999999999999999999999999999999999999
547
sage: floor(int(10^50))
548
100000000000000000000000000000000000000000000000000
549
"""
550
try:
551
return x.floor()
552
except AttributeError:
553
if isinstance(x, (int, long)):
554
return Integer(x)
555
elif isinstance(x, (float, complex)):
556
return Integer(int(math.floor(x)))
557
elif type(x).__module__ == 'numpy':
558
import numpy
559
return numpy.floor(x)
560
561
x_original = x
562
563
from sage.rings.all import RealIntervalField
564
565
# If x can be coerced into a real interval, then we should
566
# try increasing the number of bits of precision until
567
# we get the floor at each of the endpoints is the same.
568
# The precision will continue to be increased up to maximum_bits
569
# of precision at which point it will raise a value error.
570
bits = 53
571
try:
572
x_interval = RealIntervalField(bits)(x)
573
upper_floor = x_interval.upper().floor()
574
lower_floor = x_interval.lower().floor()
575
576
while upper_floor != lower_floor and bits < maximum_bits:
577
bits += 100
578
x_interval = RealIntervalField(bits)(x)
579
upper_floor = x_interval.upper().floor()
580
lower_floor = x_interval.lower().floor()
581
582
if bits < maximum_bits:
583
return lower_floor
584
else:
585
try:
586
return floor(SR(x).full_simplify())
587
except ValueError:
588
pass
589
raise ValueError, "x (= %s) requires more than %s bits of precision to compute its floor"%(x, maximum_bits)
590
591
except TypeError:
592
# If x cannot be coerced into a RealField, then
593
# it should be left as a symbolic expression.
594
return BuiltinFunction.__call__(self, SR(x_original))
595
596
def _eval_(self, x):
597
"""
598
EXAMPLES::
599
600
sage: floor(x).subs(x==7.5)
601
7
602
sage: floor(x)
603
floor(x)
604
"""
605
try:
606
return x.floor()
607
except AttributeError:
608
if isinstance(x, (int, long)):
609
return Integer(x)
610
elif isinstance(x, (float, complex)):
611
return Integer(int(math.floor(x)))
612
return None
613
614
floor = Function_floor()
615
616
class Function_gamma(GinacFunction):
617
def __init__(self):
618
r"""
619
The Gamma function. This is defined by
620
621
.. math::
622
623
\Gamma(z) = \int_0^\infty t^{z-1}e^{-t} dt
624
625
for complex input `z` with real part greater than zero, and by
626
analytic continuation on the rest of the complex plane (except
627
for negative integers, which are poles).
628
629
It is computed by various libraries within Sage, depending on
630
the input type.
631
632
EXAMPLES::
633
634
sage: from sage.functions.other import gamma1
635
sage: gamma1(CDF(0.5,14))
636
-4.05370307804e-10 - 5.77329983455e-10*I
637
sage: gamma1(CDF(I))
638
-0.154949828302 - 0.498015668118*I
639
640
Recall that `\Gamma(n)` is `n-1` factorial::
641
642
sage: gamma1(11) == factorial(10)
643
True
644
sage: gamma1(6)
645
120
646
sage: gamma1(1/2)
647
sqrt(pi)
648
sage: gamma1(-1)
649
Infinity
650
sage: gamma1(I)
651
gamma(I)
652
sage: gamma1(x/2)(x=5)
653
3/4*sqrt(pi)
654
655
sage: gamma1(float(6))
656
120.0
657
sage: gamma1(x)
658
gamma(x)
659
660
::
661
662
sage: gamma1(pi)
663
gamma(pi)
664
sage: gamma1(i)
665
gamma(I)
666
sage: gamma1(i).n()
667
-0.154949828301811 - 0.498015668118356*I
668
sage: gamma1(int(5))
669
24
670
671
::
672
673
sage: conjugate(gamma(x))
674
gamma(conjugate(x))
675
676
::
677
678
sage: plot(gamma1(x),(x,1,5))
679
680
To prevent automatic evaluation use the ``hold`` argument::
681
682
sage: gamma1(1/2,hold=True)
683
gamma(1/2)
684
685
To then evaluate again, we currently must use Maxima via
686
:meth:`sage.symbolic.expression.Expression.simplify`::
687
688
sage: gamma1(1/2,hold=True).simplify()
689
sqrt(pi)
690
691
TESTS:
692
693
We verify that we can convert this function to Maxima and
694
convert back to Sage::
695
696
sage: z = var('z')
697
sage: maxima(gamma1(z)).sage()
698
gamma(z)
699
sage: latex(gamma1(z))
700
\Gamma\left(z\right)
701
702
Test that Trac ticket 5556 is fixed::
703
704
sage: gamma1(3/4)
705
gamma(3/4)
706
707
sage: gamma1(3/4).n(100)
708
1.2254167024651776451290983034
709
710
Check that negative integer input works::
711
712
sage: (-1).gamma()
713
Infinity
714
sage: (-1.).gamma()
715
NaN
716
sage: CC(-1).gamma()
717
Infinity
718
sage: RDF(-1).gamma()
719
NaN
720
sage: CDF(-1).gamma()
721
Infinity
722
723
Check if #8297 is fixed::
724
725
sage: latex(gamma(1/4))
726
\Gamma\left(\frac{1}{4}\right)
727
728
Test pickling::
729
730
sage: loads(dumps(gamma(x)))
731
gamma(x)
732
"""
733
GinacFunction.__init__(self, "gamma", latex_name=r'\Gamma',
734
ginac_name='tgamma',
735
conversions={'mathematica':'Gamma','maple':'GAMMA'})
736
737
def __call__(self, x, prec=None, coerce=True, hold=False):
738
"""
739
Note that the ``prec`` argument is deprecated. The precision for
740
the result is deduced from the precision of the input. Convert
741
the input to a higher precision explicitly if a result with higher
742
precision is desired.::
743
744
sage: t = gamma(RealField(100)(2.5)); t
745
1.3293403881791370204736256125
746
sage: t.prec()
747
100
748
749
sage: gamma(6, prec=53)
750
doctest:...: DeprecationWarning: The prec keyword argument is deprecated. Explicitly set the precision of the input, for example gamma(RealField(300)(1)), or use the prec argument to .n() for exact inputs, e.g., gamma(1).n(300), instead.
751
120.000000000000
752
753
TESTS::
754
755
sage: gamma(pi,prec=100)
756
2.2880377953400324179595889091
757
758
sage: gamma(3/4,prec=100)
759
1.2254167024651776451290983034
760
"""
761
if prec is not None:
762
from sage.misc.misc import deprecation
763
deprecation("The prec keyword argument is deprecated. Explicitly set the precision of the input, for example gamma(RealField(300)(1)), or use the prec argument to .n() for exact inputs, e.g., gamma(1).n(300), instead.")
764
import mpmath
765
return mpmath_utils.call(mpmath.gamma, x, prec=prec)
766
767
# this is a kludge to keep
768
# sage: Q.<i> = NumberField(x^2+1)
769
# sage: gamma(i)
770
# working, since number field elements cannot be coerced into SR
771
# without specifying an explicit embedding into CC any more
772
try:
773
res = GinacFunction.__call__(self, x, coerce=coerce, hold=hold)
774
except TypeError, err:
775
# the __call__() method returns a TypeError for fast float arguments
776
# as well, we only proceed if the error message says that
777
# the arguments cannot be coerced to SR
778
if not str(err).startswith("cannot coerce"):
779
raise
780
781
from sage.misc.misc import deprecation
782
deprecation("Calling symbolic functions with arguments that cannot be coerced into symbolic expressions is deprecated.")
783
parent = RR if prec is None else RealField(prec)
784
try:
785
x = parent(x)
786
except (ValueError, TypeError):
787
x = parent.complex_field()(x)
788
res = GinacFunction.__call__(self, x, coerce=coerce, hold=hold)
789
790
return res
791
792
gamma1 = Function_gamma()
793
794
class Function_log_gamma(GinacFunction):
795
def __init__(self):
796
r"""
797
The principal branch of the logarithm of Gamma function.
798
Gamma is defined for complex input `z` with real part greater
799
than zero, and by analytic continuation on the rest of the
800
complex plane (except for negative integers, which are poles).
801
802
It is computed by the `log_gamma` function for the number type,
803
or by `lgamma` in Ginac, failing that.
804
805
EXAMPLES:
806
807
Numerical evaluation happens when appropriate, to the
808
appropriate accuracy (see #10072)::
809
810
sage: log_gamma(6)
811
log(120)
812
sage: log_gamma(6.)
813
4.78749174278205
814
sage: log_gamma(6).n()
815
4.78749174278205
816
sage: log_gamma(RealField(100)(6))
817
4.7874917427820459942477009345
818
sage: log_gamma(2.4+i)
819
-0.0308566579348816 + 0.693427705955790*I
820
sage: log_gamma(-3.1)
821
0.400311696703985
822
823
Symbolic input works (see #10075)::
824
825
sage: log_gamma(3*x)
826
log_gamma(3*x)
827
sage: log_gamma(3+i)
828
log_gamma(I + 3)
829
sage: log_gamma(3+i+x)
830
log_gamma(x + I + 3)
831
832
To get evaluation of input for which gamma
833
is negative and the ceiling is even, we must
834
explicitly make the input complex. This is
835
a known issue, see #12521::
836
837
sage: log_gamma(-2.1)
838
NaN
839
sage: log_gamma(CC(-2.1))
840
1.53171380819509 + 3.14159265358979*I
841
842
In order to prevent evaluation, use the `hold` argument;
843
to evaluate a held expression, use the `n()` numerical
844
evaluation method::
845
846
sage: log_gamma(SR(5),hold=True)
847
log_gamma(5)
848
sage: log_gamma(SR(5),hold=True).n()
849
3.17805383034795
850
851
TESTS::
852
853
sage: log_gamma(-2.1+i)
854
-1.90373724496982 - 0.901638463592247*I
855
sage: log_gamma(pari(6))
856
4.78749174278205
857
sage: log_gamma(CC(6))
858
4.78749174278205
859
sage: log_gamma(CC(-2.5))
860
-0.0562437164976740 + 3.14159265358979*I
861
862
``conjugate(log_gamma(x))==log_gamma(conjugate(x))`` unless on the branch
863
cut, which runs along the negative real axis.::
864
865
sage: conjugate(log_gamma(x))
866
conjugate(log_gamma(x))
867
sage: var('y', domain='positive')
868
y
869
sage: conjugate(log_gamma(y))
870
log_gamma(y)
871
sage: conjugate(log_gamma(y+I))
872
conjugate(log_gamma(y + I))
873
sage: log_gamma(-2)
874
+Infinity
875
sage: conjugate(log_gamma(-2))
876
+Infinity
877
"""
878
GinacFunction.__init__(self, "log_gamma", latex_name=r'\log\Gamma',
879
conversions={'mathematica':'LogGamma','maxima':'log_gamma'})
880
881
log_gamma = Function_log_gamma()
882
883
class Function_gamma_inc(BuiltinFunction):
884
def __init__(self):
885
r"""
886
The incomplete gamma function.
887
888
EXAMPLES::
889
890
sage: gamma_inc(CDF(0,1), 3)
891
0.00320857499337 + 0.0124061858119*I
892
sage: gamma_inc(RDF(1), 3)
893
0.0497870683678639
894
sage: gamma_inc(3,2)
895
gamma(3, 2)
896
sage: gamma_inc(x,0)
897
gamma(x)
898
sage: latex(gamma_inc(3,2))
899
\Gamma\left(3, 2\right)
900
sage: loads(dumps((gamma_inc(3,2))))
901
gamma(3, 2)
902
sage: i = ComplexField(30).0; gamma_inc(2, 1 + i)
903
0.70709210 - 0.42035364*I
904
sage: gamma_inc(2., 5)
905
0.0404276819945128
906
"""
907
BuiltinFunction.__init__(self, "gamma", nargs=2, latex_name=r"\Gamma",
908
conversions={'maxima':'gamma_incomplete', 'mathematica':'Gamma',
909
'maple':'GAMMA'})
910
911
def _eval_(self, x, y):
912
"""
913
EXAMPLES::
914
915
sage: gamma_inc(2.,0)
916
1.00000000000000
917
sage: gamma_inc(2,0)
918
1
919
sage: gamma_inc(1/2,2)
920
-(erf(sqrt(2)) - 1)*sqrt(pi)
921
sage: gamma_inc(1/2,1)
922
-(erf(1) - 1)*sqrt(pi)
923
sage: gamma_inc(1/2,0)
924
sqrt(pi)
925
sage: gamma_inc(x,0)
926
gamma(x)
927
sage: gamma_inc(1,2)
928
e^(-2)
929
sage: gamma_inc(0,2)
930
-Ei(-2)
931
"""
932
if not isinstance(x, Expression) and not isinstance(y, Expression) and \
933
(is_inexact(x) or is_inexact(y)):
934
x, y = coercion_model.canonical_coercion(x, y)
935
return self._evalf_(x, y, parent(x))
936
937
if y == 0:
938
return gamma(x)
939
if x == 1:
940
return exp(-y)
941
if x == 0:
942
return -Ei(-y)
943
if x == Rational(1)/2: #only for x>0
944
return sqrt(pi)*(1-erf(sqrt(y)))
945
return None
946
947
def _evalf_(self, x, y, parent=None):
948
"""
949
EXAMPLES::
950
951
sage: gamma_inc(0,2)
952
-Ei(-2)
953
sage: gamma_inc(0,2.)
954
0.0489005107080611
955
sage: gamma_inc(3,2).n()
956
1.35335283236613
957
"""
958
try:
959
return x.gamma_inc(y)
960
except AttributeError:
961
if not (is_ComplexNumber(x)):
962
if is_ComplexNumber(y):
963
C = y.parent()
964
else:
965
C = ComplexField()
966
x = C(x)
967
return x.gamma_inc(y)
968
969
# synonym.
970
incomplete_gamma = gamma_inc=Function_gamma_inc()
971
972
def gamma(a, *args, **kwds):
973
r"""
974
Gamma and incomplete gamma functions.
975
This is defined by the integral
976
977
.. math::
978
979
\Gamma(a, z) = \int_z^\infty t^{a-1}e^{-t} dt
980
981
EXAMPLES::
982
983
Recall that `\Gamma(n)` is `n-1` factorial::
984
985
sage: gamma(11) == factorial(10)
986
True
987
sage: gamma(6)
988
120
989
sage: gamma(1/2)
990
sqrt(pi)
991
sage: gamma(-4/3)
992
gamma(-4/3)
993
sage: gamma(-1)
994
Infinity
995
sage: gamma(0)
996
Infinity
997
998
::
999
1000
sage: gamma_inc(3,2)
1001
gamma(3, 2)
1002
sage: gamma_inc(x,0)
1003
gamma(x)
1004
1005
::
1006
1007
sage: gamma(5, hold=True)
1008
gamma(5)
1009
sage: gamma(x, 0, hold=True)
1010
gamma(x, 0)
1011
1012
::
1013
1014
sage: gamma(CDF(0.5,14))
1015
-4.05370307804e-10 - 5.77329983455e-10*I
1016
sage: gamma(CDF(I))
1017
-0.154949828302 - 0.498015668118*I
1018
1019
The gamma function only works with input that can be coerced to the
1020
Symbolic Ring::
1021
1022
sage: Q.<i> = NumberField(x^2+1)
1023
sage: gamma(i)
1024
doctest:...: DeprecationWarning: Calling symbolic functions with arguments that cannot be coerced into symbolic expressions is deprecated.
1025
-0.154949828301811 - 0.498015668118356*I
1026
1027
We make an exception for elements of AA or QQbar, which cannot be
1028
coerced into symbolic expressions to allow this usage::
1029
1030
sage: t = QQbar(sqrt(2)) + sqrt(3); t
1031
3.146264369941973?
1032
sage: t.parent()
1033
Algebraic Field
1034
1035
Symbolic functions convert the arguments to symbolic expressions if they
1036
are in QQbar or AA::
1037
1038
sage: gamma(QQbar(I))
1039
-0.154949828301811 - 0.498015668118356*I
1040
"""
1041
if not args:
1042
return gamma1(a, **kwds)
1043
if len(args) > 1:
1044
raise TypeError, "Symbolic function gamma takes at most 2 arguments (%s given)"%(len(args)+1)
1045
return incomplete_gamma(a,args[0],**kwds)
1046
1047
# We have to add the wrapper function manually to the symbol_table when we have
1048
# two functions with different number of arguments and the same name
1049
symbol_table['functions']['gamma'] = gamma
1050
1051
class Function_psi1(GinacFunction):
1052
def __init__(self):
1053
r"""
1054
The digamma function, `\psi(x)`, is the logarithmic derivative of the
1055
gamma function.
1056
1057
.. math::
1058
1059
\psi(x) = \frac{d}{dx} \log(\Gamma(x)) = \frac{\Gamma'(x)}{\Gamma(x)}
1060
1061
EXAMPLES::
1062
1063
sage: from sage.functions.other import psi1
1064
sage: psi1(x)
1065
psi(x)
1066
sage: psi1(x).derivative(x)
1067
psi(1, x)
1068
1069
::
1070
1071
sage: psi1(3)
1072
-euler_gamma + 3/2
1073
1074
::
1075
1076
sage: psi(.5)
1077
-1.96351002602142
1078
sage: psi(RealField(100)(.5))
1079
-1.9635100260214234794409763330
1080
1081
TESTS::
1082
1083
sage: latex(psi1(x))
1084
\psi\left(x\right)
1085
sage: loads(dumps(psi1(x)+1))
1086
psi(x) + 1
1087
1088
sage: t = psi1(x); t
1089
psi(x)
1090
sage: t.subs(x=.2)
1091
-5.28903989659219
1092
"""
1093
GinacFunction.__init__(self, "psi", nargs=1, latex_name='\psi',
1094
conversions=dict(maxima='psi[0]', mathematica='PolyGamma'))
1095
1096
class Function_psi2(GinacFunction):
1097
def __init__(self):
1098
r"""
1099
Derivatives of the digamma function `\psi(x)`. T
1100
1101
EXAMPLES::
1102
1103
sage: from sage.functions.other import psi2
1104
sage: psi2(2, x)
1105
psi(2, x)
1106
sage: psi2(2, x).derivative(x)
1107
psi(3, x)
1108
sage: n = var('n')
1109
sage: psi2(n, x).derivative(x)
1110
psi(n + 1, x)
1111
1112
::
1113
1114
sage: psi2(0, x)
1115
psi(x)
1116
sage: psi2(-1, x)
1117
log(gamma(x))
1118
sage: psi2(3, 1)
1119
1/15*pi^4
1120
1121
::
1122
1123
sage: psi2(2, .5).n()
1124
-16.8287966442343
1125
sage: psi2(2, .5).n(100)
1126
-16.828796644234319995596334261
1127
1128
Tests::
1129
1130
sage: psi2(n, x).derivative(n)
1131
Traceback (most recent call last):
1132
...
1133
RuntimeError: cannot diff psi(n,x) with respect to n
1134
1135
sage: latex(psi2(2,x))
1136
\psi\left(2, x\right)
1137
sage: loads(dumps(psi2(2,x)+1))
1138
psi(2, x) + 1
1139
"""
1140
GinacFunction.__init__(self, "psi", nargs=2, latex_name='\psi',
1141
conversions=dict(mathematica='PolyGamma'))
1142
1143
def _maxima_init_evaled_(self, n, x):
1144
"""
1145
EXAMPLES:
1146
1147
These are indirect doctests for this function.::
1148
1149
sage: from sage.functions.other import psi2
1150
sage: psi2(2, x)._maxima_()
1151
psi[2](x)
1152
sage: psi2(4, x)._maxima_()
1153
psi[4](x)
1154
"""
1155
return "psi[%s](%s)"%(n, x)
1156
1157
psi1 = Function_psi1()
1158
psi2 = Function_psi2()
1159
1160
def psi(x, *args, **kwds):
1161
r"""
1162
The digamma function, `\psi(x)`, is the logarithmic derivative of the
1163
gamma function.
1164
1165
.. math::
1166
1167
\psi(x) = \frac{d}{dx} \log(\Gamma(x)) = \frac{\Gamma'(x)}{\Gamma(x)}
1168
1169
We represent the `n`-th derivative of the digamma function with
1170
`\psi(n, x)` or `psi(n, x)`.
1171
1172
EXAMPLES::
1173
1174
sage: psi(x)
1175
psi(x)
1176
sage: psi(.5)
1177
-1.96351002602142
1178
sage: psi(3)
1179
-euler_gamma + 3/2
1180
sage: psi(1, 5)
1181
1/6*pi^2 - 205/144
1182
sage: psi(1, x)
1183
psi(1, x)
1184
sage: psi(1, x).derivative(x)
1185
psi(2, x)
1186
1187
::
1188
1189
sage: psi(3, hold=True)
1190
psi(3)
1191
sage: psi(1, 5, hold=True)
1192
psi(1, 5)
1193
1194
TESTS::
1195
1196
sage: psi(2, x, 3)
1197
Traceback (most recent call last):
1198
...
1199
TypeError: Symbolic function psi takes at most 2 arguments (3 given)
1200
"""
1201
if not args:
1202
return psi1(x, **kwds)
1203
if len(args) > 1:
1204
raise TypeError, "Symbolic function psi takes at most 2 arguments (%s given)"%(len(args)+1)
1205
return psi2(x,args[0],**kwds)
1206
1207
# We have to add the wrapper function manually to the symbol_table when we have
1208
# two functions with different number of arguments and the same name
1209
symbol_table['functions']['psi'] = psi
1210
1211
class Function_factorial(GinacFunction):
1212
def __init__(self):
1213
r"""
1214
Returns the factorial of `n`.
1215
1216
INPUT:
1217
1218
- ``n`` - any complex argument (except negative
1219
integers) or any symbolic expression
1220
1221
1222
OUTPUT: an integer or symbolic expression
1223
1224
EXAMPLES::
1225
1226
sage: x = var('x')
1227
sage: factorial(0)
1228
1
1229
sage: factorial(4)
1230
24
1231
sage: factorial(10)
1232
3628800
1233
sage: factorial(6) == 6*5*4*3*2
1234
True
1235
sage: f = factorial(x + factorial(x)); f
1236
factorial(x + factorial(x))
1237
sage: f(x=3)
1238
362880
1239
sage: factorial(x)^2
1240
factorial(x)^2
1241
1242
To prevent automatic evaluation use the ``hold`` argument::
1243
1244
sage: factorial(5,hold=True)
1245
factorial(5)
1246
1247
To then evaluate again, we currently must use Maxima via
1248
:meth:`sage.symbolic.expression.Expression.simplify`::
1249
1250
sage: factorial(5,hold=True).simplify()
1251
120
1252
1253
We can also give input other than nonnegative integers. For
1254
other nonnegative numbers, the :func:`gamma` function is used::
1255
1256
sage: factorial(1/2)
1257
1/2*sqrt(pi)
1258
sage: factorial(3/4)
1259
gamma(7/4)
1260
sage: factorial(2.3)
1261
2.68343738195577
1262
1263
But negative input always fails::
1264
1265
sage: factorial(-32)
1266
Traceback (most recent call last):
1267
...
1268
ValueError: factorial -- self = (-32) must be nonnegative
1269
1270
TESTS:
1271
1272
We verify that we can convert this function to Maxima and
1273
bring it back into Sage.::
1274
1275
sage: z = var('z')
1276
sage: factorial._maxima_init_()
1277
'factorial'
1278
sage: maxima(factorial(z))
1279
z!
1280
sage: _.sage()
1281
factorial(z)
1282
sage: k = var('k')
1283
sage: factorial(k)
1284
factorial(k)
1285
1286
sage: factorial(3.14)
1287
7.173269190187...
1288
1289
Test latex typesetting::
1290
1291
sage: latex(factorial(x))
1292
x!
1293
sage: latex(factorial(2*x))
1294
\left(2 \, x\right)!
1295
sage: latex(factorial(sin(x)))
1296
\sin\left(x\right)!
1297
sage: latex(factorial(sqrt(x+1)))
1298
\left(\sqrt{x + 1}\right)!
1299
sage: latex(factorial(sqrt(x)))
1300
\sqrt{x}!
1301
sage: latex(factorial(x^(2/3)))
1302
\left(x^{\frac{2}{3}}\right)!
1303
1304
sage: latex(factorial)
1305
{\rm factorial}
1306
1307
Test pickling::
1308
1309
sage: loads(dumps(factorial))
1310
factorial
1311
"""
1312
GinacFunction.__init__(self, "factorial", latex_name='{\\rm factorial}',
1313
conversions=dict(maxima='factorial', mathematica='Factorial'))
1314
1315
def _eval_(self, x):
1316
"""
1317
Evaluate the factorial function.
1318
1319
Note that this method overrides the eval method defined in GiNaC
1320
which calls numeric evaluation on all numeric input. We preserve
1321
exact results if the input is a rational number.
1322
1323
EXAMPLES::
1324
1325
sage: k = var('k')
1326
sage: k.factorial()
1327
factorial(k)
1328
sage: SR(1/2).factorial()
1329
1/2*sqrt(pi)
1330
sage: SR(3/4).factorial()
1331
gamma(7/4)
1332
sage: SR(5).factorial()
1333
120
1334
sage: SR(3245908723049857203948572398475r).factorial()
1335
factorial(3245908723049857203948572398475L)
1336
sage: SR(3245908723049857203948572398475).factorial()
1337
factorial(3245908723049857203948572398475)
1338
"""
1339
if isinstance(x, Rational):
1340
return gamma(x+1)
1341
elif isinstance(x, (Integer, int)) or \
1342
(not isinstance(x, Expression) and is_inexact(x)):
1343
return py_factorial_py(x)
1344
1345
return None
1346
1347
factorial = Function_factorial()
1348
1349
class Function_binomial(GinacFunction):
1350
def __init__(self):
1351
r"""
1352
Return the binomial coefficient
1353
1354
.. math::
1355
1356
\binom{x}{m} = x (x-1) \cdots (x-m+1) / m!
1357
1358
1359
which is defined for `m \in \ZZ` and any
1360
`x`. We extend this definition to include cases when
1361
`x-m` is an integer but `m` is not by
1362
1363
.. math::
1364
1365
\binom{x}{m}= \binom{x}{x-m}
1366
1367
If `m < 0`, return `0`.
1368
1369
INPUT:
1370
1371
- ``x``, ``m`` - numbers or symbolic expressions. Either ``m``
1372
or ``x-m`` must be an integer, else the output is symbolic.
1373
1374
OUTPUT: number or symbolic expression (if input is symbolic)
1375
1376
EXAMPLES::
1377
1378
sage: binomial(5,2)
1379
10
1380
sage: binomial(2,0)
1381
1
1382
sage: binomial(1/2, 0)
1383
1
1384
sage: binomial(3,-1)
1385
0
1386
sage: binomial(20,10)
1387
184756
1388
sage: binomial(-2, 5)
1389
-6
1390
sage: binomial(RealField()('2.5'), 2)
1391
1.87500000000000
1392
sage: n=var('n'); binomial(n,2)
1393
1/2*(n - 1)*n
1394
sage: n=var('n'); binomial(n,n)
1395
1
1396
sage: n=var('n'); binomial(n,n-1)
1397
n
1398
sage: binomial(2^100, 2^100)
1399
1
1400
1401
::
1402
1403
sage: k, i = var('k,i')
1404
sage: binomial(k,i)
1405
binomial(k, i)
1406
1407
We can use a ``hold`` parameter to prevent automatic evaluation,
1408
but only using method notation::
1409
1410
sage: SR(5).binomial(3, hold=True)
1411
binomial(5, 3)
1412
sage: SR(5).binomial(3, hold=True).simplify()
1413
10
1414
1415
TESTS: We verify that we can convert this function to Maxima and
1416
bring it back into Sage.
1417
1418
::
1419
1420
sage: n,k = var('n,k')
1421
sage: maxima(binomial(n,k))
1422
binomial(n,k)
1423
sage: _.sage()
1424
binomial(n, k)
1425
sage: sage.functions.other.binomial._maxima_init_() # temporary workaround until we can get symbolic binomial to import in global namespace, if that's desired
1426
'binomial'
1427
1428
Test pickling::
1429
1430
sage: loads(dumps(binomial(n,k)))
1431
binomial(n, k)
1432
"""
1433
GinacFunction.__init__(self, "binomial", nargs=2,
1434
conversions=dict(maxima='binomial', mathematica='Binomial'))
1435
1436
binomial = Function_binomial()
1437
1438
class Function_beta(GinacFunction):
1439
def __init__(self):
1440
r"""
1441
Return the beta function. This is defined by
1442
1443
.. math::
1444
1445
B(p,q) = \int_0^1 t^{p-1}(1-t)^{1-q} dt
1446
1447
for complex or symbolic input `p` and `q`.
1448
Note that the order of inputs does not matter: `B(p,q)=B(q,p)`.
1449
1450
GiNaC is used to compute `B(p,q)`. However, complex inputs
1451
are not yet handled in general. When GiNaC raises an error on
1452
such inputs, we raise a NotImplementedError.
1453
1454
If either input is 1, GiNaC returns the reciprocal of the
1455
other. In other cases, GiNaC uses one of the following
1456
formulas:
1457
1458
.. math::
1459
1460
B(p,q) = \Gamma(p)\Gamma(q)/\Gamma(p+q)
1461
1462
or
1463
1464
.. math::
1465
1466
B(p,q) = (-1)^q B(1-p-q, q).
1467
1468
1469
For numerical inputs, GiNaC uses the formula
1470
1471
.. math::
1472
1473
B(p,q) = \exp[\log\Gamma(p)+\log\Gamma(q)-\log\Gamma(p+q)]
1474
1475
1476
INPUT:
1477
1478
- ``p`` - number or symbolic expression
1479
1480
- ``q`` - number or symbolic expression
1481
1482
1483
OUTPUT: number or symbolic expression (if input is symbolic)
1484
1485
EXAMPLES::
1486
1487
sage: beta(3,2)
1488
1/12
1489
sage: beta(3,1)
1490
1/3
1491
sage: beta(1/2,1/2)
1492
beta(1/2, 1/2)
1493
sage: beta(-1,1)
1494
-1
1495
sage: beta(-1/2,-1/2)
1496
0
1497
sage: beta(x/2,3)
1498
beta(1/2*x, 3)
1499
sage: beta(.5,.5)
1500
3.14159265358979
1501
sage: beta(1,2.0+I)
1502
0.400000000000000 - 0.200000000000000*I
1503
sage: beta(3,x+I)
1504
beta(x + I, 3)
1505
1506
Note that the order of arguments does not matter::
1507
1508
sage: beta(1/2,3*x)
1509
beta(3*x, 1/2)
1510
1511
The result is symbolic if exact input is given::
1512
1513
sage: beta(2,1+5*I)
1514
beta(2, 5*I + 1)
1515
sage: beta(2, 2.)
1516
0.166666666666667
1517
sage: beta(I, 2.)
1518
-0.500000000000000 - 0.500000000000000*I
1519
sage: beta(2., 2)
1520
0.166666666666667
1521
sage: beta(2., I)
1522
-0.500000000000000 - 0.500000000000000*I
1523
1524
Test pickling::
1525
1526
sage: loads(dumps(beta))
1527
beta
1528
"""
1529
GinacFunction.__init__(self, "beta", nargs=2,
1530
conversions=dict(maxima='beta', mathematica='Beta'))
1531
1532
beta = Function_beta()
1533
1534
def _do_sqrt(x, prec=None, extend=True, all=False):
1535
r"""
1536
Used internally to compute the square root of x.
1537
1538
INPUT:
1539
1540
- ``x`` - a number
1541
1542
- ``prec`` - None (default) or a positive integer
1543
(bits of precision) If not None, then compute the square root
1544
numerically to prec bits of precision.
1545
1546
- ``extend`` - bool (default: True); this is a place
1547
holder, and is always ignored since in the symbolic ring everything
1548
has a square root.
1549
1550
- ``extend`` - bool (default: True); whether to extend
1551
the base ring to find roots. The extend parameter is ignored if
1552
prec is a positive integer.
1553
1554
- ``all`` - bool (default: False); whether to return
1555
a list of all the square roots of x.
1556
1557
1558
EXAMPLES::
1559
1560
sage: from sage.functions.other import _do_sqrt
1561
sage: _do_sqrt(3)
1562
sqrt(3)
1563
sage: _do_sqrt(3,prec=10)
1564
1.7
1565
sage: _do_sqrt(3,prec=100)
1566
1.7320508075688772935274463415
1567
sage: _do_sqrt(3,all=True)
1568
[sqrt(3), -sqrt(3)]
1569
1570
Note that the extend parameter is ignored in the symbolic ring::
1571
1572
sage: _do_sqrt(3,extend=False)
1573
sqrt(3)
1574
"""
1575
from sage.rings.all import RealField, ComplexField
1576
if prec:
1577
if x >= 0:
1578
return RealField(prec)(x).sqrt(all=all)
1579
else:
1580
return ComplexField(prec)(x).sqrt(all=all)
1581
if x == -1:
1582
from sage.symbolic.pynac import I
1583
z = I
1584
else:
1585
z = SR(x) ** one_half
1586
1587
if all:
1588
if z:
1589
return [z, -z]
1590
else:
1591
return [z]
1592
return z
1593
1594
def sqrt(x, *args, **kwds):
1595
r"""
1596
INPUT:
1597
1598
- ``x`` - a number
1599
1600
- ``prec`` - integer (default: None): if None, returns
1601
an exact square root; otherwise returns a numerical square root if
1602
necessary, to the given bits of precision.
1603
1604
- ``extend`` - bool (default: True); this is a place
1605
holder, and is always ignored or passed to the sqrt function for x,
1606
since in the symbolic ring everything has a square root.
1607
1608
- ``all`` - bool (default: False); if True, return all
1609
square roots of self, instead of just one.
1610
1611
EXAMPLES::
1612
1613
sage: sqrt(-1)
1614
I
1615
sage: sqrt(2)
1616
sqrt(2)
1617
sage: sqrt(2)^2
1618
2
1619
sage: sqrt(4)
1620
2
1621
sage: sqrt(4,all=True)
1622
[2, -2]
1623
sage: sqrt(x^2)
1624
sqrt(x^2)
1625
sage: sqrt(2).n()
1626
1.41421356237310
1627
1628
To prevent automatic evaluation, one can use the ``hold`` parameter
1629
after coercing to the symbolic ring::
1630
1631
sage: sqrt(SR(4),hold=True)
1632
sqrt(4)
1633
sage: sqrt(4,hold=True)
1634
Traceback (most recent call last):
1635
...
1636
TypeError: _do_sqrt() got an unexpected keyword argument 'hold'
1637
1638
This illustrates that the bug reported in #6171 has been fixed::
1639
1640
sage: a = 1.1
1641
sage: a.sqrt(prec=100) # this is supposed to fail
1642
Traceback (most recent call last):
1643
...
1644
TypeError: sqrt() got an unexpected keyword argument 'prec'
1645
sage: sqrt(a, prec=100)
1646
1.0488088481701515469914535137
1647
sage: sqrt(4.00, prec=250)
1648
2.0000000000000000000000000000000000000000000000000000000000000000000000000
1649
1650
One can use numpy input as well::
1651
1652
sage: import numpy
1653
sage: a = numpy.arange(2,5)
1654
sage: sqrt(a)
1655
array([ 1.41421356, 1.73205081, 2. ])
1656
"""
1657
if isinstance(x, float):
1658
return math.sqrt(x)
1659
elif type(x).__module__ == 'numpy':
1660
from numpy import sqrt
1661
return sqrt(x)
1662
try:
1663
return x.sqrt(*args, **kwds)
1664
# The following includes TypeError to catch cases where sqrt
1665
# is called with a "prec" keyword, for example, but the sqrt
1666
# method for x doesn't accept such a keyword.
1667
except (AttributeError, TypeError):
1668
pass
1669
return _do_sqrt(x, *args, **kwds)
1670
1671
# register sqrt in pynac symbol_table for conversion back from other systems
1672
register_symbol(sqrt, dict(mathematica='Sqrt'))
1673
symbol_table['functions']['sqrt'] = sqrt
1674
1675
Function_sqrt = type('deprecated_sqrt', (),
1676
{'__call__': staticmethod(sqrt),
1677
'__setstate__': lambda x, y: None})
1678
1679
class Function_arg(BuiltinFunction):
1680
def __init__(self):
1681
r"""
1682
The argument function for complex numbers.
1683
1684
EXAMPLES::
1685
1686
sage: arg(3+i)
1687
arctan(1/3)
1688
sage: arg(-1+i)
1689
3/4*pi
1690
sage: arg(2+2*i)
1691
1/4*pi
1692
sage: arg(2+x)
1693
arg(x + 2)
1694
sage: arg(2.0+i+x)
1695
arg(x + 2.00000000000000 + 1.00000000000000*I)
1696
sage: arg(-3)
1697
pi
1698
sage: arg(3)
1699
0
1700
sage: arg(0)
1701
0
1702
sage: latex(arg(x))
1703
{\rm arg}\left(x\right)
1704
sage: maxima(arg(x))
1705
atan2(0,x)
1706
sage: maxima(arg(2+i))
1707
atan(1/2)
1708
sage: maxima(arg(sqrt(2)+i))
1709
atan(1/sqrt(2))
1710
sage: arg(2+i)
1711
arctan(1/2)
1712
sage: arg(sqrt(2)+i)
1713
arg(sqrt(2) + I)
1714
sage: arg(sqrt(2)+i).simplify()
1715
arctan(1/2*sqrt(2))
1716
1717
TESTS::
1718
1719
sage: arg(0.0)
1720
0.000000000000000
1721
sage: arg(3.0)
1722
0.000000000000000
1723
sage: arg(-2.5)
1724
3.14159265358979
1725
sage: arg(2.0+3*i)
1726
0.982793723247329
1727
"""
1728
BuiltinFunction.__init__(self, "arg",
1729
conversions=dict(maxima='carg', mathematica='Arg'))
1730
1731
def _eval_(self, x):
1732
"""
1733
EXAMPLES::
1734
1735
sage: arg(3+i)
1736
arctan(1/3)
1737
sage: arg(-1+i)
1738
3/4*pi
1739
sage: arg(2+2*i)
1740
1/4*pi
1741
sage: arg(2+x)
1742
arg(x + 2)
1743
sage: arg(2.0+i+x)
1744
arg(x + 2.00000000000000 + 1.00000000000000*I)
1745
sage: arg(-3)
1746
pi
1747
sage: arg(3)
1748
0
1749
sage: arg(0)
1750
0
1751
sage: arg(sqrt(2)+i)
1752
arg(sqrt(2) + I)
1753
1754
"""
1755
if not isinstance(x,Expression): # x contains no variables
1756
if parent(x)(0)==x: #compatibility with maxima
1757
return parent(x)(0)
1758
else:
1759
if is_inexact(x): # inexact complex numbers, e.g. 2.0+i
1760
return self._evalf_(x, parent(x))
1761
else: # exact complex numbers, e.g. 2+i
1762
return arctan2(imag_part(x),real_part(x))
1763
else:
1764
# x contains variables, e.g. 2+i+y or 2.0+i+y
1765
# or x involves an expression such as sqrt(2)
1766
return None
1767
1768
def _evalf_(self, x, parent_d=None):
1769
"""
1770
EXAMPLES::
1771
1772
sage: arg(0.0)
1773
0.000000000000000
1774
sage: arg(3.0)
1775
0.000000000000000
1776
sage: arg(3.00000000000000000000000000)
1777
0.00000000000000000000000000
1778
sage: arg(3.00000000000000000000000000).prec()
1779
90
1780
sage: arg(ComplexIntervalField(90)(3)).prec()
1781
90
1782
sage: arg(ComplexIntervalField(90)(3)).parent()
1783
Real Interval Field with 90 bits of precision
1784
sage: arg(3.0r)
1785
0.000000000000000
1786
sage: arg(RDF(3))
1787
0.0
1788
sage: arg(RDF(3)).parent()
1789
Real Double Field
1790
sage: arg(-2.5)
1791
3.14159265358979
1792
sage: arg(2.0+3*i)
1793
0.982793723247329
1794
"""
1795
try:
1796
return x.arg()
1797
except AttributeError:
1798
pass
1799
# try to find a parent that support .arg()
1800
if parent_d is None:
1801
parent_d = parent(x)
1802
try:
1803
parent_d = parent_d.complex_field()
1804
except AttributeError:
1805
from sage.rings.complex_field import ComplexField
1806
try:
1807
parent_d = ComplexField(x.prec())
1808
except AttributeError:
1809
parent_d = ComplexField()
1810
1811
return parent_d(x).arg()
1812
1813
arg=Function_arg()
1814
1815
1816
############################
1817
# Real and Imaginary Parts #
1818
############################
1819
class Function_real_part(GinacFunction):
1820
def __init__(self):
1821
r"""
1822
Returns the real part of the (possibly complex) input.
1823
1824
It is possible to prevent automatic evaluation using the
1825
``hold`` parameter::
1826
1827
sage: real_part(I,hold=True)
1828
real_part(I)
1829
1830
To then evaluate again, we currently must use Maxima via
1831
:meth:`sage.symbolic.expression.Expression.simplify`::
1832
1833
sage: real_part(I,hold=True).simplify()
1834
0
1835
1836
EXAMPLES::
1837
1838
sage: z = 1+2*I
1839
sage: real(z)
1840
1
1841
sage: real(5/3)
1842
5/3
1843
sage: a = 2.5
1844
sage: real(a)
1845
2.50000000000000
1846
sage: type(real(a))
1847
<type 'sage.rings.real_mpfr.RealLiteral'>
1848
sage: real(1.0r)
1849
1.0
1850
1851
TESTS::
1852
1853
sage: loads(dumps(real_part))
1854
real_part
1855
1856
Check if #6401 is fixed::
1857
1858
sage: latex(x.real())
1859
\Re \left( x \right)
1860
1861
sage: f(x) = function('f',x)
1862
sage: latex( f(x).real())
1863
\Re \left( f\left(x\right) \right)
1864
"""
1865
GinacFunction.__init__(self, "real_part",
1866
conversions=dict(maxima='realpart'))
1867
1868
def _eval_numpy_(self, x):
1869
"""
1870
EXAMPLES::
1871
1872
sage: import numpy
1873
sage: a = numpy.array([1+2*I, -2-3*I], dtype=numpy.complex)
1874
sage: real_part(a)
1875
array([ 1., -2.])
1876
"""
1877
import numpy
1878
return numpy.real(x)
1879
1880
real = real_part = Function_real_part()
1881
1882
class Function_imag_part(GinacFunction):
1883
def __init__(self):
1884
r"""
1885
Returns the imaginary part of the (possibly complex) input.
1886
1887
It is possible to prevent automatic evaluation using the
1888
``hold`` parameter::
1889
1890
sage: imag_part(I,hold=True)
1891
imag_part(I)
1892
1893
To then evaluate again, we currently must use Maxima via
1894
:meth:`sage.symbolic.expression.Expression.simplify`::
1895
1896
sage: imag_part(I,hold=True).simplify()
1897
1
1898
1899
TESTS::
1900
1901
sage: z = 1+2*I
1902
sage: imaginary(z)
1903
2
1904
sage: imag(z)
1905
2
1906
sage: loads(dumps(imag_part))
1907
imag_part
1908
1909
Check if #6401 is fixed::
1910
1911
sage: latex(x.imag())
1912
\Im \left( x \right)
1913
1914
sage: f(x) = function('f',x)
1915
sage: latex( f(x).imag())
1916
\Im \left( f\left(x\right) \right)
1917
"""
1918
GinacFunction.__init__(self, "imag_part",
1919
conversions=dict(maxima='imagpart'))
1920
1921
def _eval_numpy_(self, x):
1922
"""
1923
EXAMPLES::
1924
1925
sage: import numpy
1926
sage: a = numpy.array([1+2*I, -2-3*I], dtype=numpy.complex)
1927
sage: imag_part(a)
1928
array([ 2., -3.])
1929
"""
1930
import numpy
1931
return numpy.imag(x)
1932
1933
imag = imag_part = imaginary = Function_imag_part()
1934
1935
1936
############################
1937
# Complex Conjugate #
1938
############################
1939
class Function_conjugate(GinacFunction):
1940
def __init__(self):
1941
r"""
1942
Returns the complex conjugate of the input.
1943
1944
It is possible to prevent automatic evaluation using the
1945
``hold`` parameter::
1946
1947
sage: conjugate(I,hold=True)
1948
conjugate(I)
1949
1950
To then evaluate again, we currently must use Maxima via
1951
:meth:`sage.symbolic.expression.Expression.simplify`::
1952
1953
sage: conjugate(I,hold=True).simplify()
1954
-I
1955
1956
TESTS::
1957
1958
sage: x,y = var('x,y')
1959
sage: x.conjugate()
1960
conjugate(x)
1961
sage: latex(conjugate(x))
1962
\overline{x}
1963
sage: f = function('f')
1964
sage: latex(f(x).conjugate())
1965
\overline{f\left(x\right)}
1966
sage: f = function('psi',x,y)
1967
sage: latex(f.conjugate())
1968
\overline{\psi\left(x, y\right)}
1969
sage: x.conjugate().conjugate()
1970
x
1971
sage: x.conjugate().operator()
1972
conjugate
1973
sage: x.conjugate().operator() == conjugate
1974
True
1975
1976
Check if #8755 is fixed::
1977
1978
sage: conjugate(sqrt(-3))
1979
conjugate(sqrt(-3))
1980
sage: conjugate(sqrt(3))
1981
sqrt(3)
1982
sage: conjugate(sqrt(x))
1983
conjugate(sqrt(x))
1984
sage: conjugate(x^2)
1985
conjugate(x)^2
1986
sage: var('y',domain='positive')
1987
y
1988
sage: conjugate(sqrt(y))
1989
sqrt(y)
1990
1991
Check if #10964 is fixed::
1992
1993
sage: z= I*sqrt(-3); z
1994
I*sqrt(-3)
1995
sage: conjugate(z)
1996
-I*conjugate(sqrt(-3))
1997
sage: var('a')
1998
a
1999
sage: conjugate(a*sqrt(-2)*sqrt(-3))
2000
conjugate(a)*conjugate(sqrt(-3))*conjugate(sqrt(-2))
2001
2002
Test pickling::
2003
2004
sage: loads(dumps(conjugate))
2005
conjugate
2006
"""
2007
GinacFunction.__init__(self, "conjugate")
2008
2009
conjugate = Function_conjugate()
2010
2011