1# cython: language_level=3
2
3"""
4Mathematical Constants
5
6Numeric, Arithmetic, or Symbolic constants like Pi, E, or Infinity.
7"""
8
9import math
10import mpmath
11import numpy
12import sympy
13
14from mathics.version import __version__  # noqa used in loading to check consistency.
15
16from mathics.builtin.base import Predefined, SympyObject
17from mathics.core.expression import (
18    MachineReal,
19    PrecisionReal,
20    Symbol,
21    strip_context,
22)
23from mathics.core.numbers import get_precision, PrecisionValueError, machine_precision
24
25
26def mp_constant(fn: str, d=None) -> mpmath.mpf:
27    """
28    Return the mpmath constant _fn_ with integer precision _d_.
29    """
30    if d is None:
31        return getattr(mpmath, fn)()
32    else:
33        # TODO: In some functions like Pi, you can
34        # ask for a certain number of digits, but the
35        # accuracy will be less than that. Figure out
36        # what's up and compensate somehow.
37        mpmath.mp.dps = int_d = int(d * 3.321928)
38        return getattr(mpmath, fn)(prec=int_d)
39
40
41def mp_convert_constant(obj, **kwargs):
42    if isinstance(obj, mpmath.ctx_mp_python._constant):
43        prec = kwargs.get("prec", None)
44        if prec is not None:
45            return sympy.Float(obj(prec=prec))
46        return sympy.Float(obj)
47    return obj
48
49
50def numpy_constant(name: str, d=None) -> float:
51    if d:
52        # by mmatera: Here I have a question:
53        # 0.0123`2 should be rounded to
54        # 0.01 or to 0.0123?
55        # (absolute versus relative accuracy)
56        val = getattr(numpy, name)
57        val = numpy.round(val, d)
58        return val
59    else:
60        return getattr(numpy, name)
61
62
63def sympy_constant(fn, d=None):
64    return getattr(sympy, fn).evalf(n=d)
65
66
67class _Constant_Common(Predefined):
68
69    attributes = ("Constant", "Protected", "ReadProtected")
70    nargs = 0
71    options = {"Method": "Automatic"}
72
73    def apply_N(self, precision, evaluation, options={}):
74        "N[%(name)s, precision_?NumericQ, OptionsPattern[%(name)s]]"
75
76        preference = self.get_option(options, "Method", evaluation).get_string_value()
77        if preference == "Automatic":
78            return self.get_constant(precision, evaluation)
79        else:
80            return self.get_constant(precision, evaluation, preference)
81
82    def apply_N2(self, evaluation, options={}):
83        "N[%(name)s, OptionsPattern[%(name)s]]"
84        return self.apply_N(None, evaluation, options)
85
86    def is_constant(self) -> bool:
87        return True
88
89    def get_constant(self, precision, evaluation, preference=None):
90        # first, determine the precision
91        machine_d = int(0.30103 * machine_precision)
92        d = None
93        if precision:
94            try:
95                d = get_precision(precision, evaluation)
96            except PrecisionValueError:
97                pass
98
99        if d is None:
100            d = machine_d
101
102        # If preference not especified, determine it
103        # from the precision.
104        if preference is None:
105            if d <= machine_d:
106                preference = "numpy"
107            else:
108                preference = "mpmath"
109        # If preference is not valid, send a message and return.
110        if not (preference in ("sympy", "numpy", "mpmath")):
111            evaluation.message(f'{preference} not in ("sympy", "numpy", "mpmath")')
112            return
113        # Try to determine the numeric value
114        value = None
115        if preference == "mpmath" and not hasattr(self, "mpmath_name"):
116            preference = "numpy"
117        elif preference == "sympy" and not hasattr(self, "sympy_name"):
118            preference = "numpy"
119
120        if preference == "numpy" and not hasattr(self, "numpy_name"):
121            if hasattr(self, "sympy_name"):
122                preference = "sympy"
123            elif hasattr(self, "mpmath_name"):
124                preference = "mpmath"
125            else:
126                preference = ""
127        if preference == "numpy":
128            value = numpy_constant(self.numpy_name)
129            if d == machine_d:
130                return MachineReal(value)
131        if preference == "sympy":
132            value = sympy_constant(self.sympy_name, d + 2)
133        if preference == "mpmath":
134            value = mp_constant(self.mpmath_name, d * 2)
135        if value:
136            return PrecisionReal(sympy.Float(str(value), d))
137        # If the value is not available, return none
138        # and keep it unevaluated.
139        return
140
141
142class MPMathConstant(_Constant_Common):
143    """Representation of a constant in mpmath, e.g. Pi, E, I, etc."""
144
145    # Subclasses should define this.
146    mpmath_name = None
147
148    mathics_to_mpmath = {}
149
150    def __init__(self, *args, **kwargs):
151        super().__init__(*args, **kwargs)
152        if self.mpmath_name is None:
153            self.mpmath_name = strip_context(self.get_name()).lower()
154        self.mathics_to_mpmath[self.__class__.__name__] = self.mpmath_name
155
156    def to_mpmath(self, args):
157        if self.mpmath_name is None or len(args) != 0:
158            return None
159        return getattr(mpmath, self.mpmath_name)
160
161
162class NumpyConstant(_Constant_Common):
163    """Representation of a constant in numpy, e.g. Pi, E, etc."""
164
165    # Subclasses should define this.
166    numpy_name = None
167
168    mathics_to_numpy = {}
169
170    def __init__(self, *args, **kwargs):
171        super().__init__(*args, **kwargs)
172        if self.numpy_name is None:
173            self.numpy_name = strip_context(self.get_name()).lower()
174        self.mathics_to_numpy[self.__class__.__name__] = self.numpy_name
175
176    def to_numpy(self, args):
177        if self.numpy_name is None or len(args) != 0:
178            return None
179        return self.get_constant()
180
181
182class SympyConstant(_Constant_Common, SympyObject):
183    """Representation of a constant in Sympy, e.g. Pi, E, I, Catalan, etc."""
184
185    # Subclasses should define this.
186    sympy_name = None
187
188    def to_sympy(self, expr=None, **kwargs):
189        if expr is None or expr.is_atom():
190            result = getattr(sympy, self.sympy_name)
191            if kwargs.get("evaluate", False):
192                result = mp_convert_constant(result, **kwargs)
193            return result
194        else:
195            # there is no "native" SymPy expression for e.g. E[x]
196            return None
197
198
199class Catalan(MPMathConstant, SympyConstant):
200    """
201    <dl>
202    <dt>'Catalan'
203        <dd>is Catalan's constant with numerical value \u2243 0.915966.
204    </dl>
205
206    >> Catalan // N
207     = 0.915965594177219
208
209    >> N[Catalan, 20]
210     = 0.91596559417721901505
211    """
212
213    mpmath_name = "catalan"
214    # numpy_name = "catalan"  ## This is not defined in numpy
215    sympy_name = "Catalan"
216
217
218class ComplexInfinity(SympyConstant):
219    """
220    <dl>
221    <dt>'ComplexInfinity'
222        <dd>represents an infinite complex quantity of undetermined direction.
223    </dl>
224
225    >> 1 / ComplexInfinity
226     = 0
227    >> ComplexInfinity * Infinity
228     = ComplexInfinity
229    >> FullForm[ComplexInfinity]
230     = DirectedInfinity[]
231
232    ## Issue689
233    #> ComplexInfinity + ComplexInfinity
234     : Indeterminate expression ComplexInfinity + ComplexInfinity encountered.
235     = Indeterminate
236    #> ComplexInfinity + Infinity
237     : Indeterminate expression ComplexInfinity + Infinity encountered.
238     = Indeterminate
239    """
240
241    sympy_name = "zoo"
242
243    rules = {
244        "ComplexInfinity": "DirectedInfinity[]",
245    }
246
247
248class Degree(MPMathConstant, NumpyConstant, SympyConstant):
249    """
250    <dl>
251      <dt>'Degree'
252      <dd>is the number of radians in one degree. It has a numerical value of \u03c0 / 180.
253    </dl>
254    >> Cos[60 Degree]
255     = 1 / 2
256
257    Degree has the value of Pi / 180
258    >> Degree == Pi / 180
259     = True
260
261    #> Cos[Degree[x]]
262     = Cos[Degree[x]]
263
264    ## Issue 274
265    #> \\[Degree] == ° == Degree
266     = True
267
268    #> N[Degree]
269     = 0.0174533
270    #> N[Degree, 30]
271     = 0.0174532925199432957692369076849
272    """
273
274    mpmath_name = "degree"
275
276    def to_sympy(self, expr=None, **kwargs):
277        if expr == Symbol("System`Degree"):
278            # return mpmath.degree
279            return sympy.pi / 180
280
281    def to_numpy(self, expr=None, **kwargs):
282        if expr == Symbol("System`Degree"):
283            # return mpmath.degree
284            return numpy.pi / 180
285
286    def apply_N(self, precision, evaluation, options={}):
287        "N[Degree, precision_, OptionsPattern[%(name)s]]"
288        try:
289            if precision:
290                d = get_precision(precision, evaluation)
291            else:
292                d = get_precision(Symbol("System`MachinePrecision"), evaluation)
293        except PrecisionValueError:
294            return
295
296        # FIXME: There are all sorts of interactions between in the trig functions,
297        # that are expected to work out right. Until we have convertion between
298        # mpmath and sympy worked out so that values can be made the to the same
299        # precision and compared. we have to not use mpmath right now.
300        # return self.get_constant(precision, evaluation, preference="mpmath")
301
302        if d is None:
303            return MachineReal(math.pi / 180)
304        else:
305            return PrecisionReal((sympy.pi / 180).n(d))
306
307
308class E(MPMathConstant, NumpyConstant, SympyConstant):
309    """
310        <dl>
311        <dt>'E'</dt>
312            <dd>is the constant \u2147 with numerical value \u2243 2.71828.
313        </dl>
314
315        >> N[E]
316         = 2.71828
317        >> N[E, 50]
318         = 2.7182818284590452353602874713526624977572470937000
319
320        #> 5. E
321         = 13.5914
322    """
323
324    mpmath_name = "e"
325    numpy_name = "e"
326    sympy_name = "E"
327
328    def apply_N(self, precision, evaluation, options={}):
329        "N[E, precision_, OptionsPattern[%(name)s]]"
330        return self.get_constant(precision, evaluation)
331
332
333class EulerGamma(MPMathConstant, NumpyConstant, SympyConstant):
334    """
335    <dl>
336      <dt>'EulerGamma'</dt>
337      <dd>is Euler's constant \u03b3 with numerial value \u2243 0.577216.
338    </dl>
339
340    >> EulerGamma // N
341     = 0.577216
342
343    >> N[EulerGamma, 40]
344     = 0.5772156649015328606065120900824024310422
345    """
346
347    mpmath_name = "euler"
348    numpy_name = "euler_gamma"
349    sympy_name = "EulerGamma"
350
351
352class Glaisher(MPMathConstant):
353    """
354    <dl>
355      <dt>'Glaisher'</dt>
356      <dd>is Glaisher's constant, with numerical value \u2243 1.28243.
357    </dl>
358
359    >> N[Glaisher]
360     = 1.28242712910062
361    >> N[Glaisher, 50]
362     = 1.2824271291006226368753425688697917277676889273250
363     # 1.2824271291006219541941391071304678916931152343750
364    """
365
366    mpmath_name = "glaisher"
367
368
369class GoldenRatio(MPMathConstant, SympyConstant):
370    """
371    <dl>
372      <dt>'GoldenRatio'
373      <dd>is the golden ratio, \u03D5 = (1+Sqrt[5])/2.
374    </dl>
375
376    >> GoldenRatio // N
377     = 1.61803398874989
378    >> N[GoldenRatio, 40]
379     = 1.618033988749894848204586834365638117720
380    """
381
382    sympy_name = "GoldenRatio"
383    mpmath_name = "phi"
384
385
386class Indeterminate(SympyConstant):
387    """
388    <dl>
389    <dt>'Indeterminate'</dt>
390        <dd>represents an indeterminate result.
391    </dl>
392
393    >> 0^0
394     : Indeterminate expression 0 ^ 0 encountered.
395     = Indeterminate
396
397    >> Tan[Indeterminate]
398     = Indeterminate
399    """
400
401    sympy_name = "nan"
402
403
404class Infinity(SympyConstant):
405    """
406    <dl>
407    <dt>'Infinity'
408        <dd>represents an infinite real quantity.
409    </dl>
410
411    >> 1 / Infinity
412     = 0
413    >> Infinity + 100
414     = Infinity
415
416    Use 'Infinity' in sum and limit calculations:
417    >> Sum[1/x^2, {x, 1, Infinity}]
418     = Pi ^ 2 / 6
419
420    #> FullForm[Infinity]
421     = DirectedInfinity[1]
422    #> (2 + 3.5*I) / Infinity
423     = 0.
424    #> Infinity + Infinity
425     = Infinity
426    #> Infinity / Infinity
427     : Indeterminate expression 0 Infinity encountered.
428     = Indeterminate
429    """
430
431    sympy_name = "oo"
432    numpy_name = "Inf"
433    mpmath_name = "inf"
434    python_equivalent = math.inf
435
436    rules = {
437        "Infinity": "DirectedInfinity[1]",
438        "MakeBoxes[Infinity, f:StandardForm|TraditionalForm]": ('"\\[Infinity]"'),
439    }
440
441
442class Khinchin(MPMathConstant):
443    """
444    <dl>
445      <dt>'Khinchin'</dt>
446      <dd>is Khinchin's constant, with numerical value \u2243 2.68545.
447    </dl>
448
449    >> N[Khinchin]
450     = 2.68545200106531
451    >> N[Khinchin, 50]
452     = 2.6854520010653064453097148354817956938203822939945
453     # = 2.6854520010653075701156922150403261184692382812500
454    """
455
456    mpmath_name = "khinchin"
457
458
459class Pi(MPMathConstant, SympyConstant):
460    """
461    <dl>
462      <dt>'Pi'</dt>
463      <dd>is the constant \u03c0.
464    </dl>
465
466    >> N[Pi]
467     = 3.14159
468
469    Force using the value given from numpy to compute Pi.
470    >> N[Pi, Method->"numpy"]
471     = 3.14159
472
473    Force using the value given from sympy to compute Pi to 3 places,
474    two places after the decimal point.
475
476    Note that sympy is the default method.
477    >> N[Pi, 3, Method->"sympy"]
478     = 3.14
479
480     >> N[Pi, 50]
481     = 3.1415926535897932384626433832795028841971693993751
482    >> Attributes[Pi]
483     = {Constant, Protected, ReadProtected}
484    """
485
486    sympy_name = "pi"
487    mpmath_name = "pi"
488    numpy_name = "pi"
489