text: # cython: language_level=3
  # -*- coding: utf-8 -*-
  
  # Note: docstring is flowed in documentation. Line breaks in the docstring will appear in the
  # printed output, so be carful not to add then mid-sentence.
  
  """
  Numeric Evaluation and Precision
  
  Support for numeric evaluation with arbitrary precision is just a proof-of-concept.
  
  Precision is not "guarded" through the evaluation process. Only integer precision is supported.
  However, things like 'N[Pi, 100]' should work as expected.
  """
  from mathics.version import __version__  # noqa used in loading to check consistency.
  import sympy
  import mpmath
  import numpy as np
  import math
  import hashlib
  import zlib
  from collections import namedtuple
  from contextlib import contextmanager
  from itertools import chain, product
  from functools import lru_cache
  
  
  from mathics.builtin.base import Builtin, Predefined
  from mathics.core.numbers import (
      dps,
      convert_int_to_digit_list,
      machine_precision,
      machine_epsilon,
      get_precision,
      PrecisionValueError,
  )
  from mathics.core.expression import (
      Complex,
      Expression,
      Integer,
      Integer0,
      MachineReal,
      Number,
      Rational,
      Real,
      String,
      Symbol,
      SymbolFalse,
      SymbolTrue,
      SymbolList,
      SymbolN,
      from_python,
  )
  from mathics.core.convert import from_sympy
  
  
  @lru_cache(maxsize=1024)
  def log_n_b(py_n, py_b) -> int:
      return int(mpmath.ceil(mpmath.log(py_n, py_b))) if py_n != 0 and py_n != 1 else 1
  
  
  class N(Builtin):
      """
      <dl>
      <dt>'N[$expr$, $prec$]'
          <dd>evaluates $expr$ numerically with a precision of $prec$ digits.
      </dl>
      >> N[Pi, 50]
       = 3.1415926535897932384626433832795028841971693993751
  
      >> N[1/7]
       = 0.142857
  
      >> N[1/7, 5]
       = 0.14286
  
      You can manually assign numerical values to symbols.
      When you do not specify a precision, 'MachinePrecision' is taken.
      >> N[a] = 10.9
       = 10.9
      >> a
       = a
  
      'N' automatically threads over expressions, except when a symbol has
       attributes 'NHoldAll', 'NHoldFirst', or 'NHoldRest'.
      >> N[a + b]
       = 10.9 + b
      >> N[a, 20]
       = a
      >> N[a, 20] = 11;
      >> N[a + b, 20]
       = 11.000000000000000000 + b
      >> N[f[a, b]]
       = f[10.9, b]
      >> SetAttributes[f, NHoldAll]
      >> N[f[a, b]]
       = f[a, b]
  
      The precision can be a pattern:
      >> N[c, p_?(#>10&)] := p
      >> N[c, 3]
       = c
      >> N[c, 11]
       = 11.000000000
  
      You can also use 'UpSet' or 'TagSet' to specify values for 'N':
      >> N[d] ^= 5;
      However, the value will not be stored in 'UpValues', but
      in 'NValues' (as for 'Set'):
      >> UpValues[d]
       = {}
      >> NValues[d]
       = {HoldPattern[N[d, MachinePrecision]] :> 5}
      >> e /: N[e] = 6;
      >> N[e]
       = 6.
  
      Values for 'N[$expr$]' must be associated with the head of $expr$:
      >> f /: N[e[f]] = 7;
       : Tag f not found or too deep for an assigned rule.
  
      You can use 'Condition':
      >> N[g[x_, y_], p_] := x + y * Pi /; x + y > 3
      >> SetAttributes[g, NHoldRest]
      >> N[g[1, 1]]
       = g[1., 1]
      >> N[g[2, 2]] // InputForm
       = 8.283185307179586
  
      The precision of the result is no higher than the precision of the input
      >> N[Exp[0.1], 100]
       = 1.10517
      >> % // Precision
       = MachinePrecision
      >> N[Exp[1/10], 100]
       = 1.105170918075647624811707826490246668224547194737518718792863289440967966747654302989143318970748654
      >> % // Precision
       = 100.
      >> N[Exp[1.0`20], 100]
       = 2.7182818284590452354
      >> % // Precision
       = 20.
  
      #> p=N[Pi,100]
       = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
      #> ToString[p]
       = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
  
      #> N[1.012345678901234567890123, 20]
       = 1.0123456789012345679
  
      #> N[I, 30]
       = 1.00000000000000000000000000000 I
  
      #> N[1.012345678901234567890123, 50]
       = 1.01234567890123456789012
      #> % // Precision
       = 24.
      """
  
      messages = {
          "precbd": ("Requested precision `1` is not a " + "machine-sized real number."),
          "preclg": (
              "Requested precision `1` is larger than $MaxPrecision. "
              + "Using current $MaxPrecision of `2` instead. "
              + "$MaxPrecision = Infinity specifies that any precision "
              + "should be allowed."
          ),
          "precsm": (
              "Requested precision `1` is smaller than "
              + "$MinPrecision. Using current $MinPrecision of "
              + "`2` instead."
          ),
      }
  
      rules = {
          "N[expr_]": "N[expr, MachinePrecision]",
      }
  
      def apply_with_prec(self, expr, prec, evaluation):
          "N[expr_, prec_]"
  
          try:
              d = get_precision(prec, evaluation)
          except PrecisionValueError:
              return
  
          if expr.get_head_name() in ("System`List", "System`Rule"):
              return Expression(
                  expr.head,
                  *[self.apply_with_prec(leaf, prec, evaluation) for leaf in expr.leaves],
              )
  
          # Special case for the Root builtin
          if expr.has_form("Root", 2):
              return from_sympy(sympy.N(expr.to_sympy(), d))
  
          if isinstance(expr, Number):
              return expr.round(d)
  
          name = expr.get_lookup_name()
          if name != "":
              nexpr = Expression(SymbolN, expr, prec)
              result = evaluation.definitions.get_value(
                  name, "System`NValues", nexpr, evaluation
              )
              if result is not None:
                  if not result.sameQ(nexpr):
                      result = Expression(SymbolN, result, prec).evaluate(evaluation)
                  return result
  
          if expr.is_atom():
              return expr
          else:
              attributes = expr.head.get_attributes(evaluation.definitions)
              if "System`NHoldAll" in attributes:
                  eval_range = ()
              elif "System`NHoldFirst" in attributes:
                  eval_range = range(1, len(expr.leaves))
              elif "System`NHoldRest" in attributes:
                  if len(expr.leaves) > 0:
                      eval_range = (0,)
                  else:
                      eval_range = ()
              else:
                  eval_range = range(len(expr.leaves))
              head = Expression(SymbolN, expr.head, prec).evaluate(evaluation)
              leaves = expr.get_mutable_leaves()
              for index in eval_range:
                  leaves[index] = Expression(SymbolN, leaves[index], prec).evaluate(
                      evaluation
                  )
              return Expression(head, *leaves)
  
  
  def _scipy_interface(integrator, options_map, mandatory=None, adapt_func=None):
      """
      This function provides a proxy for scipy.integrate
      functions, adapting the parameters.
      """
  
      def _scipy_proxy_func_filter(fun, a, b, **opts):
          native_opts = {}
          if mandatory:
              native_opts.update(mandatory)
          for opt, val in opts.items():
              native_opt = options_map.get(opt, None)
              if native_opt:
                  if native_opt[1]:
                      val = native_opt[1](val)
                  native_opts[native_opt[0]] = val
          return adapt_func(integrator(fun, a, b, **native_opts))
  
      def _scipy_proxy_func(fun, a, b, **opts):
          native_opts = {}
          if mandatory:
              native_opts.update(mandatory)
          for opt, val in opts.items():
              native_opt = options_map.get(opt, None)
              if native_opt:
                  if native_opt[1]:
                      val = native_opt[1](val)
                  native_opts[native_opt[0]] = val
          return integrator(fun, a, b, **native_opts)
  
      return _scipy_proxy_func_filter if adapt_func else _scipy_proxy_func
  
  
  def _internal_adaptative_simpsons_rule(f, a, b, **opts):
      """
      1D adaptative Simpson's rule integrator
      Adapted from https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method
         by @mmatera
  
      TODO: handle weak divergences
      """
      wsr = 1.0 / 6.0
  
      tol = opts.get("tol")
      if not tol:
          tol = 1.0e-10
  
      maxrec = opts.get("maxrec")
      if not maxrec:
          maxrec = 150
  
      def _quad_simpsons_mem(f, a, fa, b, fb):
          """Evaluates the Simpson's Rule, also returning m and f(m) to reuse"""
          m = 0.5 * (a + b)
          try:
              fm = f(m)
          except ZeroDivisionError:
              fm = None
  
          if fm is None or np.isinf(fm):
              m = m + 1e-10
              fm = f(m)
          return (m, fm, wsr * abs(b - a) * (fa + 4.0 * fm + fb))
  
      def _quad_asr(f, a, fa, b, fb, eps, whole, m, fm, maxrec):
          """
          Efficient recursive implementation of adaptive Simpson's rule.
          Function values at the start, middle, end of the intervals
          are retained.
          """
          maxrec = maxrec - 1
          try:
              left = _quad_simpsons_mem(f, a, fa, m, fm)
              lm, flm, left = left
              right = _quad_simpsons_mem(f, m, fm, b, fb)
              rm, frm, right = right
  
              delta = left + right - whole
              err = abs(delta)
              if err <= 15 * eps or maxrec == 0:
                  return (left + right + delta / 15, err)
              left = _quad_asr(f, a, fa, m, fm, 0.5 * eps, left, lm, flm, maxrec)
              right = _quad_asr(f, m, fm, b, fb, 0.5 * eps, right, rm, frm, maxrec)
              return (left[0] + right[0], left[1] + right[1])
          except Exception:
              raise
  
      def ensure_evaluation(f, x):
          try:
              val = f(x)
          except ZeroDivisionError:
              return None
          if np.isinf(val):
              return None
          return val
  
      invert_interval = False
      if a > b:
          b, a, invert_interval = a, b, True
  
      fa, fb = ensure_evaluation(f, a), ensure_evaluation(f, b)
      if fa is None:
          x = 10.0 * machine_epsilon if a == 0 else a * (1.0 + 10.0 * machine_epsilon)
          fa = ensure_evaluation(f, x)
          if fa is None:
              raise Exception(f"Function undefined around {a}. Cannot integrate")
      if fb is None:
          x = -10.0 * machine_epsilon if b == 0 else b * (1.0 - 10.0 * machine_epsilon)
          fb = ensure_evaluation(f, x)
          if fb is None:
              raise Exception(f"Function undefined around {b}. Cannot integrate")
  
      m, fm, whole = _quad_simpsons_mem(f, a, fa, b, fb)
      if invert_interval:
          return -_quad_asr(f, a, fa, b, fb, tol, whole, m, fm, maxrec)
      else:
          return _quad_asr(f, a, fa, b, fb, tol, whole, m, fm, maxrec)
  
  
  def _fubini(func, ranges, **opts):
      if not ranges:
          return 0.0
      a, b = ranges[0]
      integrator = opts["integrator"]
      tol = opts.get("tol")
      if tol is None:
          opts["tol"] = 1.0e-10
          tol = 1.0e-10
  
      if len(ranges) > 1:
  
          def subintegral(*u):
              def ff(*z):
                  return func(*(u + z))
  
              val = _fubini(ff, ranges[1:], **opts)[0]
              return val
  
          opts["tol"] = 4.0 * tol
          val = integrator(subintegral, a, b, **opts)
          return val
      else:
          val = integrator(func, a, b, **opts)
          return val
  
  
  class NIntegrate(Builtin):
      """
      <dl>
         <dt>'NIntegrate[$expr$, $interval$]'
         <dd>returns a numeric approximation to the definite integral of $expr$ with limits $interval$ and with a precision of $prec$ digits.
  
          <dt>'NIntegrate[$expr$, $interval1$, $interval2$, ...]'
          <dd>returns a numeric approximation to the multiple integral of $expr$ with limits $interval1$, $interval2$ and with a precision of $prec$ digits.
      </dl>
  
      >> NIntegrate[Exp[-x],{x,0,Infinity},Tolerance->1*^-6]
       = 1.
      >> NIntegrate[Exp[x],{x,-Infinity, 0},Tolerance->1*^-6]
       = 1.
      >> NIntegrate[Exp[-x^2/2.],{x,-Infinity, Infinity},Tolerance->1*^-6]
       = 2.50663
  
      >> Table[1./NIntegrate[x^k,{x,0,1},Tolerance->1*^-6], {k,0,6}]
       : The specified method failed to return a number. Falling back into the internal evaluator.
       = {1., 2., 3., 4., 5., 6., 7.}
  
      >> NIntegrate[1 / z, {z, -1 - I, 1 - I, 1 + I, -1 + I, -1 - I}, Tolerance->1.*^-4]
       : Integration over a complex domain is not implemented yet
       = NIntegrate[1 / z, {z, -1 - I, 1 - I, 1 + I, -1 + I, -1 - I}, Tolerance -> 0.0001]
       ## = 6.2832 I
  
      Integrate singularities with weak divergences:
      >> Table[ NIntegrate[x^(1./k-1.), {x,0,1.}, Tolerance->1*^-6], {k,1,7.} ]
       = {1., 2., 3., 4., 5., 6., 7.}
  
      Mutiple Integrals :
      >> NIntegrate[x * y,{x, 0, 1}, {y, 0, 1}]
       = 0.25
  
      """
  
      messages = {
          "bdmtd": "The Method option should be a built-in method name.",
          "inumr": (
              "The integrand `1` has evaluated to non-numerical "
              + "values for all sampling points in the region "
              + "with boundaries `2`"
          ),
          "nlim": "`1` = `2` is not a valid limit of integration.",
          "ilim": "Invalid integration variable or limit(s) in `1`.",
          "mtdfail": (
              "The specified method failed to return a "
              + "number. Falling back into the internal "
              + "evaluator."
          ),
          "cmpint": ("Integration over a complex domain is not " + "implemented yet"),
      }
  
      options = {
          "Method": '"Automatic"',
          "Tolerance": "1*^-10",
          "Accuracy": "1*^-10",
          "MaxRecursion": "10",
      }
  
      methods = {
          "Automatic": (None, False),
      }
  
      def __init__(self, *args, **kwargs):
          super().__init__(*args, **kwargs)
          self.methods["Internal"] = (_internal_adaptative_simpsons_rule, False)
          try:
              from scipy.integrate import romberg, quad, nquad
              self.methods["NQuadrature"] = (
                  _scipy_interface(
                      nquad, {}, {"full_output": 1}, lambda res: (res[0], res[1])
                  ),
                  True,
              )
              self.methods["Quadrature"] = (
                  _scipy_interface(
                      quad,
                      {
                          "tol": ("epsabs", None),
                          "maxrec": ("limit", lambda maxrec: int(2 ** maxrec)),
                      },
                      {"full_output": 1},
                      lambda res: (res[0], res[1]),
                  ),
                  False,
              )
              self.methods["Romberg"] = (
                  _scipy_interface(
                      romberg,
                      {"tol": ("tol", None), "maxrec": ("divmax", None)},
                      None,
                      lambda x: (x, np.nan),
                  ),
                  False,
              )
              self.methods["Automatic"] = self.methods["Quadrature"]
          except Exception:
              self.methods["Automatic"] = self.methods["Internal"]
              self.methods["Simpson"] = self.methods["Internal"]
  
          self.messages["bdmtd"] = (
              "The Method option should be a "
              + "built-in method name in {`"
              + "`, `".join(list(self.methods))
              + "`}. Using `Automatic`"
          )
  
      @staticmethod
      def decompose_domain(interval, evaluation):
          if interval.has_form("System`Sequence", 1, None):
              intervals = []
              for leaf in interval.leaves:
                  inner_interval = NIntegrate.decompose_domain(leaf, evaluation)
                  if inner_interval:
                      intervals.append(inner_interval)
                  else:
                      evaluation.message("ilim", leaf)
                      return None
              return intervals
  
          if interval.has_form("System`List", 3, None):
              intervals = []
              intvar = interval.leaves[0]
              if not isinstance(intvar, Symbol):
                  evaluation.message("ilim", interval)
                  return None
              boundaries = [a for a in interval.leaves[1:]]
              if any([b.get_head_name() == "System`Complex" for b in boundaries]):
                  intvar = Expression(
                      "List", intvar, Expression("Blank", Symbol("Complex"))
                  )
              for i in range(len(boundaries) - 1):
                  intervals.append((boundaries[i], boundaries[i + 1]))
              if len(intervals) > 0:
                  return (intvar, intervals)
  
          evaluation.message("ilim", interval)
          return None
  
      def apply_with_func_domain(self, func, domain, evaluation, options):
          "%(name)s[func_, domain__, OptionsPattern[%(name)s]]"
          method = options["System`Method"].evaluate(evaluation)
          method_options = {}
          if method.has_form("System`List", 2):
              method = method.leaves[0]
              method_options.update(method.leaves[1].get_option_values())
          if isinstance(method, String):
              method = method.value
          elif isinstance(method, Symbol):
              method = method.get_name()
          else:
              evaluation.message("NIntegrate", "bdmtd", method)
              return
          tolerance = options["System`Tolerance"].evaluate(evaluation)
          tolerance = float(tolerance.value)
          accuracy = options["System`Accuracy"].evaluate(evaluation)
          accuracy = accuracy.value
          maxrecursion = options["System`MaxRecursion"].evaluate(evaluation)
          maxrecursion = maxrecursion.value
          nintegrate_method = self.methods.get(method, None)
          if nintegrate_method is None:
              evaluation.message("NIntegrate", "bdmtd", method)
              nintegrate_method = self.methods.get("Automatic")
          if type(nintegrate_method) is tuple:
              nintegrate_method, is_multidimensional = nintegrate_method
          else:
              is_multidimensional = False
  
          domain = self.decompose_domain(domain, evaluation)
          if not domain:
              return
          if not isinstance(domain, list):
              domain = [domain]
  
          coords = [axis[0] for axis in domain]
          # If any of the points in the integration domain is complex,
          # stop the evaluation...
          if any([c.get_head_name() == "System`List" for c in coords]):
              evaluation.message("NIntegrate", "cmpint")
              return
  
          intvars = Expression(SymbolList, *coords)
          integrand = Expression("Compile", intvars, func).evaluate(evaluation)
  
          if len(integrand.leaves) >= 3:
              integrand = integrand.leaves[2].cfunc
          else:
              evaluation.message("inumer", func, domain)
              return
          results = []
          for subdomain in product(*[axis[1] for axis in domain]):
              # On each subdomain, check if the region is bounded.
              # If not, implement a coordinate map
              func2 = integrand
              subdomain2 = []
              coordtransform = []
              nulldomain = False
              for i, r in enumerate(subdomain):
                  a = r[0].evaluate(evaluation)
                  b = r[1].evaluate(evaluation)
                  if a == b:
                      nulldomain = True
                      break
                  elif a.get_head_name() == "System`DirectedInfinity":
                      if b.get_head_name() == "System`DirectedInfinity":
                          a = a.to_python()
                          b = b.to_python()
                          le = 1 - machine_epsilon
                          if a == b:
                              nulldomain = True
                              break
                          elif a < b:
                              subdomain2.append([-le, le])
                          else:
                              subdomain2.append([le, -le])
                          coordtransform.append(
                              (np.arctanh, lambda u: 1.0 / (1.0 - u ** 2))
                          )
                      else:
                          if not b.is_numeric():
                              evaluation.message("nlim", coords[i], b)
                              return
                          z = a.leaves[0].value
                          b = b.value
                          subdomain2.append([machine_epsilon, 1.0])
                          coordtransform.append(
                              (lambda u: b - z + z / u, lambda u: -z * u ** (-2.0))
                          )
                  elif b.get_head_name() == "System`DirectedInfinity":
                      if not a.is_numeric():
                          evaluation.message("nlim", coords[i], a)
                          return
                      a = a.value
                      z = b.leaves[0].value
                      subdomain2.append([machine_epsilon, 1.0])
                      coordtransform.append(
                          (lambda u: a - z + z / u, lambda u: z * u ** (-2.0))
                      )
                  elif a.is_numeric() and b.is_numeric():
                      a = Expression(SymbolN, a).evaluate(evaluation).value
                      b = Expression(SymbolN, b).evaluate(evaluation).value
                      subdomain2.append([a, b])
                      coordtransform.append(None)
                  else:
                      for x in (a, b):
                          if not x.is_numeric():
                              evaluation.message("nlim", coords[i], x)
                      return
  
              if nulldomain:
                  continue
  
              if any(coordtransform):
                  func2 = lambda *u: (
                      integrand(
                          *[
                              x[0](u[i]) if x else u[i]
                              for i, x in enumerate(coordtransform)
                          ]
                      )
                      * np.prod(
                          [jac[1](u[i]) for i, jac in enumerate(coordtransform) if jac]
                      )
                  )
              opts = {
                  "acur": accuracy,
                  "tol": tolerance,
                  "maxrec": maxrecursion,
              }
              opts.update(method_options)
              try:
                  if len(subdomain2) > 1:
                      if is_multidimensional:
                          nintegrate_method(func2, subdomain2, **opts)
                      else:
                          val = _fubini(
                              func2, subdomain2, integrator=nintegrate_method, **opts
                          )
                  else:
                      val = nintegrate_method(func2, *(subdomain2[0]), **opts)
              except Exception:
                  val = None
  
              if val is None:
                  evaluation.message("NIntegrate", "mtdfail")
                  if len(subdomain2) > 1:
                      val = _fubini(
                          func2,
                          subdomain2,
                          integrator=_internal_adaptative_simpsons_rule,
                          **opts,
                      )
                  else:
                      val = _internal_adaptative_simpsons_rule(
                          func2, *(subdomain2[0]), **opts
                      )
              results.append(val)
  
          result = sum([r[0] for r in results])
          return from_python(result)
  
  
  class MachinePrecision(Predefined):
      """
      <dl>
      <dt>'MachinePrecision'
          <dd>represents the precision of machine precision numbers.
      </dl>
  
      >> N[MachinePrecision]
       = 15.9546
      >> N[MachinePrecision, 30]
       = 15.9545897701910033463281614204
  
      #> N[E, MachinePrecision]
       = 2.71828
  
      #> Round[MachinePrecision]
       = 16
      """
  
      rules = {
          "N[MachinePrecision, prec_]": ("N[Log[10, 2] * %i, prec]" % machine_precision),
      }
  
  
  class MachineEpsilon_(Predefined):
      """
      <dl>
      <dt>'$MachineEpsilon'
          <dd>is the distance between '1.0' and the next
              nearest representable machine-precision number.
      </dl>
  
      >> $MachineEpsilon
       = 2.22045*^-16
  
      >> x = 1.0 + {0.4, 0.5, 0.6} $MachineEpsilon;
      >> x - 1
       = {0., 0., 2.22045*^-16}
      """
  
      name = "$MachineEpsilon"
  
      def evaluate(self, evaluation):
          return MachineReal(machine_epsilon)
  
  
  class MachinePrecision_(Predefined):
      """
      <dl>
      <dt>'$MachinePrecision'
          <dd>is the number of decimal digits of precision for
              machine-precision numbers.
      </dl>
  
      >> $MachinePrecision
       = 15.9546
      """
  
      name = "$MachinePrecision"
  
      rules = {
          "$MachinePrecision": "N[MachinePrecision]",
      }
  
  
  class Precision(Builtin):
      """
      <dl>
      <dt>'Precision[$expr$]'
          <dd>examines the number of significant digits of $expr$.
      </dl>
      This is rather a proof-of-concept than a full implementation.
      Precision of compound expression is not supported yet.
      >> Precision[1]
       = Infinity
      >> Precision[1/2]
       = Infinity
      >> Precision[0.5]
       = MachinePrecision
  
      #> Precision[0.0]
       = MachinePrecision
      #> Precision[0.000000000000000000000000000000000000]
       = 0.
      #> Precision[-0.0]
       = MachinePrecision
      #> Precision[-0.000000000000000000000000000000000000]
       = 0.
  
      #> 1.0000000000000000 // Precision
       = MachinePrecision
      #> 1.00000000000000000 // Precision
       = 17.
  
      #> 0.4 + 2.4 I // Precision
       = MachinePrecision
      #> Precision[2 + 3 I]
       = Infinity
  
      #> Precision["abc"]
       = Infinity
      """
  
      rules = {
          "Precision[z_?MachineNumberQ]": "MachinePrecision",
      }
  
      def apply(self, z, evaluation):
          "Precision[z_]"
  
          if not z.is_inexact():
              return Symbol("Infinity")
          elif z.to_sympy().is_zero:
              return Real(0)
          else:
              return Real(dps(z.get_precision()))
  
  
  class MinPrecision(Builtin):
      """
      <dl>
      <dt>'$MinPrecision'
        <dd>represents the minimum number of digits of precision
            permitted in abitrary-precision numbers.
      </dl>
  
      >> $MinPrecision
       = 0
  
      >> $MinPrecision = 10;
  
      >> N[Pi, 9]
       : Requested precision 9 is smaller than $MinPrecision. Using current $MinPrecision of 10. instead.
       = 3.141592654
  
      #> N[Pi, 10]
       = 3.141592654
  
      #> $MinPrecision = x
       : Cannot set $MinPrecision to x; value must be a non-negative number.
       = x
      #> $MinPrecision = -Infinity
       : Cannot set $MinPrecision to -Infinity; value must be a non-negative number.
       = -Infinity
      #> $MinPrecision = -1
       : Cannot set $MinPrecision to -1; value must be a non-negative number.
       = -1
      #> $MinPrecision = 0;
  
      #> $MaxPrecision = 10;
      #> $MinPrecision = 15
       : Cannot set $MinPrecision such that $MaxPrecision < $MinPrecision.
       = 15
      #> $MinPrecision
       = 0
      #> $MaxPrecision = Infinity;
      """
  
      name = "$MinPrecision"
      rules = {
          "$MinPrecision": "0",
      }
  
      messages = {
          "precset": "Cannot set `1` to `2`; value must be a non-negative number.",
          "preccon": "Cannot set `1` such that $MaxPrecision < $MinPrecision.",
      }
  
  
  class MaxPrecision(Predefined):
      """
      <dl>
      <dt>'$MaxPrecision'
        <dd>represents the maximum number of digits of precision
            permitted in abitrary-precision numbers.
      </dl>
  
      >> $MaxPrecision
       = Infinity
  
      >> $MaxPrecision = 10;
  
      >> N[Pi, 11]
       : Requested precision 11 is larger than $MaxPrecision. Using current $MaxPrecision of 10. instead. $MaxPrecision = Infinity specifies that any precision should be allowed.
       = 3.141592654
  
      #> N[Pi, 10]
       = 3.141592654
  
      #> $MaxPrecision = x
       : Cannot set $MaxPrecision to x; value must be a positive number or Infinity.
       = x
      #> $MaxPrecision = -Infinity
       : Cannot set $MaxPrecision to -Infinity; value must be a positive number or Infinity.
       = -Infinity
      #> $MaxPrecision = 0
       : Cannot set $MaxPrecision to 0; value must be a positive number or Infinity.
       = 0
      #> $MaxPrecision = Infinity;
  
      #> $MinPrecision = 15;
      #> $MaxPrecision = 10
       : Cannot set $MaxPrecision such that $MaxPrecision < $MinPrecision.
       = 10
      #> $MaxPrecision
       = Infinity
      #> $MinPrecision = 0;
      """
  
      name = "$MaxPrecision"
  
      rules = {
          "$MaxPrecision": "Infinity",
      }
  
      messages = {
          "precset": "Cannot set `1` to `2`; value must be a positive number or Infinity.",
          "preccon": "Cannot set `1` such that $MaxPrecision < $MinPrecision.",
      }
  
  
  class Round(Builtin):
      """
      <dl>
      <dt>'Round[$expr$]'
          <dd>rounds $expr$ to the nearest integer.
      <dt>'Round[$expr$, $k$]'
          <dd>rounds $expr$ to the closest multiple of $k$.
      </dl>
  
      >> Round[10.6]
       = 11
      >> Round[0.06, 0.1]
       = 0.1
      >> Round[0.04, 0.1]
       = 0.
  
      Constants can be rounded too
      >> Round[Pi, .5]
       = 3.
      >> Round[Pi^2]
       = 10
  
      Round to exact value
      >> Round[2.6, 1/3]
       = 8 / 3
      >> Round[10, Pi]
       = 3 Pi
  
      Round complex numbers
      >> Round[6/(2 + 3 I)]
       = 1 - I
      >> Round[1 + 2 I, 2 I]
       = 2 I
  
      Round Negative numbers too
      >> Round[-1.4]
       = -1
  
      Expressions other than numbers remain unevaluated:
      >> Round[x]
       = Round[x]
      >> Round[1.5, k]
       = Round[1.5, k]
      """
  
      attributes = ("Listable", "NumericFunction")
  
      rules = {
          "Round[expr_?NumericQ]": "Round[Re[expr], 1] + I * Round[Im[expr], 1]",
          "Round[expr_Complex, k_?RealNumberQ]": (
              "Round[Re[expr], k] + I * Round[Im[expr], k]"
          ),
      }
  
      def apply(self, expr, k, evaluation):
          "Round[expr_?NumericQ, k_?NumericQ]"
  
          n = Expression("Divide", expr, k).round_to_float(
              evaluation, permit_complex=True
          )
          if n is None:
              return
          elif isinstance(n, complex):
              n = round(n.real)
          else:
              n = round(n)
          n = int(n)
          return Expression("Times", Integer(n), k)
  
  
  class Rationalize(Builtin):
      """
      <dl>
      <dt>'Rationalize[$x$]'
          <dd>converts a real number $x$ to a nearby rational number.
      <dt>'Rationalize[$x$, $dx$]'
          <dd>finds the rational number within $dx$ of $x$ with the smallest denominator.
      </dl>
  
      >> Rationalize[2.2]
      = 11 / 5
  
      Not all numbers can be well approximated.
      >> Rationalize[N[Pi]]
       = 3.14159
  
      Find the exact rational representation of 'N[Pi]'
      >> Rationalize[N[Pi], 0]
       = 245850922 / 78256779
  
      #> Rationalize[1.6 + 0.8 I]
       = 8 / 5 + 4 I / 5
  
      #> Rationalize[N[Pi] + 0.8 I, 1*^-6]
       = 355 / 113 + 4 I / 5
  
      #> Rationalize[N[Pi] + 0.8 I, x]
       : Tolerance specification x must be a non-negative number.
       = Rationalize[3.14159 + 0.8 I, x]
  
      #> Rationalize[N[Pi] + 0.8 I, -1]
       : Tolerance specification -1 must be a non-negative number.
       = Rationalize[3.14159 + 0.8 I, -1]
  
      #> Rationalize[N[Pi] + 0.8 I, 0]
       = 245850922 / 78256779 + 4 I / 5
  
      #> Rationalize[17 / 7]
       = 17 / 7
  
      #> Rationalize[x]
       = x
  
      #> Table[Rationalize[E, 0.1^n], {n, 1, 10}]
       = {8 / 3, 19 / 7, 87 / 32, 193 / 71, 1071 / 394, 2721 / 1001, 15062 / 5541, 23225 / 8544, 49171 / 18089, 419314 / 154257}
  
      #> Rationalize[x, y]
       : Tolerance specification y must be a non-negative number.
       = Rationalize[x, y]
      """
  
      messages = {
          "tolnn": "Tolerance specification `1` must be a non-negative number.",
      }
  
      rules = {
          "Rationalize[z_Complex]": "Rationalize[Re[z]] + I Rationalize[Im[z]]",
          "Rationalize[z_Complex, dx_?Internal`RealValuedNumberQ]/;dx >= 0": "Rationalize[Re[z], dx] + I Rationalize[Im[z], dx]",
      }
  
      def apply(self, x, evaluation):
          "Rationalize[x_]"
  
          py_x = x.to_sympy()
          if py_x is None or (not py_x.is_number) or (not py_x.is_real):
              return x
          return from_sympy(self.find_approximant(py_x))
  
      @staticmethod
      def find_approximant(x):
          c = 1e-4
          it = sympy.ntheory.continued_fraction_convergents(
              sympy.ntheory.continued_fraction_iterator(x)
          )
          for i in it:
              p, q = i.as_numer_denom()
              tol = c / q ** 2
              if abs(i - x) <= tol:
                  return i
              if tol < machine_epsilon:
                  break
          return x
  
      @staticmethod
      def find_exact(x):
          p, q = x.as_numer_denom()
          it = sympy.ntheory.continued_fraction_convergents(
              sympy.ntheory.continued_fraction_iterator(x)
          )
          for i in it:
              p, q = i.as_numer_denom()
              if abs(x - i) < machine_epsilon:
                  return i
  
      def apply_dx(self, x, dx, evaluation):
          "Rationalize[x_, dx_]"
          py_x = x.to_sympy()
          if py_x is None:
              return x
          py_dx = dx.to_sympy()
          if (
              py_dx is None
              or (not py_dx.is_number)
              or (not py_dx.is_real)
              or py_dx.is_negative
          ):
              return evaluation.message("Rationalize", "tolnn", dx)
          elif py_dx == 0:
              return from_sympy(self.find_exact(py_x))
          a = self.approx_interval_continued_fraction(py_x - py_dx, py_x + py_dx)
          sym_x = sympy.ntheory.continued_fraction_reduce(a)
          return Rational(sym_x)
  
      @staticmethod
      def approx_interval_continued_fraction(xmin, xmax):
          result = []
          a_gen = sympy.ntheory.continued_fraction_iterator(xmin)
          b_gen = sympy.ntheory.continued_fraction_iterator(xmax)
          while True:
              a, b = next(a_gen), next(b_gen)
              if a == b:
                  result.append(a)
              else:
                  result.append(min(a, b) + 1)
                  break
          return result
  
  
  def chop(expr, delta=10.0 ** (-10.0)):
      if isinstance(expr, Real):
          if expr.is_nan(expr):
              return expr
          if -delta < expr.get_float_value() < delta:
              return Integer0
      elif isinstance(expr, Complex) and expr.is_inexact():
          real, imag = expr.real, expr.imag
          if -delta < real.get_float_value() < delta:
              real = Integer0
          if -delta < imag.get_float_value() < delta:
              imag = Integer0
          return Complex(real, imag)
      elif isinstance(expr, Expression):
          return Expression(chop(expr.head), *[chop(leaf) for leaf in expr.leaves])
      return expr
  
  
  class Chop(Builtin):
      """
      <dl>
      <dt>'Chop[$expr$]'
          <dd>replaces floating point numbers close to 0 by 0.
      <dt>'Chop[$expr$, $delta$]'
          <dd>uses a tolerance of $delta$. The default tolerance is '10^-10'.
      </dl>
  
      >> Chop[10.0 ^ -16]
       = 0
      >> Chop[10.0 ^ -9]
       = 1.*^-9
      >> Chop[10 ^ -11 I]
       = I / 100000000000
      >> Chop[0. + 10 ^ -11 I]
       = 0
      """
  
      messages = {
          "tolnn": "Tolerance specification a must be a non-negative number.",
      }
  
      rules = {
          "Chop[expr_]": "Chop[expr, 10^-10]",
      }
  
      def apply(self, expr, delta, evaluation):
          "Chop[expr_, delta_:(10^-10)]"
  
          delta = delta.round_to_float(evaluation)
          if delta is None or delta < 0:
              return evaluation.message("Chop", "tolnn")
  
          return chop(expr, delta=delta)
  
  
  class NumericQ(Builtin):
      """
      <dl>
      <dt>'NumericQ[$expr$]'
          <dd>tests whether $expr$ represents a numeric quantity.
      </dl>
  
      >> NumericQ[2]
       = True
      >> NumericQ[Sqrt[Pi]]
       = True
      >> NumberQ[Sqrt[Pi]]
       = False
      """
  
      def apply(self, expr, evaluation):
          "NumericQ[expr_]"
  
          def test(expr):
              if isinstance(expr, Expression):
                  attr = evaluation.definitions.get_attributes(expr.head.get_name())
                  return "System`NumericFunction" in attr and all(
                      test(leaf) for leaf in expr.leaves
                  )
              else:
                  return expr.is_numeric()
  
          return SymbolTrue if test(expr) else SymbolFalse
  
  
  class RealValuedNumericQ(Builtin):
      # No docstring since this is internal and it will mess up documentation.
      # FIXME: Perhaps in future we will have a more explicite way to indicate not
      # to add something to the docs.
      context = "Internal`"
  
      rules = {
          "Internal`RealValuedNumericQ[x_]": "Head[N[x]] === Real",
      }
  
  
  class RealValuedNumberQ(Builtin):
      # No docstring since this is internal and it will mess up documentation.
      # FIXME: Perhaps in future we will have a more explicite way to indicate not
      # to add something to the docs.
      context = "Internal`"
  
      rules = {
          "Internal`RealValuedNumberQ[x_Real]": "True",
          "Internal`RealValuedNumberQ[x_Integer]": "True",
          "Internal`RealValuedNumberQ[x_Rational]": "True",
          "Internal`RealValuedNumberQ[x_]": "False",
      }
  
  
  class IntegerDigits(Builtin):
      """
      <dl>
      <dt>'IntegerDigits[$n$]'
          <dd>returns a list of the base-10 digits in the integer $n$.
      <dt>'IntegerDigits[$n$, $base$]'
          <dd>returns a list of the base-$base$ digits in $n$.
      <dt>'IntegerDigits[$n$, $base$, $length$]'
          <dd>returns a list of length $length$, truncating or padding
          with zeroes on the left as necessary.
      </dl>
  
      >> IntegerDigits[76543]
       = {7, 6, 5, 4, 3}
  
      The sign of $n$ is discarded:
      >> IntegerDigits[-76543]
       = {7, 6, 5, 4, 3}
  
      >> IntegerDigits[15, 16]
       = {15}
      >> IntegerDigits[1234, 16]
       = {4, 13, 2}
      >> IntegerDigits[1234, 10, 5]
       = {0, 1, 2, 3, 4}
  
      #> IntegerDigits[1000, 10]
       = {1, 0, 0, 0}
  
      #> IntegerDigits[0]
       = {0}
      """
  
      attributes = ("Listable",)
  
      messages = {
          "int": "Integer expected at position 1 in `1`",
          "ibase": "Base `1` is not an integer greater than 1.",
      }
  
      rules = {
          "IntegerDigits[n_]": "IntegerDigits[n, 10]",
      }
  
      def apply_len(self, n, base, length, evaluation):
          "IntegerDigits[n_, base_, length_]"
  
          if not (isinstance(length, Integer) and length.get_int_value() >= 0):
              return evaluation.message("IntegerDigits", "intnn")
  
          return self.apply(n, base, evaluation, nr_elements=length.get_int_value())
  
      def apply(self, n, base, evaluation, nr_elements=None):
          "IntegerDigits[n_, base_]"
  
          if not (isinstance(n, Integer)):
              return evaluation.message(
                  "IntegerDigits", "int", Expression("IntegerDigits", n, base)
              )
  
          if not (isinstance(base, Integer) and base.get_int_value() > 1):
              return evaluation.message("IntegerDigits", "ibase", base)
  
          if nr_elements == 0:
              # trivial case: we don't want any digits
              return Expression(SymbolList)
  
          digits = convert_int_to_digit_list(n.get_int_value(), base.get_int_value())
  
          if nr_elements is not None:
              if len(digits) >= nr_elements:
                  # Truncate, preserving the digits on the right
                  digits = digits[-nr_elements:]
              else:
                  # Pad with zeroes
                  digits = [0] * (nr_elements - len(digits)) + digits
  
          return Expression(SymbolList, *digits)
  
  
  def check_finite_decimal(denominator):
      # The rational number is finite decimal if the denominator has form 2^a * 5^b
      while denominator % 5 == 0:
          denominator = denominator / 5
  
      while denominator % 2 == 0:
          denominator = denominator / 2
  
      return True if denominator == 1 else False
  
  
  def convert_repeating_decimal(numerator, denominator, base):
      head = [x for x in str(numerator // denominator)]
      tails = []
      subresults = [numerator % denominator]
      numerator %= denominator
  
      while numerator != 0:  # only rational input can go to this case
          numerator *= base
          result_digit, numerator = divmod(numerator, denominator)
          tails.append(str(result_digit))
          if numerator not in subresults:
              subresults.append(numerator)
          else:
              break
  
      for i in range(len(head) - 1, -1, -1):
          j = len(tails) - 1
          if head[i] != tails[j]:
              break
          else:
              del tails[j]
              tails.insert(0, head[i])
              del head[i]
              j = j - 1
  
      # truncate all leading 0's
      if all(elem == "0" for elem in head):
          for i in range(0, len(tails)):
              if tails[0] == "0":
                  tails = tails[1:] + [str(0)]
              else:
                  break
      return (head, tails)
  
  
  def convert_float_base(x, base, precision=10):
  
      length_of_int = 0 if x == 0 else int(mpmath.log(x, base))
      # iexps = list(range(length_of_int, -1, -1))
  
      def convert_int(x, base, exponents):
          out = []
          for e in range(0, exponents + 1):
              d = x % base
              out.append(d)
              x = x / base
              if x == 0:
                  break
          out.reverse()
          return out
  
      def convert_float(x, base, exponents):
          out = []
          for e in range(0, exponents):
              d = int(x * base)
              out.append(d)
              x = (x * base) - d
              if x == 0:
                  break
          return out
  
      int_part = convert_int(int(x), base, length_of_int)
      if isinstance(x, (float, sympy.Float)):
          # fexps = list(range(-1, -int(precision + 1), -1))
          real_part = convert_float(x - int(x), base, precision + 1)
          return int_part + real_part
      elif isinstance(x, int):
          return int_part
      else:
          raise TypeError(x)
  
  
  class RealDigits(Builtin):
      """
      <dl>
      <dt>'RealDigits[$n$]'
          <dd>returns the decimal representation of the real number $n$ as list of digits, together with the number of digits that are to the left of the decimal point.
  
      <dt>'RealDigits[$n$, $b$]'
          <dd>returns a list of base_$b$ representation of the real number $n$.
  
      <dt>'RealDigits[$n$, $b$, $len$]'
          <dd>returns a list of $len$ digits.
  
      <dt>'RealDigits[$n$, $b$, $len$, $p$]'
          <dd>return $len$ digits starting with the coefficient of $b$^$p$
      </dl>
  
      Return the list of digits and exponent:
      >> RealDigits[123.55555]
       = {{1, 2, 3, 5, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0}, 3}
  
      Return an explicit recurring decimal form:
      >> RealDigits[19 / 7]
       = {{2, {7, 1, 4, 2, 8, 5}}, 1}
  
      The 10000th digit of  is an 8:
      >> RealDigits[Pi, 10, 1, -10000]
      = {{8}, -9999}
  
      20 digits starting with the coefficient of 10^-5:
      >> RealDigits[Pi, 10, 20, -5]
       = {{9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3, 8, 4, 6, 2, 6, 4, 3}, -4}
  
      RealDigits gives Indeterminate if more digits than the precision are requested:
      >> RealDigits[123.45, 10, 18]
       = {{1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, Indeterminate, Indeterminate}, 3}
  
      #> RealDigits[-1.25, -1]
       : Base -1 is not a real number greater than 1.
       = RealDigits[-1.25, -1]
  
      Return 25 digits of in base 10:
      >> RealDigits[Pi, 10, 25]
       = {{3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3, 8, 4, 6, 2, 6, 4, 3}, 1}
  
      #> RealDigits[-Pi]
       : The number of digits to return cannot be determined.
       = RealDigits[-Pi]
  
      #> RealDigits[I, 7]
       : The value I is not a real number.
      = RealDigits[I, 7]
  
      #> RealDigits[Pi]
       : The number of digits to return cannot be determined.
       = RealDigits[Pi]
  
      #> RealDigits[3 + 4 I]
       : The value 3 + 4 I is not a real number.
       = RealDigits[3 + 4 I]
  
  
      #> RealDigits[3.14, 10, 1.5]
       : Non-negative machine-sized integer expected at position 3 in RealDigits[3.14, 10, 1.5].
       = RealDigits[3.14, 10, 1.5]
  
      #> RealDigits[3.14, 10, 1, 1.5]
       : Machine-sized integer expected at position 4 in RealDigits[3.14, 10, 1, 1.5].
       = RealDigits[3.14, 10, 1, 1.5]
  
      """
  
      attributes = ("Listable",)
  
      messages = {
          "realx": "The value `1` is not a real number.",
          "ndig": "The number of digits to return cannot be determined.",
          "rbase": "Base `1` is not a real number greater than 1.",
          "intnm": "Non-negative machine-sized integer expected at position 3 in `1`.",
          "intm": "Machine-sized integer expected at position 4 in `1`.",
      }
  
      def apply_complex(self, n, var, evaluation):
          "%(name)s[n_Complex, var___]"
          return evaluation.message("RealDigits", "realx", n)
  
      def apply_rational_with_base(self, n, b, evaluation):
          "%(name)s[n_Rational, b_Integer]"
          # expr = Expression("RealDigits", n)
          py_n = abs(n.value)
          py_b = b.get_int_value()
          if check_finite_decimal(n.denominator().get_int_value()) and not py_b % 2:
              return self.apply_with_base(n, b, evaluation)
          else:
              exp = int(mpmath.ceil(mpmath.log(py_n, py_b)))
              (head, tails) = convert_repeating_decimal(
                  py_n.as_numer_denom()[0], py_n.as_numer_denom()[1], py_b
              )
  
              leaves = []
              for x in head:
                  if not x == "0":
                      leaves.append(from_python(int(x)))
              leaves.append(from_python(tails))
              list_str = Expression(SymbolList, *leaves)
          return Expression(SymbolList, list_str, exp)
  
      def apply_rational_without_base(self, n, evaluation):
          "%(name)s[n_Rational]"
  
          return self.apply_rational_with_base(n, from_python(10), evaluation)
  
      def apply(self, n, evaluation):
          "%(name)s[n_]"
  
          # Handling the testcases that throw the error message and return the ouput that doesn't include `base` argument
          if isinstance(n, Symbol) and n.name.startswith("System`"):
              return evaluation.message("RealDigits", "ndig", n)
  
          if n.is_numeric():
              return self.apply_with_base(n, from_python(10), evaluation)
  
      def apply_with_base(self, n, b, evaluation, nr_elements=None, pos=None):
          "%(name)s[n_?NumericQ, b_Integer]"
  
          expr = Expression("RealDigits", n)
          rational_no = (
              True if isinstance(n, Rational) else False
          )  # it is used for checking whether the input n is a rational or not
          py_b = b.get_int_value()
          if isinstance(n, (Expression, Symbol, Rational)):
              pos_len = abs(pos) + 1 if pos is not None and pos < 0 else 1
              if nr_elements is not None:
                  n = Expression(
                      "N", n, int(mpmath.log(py_b ** (nr_elements + pos_len), 10)) + 1
                  ).evaluate(evaluation)
              else:
                  if rational_no:
                      n = Expression(SymbolN, n).evaluate(evaluation)
                  else:
                      return evaluation.message("RealDigits", "ndig", expr)
          py_n = abs(n.value)
  
          if not py_b > 1:
              return evaluation.message("RealDigits", "rbase", py_b)
  
          if isinstance(py_n, complex):
              return evaluation.message("RealDigits", "realx", expr)
  
          if isinstance(n, Integer):
              display_len = (
                  int(mpmath.floor(mpmath.log(py_n, py_b)))
                  if py_n != 0 and py_n != 1
                  else 1
              )
          else:
              display_len = int(
                  Expression(
                      "N",
                      Expression(
                          "Round",
                          Expression(
                              "Divide",
                              Expression("Precision", py_n),
                              Expression("Log", 10, py_b),
                          ),
                      ),
                  )
                  .evaluate(evaluation)
                  .to_python()
              )
  
          exp = log_n_b(py_n, py_b)
  
          if py_n == 0 and nr_elements is not None:
              exp = 0
  
          digits = []
          if not py_b == 10:
              digits = convert_float_base(py_n, py_b, display_len - exp)
              # truncate all the leading 0's
              i = 0
              while digits and digits[i] == 0:
                  i += 1
              digits = digits[i:]
  
              if not isinstance(n, Integer):
                  if len(digits) > display_len:
                      digits = digits[: display_len - 1]
          else:
              # drop any leading zeroes
              for x in str(py_n):
                  if x != "." and (digits or x != "0"):
                      digits.append(x)
  
          if pos is not None:
              temp = exp
              exp = pos + 1
              move = temp - 1 - pos
              if move <= 0:
                  digits = [0] * abs(move) + digits
              else:
                  digits = digits[abs(move) :]
                  display_len = display_len - move
  
          leaves = []
          for x in digits:
              if x == "e" or x == "E":
                  break
              # Convert to Mathics' list format
              leaves.append(from_python(int(x)))
  
          if not rational_no:
              while len(leaves) < display_len:
                  leaves.append(from_python(0))
  
          if nr_elements is not None:
              # display_len == nr_elements
              if len(leaves) >= nr_elements:
                  # Truncate, preserving the digits on the right
                  leaves = leaves[:nr_elements]
              else:
                  if isinstance(n, Integer):
                      while len(leaves) < nr_elements:
                          leaves.append(from_python(0))
                  else:
                      # Adding Indeterminate if the length is greater than the precision
                      while len(leaves) < nr_elements:
                          leaves.append(from_python(Symbol("Indeterminate")))
          list_str = Expression(SymbolList, *leaves)
          return Expression(SymbolList, list_str, exp)
  
      def apply_with_base_and_length(self, n, b, length, evaluation, pos=None):
          "%(name)s[n_?NumericQ, b_Integer, length_]"
          leaves = []
          if pos is not None:
              leaves.append(from_python(pos))
          expr = Expression("RealDigits", n, b, length, *leaves)
          if not (isinstance(length, Integer) and length.get_int_value() >= 0):
              return evaluation.message("RealDigits", "intnm", expr)
  
          return self.apply_with_base(
              n, b, evaluation, nr_elements=length.get_int_value(), pos=pos
          )
  
      def apply_with_base_length_and_precision(self, n, b, length, p, evaluation):
          "%(name)s[n_?NumericQ, b_Integer, length_, p_]"
          if not isinstance(p, Integer):
              return evaluation.message(
                  "RealDigits", "intm", Expression("RealDigits", n, b, length, p)
              )
  
          return self.apply_with_base_and_length(
              n, b, length, evaluation, pos=p.get_int_value()
          )
  
  
  class _ZLibHash:  # make zlib hashes behave as if they were from hashlib
      def __init__(self, fn):
          self._bytes = b""
          self._fn = fn
  
      def update(self, bytes):
          self._bytes += bytes
  
      def hexdigest(self):
          return format(self._fn(self._bytes), "x")
  
  
  class Hash(Builtin):
      """
      <dl>
      <dt>'Hash[$expr$]'
        <dd>returns an integer hash for the given $expr$.
      <dt>'Hash[$expr$, $type$]'
        <dd>returns an integer hash of the specified $type$ for the given $expr$.</dd>
        <dd>The types supported are "MD5", "Adler32", "CRC32", "SHA", "SHA224", "SHA256", "SHA384", and "SHA512".</dd>
      <dt>'Hash[$expr$, $type$, $format$]'
        <dd>Returns the hash in the specified format.</dd>
      </dl>
  
      > Hash["The Adventures of Huckleberry Finn"]
      = 213425047836523694663619736686226550816
  
      > Hash["The Adventures of Huckleberry Finn", "SHA256"]
      = 95092649594590384288057183408609254918934351811669818342876362244564858646638
  
      > Hash[1/3]
      = 56073172797010645108327809727054836008
  
      > Hash[{a, b, {c, {d, e, f}}}]
      = 135682164776235407777080772547528225284
  
      > Hash[SomeHead[3.1415]]
      = 58042316473471877315442015469706095084
  
      >> Hash[{a, b, c}, "xyzstr"]
       = Hash[{a, b, c}, xyzstr, Integer]
      """
  
      rules = {
          "Hash[expr_]": 'Hash[expr, "MD5", "Integer"]',
          "Hash[expr_, type_String]": 'Hash[expr, type, "Integer"]',
      }
  
      attributes = ("Protected", "ReadProtected")
  
      # FIXME md2
      _supported_hashes = {
          "Adler32": lambda: _ZLibHash(zlib.adler32),
          "CRC32": lambda: _ZLibHash(zlib.crc32),
          "MD5": hashlib.md5,
          "SHA": hashlib.sha1,
          "SHA224": hashlib.sha224,
          "SHA256": hashlib.sha256,
          "SHA384": hashlib.sha384,
          "SHA512": hashlib.sha512,
      }
  
      @staticmethod
      def compute(user_hash, py_hashtype, py_format):
          hash_func = Hash._supported_hashes.get(py_hashtype)
          if hash_func is None:  # unknown hash function?
              return  # in order to return original Expression
          h = hash_func()
          user_hash(h.update)
          res = h.hexdigest()
          if py_format in ("HexString", "HexStringLittleEndian"):
              return from_python(res)
          res = int(res, 16)
          if py_format == "DecimalString":
              return from_python(str(res))
          elif py_format == "ByteArray":
              return from_python(bytearray(res))
          return from_python(res)
  
      def apply(self, expr, hashtype, outformat, evaluation):
          "Hash[expr_, hashtype_String, outformat_String]"
          return Hash.compute(
              expr.user_hash, hashtype.get_string_value(), outformat.get_string_value()
          )
  
  
  class TypeEscalation(Exception):
      def __init__(self, mode):
          self.mode = mode
  
  
  class Fold(object):
      # allows inherited classes to specify a single algorithm implementation that
      # can be called with machine precision, arbitrary precision or symbolically.
  
      ComputationFunctions = namedtuple("ComputationFunctions", ("sin", "cos"))
  
      FLOAT = 0
      MPMATH = 1
      SYMBOLIC = 2
  
      math = {
          FLOAT: ComputationFunctions(
              cos=math.cos,
              sin=math.sin,
          ),
          MPMATH: ComputationFunctions(
              cos=mpmath.cos,
              sin=mpmath.sin,
          ),
          SYMBOLIC: ComputationFunctions(
              cos=lambda x: Expression("Cos", x),
              sin=lambda x: Expression("Sin", x),
          ),
      }
  
      operands = {
          FLOAT: lambda x: None if x is None else x.round_to_float(),
          MPMATH: lambda x: None if x is None else x.to_mpmath(),
          SYMBOLIC: lambda x: x,
      }
  
      def _operands(self, state, steps):
          raise NotImplementedError
  
      def _fold(self, state, steps, math):
          raise NotImplementedError
  
      def _spans(self, operands):
          spans = {}
          k = 0
          j = 0
  
          for mode in (self.FLOAT, self.MPMATH):
              for i, operand in enumerate(operands[k:]):
                  if operand[0] > mode:
                      break
                  j = i + k + 1
  
              if k == 0 and j == 1:  # only init state? then ignore.
                  j = 0
  
              spans[mode] = slice(k, j)
              k = j
  
          spans[self.SYMBOLIC] = slice(k, len(operands))
  
          return spans
  
      def fold(self, x, ll):
          # computes fold(x, ll) with the internal _fold function. will start
          # its evaluation machine precision, and will escalate to arbitrary
          # precision if or symbolical evaluation only if necessary. folded
          # items already computed are carried over to new evaluation modes.
  
          yield x  # initial state
  
          init = None
          operands = list(self._operands(x, ll))
          spans = self._spans(operands)
  
          for mode in (self.FLOAT, self.MPMATH, self.SYMBOLIC):
              s_operands = [y[1:] for y in operands[spans[mode]]]
  
              if not s_operands:
                  continue
  
              if mode == self.MPMATH:
                  from mathics.core.numbers import min_prec
  
                  precision = min_prec(*[t for t in chain(*s_operands) if t is not None])
                  working_precision = mpmath.workprec
              else:
  
                  @contextmanager
                  def working_precision(_):
                      yield
  
                  precision = None
  
              if mode == self.FLOAT:
  
                  def out(z):
                      return Real(z)
  
              elif mode == self.MPMATH:
  
                  def out(z):
                      return Real(z, precision)
  
              else:
  
                  def out(z):
                      return z
  
              as_operand = self.operands.get(mode)
  
              def converted_operands():
                  for y in s_operands:
                      yield tuple(as_operand(t) for t in y)
  
              with working_precision(precision):
                  c_operands = converted_operands()
  
                  if init is not None:
                      c_init = tuple(
                          (None if t is None else as_operand(from_python(t)))
                          for t in init
                      )
                  else:
                      c_init = next(c_operands)
                      init = tuple((None if t is None else out(t)) for t in c_init)
  
                  generator = self._fold(c_init, c_operands, self.math.get(mode))
  
                  for y in generator:
                      y = tuple(out(t) for t in y)
                      yield y
                      init = y
  