"""Centers of the hyperbolic components of the Mandelbrot set -- numberdb.org/T392. Let f_c(z) = z^2 + c, set F_0(c) = 0 and F_{n+1}(c) = F_n(c)^2 + c. This generator fills T392 with the roots of the exact-period Gleason polynomial G_n(c) = F_n(c) / prod_{d|n, d 0 def _intervals_overlap(first, second): return not (first.upper() < second.lower() or second.upper() < first.lower()) def _conjugate_imaginary_intervals(first, second): return _interval_contains_zero(first.imag() + second.imag()) def _absolute_diameter(value): if value is None: return ZZ(0) diameter = getattr(value, "absolute_diameter", None) if diameter is None: return ZZ(0) return abs(diameter()) def _real_projection_clusters(roots): clusters = [] current = [] current_upper = None for root in roots: real = root.real() if not current: current = [root] current_upper = real.upper() continue if real.lower() <= current_upper: current.append(root) if real.upper() > current_upper: current_upper = real.upper() else: clusters.append(current) current = [root] current_upper = real.upper() if current: clusters.append(current) return clusters def _check_real_projection_order(period, roots): """Prove no non-conjugate roots share a real part.""" for cluster in _real_projection_clusters(roots): if len(cluster) == 1: continue if len(cluster) != 2: raise ArithmeticError( "period %d has %d roots with overlapping real intervals" % ( period, len(cluster)) ) first, second = cluster first_imag = first.imag() second_imag = second.imag() opposite_signs = ( _is_negative_interval(first_imag) and _is_positive_interval(second_imag) ) or ( _is_positive_interval(first_imag) and _is_negative_interval(second_imag) ) if not opposite_signs or not _conjugate_imaginary_intervals(first, second): raise ArithmeticError( "period %d has a non-conjugate real-part tie" % period ) def centers(period): """Return the sorted centers for one exact period.""" period = int(period) if period in _CENTER_CACHE: return _CENTER_CACHE[period] polynomial = gleason_polynomial(period) isolated = complex_roots(polynomial, min_prec=ROOT_MIN_PREC) roots = [root for root, multiplicity in isolated] if any(multiplicity != 1 for _root, multiplicity in isolated): raise ArithmeticError("period %d has a multiple root" % period) roots = sorted(roots, key=_interval_sort_key) if len(roots) != polynomial.degree(): raise ArithmeticError("period %d root count mismatch" % period) _check_real_projection_order(period, roots) real_count = polynomial.number_of_real_roots() isolated_real_count = sum( 1 for root in roots if _imaginary_part_contains_zero(root) ) if isolated_real_count != real_count: raise ArithmeticError( "period %d has %d certified real roots, expected %d" % ( period, isolated_real_count, real_count) ) records = [] for index, root in enumerate(roots, start=1): root_is_real = _imaginary_part_contains_zero(root) if period == 1: value = ZZ(0) elif period == 2: value = ZZ(-1) elif root_is_real: value = root.real() else: value = root record = {"number": value} equals = EQUALS.get((period, index)) if equals: record["equals"] = equals records.append(record) _CENTER_CACHE[period] = records return records class MandelbrotHyperbolicComponentCenters(numberdb.Generator): """Generator for T392, centers of Mandelbrot hyperbolic components.""" table = TABLE parameters = ("n", "i") type = "C" digits = DIGITS rigour = "proven" def enumerate(self): for period in range(1, MAX_PERIOD + 1): for index in range(1, len(centers(period)) + 1): yield {"n": str(period), "i": str(index)} def value(self, params, digits): period = ZZ(params["n"]) index = ZZ(params["i"]) return centers(period)[int(index) - 1] def run_integrity_checks(): fingerprints = { 1: c, 2: c + 1, 3: c ** 3 + 2 * c ** 2 + c + 1, 4: c ** 6 + 3 * c ** 5 + 3 * c ** 4 + 3 * c ** 3 + 2 * c ** 2 + 1, } for period, expected in fingerprints.items(): if gleason_polynomial(period) != expected: raise ArithmeticError("G_%d does not match the convention fingerprint" % period) total = 0 real_counts = [] worst_diameter = ZZ(0) worst_at = None for period in range(1, MAX_PERIOD + 1): polynomial = gleason_polynomial(period) if polynomial.degree() != expected_degree(period): raise ArithmeticError("G_%d has degree %d, expected %d" % ( period, polynomial.degree(), expected_degree(period))) real_counts.append(int(polynomial.number_of_real_roots())) for index, record in enumerate(centers(period), start=1): total += 1 value = record["number"] real = value.real() if hasattr(value, "real") else value imag = value.imag() if hasattr(value, "imag") else None diameter = max(_absolute_diameter(real), _absolute_diameter(imag)) if diameter > worst_diameter: worst_diameter = diameter worst_at = (period, index) expected_reals = [1, 1, 1, 2, 3, 5, 9, 16, 28] if real_counts != expected_reals: raise ArithmeticError("real root counts %s, expected %s" % ( real_counts, expected_reals)) print("integrity checks passed for %d centers through period %d" % ( total, MAX_PERIOD)) print("real root counts: %s" % ",".join(str(n) for n in real_counts)) print("worst root diameter is %s at n=%s, i=%s" % ( worst_diameter, worst_at[0], worst_at[1])) print("no non-conjugate real-part ties occur through n=%d" % MAX_PERIOD) def fill_draft_once(generator, message): """Fill a fresh prose draft without the client's empty upsert probe.""" from numberdb._generate import ( _check_precision, _check_rigour, _producer, _run_name, _source_files, ) from numberdb._write import Entries, attach, submit_entries, to_text table = generator.table run = _run_name(generator) entries = Entries(*generator.parameters) for params in generator.enumerate(): params = dict(params) wanted = generator.digits_for(params) entry = generator._entry(params, wanted) value = entry["number"] identity = ",".join(str(params[name]) for name in generator.parameters) _check_rigour(generator, table, identity, value) written = to_text(value, wanted, generator.format) _check_precision(table, identity, written, wanted, lowering=False) record = dict(entry) record.pop("digits", None) entries.add(**params, **record, digits=wanted) answer = submit_entries( table, entries, message=message, produced_by=_producer(generator, os.environ.get("NUMBERDB_ASSISTED_BY", "")), upsert=False, run=run, rigour=generator.rigour, ) for name, body in sorted(_source_files(generator).items()): attach(table, name, body, run=run, message=message, rigour=generator.rigour) return answer if __name__ == "__main__": _key_from_stdin() generator = MandelbrotHyperbolicComponentCenters() run_integrity_checks() if "--publish" in sys.argv or os.environ.get("NUMBERDB_PUBLISH") == "1": print(fill_draft_once( generator, message="fill Mandelbrot hyperbolic component centers", )) else: report = generator.verify(sample=None) print(report) sys.exit(0 if report.ok else 1)