1# cython: language_level=3 2# -*- coding: utf-8 -*- 3 4 5""" 6Numeric Evaluation 7 8Support for numeric evaluation with arbitrary precision is just a 9proof-of-concept. 10Precision is not "guarded" through the evaluation process. Only 11integer precision is supported. 12However, things like 'N[Pi, 100]' should work as expected. 13""" 14from mathics.version import __version__ # noqa used in loading to check consistency. 15import sympy 16import mpmath 17import numpy as np 18import math 19import hashlib 20import zlib 21from collections import namedtuple 22from contextlib import contextmanager 23from itertools import chain, product 24from functools import lru_cache 25 26 27from mathics.builtin.base import Builtin, Predefined 28from mathics.core.numbers import ( 29 dps, 30 convert_int_to_digit_list, 31 machine_precision, 32 machine_epsilon, 33 get_precision, 34 PrecisionValueError, 35) 36from mathics.core.expression import ( 37 Complex, 38 Expression, 39 Integer, 40 Integer0, 41 MachineReal, 42 Number, 43 Rational, 44 Real, 45 String, 46 Symbol, 47 SymbolFalse, 48 SymbolTrue, 49 SymbolList, 50 SymbolN, 51 from_python, 52) 53from mathics.core.convert import from_sympy 54 55 56@lru_cache(maxsize=1024) 57def log_n_b(py_n, py_b) -> int: 58 return int(mpmath.ceil(mpmath.log(py_n, py_b))) if py_n != 0 and py_n != 1 else 1 59 60 61class N(Builtin): 62 """ 63 <dl> 64 <dt>'N[$expr$, $prec$]' 65 <dd>evaluates $expr$ numerically with a precision of $prec$ digits. 66 </dl> 67 >> N[Pi, 50] 68 = 3.1415926535897932384626433832795028841971693993751 69 70 >> N[1/7] 71 = 0.142857 72 73 >> N[1/7, 5] 74 = 0.14286 75 76 You can manually assign numerical values to symbols. 77 When you do not specify a precision, 'MachinePrecision' is taken. 78 >> N[a] = 10.9 79 = 10.9 80 >> a 81 = a 82 83 'N' automatically threads over expressions, except when a symbol has 84 attributes 'NHoldAll', 'NHoldFirst', or 'NHoldRest'. 85 >> N[a + b] 86 = 10.9 + b 87 >> N[a, 20] 88 = a 89 >> N[a, 20] = 11; 90 >> N[a + b, 20] 91 = 11.000000000000000000 + b 92 >> N[f[a, b]] 93 = f[10.9, b] 94 >> SetAttributes[f, NHoldAll] 95 >> N[f[a, b]] 96 = f[a, b] 97 98 The precision can be a pattern: 99 >> N[c, p_?(#>10&)] := p 100 >> N[c, 3] 101 = c 102 >> N[c, 11] 103 = 11.000000000 104 105 You can also use 'UpSet' or 'TagSet' to specify values for 'N': 106 >> N[d] ^= 5; 107 However, the value will not be stored in 'UpValues', but 108 in 'NValues' (as for 'Set'): 109 >> UpValues[d] 110 = {} 111 >> NValues[d] 112 = {HoldPattern[N[d, MachinePrecision]] :> 5} 113 >> e /: N[e] = 6; 114 >> N[e] 115 = 6. 116 117 Values for 'N[$expr$]' must be associated with the head of $expr$: 118 >> f /: N[e[f]] = 7; 119 : Tag f not found or too deep for an assigned rule. 120 121 You can use 'Condition': 122 >> N[g[x_, y_], p_] := x + y * Pi /; x + y > 3 123 >> SetAttributes[g, NHoldRest] 124 >> N[g[1, 1]] 125 = g[1., 1] 126 >> N[g[2, 2]] // InputForm 127 = 8.283185307179586 128 129 The precision of the result is no higher than the precision of the input 130 >> N[Exp[0.1], 100] 131 = 1.10517 132 >> % // Precision 133 = MachinePrecision 134 >> N[Exp[1/10], 100] 135 = 1.105170918075647624811707826490246668224547194737518718792863289440967966747654302989143318970748654 136 >> % // Precision 137 = 100. 138 >> N[Exp[1.0`20], 100] 139 = 2.7182818284590452354 140 >> % // Precision 141 = 20. 142 143 #> p=N[Pi,100] 144 = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068 145 #> ToString[p] 146 = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068 147 #> 3.14159 * "a string" 148 = 3.14159 a string 149 150 #> N[Pi, Pi] 151 = 3.14 152 153 #> N[1/9, 30] 154 = 0.111111111111111111111111111111 155 #> Precision[%] 156 = 30. 157 158 #> N[1.5, 30] 159 = 1.5 160 #> Precision[%] 161 = MachinePrecision 162 #> N[1.5, 5] 163 = 1.5 164 #> Precision[%] 165 = MachinePrecision 166 167 #> {N[x], N[x, 30], N["abc"], N["abc", 30]} 168 = {x, x, abc, abc} 169 170 #> N[I, 30] 171 = 1.00000000000000000000000000000 I 172 173 #> N[1.01234567890123456789] 174 = 1.01235 175 #> N[1.012345678901234567890123, 20] 176 = 1.0123456789012345679 177 #> N[1.012345678901234567890123, 5] 178 = 1.0123 179 #> % // Precision 180 = 5. 181 #> N[1.012345678901234567890123, 50] 182 = 1.01234567890123456789012 183 #> % // Precision 184 = 24. 185 186 #> N[1.01234567890123456789`] 187 = 1.01235 188 #> N[1.01234567890123456789`, 20] 189 = 1.01235 190 #> % // Precision 191 = MachinePrecision 192 #> N[1.01234567890123456789`, 2] 193 = 1.01235 194 #> % // Precision 195 = MachinePrecision 196 """ 197 198 messages = { 199 "precbd": ("Requested precision `1` is not a " + "machine-sized real number."), 200 "preclg": ( 201 "Requested precision `1` is larger than $MaxPrecision. " 202 + "Using current $MaxPrecision of `2` instead. " 203 + "$MaxPrecision = Infinity specifies that any precision " 204 + "should be allowed." 205 ), 206 "precsm": ( 207 "Requested precision `1` is smaller than " 208 + "$MinPrecision. Using current $MinPrecision of " 209 + "`2` instead." 210 ), 211 } 212 213 rules = { 214 "N[expr_]": "N[expr, MachinePrecision]", 215 } 216 217 def apply_with_prec(self, expr, prec, evaluation): 218 "N[expr_, prec_]" 219 220 try: 221 d = get_precision(prec, evaluation) 222 except PrecisionValueError: 223 return 224 225 if expr.get_head_name() in ("System`List", "System`Rule"): 226 return Expression( 227 expr.head, 228 *[self.apply_with_prec(leaf, prec, evaluation) for leaf in expr.leaves], 229 ) 230 231 # Special case for the Root builtin 232 if expr.has_form("Root", 2): 233 return from_sympy(sympy.N(expr.to_sympy(), d)) 234 235 if isinstance(expr, Number): 236 return expr.round(d) 237 238 name = expr.get_lookup_name() 239 if name != "": 240 nexpr = Expression(SymbolN, expr, prec) 241 result = evaluation.definitions.get_value( 242 name, "System`NValues", nexpr, evaluation 243 ) 244 if result is not None: 245 if not result.sameQ(nexpr): 246 result = Expression(SymbolN, result, prec).evaluate(evaluation) 247 return result 248 249 if expr.is_atom(): 250 return expr 251 else: 252 attributes = expr.head.get_attributes(evaluation.definitions) 253 if "System`NHoldAll" in attributes: 254 eval_range = () 255 elif "System`NHoldFirst" in attributes: 256 eval_range = range(1, len(expr.leaves)) 257 elif "System`NHoldRest" in attributes: 258 if len(expr.leaves) > 0: 259 eval_range = (0,) 260 else: 261 eval_range = () 262 else: 263 eval_range = range(len(expr.leaves)) 264 head = Expression(SymbolN, expr.head, prec).evaluate(evaluation) 265 leaves = expr.get_mutable_leaves() 266 for index in eval_range: 267 leaves[index] = Expression(SymbolN, leaves[index], prec).evaluate( 268 evaluation 269 ) 270 return Expression(head, *leaves) 271 272 273def _scipy_interface(integrator, options_map, mandatory=None, adapt_func=None): 274 """ 275 This function provides a proxy for scipy.integrate 276 functions, adapting the parameters. 277 """ 278 279 def _scipy_proxy_func_filter(fun, a, b, **opts): 280 native_opts = {} 281 if mandatory: 282 native_opts.update(mandatory) 283 for opt, val in opts.items(): 284 native_opt = options_map.get(opt, None) 285 if native_opt: 286 if native_opt[1]: 287 val = native_opt[1](val) 288 native_opts[native_opt[0]] = val 289 return adapt_func(integrator(fun, a, b, **native_opts)) 290 291 def _scipy_proxy_func(fun, a, b, **opts): 292 native_opts = {} 293 if mandatory: 294 native_opts.update(mandatory) 295 for opt, val in opts.items(): 296 native_opt = options_map.get(opt, None) 297 if native_opt: 298 if native_opt[1]: 299 val = native_opt[1](val) 300 native_opts[native_opt[0]] = val 301 return integrator(fun, a, b, **native_opts) 302 303 return _scipy_proxy_func_filter if adapt_func else _scipy_proxy_func 304 305 306def _internal_adaptative_simpsons_rule(f, a, b, **opts): 307 """ 308 1D adaptative Simpson's rule integrator 309 Adapted from https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method 310 by @mmatera 311 312 TODO: handle weak divergences 313 """ 314 wsr = 1.0 / 6.0 315 316 tol = opts.get("tol") 317 if not tol: 318 tol = 1.0e-10 319 320 maxrec = opts.get("maxrec") 321 if not maxrec: 322 maxrec = 150 323 324 def _quad_simpsons_mem(f, a, fa, b, fb): 325 """Evaluates the Simpson's Rule, also returning m and f(m) to reuse""" 326 m = 0.5 * (a + b) 327 try: 328 fm = f(m) 329 except ZeroDivisionError: 330 fm = None 331 332 if fm is None or np.isinf(fm): 333 m = m + 1e-10 334 fm = f(m) 335 return (m, fm, wsr * abs(b - a) * (fa + 4.0 * fm + fb)) 336 337 def _quad_asr(f, a, fa, b, fb, eps, whole, m, fm, maxrec): 338 """ 339 Efficient recursive implementation of adaptive Simpson's rule. 340 Function values at the start, middle, end of the intervals 341 are retained. 342 """ 343 maxrec = maxrec - 1 344 try: 345 left = _quad_simpsons_mem(f, a, fa, m, fm) 346 lm, flm, left = left 347 right = _quad_simpsons_mem(f, m, fm, b, fb) 348 rm, frm, right = right 349 350 delta = left + right - whole 351 err = abs(delta) 352 if err <= 15 * eps or maxrec == 0: 353 return (left + right + delta / 15, err) 354 left = _quad_asr(f, a, fa, m, fm, 0.5 * eps, left, lm, flm, maxrec) 355 right = _quad_asr(f, m, fm, b, fb, 0.5 * eps, right, rm, frm, maxrec) 356 return (left[0] + right[0], left[1] + right[1]) 357 except Exception: 358 raise 359 360 def ensure_evaluation(f, x): 361 try: 362 val = f(x) 363 except ZeroDivisionError: 364 return None 365 if np.isinf(val): 366 return None 367 return val 368 369 invert_interval = False 370 if a > b: 371 b, a, invert_interval = a, b, True 372 373 fa, fb = ensure_evaluation(f, a), ensure_evaluation(f, b) 374 if fa is None: 375 x = 10.0 * machine_epsilon if a == 0 else a * (1.0 + 10.0 * machine_epsilon) 376 fa = ensure_evaluation(f, x) 377 if fa is None: 378 raise Exception(f"Function undefined around {a}. Cannot integrate") 379 if fb is None: 380 x = -10.0 * machine_epsilon if b == 0 else b * (1.0 - 10.0 * machine_epsilon) 381 fb = ensure_evaluation(f, x) 382 if fb is None: 383 raise Exception(f"Function undefined around {b}. Cannot integrate") 384 385 m, fm, whole = _quad_simpsons_mem(f, a, fa, b, fb) 386 if invert_interval: 387 return -_quad_asr(f, a, fa, b, fb, tol, whole, m, fm, maxrec) 388 else: 389 return _quad_asr(f, a, fa, b, fb, tol, whole, m, fm, maxrec) 390 391 392def _fubini(func, ranges, **opts): 393 if not ranges: 394 return 0.0 395 a, b = ranges[0] 396 integrator = opts["integrator"] 397 tol = opts.get("tol") 398 if tol is None: 399 opts["tol"] = 1.0e-10 400 tol = 1.0e-10 401 402 if len(ranges) > 1: 403 404 def subintegral(*u): 405 def ff(*z): 406 return func(*(u + z)) 407 408 val = _fubini(ff, ranges[1:], **opts)[0] 409 return val 410 411 opts["tol"] = 4.0 * tol 412 val = integrator(subintegral, a, b, **opts) 413 return val 414 else: 415 val = integrator(func, a, b, **opts) 416 return val 417 418 419class NIntegrate(Builtin): 420 """ 421 <dl> 422 <dt>'NIntegrate[$expr$, $interval$]' 423 <dd>returns a numeric approximation to the definite integral of $expr$ with limits $interval$ and with a precision of $prec$ digits. 424 425 <dt>'NIntegrate[$expr$, $interval1$, $interval2$, ...]' 426 <dd>returns a numeric approximation to the multiple integral of $expr$ with limits $interval1$, $interval2$ and with a precision of $prec$ digits. 427 </dl> 428 429 >> NIntegrate[Exp[-x],{x,0,Infinity},Tolerance->1*^-6] 430 = 1. 431 >> NIntegrate[Exp[x],{x,-Infinity, 0},Tolerance->1*^-6] 432 = 1. 433 >> NIntegrate[Exp[-x^2/2.],{x,-Infinity, Infinity},Tolerance->1*^-6] 434 = 2.50663 435 436 >> Table[1./NIntegrate[x^k,{x,0,1},Tolerance->1*^-6], {k,0,6}] 437 : The specified method failed to return a number. Falling back into the internal evaluator. 438 = {1., 2., 3., 4., 5., 6., 7.} 439 440 >> NIntegrate[1 / z, {z, -1 - I, 1 - I, 1 + I, -1 + I, -1 - I}, Tolerance->1.*^-4] 441 : Integration over a complex domain is not implemented yet 442 = NIntegrate[1 / z, {z, -1 - I, 1 - I, 1 + I, -1 + I, -1 - I}, Tolerance -> 0.0001] 443 ## = 6.2832 I 444 445 Integrate singularities with weak divergences: 446 >> Table[ NIntegrate[x^(1./k-1.), {x,0,1.}, Tolerance->1*^-6], {k,1,7.} ] 447 = {1., 2., 3., 4., 5., 6., 7.} 448 449 Mutiple Integrals : 450 >> NIntegrate[x * y,{x, 0, 1}, {y, 0, 1}] 451 = 0.25 452 453 """ 454 455 messages = { 456 "bdmtd": "The Method option should be a built-in method name.", 457 "inumr": ( 458 "The integrand `1` has evaluated to non-numerical " 459 + "values for all sampling points in the region " 460 + "with boundaries `2`" 461 ), 462 "nlim": "`1` = `2` is not a valid limit of integration.", 463 "ilim": "Invalid integration variable or limit(s) in `1`.", 464 "mtdfail": ( 465 "The specified method failed to return a " 466 + "number. Falling back into the internal " 467 + "evaluator." 468 ), 469 "cmpint": ("Integration over a complex domain is not " + "implemented yet"), 470 } 471 472 options = { 473 "Method": '"Automatic"', 474 "Tolerance": "1*^-10", 475 "Accuracy": "1*^-10", 476 "MaxRecursion": "10", 477 } 478 479 methods = { 480 "Automatic": (None, False), 481 } 482 483 def __init__(self, *args, **kwargs): 484 super().__init__(*args, **kwargs) 485 self.methods["Internal"] = (_internal_adaptative_simpsons_rule, False) 486 try: 487 from scipy.integrate import romberg, quad, nquad 488 489 self.methods["NQuadrature"] = ( 490 _scipy_interface( 491 nquad, {}, {"full_output": 1}, lambda res: (res[0], res[1]) 492 ), 493 True, 494 ) 495 self.methods["Quadrature"] = ( 496 _scipy_interface( 497 quad, 498 { 499 "tol": ("epsabs", None), 500 "maxrec": ("limit", lambda maxrec: int(2 ** maxrec)), 501 }, 502 {"full_output": 1}, 503 lambda res: (res[0], res[1]), 504 ), 505 False, 506 ) 507 self.methods["Romberg"] = ( 508 _scipy_interface( 509 romberg, 510 {"tol": ("tol", None), "maxrec": ("divmax", None)}, 511 None, 512 lambda x: (x, np.nan), 513 ), 514 False, 515 ) 516 self.methods["Automatic"] = self.methods["Quadrature"] 517 except Exception: 518 self.methods["Automatic"] = self.methods["Internal"] 519 self.methods["Simpson"] = self.methods["Internal"] 520 521 self.messages["bdmtd"] = ( 522 "The Method option should be a " 523 + "built-in method name in {`" 524 + "`, `".join(list(self.methods)) 525 + "`}. Using `Automatic`" 526 ) 527 528 @staticmethod 529 def decompose_domain(interval, evaluation): 530 if interval.has_form("System`Sequence", 1, None): 531 intervals = [] 532 for leaf in interval.leaves: 533 inner_interval = NIntegrate.decompose_domain(leaf, evaluation) 534 if inner_interval: 535 intervals.append(inner_interval) 536 else: 537 evaluation.message("ilim", leaf) 538 return None 539 return intervals 540 541 if interval.has_form("System`List", 3, None): 542 intervals = [] 543 intvar = interval.leaves[0] 544 if not isinstance(intvar, Symbol): 545 evaluation.message("ilim", interval) 546 return None 547 boundaries = [a for a in interval.leaves[1:]] 548 if any([b.get_head_name() == "System`Complex" for b in boundaries]): 549 intvar = Expression( 550 "List", intvar, Expression("Blank", Symbol("Complex")) 551 ) 552 for i in range(len(boundaries) - 1): 553 intervals.append((boundaries[i], boundaries[i + 1])) 554 if len(intervals) > 0: 555 return (intvar, intervals) 556 557 evaluation.message("ilim", interval) 558 return None 559 560 def apply_with_func_domain(self, func, domain, evaluation, options): 561 "%(name)s[func_, domain__, OptionsPattern[%(name)s]]" 562 method = options["System`Method"].evaluate(evaluation) 563 method_options = {} 564 if method.has_form("System`List", 2): 565 method = method.leaves[0] 566 method_options.update(method.leaves[1].get_option_values()) 567 if isinstance(method, String): 568 method = method.value 569 elif isinstance(method, Symbol): 570 method = method.get_name() 571 else: 572 evaluation.message("NIntegrate", "bdmtd", method) 573 return 574 tolerance = options["System`Tolerance"].evaluate(evaluation) 575 tolerance = float(tolerance.value) 576 accuracy = options["System`Accuracy"].evaluate(evaluation) 577 accuracy = accuracy.value 578 maxrecursion = options["System`MaxRecursion"].evaluate(evaluation) 579 maxrecursion = maxrecursion.value 580 nintegrate_method = self.methods.get(method, None) 581 if nintegrate_method is None: 582 evaluation.message("NIntegrate", "bdmtd", method) 583 nintegrate_method = self.methods.get("Automatic") 584 if type(nintegrate_method) is tuple: 585 nintegrate_method, is_multidimensional = nintegrate_method 586 else: 587 is_multidimensional = False 588 589 domain = self.decompose_domain(domain, evaluation) 590 if not domain: 591 return 592 if not isinstance(domain, list): 593 domain = [domain] 594 595 coords = [axis[0] for axis in domain] 596 # If any of the points in the integration domain is complex, 597 # stop the evaluation... 598 if any([c.get_head_name() == "System`List" for c in coords]): 599 evaluation.message("NIntegrate", "cmpint") 600 return 601 602 intvars = Expression(SymbolList, *coords) 603 integrand = Expression("Compile", intvars, func).evaluate(evaluation) 604 605 if len(integrand.leaves) >= 3: 606 integrand = integrand.leaves[2].cfunc 607 else: 608 evaluation.message("inumer", func, domain) 609 return 610 results = [] 611 for subdomain in product(*[axis[1] for axis in domain]): 612 # On each subdomain, check if the region is bounded. 613 # If not, implement a coordinate map 614 func2 = integrand 615 subdomain2 = [] 616 coordtransform = [] 617 nulldomain = False 618 for i, r in enumerate(subdomain): 619 a = r[0].evaluate(evaluation) 620 b = r[1].evaluate(evaluation) 621 if a == b: 622 nulldomain = True 623 break 624 elif a.get_head_name() == "System`DirectedInfinity": 625 if b.get_head_name() == "System`DirectedInfinity": 626 a = a.to_python() 627 b = b.to_python() 628 le = 1 - machine_epsilon 629 if a == b: 630 nulldomain = True 631 break 632 elif a < b: 633 subdomain2.append([-le, le]) 634 else: 635 subdomain2.append([le, -le]) 636 coordtransform.append( 637 (np.arctanh, lambda u: 1.0 / (1.0 - u ** 2)) 638 ) 639 else: 640 if not b.is_numeric(): 641 evaluation.message("nlim", coords[i], b) 642 return 643 z = a.leaves[0].value 644 b = b.value 645 subdomain2.append([machine_epsilon, 1.0]) 646 coordtransform.append( 647 (lambda u: b - z + z / u, lambda u: -z * u ** (-2.0)) 648 ) 649 elif b.get_head_name() == "System`DirectedInfinity": 650 if not a.is_numeric(): 651 evaluation.message("nlim", coords[i], a) 652 return 653 a = a.value 654 z = b.leaves[0].value 655 subdomain2.append([machine_epsilon, 1.0]) 656 coordtransform.append( 657 (lambda u: a - z + z / u, lambda u: z * u ** (-2.0)) 658 ) 659 elif a.is_numeric() and b.is_numeric(): 660 a = Expression(SymbolN, a).evaluate(evaluation).value 661 b = Expression(SymbolN, b).evaluate(evaluation).value 662 subdomain2.append([a, b]) 663 coordtransform.append(None) 664 else: 665 for x in (a, b): 666 if not x.is_numeric(): 667 evaluation.message("nlim", coords[i], x) 668 return 669 670 if nulldomain: 671 continue 672 673 if any(coordtransform): 674 func2 = lambda *u: ( 675 integrand( 676 *[ 677 x[0](u[i]) if x else u[i] 678 for i, x in enumerate(coordtransform) 679 ] 680 ) 681 * np.prod( 682 [jac[1](u[i]) for i, jac in enumerate(coordtransform) if jac] 683 ) 684 ) 685 opts = { 686 "acur": accuracy, 687 "tol": tolerance, 688 "maxrec": maxrecursion, 689 } 690 opts.update(method_options) 691 try: 692 if len(subdomain2) > 1: 693 if is_multidimensional: 694 nintegrate_method(func2, subdomain2, **opts) 695 else: 696 val = _fubini( 697 func2, subdomain2, integrator=nintegrate_method, **opts 698 ) 699 else: 700 val = nintegrate_method(func2, *(subdomain2[0]), **opts) 701 except Exception: 702 val = None 703 704 if val is None: 705 evaluation.message("NIntegrate", "mtdfail") 706 if len(subdomain2) > 1: 707 val = _fubini( 708 func2, 709 subdomain2, 710 integrator=_internal_adaptative_simpsons_rule, 711 **opts, 712 ) 713 else: 714 val = _internal_adaptative_simpsons_rule( 715 func2, *(subdomain2[0]), **opts 716 ) 717 results.append(val) 718 719 result = sum([r[0] for r in results]) 720 # error = sum([r[1] for r in results]) -> use it when accuracy 721 # be implemented... 722 return from_python(result) 723 724 725class MachinePrecision(Predefined): 726 """ 727 <dl> 728 <dt>'MachinePrecision' 729 <dd>represents the precision of machine precision numbers. 730 </dl> 731 732 >> N[MachinePrecision] 733 = 15.9546 734 >> N[MachinePrecision, 30] 735 = 15.9545897701910033463281614204 736 737 #> N[E, MachinePrecision] 738 = 2.71828 739 740 #> Round[MachinePrecision] 741 = 16 742 """ 743 744 rules = { 745 "N[MachinePrecision, prec_]": ("N[Log[10, 2] * %i, prec]" % machine_precision), 746 } 747 748 749class MachineEpsilon_(Predefined): 750 """ 751 <dl> 752 <dt>'$MachineEpsilon' 753 <dd>is the distance between '1.0' and the next 754 nearest representable machine-precision number. 755 </dl> 756 757 >> $MachineEpsilon 758 = 2.22045*^-16 759 760 >> x = 1.0 + {0.4, 0.5, 0.6} $MachineEpsilon; 761 >> x - 1 762 = {0., 0., 2.22045*^-16} 763 """ 764 765 name = "$MachineEpsilon" 766 767 def evaluate(self, evaluation): 768 return MachineReal(machine_epsilon) 769 770 771class MachinePrecision_(Predefined): 772 """ 773 <dl> 774 <dt>'$MachinePrecision' 775 <dd>is the number of decimal digits of precision for 776 machine-precision numbers. 777 </dl> 778 779 >> $MachinePrecision 780 = 15.9546 781 """ 782 783 name = "$MachinePrecision" 784 785 rules = { 786 "$MachinePrecision": "N[MachinePrecision]", 787 } 788 789 790class Precision(Builtin): 791 """ 792 <dl> 793 <dt>'Precision[$expr$]' 794 <dd>examines the number of significant digits of $expr$. 795 </dl> 796 This is rather a proof-of-concept than a full implementation. 797 Precision of compound expression is not supported yet. 798 >> Precision[1] 799 = Infinity 800 >> Precision[1/2] 801 = Infinity 802 >> Precision[0.5] 803 = MachinePrecision 804 805 #> Precision[0.0] 806 = MachinePrecision 807 #> Precision[0.000000000000000000000000000000000000] 808 = 0. 809 #> Precision[-0.0] 810 = MachinePrecision 811 #> Precision[-0.000000000000000000000000000000000000] 812 = 0. 813 814 #> 1.0000000000000000 // Precision 815 = MachinePrecision 816 #> 1.00000000000000000 // Precision 817 = 17. 818 819 #> 0.4 + 2.4 I // Precision 820 = MachinePrecision 821 #> Precision[2 + 3 I] 822 = Infinity 823 824 #> Precision["abc"] 825 = Infinity 826 """ 827 828 rules = { 829 "Precision[z_?MachineNumberQ]": "MachinePrecision", 830 } 831 832 def apply(self, z, evaluation): 833 "Precision[z_]" 834 835 if not z.is_inexact(): 836 return Symbol("Infinity") 837 elif z.to_sympy().is_zero: 838 return Real(0) 839 else: 840 return Real(dps(z.get_precision())) 841 842 843class MinPrecision(Builtin): 844 """ 845 <dl> 846 <dt>'$MinPrecision' 847 <dd>represents the minimum number of digits of precision 848 permitted in abitrary-precision numbers. 849 </dl> 850 851 >> $MinPrecision 852 = 0 853 854 >> $MinPrecision = 10; 855 856 >> N[Pi, 9] 857 : Requested precision 9 is smaller than $MinPrecision. Using current $MinPrecision of 10. instead. 858 = 3.141592654 859 860 #> N[Pi, 10] 861 = 3.141592654 862 863 #> $MinPrecision = x 864 : Cannot set $MinPrecision to x; value must be a non-negative number. 865 = x 866 #> $MinPrecision = -Infinity 867 : Cannot set $MinPrecision to -Infinity; value must be a non-negative number. 868 = -Infinity 869 #> $MinPrecision = -1 870 : Cannot set $MinPrecision to -1; value must be a non-negative number. 871 = -1 872 #> $MinPrecision = 0; 873 874 #> $MaxPrecision = 10; 875 #> $MinPrecision = 15 876 : Cannot set $MinPrecision such that $MaxPrecision < $MinPrecision. 877 = 15 878 #> $MinPrecision 879 = 0 880 #> $MaxPrecision = Infinity; 881 """ 882 883 name = "$MinPrecision" 884 rules = { 885 "$MinPrecision": "0", 886 } 887 888 messages = { 889 "precset": "Cannot set `1` to `2`; value must be a non-negative number.", 890 "preccon": "Cannot set `1` such that $MaxPrecision < $MinPrecision.", 891 } 892 893 894class MaxPrecision(Predefined): 895 """ 896 <dl> 897 <dt>'$MaxPrecision' 898 <dd>represents the maximum number of digits of precision 899 permitted in abitrary-precision numbers. 900 </dl> 901 902 >> $MaxPrecision 903 = Infinity 904 905 >> $MaxPrecision = 10; 906 907 >> N[Pi, 11] 908 : Requested precision 11 is larger than $MaxPrecision. Using current $MaxPrecision of 10. instead. $MaxPrecision = Infinity specifies that any precision should be allowed. 909 = 3.141592654 910 911 #> N[Pi, 10] 912 = 3.141592654 913 914 #> $MaxPrecision = x 915 : Cannot set $MaxPrecision to x; value must be a positive number or Infinity. 916 = x 917 #> $MaxPrecision = -Infinity 918 : Cannot set $MaxPrecision to -Infinity; value must be a positive number or Infinity. 919 = -Infinity 920 #> $MaxPrecision = 0 921 : Cannot set $MaxPrecision to 0; value must be a positive number or Infinity. 922 = 0 923 #> $MaxPrecision = Infinity; 924 925 #> $MinPrecision = 15; 926 #> $MaxPrecision = 10 927 : Cannot set $MaxPrecision such that $MaxPrecision < $MinPrecision. 928 = 10 929 #> $MaxPrecision 930 = Infinity 931 #> $MinPrecision = 0; 932 """ 933 934 name = "$MaxPrecision" 935 936 rules = { 937 "$MaxPrecision": "Infinity", 938 } 939 940 messages = { 941 "precset": "Cannot set `1` to `2`; value must be a positive number or Infinity.", 942 "preccon": "Cannot set `1` such that $MaxPrecision < $MinPrecision.", 943 } 944 945 946class Round(Builtin): 947 """ 948 <dl> 949 <dt>'Round[$expr$]' 950 <dd>rounds $expr$ to the nearest integer. 951 <dt>'Round[$expr$, $k$]' 952 <dd>rounds $expr$ to the closest multiple of $k$. 953 </dl> 954 955 >> Round[10.6] 956 = 11 957 >> Round[0.06, 0.1] 958 = 0.1 959 >> Round[0.04, 0.1] 960 = 0. 961 962 Constants can be rounded too 963 >> Round[Pi, .5] 964 = 3. 965 >> Round[Pi^2] 966 = 10 967 968 Round to exact value 969 >> Round[2.6, 1/3] 970 = 8 / 3 971 >> Round[10, Pi] 972 = 3 Pi 973 974 Round complex numbers 975 >> Round[6/(2 + 3 I)] 976 = 1 - I 977 >> Round[1 + 2 I, 2 I] 978 = 2 I 979 980 Round Negative numbers too 981 >> Round[-1.4] 982 = -1 983 984 Expressions other than numbers remain unevaluated: 985 >> Round[x] 986 = Round[x] 987 >> Round[1.5, k] 988 = Round[1.5, k] 989 """ 990 991 attributes = ("Listable", "NumericFunction") 992 993 rules = { 994 "Round[expr_?NumericQ]": "Round[Re[expr], 1] + I * Round[Im[expr], 1]", 995 "Round[expr_Complex, k_?RealNumberQ]": ( 996 "Round[Re[expr], k] + I * Round[Im[expr], k]" 997 ), 998 } 999 1000 def apply(self, expr, k, evaluation): 1001 "Round[expr_?NumericQ, k_?NumericQ]" 1002 1003 n = Expression("Divide", expr, k).round_to_float( 1004 evaluation, permit_complex=True 1005 ) 1006 if n is None: 1007 return 1008 elif isinstance(n, complex): 1009 n = round(n.real) 1010 else: 1011 n = round(n) 1012 n = int(n) 1013 return Expression("Times", Integer(n), k) 1014 1015 1016class Rationalize(Builtin): 1017 """ 1018 <dl> 1019 <dt>'Rationalize[$x$]' 1020 <dd>converts a real number $x$ to a nearby rational number. 1021 <dt>'Rationalize[$x$, $dx$]' 1022 <dd>finds the rational number within $dx$ of $x$ with the smallest denominator. 1023 </dl> 1024 1025 >> Rationalize[2.2] 1026 = 11 / 5 1027 1028 Not all numbers can be well approximated. 1029 >> Rationalize[N[Pi]] 1030 = 3.14159 1031 1032 Find the exact rational representation of 'N[Pi]' 1033 >> Rationalize[N[Pi], 0] 1034 = 245850922 / 78256779 1035 1036 #> Rationalize[1.6 + 0.8 I] 1037 = 8 / 5 + 4 I / 5 1038 1039 #> Rationalize[N[Pi] + 0.8 I, 1*^-6] 1040 = 355 / 113 + 4 I / 5 1041 1042 #> Rationalize[N[Pi] + 0.8 I, x] 1043 : Tolerance specification x must be a non-negative number. 1044 = Rationalize[3.14159 + 0.8 I, x] 1045 1046 #> Rationalize[N[Pi] + 0.8 I, -1] 1047 : Tolerance specification -1 must be a non-negative number. 1048 = Rationalize[3.14159 + 0.8 I, -1] 1049 1050 #> Rationalize[N[Pi] + 0.8 I, 0] 1051 = 245850922 / 78256779 + 4 I / 5 1052 1053 #> Rationalize[17 / 7] 1054 = 17 / 7 1055 1056 #> Rationalize[x] 1057 = x 1058 1059 #> Table[Rationalize[E, 0.1^n], {n, 1, 10}] 1060 = {8 / 3, 19 / 7, 87 / 32, 193 / 71, 1071 / 394, 2721 / 1001, 15062 / 5541, 23225 / 8544, 49171 / 18089, 419314 / 154257} 1061 1062 #> Rationalize[x, y] 1063 : Tolerance specification y must be a non-negative number. 1064 = Rationalize[x, y] 1065 """ 1066 1067 messages = { 1068 "tolnn": "Tolerance specification `1` must be a non-negative number.", 1069 } 1070 1071 rules = { 1072 "Rationalize[z_Complex]": "Rationalize[Re[z]] + I Rationalize[Im[z]]", 1073 "Rationalize[z_Complex, dx_?Internal`RealValuedNumberQ]/;dx >= 0": "Rationalize[Re[z], dx] + I Rationalize[Im[z], dx]", 1074 } 1075 1076 def apply(self, x, evaluation): 1077 "Rationalize[x_]" 1078 1079 py_x = x.to_sympy() 1080 if py_x is None or (not py_x.is_number) or (not py_x.is_real): 1081 return x 1082 return from_sympy(self.find_approximant(py_x)) 1083 1084 @staticmethod 1085 def find_approximant(x): 1086 c = 1e-4 1087 it = sympy.ntheory.continued_fraction_convergents( 1088 sympy.ntheory.continued_fraction_iterator(x) 1089 ) 1090 for i in it: 1091 p, q = i.as_numer_denom() 1092 tol = c / q ** 2 1093 if abs(i - x) <= tol: 1094 return i 1095 if tol < machine_epsilon: 1096 break 1097 return x 1098 1099 @staticmethod 1100 def find_exact(x): 1101 p, q = x.as_numer_denom() 1102 it = sympy.ntheory.continued_fraction_convergents( 1103 sympy.ntheory.continued_fraction_iterator(x) 1104 ) 1105 for i in it: 1106 p, q = i.as_numer_denom() 1107 if abs(x - i) < machine_epsilon: 1108 return i 1109 1110 def apply_dx(self, x, dx, evaluation): 1111 "Rationalize[x_, dx_]" 1112 py_x = x.to_sympy() 1113 if py_x is None: 1114 return x 1115 py_dx = dx.to_sympy() 1116 if ( 1117 py_dx is None 1118 or (not py_dx.is_number) 1119 or (not py_dx.is_real) 1120 or py_dx.is_negative 1121 ): 1122 return evaluation.message("Rationalize", "tolnn", dx) 1123 elif py_dx == 0: 1124 return from_sympy(self.find_exact(py_x)) 1125 a = self.approx_interval_continued_fraction(py_x - py_dx, py_x + py_dx) 1126 sym_x = sympy.ntheory.continued_fraction_reduce(a) 1127 return Rational(sym_x) 1128 1129 @staticmethod 1130 def approx_interval_continued_fraction(xmin, xmax): 1131 result = [] 1132 a_gen = sympy.ntheory.continued_fraction_iterator(xmin) 1133 b_gen = sympy.ntheory.continued_fraction_iterator(xmax) 1134 while True: 1135 a, b = next(a_gen), next(b_gen) 1136 if a == b: 1137 result.append(a) 1138 else: 1139 result.append(min(a, b) + 1) 1140 break 1141 return result 1142 1143 1144def chop(expr, delta=10.0 ** (-10.0)): 1145 if isinstance(expr, Real): 1146 if -delta < expr.get_float_value() < delta: 1147 return Integer0 1148 elif isinstance(expr, Complex) and expr.is_inexact(): 1149 real, imag = expr.real, expr.imag 1150 if -delta < real.get_float_value() < delta: 1151 real = Integer0 1152 if -delta < imag.get_float_value() < delta: 1153 imag = Integer0 1154 return Complex(real, imag) 1155 elif isinstance(expr, Expression): 1156 return Expression(chop(expr.head), *[chop(leaf) for leaf in expr.leaves]) 1157 return expr 1158 1159 1160class Chop(Builtin): 1161 """ 1162 <dl> 1163 <dt>'Chop[$expr$]' 1164 <dd>replaces floating point numbers close to 0 by 0. 1165 <dt>'Chop[$expr$, $delta$]' 1166 <dd>uses a tolerance of $delta$. The default tolerance is '10^-10'. 1167 </dl> 1168 1169 >> Chop[10.0 ^ -16] 1170 = 0 1171 >> Chop[10.0 ^ -9] 1172 = 1.*^-9 1173 >> Chop[10 ^ -11 I] 1174 = I / 100000000000 1175 >> Chop[0. + 10 ^ -11 I] 1176 = 0 1177 """ 1178 1179 messages = { 1180 "tolnn": "Tolerance specification a must be a non-negative number.", 1181 } 1182 1183 rules = { 1184 "Chop[expr_]": "Chop[expr, 10^-10]", 1185 } 1186 1187 def apply(self, expr, delta, evaluation): 1188 "Chop[expr_, delta_:(10^-10)]" 1189 1190 delta = delta.round_to_float(evaluation) 1191 if delta is None or delta < 0: 1192 return evaluation.message("Chop", "tolnn") 1193 1194 return chop(expr, delta=delta) 1195 1196 1197class NumericQ(Builtin): 1198 """ 1199 <dl> 1200 <dt>'NumericQ[$expr$]' 1201 <dd>tests whether $expr$ represents a numeric quantity. 1202 </dl> 1203 1204 >> NumericQ[2] 1205 = True 1206 >> NumericQ[Sqrt[Pi]] 1207 = True 1208 >> NumberQ[Sqrt[Pi]] 1209 = False 1210 """ 1211 1212 def apply(self, expr, evaluation): 1213 "NumericQ[expr_]" 1214 1215 def test(expr): 1216 if isinstance(expr, Expression): 1217 attr = evaluation.definitions.get_attributes(expr.head.get_name()) 1218 return "System`NumericFunction" in attr and all( 1219 test(leaf) for leaf in expr.leaves 1220 ) 1221 else: 1222 return expr.is_numeric() 1223 1224 return SymbolTrue if test(expr) else SymbolFalse 1225 1226 1227class RealValuedNumericQ(Builtin): 1228 """ 1229 #> Internal`RealValuedNumericQ /@ {1, N[Pi], 1/2, Sin[1.], Pi, 3/4, aa, I} 1230 = {True, True, True, True, True, True, False, False} 1231 """ 1232 1233 context = "Internal`" 1234 1235 rules = { 1236 "Internal`RealValuedNumericQ[x_]": "Head[N[x]] === Real", 1237 } 1238 1239 1240class RealValuedNumberQ(Builtin): 1241 """ 1242 #> Internal`RealValuedNumberQ /@ {1, N[Pi], 1/2, Sin[1.], Pi, 3/4, aa, I} 1243 = {True, True, True, True, False, True, False, False} 1244 """ 1245 1246 context = "Internal`" 1247 1248 rules = { 1249 "Internal`RealValuedNumberQ[x_Real]": "True", 1250 "Internal`RealValuedNumberQ[x_Integer]": "True", 1251 "Internal`RealValuedNumberQ[x_Rational]": "True", 1252 "Internal`RealValuedNumberQ[x_]": "False", 1253 } 1254 1255 1256class IntegerDigits(Builtin): 1257 """ 1258 <dl> 1259 <dt>'IntegerDigits[$n$]' 1260 <dd>returns a list of the base-10 digits in the integer $n$. 1261 <dt>'IntegerDigits[$n$, $base$]' 1262 <dd>returns a list of the base-$base$ digits in $n$. 1263 <dt>'IntegerDigits[$n$, $base$, $length$]' 1264 <dd>returns a list of length $length$, truncating or padding 1265 with zeroes on the left as necessary. 1266 </dl> 1267 1268 >> IntegerDigits[76543] 1269 = {7, 6, 5, 4, 3} 1270 1271 The sign of $n$ is discarded: 1272 >> IntegerDigits[-76543] 1273 = {7, 6, 5, 4, 3} 1274 1275 >> IntegerDigits[15, 16] 1276 = {15} 1277 >> IntegerDigits[1234, 16] 1278 = {4, 13, 2} 1279 >> IntegerDigits[1234, 10, 5] 1280 = {0, 1, 2, 3, 4} 1281 1282 #> IntegerDigits[1000, 10] 1283 = {1, 0, 0, 0} 1284 1285 #> IntegerDigits[0] 1286 = {0} 1287 """ 1288 1289 attributes = ("Listable",) 1290 1291 messages = { 1292 "int": "Integer expected at position 1 in `1`", 1293 "ibase": "Base `1` is not an integer greater than 1.", 1294 } 1295 1296 rules = { 1297 "IntegerDigits[n_]": "IntegerDigits[n, 10]", 1298 } 1299 1300 def apply_len(self, n, base, length, evaluation): 1301 "IntegerDigits[n_, base_, length_]" 1302 1303 if not (isinstance(length, Integer) and length.get_int_value() >= 0): 1304 return evaluation.message("IntegerDigits", "intnn") 1305 1306 return self.apply(n, base, evaluation, nr_elements=length.get_int_value()) 1307 1308 def apply(self, n, base, evaluation, nr_elements=None): 1309 "IntegerDigits[n_, base_]" 1310 1311 if not (isinstance(n, Integer)): 1312 return evaluation.message( 1313 "IntegerDigits", "int", Expression("IntegerDigits", n, base) 1314 ) 1315 1316 if not (isinstance(base, Integer) and base.get_int_value() > 1): 1317 return evaluation.message("IntegerDigits", "ibase", base) 1318 1319 if nr_elements == 0: 1320 # trivial case: we don't want any digits 1321 return Expression(SymbolList) 1322 1323 digits = convert_int_to_digit_list(n.get_int_value(), base.get_int_value()) 1324 1325 if nr_elements is not None: 1326 if len(digits) >= nr_elements: 1327 # Truncate, preserving the digits on the right 1328 digits = digits[-nr_elements:] 1329 else: 1330 # Pad with zeroes 1331 digits = [0] * (nr_elements - len(digits)) + digits 1332 1333 return Expression(SymbolList, *digits) 1334 1335 1336def check_finite_decimal(denominator): 1337 # The rational number is finite decimal if the denominator has form 2^a * 5^b 1338 while denominator % 5 == 0: 1339 denominator = denominator / 5 1340 1341 while denominator % 2 == 0: 1342 denominator = denominator / 2 1343 1344 return True if denominator == 1 else False 1345 1346 1347def convert_repeating_decimal(numerator, denominator, base): 1348 head = [x for x in str(numerator // denominator)] 1349 tails = [] 1350 subresults = [numerator % denominator] 1351 numerator %= denominator 1352 1353 while numerator != 0: # only rational input can go to this case 1354 numerator *= base 1355 result_digit, numerator = divmod(numerator, denominator) 1356 tails.append(str(result_digit)) 1357 if numerator not in subresults: 1358 subresults.append(numerator) 1359 else: 1360 break 1361 1362 for i in range(len(head) - 1, -1, -1): 1363 j = len(tails) - 1 1364 if head[i] != tails[j]: 1365 break 1366 else: 1367 del tails[j] 1368 tails.insert(0, head[i]) 1369 del head[i] 1370 j = j - 1 1371 1372 # truncate all leading 0's 1373 if all(elem == "0" for elem in head): 1374 for i in range(0, len(tails)): 1375 if tails[0] == "0": 1376 tails = tails[1:] + [str(0)] 1377 else: 1378 break 1379 return (head, tails) 1380 1381 1382def convert_float_base(x, base, precision=10): 1383 1384 length_of_int = 0 if x == 0 else int(mpmath.log(x, base)) 1385 # iexps = list(range(length_of_int, -1, -1)) 1386 1387 def convert_int(x, base, exponents): 1388 out = [] 1389 for e in range(0, exponents + 1): 1390 d = x % base 1391 out.append(d) 1392 x = x / base 1393 if x == 0: 1394 break 1395 out.reverse() 1396 return out 1397 1398 def convert_float(x, base, exponents): 1399 out = [] 1400 for e in range(0, exponents): 1401 d = int(x * base) 1402 out.append(d) 1403 x = (x * base) - d 1404 if x == 0: 1405 break 1406 return out 1407 1408 int_part = convert_int(int(x), base, length_of_int) 1409 if isinstance(x, (float, sympy.Float)): 1410 # fexps = list(range(-1, -int(precision + 1), -1)) 1411 real_part = convert_float(x - int(x), base, precision + 1) 1412 return int_part + real_part 1413 elif isinstance(x, int): 1414 return int_part 1415 else: 1416 raise TypeError(x) 1417 1418 1419class RealDigits(Builtin): 1420 """ 1421 <dl> 1422 <dt>'RealDigits[$n$]' 1423 <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. 1424 1425 <dt>'RealDigits[$n$, $b$]' 1426 <dd>returns a list of base_$b$ representation of the real number $n$. 1427 1428 <dt>'RealDigits[$n$, $b$, $len$]' 1429 <dd>returns a list of $len$ digits. 1430 1431 <dt>'RealDigits[$n$, $b$, $len$, $p$]' 1432 <dd>return $len$ digits starting with the coefficient of $b$^$p$ 1433 </dl> 1434 1435 Return the list of digits and exponent: 1436 >> RealDigits[123.55555] 1437 = {{1, 2, 3, 5, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0}, 3} 1438 1439 >> RealDigits[0.000012355555] 1440 = {{1, 2, 3, 5, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0}, -4} 1441 1442 >> RealDigits[-123.55555] 1443 = {{1, 2, 3, 5, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0}, 3} 1444 1445 #> RealDigits[0.004] 1446 = {{4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, -2} 1447 1448 #> RealDigits[-1.25, -1] 1449 : Base -1 is not a real number greater than 1. 1450 = RealDigits[-1.25, -1] 1451 1452 Return 25 digits of in base 10: 1453 >> RealDigits[Pi, 10, 25] 1454 = {{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} 1455 1456 #> RealDigits[19 / 7, 10, 25] 1457 = {{2, 7, 1, 4, 2, 8, 5, 7, 1, 4, 2, 8, 5, 7, 1, 4, 2, 8, 5, 7, 1, 4, 2, 8, 5}, 1} 1458 1459 Return an explicit recurring decimal form: 1460 >> RealDigits[19 / 7] 1461 = {{2, {7, 1, 4, 2, 8, 5}}, 1} 1462 1463 #> RealDigits[100 / 21] 1464 = {{{4, 7, 6, 1, 9, 0}}, 1} 1465 1466 #> RealDigits[1.234, 2, 15] 1467 = {{1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1}, 1} 1468 1469 20 digits starting with the coefficient of 10^-5: 1470 >> RealDigits[Pi, 10, 20, -5] 1471 = {{9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3, 8, 4, 6, 2, 6, 4, 3}, -4} 1472 1473 #> RealDigits[Pi, 10, 20, 5] 1474 = {{0, 0, 0, 0, 0, 3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9}, 6} 1475 1476 The 10000th digit of is an 8: 1477 >> RealDigits[Pi, 10, 1, -10000] 1478 = {{8}, -9999} 1479 1480 #> RealDigits[Pi] 1481 : The number of digits to return cannot be determined. 1482 = RealDigits[Pi] 1483 1484 #> RealDigits[20 / 3] 1485 = {{{6}}, 1} 1486 1487 #> RealDigits[3 / 4] 1488 = {{7, 5}, 0} 1489 1490 #> RealDigits[23 / 4] 1491 = {{5, 7, 5}, 1} 1492 1493 #> RealDigits[3 + 4 I] 1494 : The value 3 + 4 I is not a real number. 1495 = RealDigits[3 + 4 I] 1496 1497 #> RealDigits[abc] 1498 = RealDigits[abc] 1499 1500 #> RealDigits[abc, 2] 1501 = RealDigits[abc, 2] 1502 1503 #> RealDigits[45] 1504 = {{4, 5}, 2} 1505 1506 #> RealDigits[{3.14, 4.5}] 1507 = {{{3, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1}, {{4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1}} 1508 1509 #> RealDigits[123.45, 40] 1510 = {{3, 3, 18, 0, 0, 0, 0, 0, 0, 0}, 2} 1511 1512 #> RealDigits[0.00012345, 2] 1513 = {{1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0}, -12} 1514 1515 #> RealDigits[12345, 2, 4] 1516 = {{1, 1, 0, 0}, 14} 1517 1518 #> RealDigits[123.45, 2, 15] 1519 = {{1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1}, 7} 1520 1521 RealDigits gives Indeterminate if more digits than the precision are requested: 1522 >> RealDigits[123.45, 10, 18] 1523 = {{1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, Indeterminate, Indeterminate}, 3} 1524 1525 #> RealDigits[0.000012345, 2] 1526 = {{1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1}, -16} 1527 1528 #> RealDigits[3.14, 10, 1.5] 1529 : Non-negative machine-sized integer expected at position 3 in RealDigits[3.14, 10, 1.5]. 1530 = RealDigits[3.14, 10, 1.5] 1531 1532 #> RealDigits[3.14, 10, 1, 1.5] 1533 : Machine-sized integer expected at position 4 in RealDigits[3.14, 10, 1, 1.5]. 1534 = RealDigits[3.14, 10, 1, 1.5] 1535 1536 #> RealDigits[Pi, 10, 20, -5] 1537 = {{9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3, 8, 4, 6, 2, 6, 4, 3}, -4} 1538 1539 #> RealDigits[305.0123, 10, 17, 0] 1540 = {{5, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, Indeterminate, Indeterminate, Indeterminate}, 1} 1541 1542 #> RealDigits[220, 140] 1543 = {{1, 80}, 2} 1544 1545 # #> RealDigits[Sqrt[3], 10, 50] 1546 # = {{1, 7, 3, 2, 0, 5, 0, 8, 0, 7, 5, 6, 8, 8, 7, 7, 2, 9, 3, 5, 2, 7, 4, 4, 6, 3, 4, 1, 5, 0, 5, 8, 7, 2, 3, 6, 6, 9, 4, 2, 8, 0, 5, 2, 5, 3, 8, 1, 0, 3}, 1} 1547 1548 #> RealDigits[0] 1549 = {{0}, 1} 1550 1551 #> RealDigits[1] 1552 = {{1}, 1} 1553 1554 #> RealDigits[0, 10, 5] 1555 = {{0, 0, 0, 0, 0}, 0} 1556 1557 #> RealDigits[11/23] 1558 = {{{4, 7, 8, 2, 6, 0, 8, 6, 9, 5, 6, 5, 2, 1, 7, 3, 9, 1, 3, 0, 4, 3}}, 0} 1559 1560 #> RealDigits[1/97] 1561 = {{{1, 0, 3, 0, 9, 2, 7, 8, 3, 5, 0, 5, 1, 5, 4, 6, 3, 9, 1, 7, 5, 2, 5, 7, 7, 3, 1, 9, 5, 8, 7, 6, 2, 8, 8, 6, 5, 9, 7, 9, 3, 8, 1, 4, 4, 3, 2, 9, 8, 9, 6, 9, 0, 7, 2, 1, 6, 4, 9, 4, 8, 4, 5, 3, 6, 0, 8, 2, 4, 7, 4, 2, 2, 6, 8, 0, 4, 1, 2, 3, 7, 1, 1, 3, 4, 0, 2, 0, 6, 1, 8, 5, 5, 6, 7, 0}}, -1} 1562 1563 #> RealDigits[1/97, 2] 1564 = {{{1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0}}, -6} 1565 1566 #> RealDigits[1/197, 260, 5] 1567 = {{1, 83, 38, 71, 69}, 0} 1568 1569 #> RealDigits[1/197, 260, 5, -6] 1570 = {{246, 208, 137, 67, 80}, -5} 1571 1572 #> RealDigits[Pi, 260, 20] 1573 = {{3, 36, 211, 172, 124, 173, 210, 42, 162, 76, 23, 206, 122, 187, 23, 245, 241, 225, 254, 98}, 1} 1574 1575 #> RealDigits[Pi, 260, 5] 1576 = {{3, 36, 211, 172, 124}, 1} 1577 1578 #> RealDigits[1/3] 1579 = {{{3}}, 0} 1580 1581 #> RealDigits[1/2, 7] 1582 = {{{3}}, 0} 1583 1584 #> RealDigits[3/2, 7] 1585 = {{1, {3}}, 1} 1586 1587 #> RealDigits[-3/2, 7] 1588 = {{1, {3}}, 1} 1589 1590 #> RealDigits[3/2, 6] 1591 = {{1, 3}, 1} 1592 1593 #> RealDigits[1, 7, 5] 1594 = {{1, 0, 0, 0, 0}, 1} 1595 1596 #> RealDigits[I, 7] 1597 : The value I is not a real number. 1598 = RealDigits[I, 7] 1599 1600 #> RealDigits[-Pi] 1601 : The number of digits to return cannot be determined. 1602 = RealDigits[-Pi] 1603 1604 #> RealDigits[Round[x + y]] 1605 = RealDigits[Round[x + y]] 1606 1607 """ 1608 1609 attributes = ("Listable",) 1610 1611 messages = { 1612 "realx": "The value `1` is not a real number.", 1613 "ndig": "The number of digits to return cannot be determined.", 1614 "rbase": "Base `1` is not a real number greater than 1.", 1615 "intnm": "Non-negative machine-sized integer expected at position 3 in `1`.", 1616 "intm": "Machine-sized integer expected at position 4 in `1`.", 1617 } 1618 1619 def apply_complex(self, n, var, evaluation): 1620 "%(name)s[n_Complex, var___]" 1621 return evaluation.message("RealDigits", "realx", n) 1622 1623 def apply_rational_with_base(self, n, b, evaluation): 1624 "%(name)s[n_Rational, b_Integer]" 1625 # expr = Expression("RealDigits", n) 1626 py_n = abs(n.value) 1627 py_b = b.get_int_value() 1628 if check_finite_decimal(n.denominator().get_int_value()) and not py_b % 2: 1629 return self.apply_with_base(n, b, evaluation) 1630 else: 1631 exp = int(mpmath.ceil(mpmath.log(py_n, py_b))) 1632 (head, tails) = convert_repeating_decimal( 1633 py_n.as_numer_denom()[0], py_n.as_numer_denom()[1], py_b 1634 ) 1635 1636 leaves = [] 1637 for x in head: 1638 if not x == "0": 1639 leaves.append(from_python(int(x))) 1640 leaves.append(from_python(tails)) 1641 list_str = Expression(SymbolList, *leaves) 1642 return Expression(SymbolList, list_str, exp) 1643 1644 def apply_rational_without_base(self, n, evaluation): 1645 "%(name)s[n_Rational]" 1646 1647 return self.apply_rational_with_base(n, from_python(10), evaluation) 1648 1649 def apply(self, n, evaluation): 1650 "%(name)s[n_]" 1651 1652 # Handling the testcases that throw the error message and return the ouput that doesn't include `base` argument 1653 if isinstance(n, Symbol) and n.name.startswith("System`"): 1654 return evaluation.message("RealDigits", "ndig", n) 1655 1656 if n.is_numeric(): 1657 return self.apply_with_base(n, from_python(10), evaluation) 1658 1659 def apply_with_base(self, n, b, evaluation, nr_elements=None, pos=None): 1660 "%(name)s[n_?NumericQ, b_Integer]" 1661 1662 expr = Expression("RealDigits", n) 1663 rational_no = ( 1664 True if isinstance(n, Rational) else False 1665 ) # it is used for checking whether the input n is a rational or not 1666 py_b = b.get_int_value() 1667 if isinstance(n, (Expression, Symbol, Rational)): 1668 pos_len = abs(pos) + 1 if pos is not None and pos < 0 else 1 1669 if nr_elements is not None: 1670 n = Expression( 1671 "N", n, int(mpmath.log(py_b ** (nr_elements + pos_len), 10)) + 1 1672 ).evaluate(evaluation) 1673 else: 1674 if rational_no: 1675 n = Expression(SymbolN, n).evaluate(evaluation) 1676 else: 1677 return evaluation.message("RealDigits", "ndig", expr) 1678 py_n = abs(n.value) 1679 1680 if not py_b > 1: 1681 return evaluation.message("RealDigits", "rbase", py_b) 1682 1683 if isinstance(py_n, complex): 1684 return evaluation.message("RealDigits", "realx", expr) 1685 1686 if isinstance(n, Integer): 1687 display_len = ( 1688 int(mpmath.floor(mpmath.log(py_n, py_b))) 1689 if py_n != 0 and py_n != 1 1690 else 1 1691 ) 1692 else: 1693 display_len = int( 1694 Expression( 1695 "N", 1696 Expression( 1697 "Round", 1698 Expression( 1699 "Divide", 1700 Expression("Precision", py_n), 1701 Expression("Log", 10, py_b), 1702 ), 1703 ), 1704 ) 1705 .evaluate(evaluation) 1706 .to_python() 1707 ) 1708 1709 exp = log_n_b(py_n, py_b) 1710 1711 if py_n == 0 and nr_elements is not None: 1712 exp = 0 1713 1714 digits = [] 1715 if not py_b == 10: 1716 digits = convert_float_base(py_n, py_b, display_len - exp) 1717 # truncate all the leading 0's 1718 i = 0 1719 while digits and digits[i] == 0: 1720 i += 1 1721 digits = digits[i:] 1722 1723 if not isinstance(n, Integer): 1724 if len(digits) > display_len: 1725 digits = digits[: display_len - 1] 1726 else: 1727 # drop any leading zeroes 1728 for x in str(py_n): 1729 if x != "." and (digits or x != "0"): 1730 digits.append(x) 1731 1732 if pos is not None: 1733 temp = exp 1734 exp = pos + 1 1735 move = temp - 1 - pos 1736 if move <= 0: 1737 digits = [0] * abs(move) + digits 1738 else: 1739 digits = digits[abs(move) :] 1740 display_len = display_len - move 1741 1742 leaves = [] 1743 for x in digits: 1744 if x == "e" or x == "E": 1745 break 1746 # Convert to Mathics' list format 1747 leaves.append(from_python(int(x))) 1748 1749 if not rational_no: 1750 while len(leaves) < display_len: 1751 leaves.append(from_python(0)) 1752 1753 if nr_elements is not None: 1754 # display_len == nr_elements 1755 if len(leaves) >= nr_elements: 1756 # Truncate, preserving the digits on the right 1757 leaves = leaves[:nr_elements] 1758 else: 1759 if isinstance(n, Integer): 1760 while len(leaves) < nr_elements: 1761 leaves.append(from_python(0)) 1762 else: 1763 # Adding Indeterminate if the length is greater than the precision 1764 while len(leaves) < nr_elements: 1765 leaves.append(from_python(Symbol("Indeterminate"))) 1766 list_str = Expression(SymbolList, *leaves) 1767 return Expression(SymbolList, list_str, exp) 1768 1769 def apply_with_base_and_length(self, n, b, length, evaluation, pos=None): 1770 "%(name)s[n_?NumericQ, b_Integer, length_]" 1771 leaves = [] 1772 if pos is not None: 1773 leaves.append(from_python(pos)) 1774 expr = Expression("RealDigits", n, b, length, *leaves) 1775 if not (isinstance(length, Integer) and length.get_int_value() >= 0): 1776 return evaluation.message("RealDigits", "intnm", expr) 1777 1778 return self.apply_with_base( 1779 n, b, evaluation, nr_elements=length.get_int_value(), pos=pos 1780 ) 1781 1782 def apply_with_base_length_and_precision(self, n, b, length, p, evaluation): 1783 "%(name)s[n_?NumericQ, b_Integer, length_, p_]" 1784 if not isinstance(p, Integer): 1785 return evaluation.message( 1786 "RealDigits", "intm", Expression("RealDigits", n, b, length, p) 1787 ) 1788 1789 return self.apply_with_base_and_length( 1790 n, b, length, evaluation, pos=p.get_int_value() 1791 ) 1792 1793 1794class _ZLibHash: # make zlib hashes behave as if they were from hashlib 1795 def __init__(self, fn): 1796 self._bytes = b"" 1797 self._fn = fn 1798 1799 def update(self, bytes): 1800 self._bytes += bytes 1801 1802 def hexdigest(self): 1803 return format(self._fn(self._bytes), "x") 1804 1805 1806class Hash(Builtin): 1807 """ 1808 <dl> 1809 <dt>'Hash[$expr$]' 1810 <dd>returns an integer hash for the given $expr$. 1811 <dt>'Hash[$expr$, $type$]' 1812 <dd>returns an integer hash of the specified $type$ for the given $expr$.</dd> 1813 <dd>The types supported are "MD5", "Adler32", "CRC32", "SHA", "SHA224", "SHA256", "SHA384", and "SHA512".</dd> 1814 <dt>'Hash[$expr$, $type$, $format$]' 1815 <dd>Returns the hash in the specified format.</dd> 1816 </dl> 1817 1818 > Hash["The Adventures of Huckleberry Finn"] 1819 = 213425047836523694663619736686226550816 1820 1821 > Hash["The Adventures of Huckleberry Finn", "SHA256"] 1822 = 95092649594590384288057183408609254918934351811669818342876362244564858646638 1823 1824 > Hash[1/3] 1825 = 56073172797010645108327809727054836008 1826 1827 > Hash[{a, b, {c, {d, e, f}}}] 1828 = 135682164776235407777080772547528225284 1829 1830 > Hash[SomeHead[3.1415]] 1831 = 58042316473471877315442015469706095084 1832 1833 >> Hash[{a, b, c}, "xyzstr"] 1834 = Hash[{a, b, c}, xyzstr, Integer] 1835 """ 1836 1837 rules = { 1838 "Hash[expr_]": 'Hash[expr, "MD5", "Integer"]', 1839 "Hash[expr_, type_String]": 'Hash[expr, type, "Integer"]', 1840 } 1841 1842 attributes = ("Protected", "ReadProtected") 1843 1844 # FIXME md2 1845 _supported_hashes = { 1846 "Adler32": lambda: _ZLibHash(zlib.adler32), 1847 "CRC32": lambda: _ZLibHash(zlib.crc32), 1848 "MD5": hashlib.md5, 1849 "SHA": hashlib.sha1, 1850 "SHA224": hashlib.sha224, 1851 "SHA256": hashlib.sha256, 1852 "SHA384": hashlib.sha384, 1853 "SHA512": hashlib.sha512, 1854 } 1855 1856 @staticmethod 1857 def compute(user_hash, py_hashtype, py_format): 1858 hash_func = Hash._supported_hashes.get(py_hashtype) 1859 if hash_func is None: # unknown hash function? 1860 return # in order to return original Expression 1861 h = hash_func() 1862 user_hash(h.update) 1863 res = h.hexdigest() 1864 if py_format in ("HexString", "HexStringLittleEndian"): 1865 return from_python(res) 1866 res = int(res, 16) 1867 if py_format == "DecimalString": 1868 return from_python(str(res)) 1869 elif py_format == "ByteArray": 1870 return from_python(bytearray(res)) 1871 return from_python(res) 1872 1873 def apply(self, expr, hashtype, outformat, evaluation): 1874 "Hash[expr_, hashtype_String, outformat_String]" 1875 return Hash.compute( 1876 expr.user_hash, hashtype.get_string_value(), outformat.get_string_value() 1877 ) 1878 1879 1880class TypeEscalation(Exception): 1881 def __init__(self, mode): 1882 self.mode = mode 1883 1884 1885class Fold(object): 1886 # allows inherited classes to specify a single algorithm implementation that 1887 # can be called with machine precision, arbitrary precision or symbolically. 1888 1889 ComputationFunctions = namedtuple("ComputationFunctions", ("sin", "cos")) 1890 1891 FLOAT = 0 1892 MPMATH = 1 1893 SYMBOLIC = 2 1894 1895 math = { 1896 FLOAT: ComputationFunctions(cos=math.cos, sin=math.sin,), 1897 MPMATH: ComputationFunctions(cos=mpmath.cos, sin=mpmath.sin,), 1898 SYMBOLIC: ComputationFunctions( 1899 cos=lambda x: Expression("Cos", x), sin=lambda x: Expression("Sin", x), 1900 ), 1901 } 1902 1903 operands = { 1904 FLOAT: lambda x: None if x is None else x.round_to_float(), 1905 MPMATH: lambda x: None if x is None else x.to_mpmath(), 1906 SYMBOLIC: lambda x: x, 1907 } 1908 1909 def _operands(self, state, steps): 1910 raise NotImplementedError 1911 1912 def _fold(self, state, steps, math): 1913 raise NotImplementedError 1914 1915 def _spans(self, operands): 1916 spans = {} 1917 k = 0 1918 j = 0 1919 1920 for mode in (self.FLOAT, self.MPMATH): 1921 for i, operand in enumerate(operands[k:]): 1922 if operand[0] > mode: 1923 break 1924 j = i + k + 1 1925 1926 if k == 0 and j == 1: # only init state? then ignore. 1927 j = 0 1928 1929 spans[mode] = slice(k, j) 1930 k = j 1931 1932 spans[self.SYMBOLIC] = slice(k, len(operands)) 1933 1934 return spans 1935 1936 def fold(self, x, ll): 1937 # computes fold(x, ll) with the internal _fold function. will start 1938 # its evaluation machine precision, and will escalate to arbitrary 1939 # precision if or symbolical evaluation only if necessary. folded 1940 # items already computed are carried over to new evaluation modes. 1941 1942 yield x # initial state 1943 1944 init = None 1945 operands = list(self._operands(x, ll)) 1946 spans = self._spans(operands) 1947 1948 for mode in (self.FLOAT, self.MPMATH, self.SYMBOLIC): 1949 s_operands = [y[1:] for y in operands[spans[mode]]] 1950 1951 if not s_operands: 1952 continue 1953 1954 if mode == self.MPMATH: 1955 from mathics.core.numbers import min_prec 1956 1957 precision = min_prec(*[t for t in chain(*s_operands) if t is not None]) 1958 working_precision = mpmath.workprec 1959 else: 1960 1961 @contextmanager 1962 def working_precision(_): 1963 yield 1964 1965 precision = None 1966 1967 if mode == self.FLOAT: 1968 1969 def out(z): 1970 return Real(z) 1971 1972 elif mode == self.MPMATH: 1973 1974 def out(z): 1975 return Real(z, precision) 1976 1977 else: 1978 1979 def out(z): 1980 return z 1981 1982 as_operand = self.operands.get(mode) 1983 1984 def converted_operands(): 1985 for y in s_operands: 1986 yield tuple(as_operand(t) for t in y) 1987 1988 with working_precision(precision): 1989 c_operands = converted_operands() 1990 1991 if init is not None: 1992 c_init = tuple( 1993 (None if t is None else as_operand(from_python(t))) 1994 for t in init 1995 ) 1996 else: 1997 c_init = next(c_operands) 1998 init = tuple((None if t is None else out(t)) for t in c_init) 1999 2000 generator = self._fold(c_init, c_operands, self.math.get(mode)) 2001 2002 for y in generator: 2003 y = tuple(out(t) for t in y) 2004 yield y 2005 init = y 2006