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