Bernoulli doubling map

An ensemble of the doubling map x' = 2x (mod 1), iterated exactly on fractions with an odd denominator.

• Name: doubling_map
• Field: mathematics/dynamical_systems
• Formalism: discrete map
• References: Bertin 2021

Factory

@register_model(
    field='mathematics/dynamical_systems',
    tags=('chaos', 'discrete map', 'bernoulli shift', 'invariant density'),
    references=('Bertin 2021',),
    defaults={'until': 50, 'dt': 1},
)
def doubling_map(
    n_points: int = 100,
    *,
    x: float = 0.4,
    spread: float = 0.01,
    q: float = DOUBLING_DENOMINATOR,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return an ensemble of the doubling map `x' = 2x (mod 1)`, iterated exactly on fractions with an odd denominator.

    In floating point `2x (mod 1)` shifts out one binary digit per step, and every orbit reaches 0 within 53 steps.
    The model keeps every point on the fractions `k / q` of an odd denominator `q` and applies the map to them exactly,
    `k' = 2k (mod q)`, through the update `x' = (2 round(q x) mod q) / q`. The orbits are the exact orbits of these
    fractions. The default `q` is a prime below `2^50` of which 2 is a primitive root, so every orbit except that of 0
    has the period `q - 1`, about `10^15`. The points start in a narrow interval and spread until their histogram
    matches the uniform invariant density.

    Args:
        n_points:
            The number of points of the ensemble.

        x:
            The left end of the interval of initial points, which are rounded to the nearest fraction `k / q`.

        spread:
            The width of the interval of initial points.

        q:
            The odd denominator of the fractions the points stay on, an integer below `2^51`.

        random_state:
            The seed or generator of the initial points.

    Returns:
        system:
            The model, with the states `points.x`.

    Examples:
        >>> from skcomplex.models.mathematics.dynamical_systems import doubling_map
        >>> float(doubling_map(10, q=101.0).elements[0].values['x'][0] * 101)
        41.0
    """
    start = np.floor(q * _draw_cloud(n_points, x, spread, random_state) + 1 / 2) / q
    return System(
        elements=[GroupElement('points', element=DoublingMapElement, size=n_points, values={'x': start, 'q': q})],
    )

Components

class DoublingMapElement(Element):
    """Doubling map `x' = 2 x (mod 1)` on the grid of the fractions of denominator `q`.

    Args:
        label:
            The label of the element.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        x:
            The initial value, between 0 and 1.

        q:
            The denominator of the grid.

    Examples:
        >>> from skcomplex.models.mathematics.dynamical_systems import DoublingMapElement
        >>> DoublingMapElement('point', q=1024.0).q
        1024.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        x: Annotated[object, StateSpec()] = 0.4,
        q: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = DOUBLING_DENOMINATOR,
    ) -> None:
        super().__init__(label, visual=visual)
        self.x = x
        self.q = q

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the next value, rounded to the grid.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the element.
        """
        return [UpdateLaw(s.x, mod(2 * floor(s.q * s.x + 1 / 2), s.q) / s.q)]

Equations

points.x[i]: x_points(t + 1) == mod(2*floor(q_points*x_points(t) + 0.5), q_points)/q_points