File Coverage

blib/lib/Math/DifferenceSet/Planar.pm

Criterion	Covered	Total	%
statement	609	609	100.0
branch	256	262	97.7
condition	66	72	91.6
subroutine	87	87	100.0
pod	48	48	100.0
total	1066	1078	98.8

package Math::DifferenceSet::Planar;

use strict;
use warnings;
use Carp qw(croak);
use Math::DifferenceSet::Planar::Data;
use Math::BigInt try => 'GMP';
use Math::Prime::Util qw(
    is_power is_prime_power euler_phi factor_exp gcd
    mulmod addmod invmod powmod divmod
);

# Math::DifferenceSet::Planar=ARRAY(...)
# ........... index ...........     # ........... value ...........
use constant _F_ORDER      =>  0;
use constant _F_BASE       =>  1;
use constant _F_EXPONENT   =>  2;
use constant _F_MODULUS    =>  3;
use constant _F_N_PLANES   =>  4;   # number of distinct planes this size
use constant _F_ELEMENTS   =>  5;   # elements arrayref
use constant _F_ROTATORS   =>  6;   # rotators arrayref, initially empty
use constant _F_INDEX_MIN  =>  7;   # index of smallest element
use constant _F_GENERATION =>  8;   # database generation
use constant _F_LOG        =>  9;   # plane logarithm value, may be undef
use constant _F_ZETA       => 10;   # "zeta" value, initially undef
use constant _F_ETA        => 11;   # "eta" value, initially undef
use constant _F_PEAK       => 12;   # peak elements arrayref, initially undef
use constant _NFIELDS      => 13;

our $VERSION = '0.018';

our $_LOG_MAX_ORDER  = 22.1807;         # limit for integer exponentiation
our $_MAX_ENUM_COUNT = 32768;           # limit for stored rotator set size
our $_MAX_MEMO_COUNT = 4096;            # limit for memoized values
our $_DEFAULT_DEPTH  = 1024;            # default check_elements() depth
our $_USE_SPACES_DB  = 1;               # enable looking up rotators
our $_MAX_FMT_COUNT  = 5;               # elements printed in messages

my $current_data = undef;               # current M::D::P::Data object
my $generation   = 0;                   # incremented on $current_data update

my %memo_n_planes = ();                 # memoized n_planes values

# ----- private subroutines -----

# calculate powers of p (mod m)
sub _multipliers {
    my ($base, $exponent, $modulus) = @_;
    my @mult = (1);
    my $n = 3 * $exponent;
    my $x = $base;
    while (@mult < $n && $x != 1) {
        push @mult, $x;
        $x = mulmod($x, $base, $modulus);
    }
    push @mult, q[...] if $x != 1;
    die "assertion failed: not $n elements: @mult" if @mult != $n;
    return @mult;
}

# return complete rotator base if known or small enough, otherwise undef
#
# If structure is known:
#
# rotators -> +-----------+
#             | radices ----------------------------> +-------+
#             +-----------+                           |  r_1  |
#             | depths  -----------------> +-------+  +-------+
#             +-----------+                |  d_1  |  |  r_2  |
#             | inverses -----> +-------+  +-------+  +-------+
#             +-----------+     |  i_1  |  |  d_2  |  |  ...  |
#                               +-------+  +-------+  +-------+
#                               |  i_2  |  |  ...  |  |  r_n  |
#                               +-------+  +-------+  +-------+
#                               |  ...  |  |  d_n  |
#                               +-------+  +-------+
#                               |  i_n  |
#                               +-------+
#
#                         n
#                        ___      j
#   R               _    | |   r   k             0  <  j   <  d
#    j j ...j       =    | |    k      (mod M),     =   k      k
#     1 2    n
#                       k = 1
#
#         d  - 1   _
#   i  r   k       =   1  (mod M),   d   >  d   >  ...  >  d   >  2
#    k  k                             1  =   2  =       =   n  =
#
#
# Otherwise, if number of rotators is small:
#
# rotators -> +-------+-------+--   --+-------+
#             |  R_1  |  R_2  |  ...  |  R_N  |
#             +-------+-------+--   --+-------+
#
sub _rotators {
    my ($this) = @_;
    my (           $base,   $exponent,   $order,   $modulus,   $rotators) =
        @{$this}[_F_BASE, _F_EXPONENT, _F_ORDER, _F_MODULUS, _F_ROTATORS];
    return $rotators if @{$rotators};
    my $space = $_USE_SPACES_DB && $this->_data->get_space($order);
    if ($space) {
        my ($radices, $depths) = $space->rotator_space;
        my $inverses = [
            map {
                invmod(
                    powmod($radices->[$_], $depths->[$_] - 1, $modulus),
                    $modulus
                )
            } 0 .. $#{$radices}
        ];
        $rotators = $this->[_F_ROTATORS] = [$radices, $depths, $inverses];
        return $rotators;
    }
    return undef if $this->[_F_N_PLANES] > $_MAX_ENUM_COUNT;
    my @mult = _multipliers($base, $exponent, $modulus);
    my $sieve = '1' x $modulus;
    substr($sieve, $_, 1) = '0' for @mult;
    @{$rotators} = (1);
    for (my $x = 2; $x < $modulus ; ++$x) {
        if (substr $sieve, $x, 1) {
            if (0 == $modulus % $x) {
                for (my $i = $x; $i < $modulus; $i += $x) {
                    substr($sieve, $i, 1) = '0';
                }
                next;
            }
            substr($sieve, $_, 1) = '0'
                for map { mulmod($_, $x, $modulus) } @mult;
            push @{$rotators}, $x;
        }
    }
    return $rotators;
}

# iterative rotator base generator, slow, but memory efficient
sub _sequential_rotators {
    my ($this) = @_;
    my (           $base,   $exponent,   $modulus,   $n_planes) =
        @{$this}[_F_BASE, _F_EXPONENT, _F_MODULUS, _F_N_PLANES];
    my @mult = _multipliers($base, $exponent, $modulus);
    shift @mult;
    my @pf = map { $_->[0] } factor_exp($modulus);
    pop @pf if $pf[-1] == $modulus;
    my $mx = 0;
    my $x  = 0;
    return sub {
        return 0 if $mx >= $n_planes;
        ELEMENT:
        while (1) {
            ++$x;
            foreach my $p (@pf) {
                next ELEMENT if !($x % $p);
            }
            foreach my $e (@mult) {
                next ELEMENT if mulmod($x, $e, $modulus) < $x;
            }
            ++$mx;
            return $x;
        }
    };
}

# structured rotator base iterator, time and space efficient
sub _structured_rotators {
    my ($this) = @_;
    my $modulus = $this->[_F_MODULUS];
    my ($radices, $depths, $inverses) = @{ $this->[_F_ROTATORS] };
    my @index = (0) x @{$radices};
    my $next = 1;
    return sub {
        return 0 if !$next;
        my $element = $next;
        my $i = 0;
        while ($i < @index) {
            if (++$index[$i] < $depths->[$i]) {
                $next = mulmod($next, $radices->[$i], $modulus);
                return $element;
            }
        }
        continue {
            $index[$i] = 0;
            $next = mulmod($next, $inverses->[$i], $modulus);
            ++$i;
        }
        $next = 0;
        return $element;
    };
}

sub _space_description {
    my ($spc) = @_;
    my $order       = $spc->order;
    my $mul_radix   = $spc->mul_radix;
    my $mul_depth   = $spc->mul_depth;
    my ($radices, $depths) = $spc->rotator_space;
    my @space = map {; "$radices->[$_]^$depths->[$_]"} 0 .. $#{$radices};
    return "$order: $mul_radix^$mul_depth [@space]";
}

# integer exponentiation
sub _pow {
    my ($base, $exponent) = @_;
    return 0 if $base <= 0 || $exponent < 0;
    return Math::BigInt->new($base)->bpow($exponent)
        if log(0+$base) * $exponent > $_LOG_MAX_ORDER;
    my $power = 1;
    while ($exponent) {
        $power *= $base if 1 & $exponent;
        $exponent >>= 1 and $base *= $base;
    }
    return $power;
}

#      0  1  2  3  4  5  6  7  8  9  10
#      3  7 14 25 28 29 31 41 61 99 103
#      7  8  9 10  0  1  2  3  4  5   6
# re-arrange elements of a planar difference set in conventional order
sub _sort_elements {
    my $modulus  = shift;
    my @elements = sort { $a <=> $b || croak "duplicate element: $a" } @_;
    my $lo = $elements[-1] - $modulus;
    my $hi = $elements[0];
    my $mx = 0;
    while ($hi - $lo != 1 && ++$mx < @elements) {
        ($lo, $hi) = ($hi, $elements[$mx]);
    }
    croak "delta 1 elements missing" if $mx >= @elements;
    $mx += @elements if !$mx--;
    push @elements, splice @elements, 0, $mx if $mx;
    return (\@elements, $mx? @elements - $mx: 0);
}

#      0  1  2  3  4  5  6  7
#     28 29 31 41  3  7 14 25
#      L        X           H
#                  L  X     H
#                 LX  H
# get index of smallest element (using bisection)
sub _index_min {
    my ($elements) = @_;
    my ($lx, $hx) = (0, $#$elements);
    my $he = $elements->[$hx];
    while ($lx < $hx) {
        my $x = ($lx + $hx) >> 1;
        my $e = $elements->[$x];
        if ($e < $he) {
            $hx = $x;
            $he = $e;
        }
        else {
            $lx = $x + 1;
        }
    }
    return $hx;
}

# return data connection, creating it if not yet open
sub _data {
    if (!defined $current_data) {
        ++$generation;
        $current_data = Math::DifferenceSet::Planar::Data->new;
    }
    return $current_data;
}

# compose a printable abbreviated description of a set
sub _fmt_set {
    my ($this) = @_;
    my @e = $this->[_F_ORDER] < $_MAX_FMT_COUNT?
        @{$this->[_F_ELEMENTS]}:
        (@{$this->[_F_ELEMENTS]}[0 .. $_MAX_FMT_COUNT-1], q[...]);
    return '{' . join(q[,], @e) . '}';
}

# identify a plane by its rotator value with respect to a given plane
# (both arguments zeta-canonized, second argument with known log)
sub _log {
    my ($this, $that) = @_;
    my $modulus = $this->[_F_MODULUS];
    my %t = ();
    foreach my $e (@{$this->[_F_ELEMENTS]}) {
        $t{$e} = 1;
    }
    my $r = 0;
    foreach my $e (@{$that->[_F_ELEMENTS]}) {
        $r = invmod($e, $modulus);
        last if $r;
    }
    my $factor = 0;
    if ($r) {
        ELEM:
        foreach my $o (@{$this->[_F_ELEMENTS]}) {
            next if !$o;
            my $ro = mulmod($r, $o, $modulus);
            foreach my $e (@{$that->[_F_ELEMENTS]}) {
                next ELEM if !exists $t{ mulmod($e, $ro, $modulus) };
            }
            $factor = $ro;
            last;
        }
    }
    else {
        my $ri = $this->iterate_rotators;
        ROT:
        while (my $ro = $ri->()) {
            foreach my $e (@{$that->[_F_ELEMENTS]}) {
                next ROT if !exists $t{ mulmod($e, $ro, $modulus) };
            }
            $factor = $ro;
            last;
        }
    }
    if ($factor) {
        my $log = $this->[_F_LOG] =
            mulmod($that->[_F_LOG], $factor, $modulus);
        return $log;
    }
    croak 'unaligned sets: ', _fmt_set($this), ' versus ', _fmt_set($that);
}

# $factor = _find_factor($ds1, $ds2);
sub _find_factor {
    my ($this, $that) = @_;
    my $order = $this->order;
    croak 'sets of same size expected' if $order != $that->order;
    my $log_this = $this->[_F_GENERATION] == $generation && $this->[_F_LOG];
    my $log_that = $this->[_F_GENERATION] == $generation && $that->[_F_LOG];
    $this->[_F_GENERATION] = $that->[_F_GENERATION] = $generation;
    if (!$log_this) {
        my $r1 = $this->zeta_canonize;
        my $r2 = $that->zeta_canonize;
        if (!$log_that) {
            my $r3 = Math::DifferenceSet::Planar->new($order)->zeta_canonize;
            $log_that = $that->[_F_LOG] = _log($r2, $r3);
        }
        $log_this = $this->[_F_LOG] =
            $this == $that? $log_that: _log($r1, $r2);
    }
    elsif (!$log_that) {
        my $r1 = $this->zeta_canonize;
        my $r2 = $that->zeta_canonize;
        $log_that = $that->[_F_LOG] = _log($r2, $r1);
    }
    return divmod($log_that, $log_this, $this->modulus);
}

# translation amount between a multiple of a set and another set
sub _delta_f {
    my ($this, $factor, $that) = @_;
    my $modulus = $this->modulus;
    my ($x)     = $this->find_delta( invmod($factor, $modulus) );
    my $s       = $that->[_F_ELEMENTS]->[0];
    return addmod($s, -mulmod($x, $factor, $modulus), $modulus);
}

# ----- class methods -----

sub list_databases {
    return Math::DifferenceSet::Planar::Data->list_databases;
}

sub set_database {
    my $class = shift;
    ++$generation;
    $current_data = Math::DifferenceSet::Planar::Data->new(@_);
    return $class->available_count;
}

# print "ok" if Math::DifferenceSet::Planar->available(9);
# print "ok" if Math::DifferenceSet::Planar->available(3, 2);
sub available {
    my ($class, $base, $exponent) = @_;
    my $order = defined($exponent)? _pow($base, $exponent): $base;
    return undef if !$order || $order > $class->_data->max_order;
    my $pds   = $class->_data->get($order, 'base');
    return !!$pds && (!defined($exponent) || $base == $pds->base);
}

# print "ok" if Math::DifferenceSet::Planar->known_space(9);
sub known_space {
    my ($class, $order) = @_;
    return 0 if $order <= 0 || $order > $class->_data->sp_max_order;
    my $spc = $class->_data->get_space($order);
    return 0 if !$spc;
    my ($rad) = $spc->rotator_space;
    return 0 + @{$rad};
}

# $desc = Math::DifferenceSet::Planar->known_space_desc(9);
sub known_space_desc {
    my ($class, $order) = @_;
    return undef if $order <= 0 || $order > $class->_data->sp_max_order;
    my $spc = $class->_data->get_space($order);
    return undef if !$spc;
    return _space_description($spc);
}

# $ds = Math::DifferenceSet::Planar->new(9);
# $ds = Math::DifferenceSet::Planar->new(3, 2);
sub new {
    my ($class, $base, $exponent) = @_;
    my $order = defined($exponent)? _pow($base, $exponent): $base;
    my $pds   = $order && $class->_data->get($order);
    if (!$pds || defined($exponent) && $base != $pds->base) {
        my $key = defined($exponent)? "$base, $exponent": $order;
        croak "PDS($key) not available";
    }
    return bless [
        $pds->order,
        $pds->base,
        $pds->exponent,
        $pds->modulus,
        $pds->n_planes,
        $pds->elements,
        [],                     # rotators
        0,                      # index_min
        $generation,
        1,                      # log
    ], $class;
}

# $ds = Math::DifferenceSet::Planar->from_elements_fast(
#   0, 1, 3, 9, 27, 49, 56, 61, 77, 81
# );
sub from_elements_fast {
    my $class    = shift;
    my $order    = $#_;
    my ($base, $exponent);
    $exponent    = is_prime_power($order, \$base)
        or croak "this implementation cannot handle order $order";
    my $modulus  = ($order + 1) * $order + 1;
    if (grep { $_ < 0 || $modulus <= $_ } @_) {
        my $max = $modulus - 1;
        croak "element values inside range 0..$max expected";
    }
    my ($elements, $index_min) = _sort_elements($modulus, @_);
    my $n_planes;
    if (!exists $memo_n_planes{$order}) {
        $n_planes = $memo_n_planes{$order} =
            euler_phi($modulus) / (3 * $exponent);
        %memo_n_planes = ($order => $n_planes)
            if $_MAX_MEMO_COUNT < keys %memo_n_planes;
    }
    else {
        $n_planes = $memo_n_planes{$order};
    }
    return bless [
        $order,
        $base,
        $exponent,
        $modulus,
        $n_planes,
        $elements,
        [],                     # rotators
        $index_min,
        $generation,
    ], $class;
}

# $ds = Math::DifferenceSet::Planar->from_elements(
#   0, 1, 3, 9, 27, 49, 56, 61, 77, 81
# );
sub from_elements {
    my ($class, @elements) = @_;
    my $this = $class->from_elements_fast(@elements);
    my $ref  = $class->new(@elements - 1);
    eval { _find_factor($ref, $this) } or
        croak "apparently not a planar difference set: ", _fmt_set($this);
    return $this;
}

# $bool = Math::DifferenceSet::Planar->verify_elements(
#   0, 1, 3, 9, 27, 49, 56, 61, 77, 81
# );
sub verify_elements {
    my ($class, @elements) = @_;
    my $order   = $#elements;
    return undef if $order <= 1;
    my $modulus = ($order + 1) * $order + 1;
    my $median  = ($modulus - 1) / 2;
    my $seen    = '0' x $median;
    foreach my $r1 (@elements) {
        return undef if $r1 < 0 || $modulus <= $r1 || $r1 != int $r1;
        foreach my $r2 (@elements) {
            last if $r1 == $r2;
            my $d = $r1 < $r2? $r2 - $r1: $modulus + $r2 - $r1;
            $d = $modulus - $d if $d > $median;
            return !1 if substr($seen, $d-1, 1)++;
        }
    }
    return $median == $seen =~ tr/1//;
}

# $bool = Math::DifferenceSet::Planar->check_elements(
#   [0, 1, 3, 9, 27, 49, 56, 61, 77, 81], 999
# );
sub check_elements {
    my ($class, $elem_listref, $depth, $factor) = @_;
    my $order        = $#{$elem_listref};
    return undef if $order <= 1;
    my $modulus      = ($order + 1) * $order + 1;
    my $median       = ($modulus - 1) / 2;
    $depth           = $median <= $_DEFAULT_DEPTH? $median: $_DEFAULT_DEPTH
        if !$depth;
    my $limit        = $median <= $depth? $median: $depth;
    my $seen         = '1' . ('0' x $limit);
    foreach my $e (@{$elem_listref}) {
        return undef if $e < 0 || $modulus <= $e || $e != int $e;
    }
    my @elements  = sort { $a <=> $b } @{$elem_listref};
    I1:
    foreach my $i1 (0 .. $order) {
        my $r1 = $elements[$i1];
        my $i2 = $i1 - 1;
        while ($i2 >= 0) {
            my $r2 = $elements[$i2];
            my $d = $r1 - $r2;
            next I1 if $d > $limit;
            return !1 if substr($seen, $d, 1)++;
            --$i2;
        }
        $i2 = $order;
        while ($i2 > $i1) {
            my $r2 = $elements[$i2];
            my $d = $modulus + $r1 - $r2;
            next I1 if $d > $limit;
            return !1 if substr($seen, $d, 1)++;
            --$i2;
        }
    }
    return !1 if 0 <= index $seen, '0';
    return 2 if $median <= $depth;
    if (!defined $factor) {
        $factor = $order if !is_power($order, 0, \$factor);
    }
    if (1 < $factor) {
        my $mx = $order;
        for (my $i = 0; $i < $order; ++$i) {
            $mx = $i, last if 1 == $elements[$i+1] - $elements[$i];
        }
        push @elements, splice @elements, 0, $mx if $mx;
        my ($multiple) = eval { _sort_elements(
            $modulus,
            map { mulmod($_, $factor, $modulus) } @elements
        )};
        return 0 if !defined $multiple;
        my $d1 = $elements[0] - $multiple->[0];
        my $d2 = $d1 >= 0? $d1 - $modulus: $d1 + $modulus;
        ($d1, $d2) = ($d2, $d1) if abs($d2) > abs($d1);         # optimization
        for (my $i = 1; $i <= $order; ++$i) {
            my $d = $elements[$i] - $multiple->[$i];
            return 0 if $d != $d1 && $d != $d2;
        }
    }
    return 1;
}

# $it1 = Math::DifferenceSet::Planar->iterate_available_sets;
# $it2 = Math::DifferenceSet::Planar->iterate_available_sets(10, 20);
# while (my $ds = $it2->()) {
#   ...
# }
sub iterate_available_sets {
    my ($class, @minmax) = @_;
    my $dit = $class->_data->iterate(@minmax);
    return sub {
        my $pds = $dit->();
        return undef if !$pds;
        return bless [
            $pds->order,
            $pds->base,
            $pds->exponent,
            $pds->modulus,
            $pds->n_planes,
            $pds->elements,
            [],                 # rotators
            0,                  # index_min
            $generation,
            1,                  # log
        ], $class;
    };
}

# $min = Math::DifferenceSet::Planar->available_min_order;
sub available_min_order   { $_[0]->_data->min_order }

# $max = Math::DifferenceSet::Planar->available_max_order;
sub available_max_order   { $_[0]->_data->max_order }

# $count = Math::DifferenceSet::Planar->available_count;
sub available_count       { $_[0]->_data->count }

# $min = Math::DifferenceSet::Planar->known_space_min_order;
sub known_space_min_order { $_[0]->_data->sp_min_order }

# $max = Math::DifferenceSet::Planar->known_space_max_order;
sub known_space_max_order { $_[0]->_data->sp_max_order }

# $count = Math::DifferenceSet::Planar->known_space_count;
sub known_space_count     { $_[0]->_data->sp_count }

# $it3 = Math::DifferenceSet::Planar->iterate_known_spaces;
# $it3 = Math::DifferenceSet::Planar->iterate_known_spaces(10,20);
# while (my $spc = $it3->()) {
#   print "$spc\n";
# }
sub iterate_known_spaces {
    my ($class, @minmax) = @_;
    my $dit = $class->_data->iterate_spaces(@minmax);
    return sub {
        my $spc = $dit->();
        return undef if !$spc;
        return _space_description($spc);
    };
}

# ----- object methods -----

# $o  = $ds->order;
# $p  = $ds->order_base;
# $n  = $ds->order_exponent;
# $m  = $ds->modulus;
# $np = $ds->n_planes;
# @e  = $ds->elements;
# $e0 = $ds->element(0);
sub order          {    $_[0]->[_F_ORDER   ]            }
sub order_base     {    $_[0]->[_F_BASE    ]            }
sub order_exponent {    $_[0]->[_F_EXPONENT]            }
sub modulus        {    $_[0]->[_F_MODULUS ]            }
sub n_planes       {    $_[0]->[_F_N_PLANES]            }
sub elements       { @{ $_[0]->[_F_ELEMENTS]          } }
sub element        {    $_[0]->[_F_ELEMENTS]->[$_[1]]   }
sub start_element  {    $_[0]->[_F_ELEMENTS]->[   0 ]   }

sub min_element { $_[0]->[_F_ELEMENTS]->[$_[0]->[_F_INDEX_MIN]  ] }
sub max_element { $_[0]->[_F_ELEMENTS]->[$_[0]->[_F_INDEX_MIN]-1] }

# @e  = $ds->elements_sorted;
sub elements_sorted {
    my ($this) = @_;
    my $elements = $this->[_F_ELEMENTS];
    return @$elements if !wantarray;
    my $index_min = $this->[_F_INDEX_MIN];
    return @{$elements}[$index_min .. $#$elements, 0 .. $index_min - 1];
}

# $ds1 = $ds->translate(1);
sub translate {
    my ($this, $delta) = @_;
    my $modulus = $this->[_F_MODULUS];
    $delta %= $modulus unless -$modulus < $delta && $delta < $modulus;
    return $this if !$delta;
    my @elements =
        map { addmod($_, $delta, $modulus) } @{$this->[_F_ELEMENTS]};
    my $that = bless [@{$this}], ref $this;
    $that->[_F_ELEMENTS]  = \@elements;
    $that->[_F_INDEX_MIN] =
        $elements[0] < $elements[-1]? 0: _index_min(\@elements);
    if (defined (my $z = $that->[_F_ZETA])) {
        $that->[_F_ZETA] = addmod($z,
            mulmod($delta, $this->[_F_ORDER]-1, $modulus), $modulus);
    }
    if (defined (my $e = $that->[_F_ETA])) {
        $that->[_F_ETA] = addmod($e,
            mulmod($delta, $this->[_F_BASE]-1, $modulus), $modulus);
    }
    $this->[_F_PEAK] = undef;
    return $that;
}

# $ds2 = $ds->canonize;
sub canonize { $_[0]->translate(- $_[0]->[_F_ELEMENTS]->[0]) }

# $ds2 = $ds->gap_canonize;
sub gap_canonize {
    my ($this) = @_;
    my $delta = ($this->largest_gap)[1];
    return $this->translate(-$delta);
}

# $ds2 = $ds->zeta_canonize;
sub zeta_canonize {
    my ($this) = @_;
    my $zeta    = $this->zeta;
    my $order   = $this->order;
    my $modulus = $this->modulus;
    if ( ($order - 1) % 3 ) {
        return $this if !$zeta;
        my $delta = divmod($zeta, $order - 1, $modulus);
        return $this->translate(-$delta);
    }
    return $this if !$zeta && !$this->contains(0);
    $modulus /= 3;
    my $delta = divmod($zeta, $order - 1, $modulus);
    $delta += $modulus while $this->contains($delta);   # 0..2 iterations
    return $this->translate(-$delta);
}

# $it  = $ds->iterate_rotators;
# while (my $m = $it->()) {
#   ...
# }
sub iterate_rotators {
    my ($this) = @_;
    my $rotators = $this->_rotators;
    return $this->_sequential_rotators if !$rotators;
    return $this->_structured_rotators if 'ARRAY' eq ref $rotators->[0];
    my $mx = 0;
    return sub { $mx < @{$rotators}? $rotators->[$mx++]: 0 };
}

# $it = $ds->iterate_planes;
# while (my $ds = $it->()) {
#   ...
# }
sub iterate_planes {
    my ($this) = @_;
    my $r_it = $this->iterate_rotators;
    return sub {
        my $r = $r_it->();
        return $r? $this->multiply($r)->canonize: undef;
    };
}

# @pm = $ds->multipliers;
sub multipliers {
    my ($this) = @_;
    my (          $base,  $exponent,  $modulus) =
        @{$this}[_F_BASE, _F_EXPONENT, _F_MODULUS];
    return 3 * $exponent if !wantarray;
    return _multipliers($base, $exponent, $modulus);
}

# $ds3 = $ds->multiply($m);
sub multiply {
    my ($this, $factor) = @_;
    my $modulus = $this->[_F_MODULUS];
    $factor %= $modulus;
    return $this if 1 == $factor;
    croak "$_[1]: factor is not coprime to modulus"
        if gcd($modulus, $factor) != 1;
    my ($elements, $index_min) = _sort_elements(
        $modulus,
        map { mulmod($_, $factor, $modulus) } @{$this->[_F_ELEMENTS]}
    );
    my $that = bless [@{$this}], ref $this;
    $that->[_F_ELEMENTS ] = $elements;
    $that->[_F_INDEX_MIN] = $index_min;
    if (defined(my $log = $that->[_F_LOG])) {
        $that->[_F_LOG] = mulmod($log, $factor, $modulus);
    }
    $that->[_F_ZETA] &&= mulmod($that->[_F_ZETA], $factor, $modulus);
    $that->[_F_ETA]  &&= mulmod($that->[_F_ETA],  $factor, $modulus);
    $this->[_F_PEAK] = undef;
    return $that;
}

# ($e1, $e2) = $ds->find_delta($delta);
sub find_delta {
    my ($this, $delta) = @_;
    my $order    = $this->order;
    my $modulus  = $this->modulus;
    my $elements = $this->[_F_ELEMENTS];
    my $de = $delta % $modulus;
    my $up = $de + $de < $modulus;
    $de = $modulus - $de if !$up;
    my ($lx, $ux, $c) = (0, 0, 0);
    my ($le, $ue) = @{$elements}[0, 0];
    my $bogus = 0;
    while ($c != $de) {
        if ($c < $de) {
            if(++$ux > $order) {
                $ux = 0;
            }
            $ue = $elements->[$ux];
            $bogus = 1, last if $ux == $lx;
        }
        else {
            $bogus = 1, last if ++$lx > $order;
            $le = $elements->[$lx];
        }
        $c = $ue < $le? $modulus + $ue - $le: $ue - $le;
    }
    croak "bogus set: delta not found: $delta (mod $modulus)" if $bogus;
    return $up? ($le, $ue): ($ue, $le);
}

# ($e1, $e2) = $ds->peak_elements
sub peak_elements {
    my ($this) = @_;
    my $peak = $this->[_F_PEAK];
    if (!defined $peak) {
        $peak = [$this->find_delta($this->modulus >> 1)];
        $this->[_F_PEAK] = $peak;
    }
    return @{$peak};
}

# ($e1, $e2, $delta) = $ds->largest_gap;
sub largest_gap {
    my ($this) = @_;
    my $modulus = $this->modulus;
    my $p_max;
    my $e_max;
    my $d_max = 0;
    my $e_pre = $this->element($this->order);
    foreach my $e ($this->elements) {
        my $d = $e - $e_pre;
        $d += $modulus if $d < 0;
        if ($d > $d_max) {
            $d_max = $d;
            $p_max = $e_pre;
            $e_max = $e;
        }
        $e_pre = $e;
    }
    return ($p_max, $e_max, $d_max);
}

# $e = $ds->zeta
sub zeta {
    my ($this) = @_;
    my $zeta = $this->[_F_ZETA];
    if (!defined $zeta) {
        my $order   = $this->[_F_ORDER];
        my $modulus = $this->[_F_MODULUS];
        my $start   = $this->[_F_ELEMENTS]->[0];
        (undef, my $x) = $this->find_delta($order + 1);
        $zeta = $this->[_F_ZETA] =
            addmod(mulmod($x, $order, $modulus), -$start, $modulus);
    }
    return $zeta;
}

# $e = $ds->eta
sub eta {
    my ($this) = @_;
    return $this->zeta if $this->[_F_EXPONENT] == 1;
    my $eta = $this->[_F_ETA];
    if (!defined $eta) {
        my $p   = $this->[_F_BASE];
        my $m   = $this->[_F_MODULUS];
        my $s   = $this->[_F_ELEMENTS]->[0];
        my ($x) = $this->find_delta( invmod($p, $m) );
        $eta = $this->[_F_ETA] = addmod(mulmod($x, $p, $m), -$s, $m);
    }
    return $eta;
}

# $bool = $ds->contains($e)
sub contains {
    my ($this, $elem) = @_;
    my $elements  = $this->[_F_ELEMENTS];
    my $lx        = $this->[_F_INDEX_MIN];
    my $hx        = $#{$elements};
    ($lx, $hx) = (0, $lx - 1) if $lx && $elements->[0] <= $elem;
    while ($lx <= $hx) {
        my $x = ($lx + $hx) >> 1;
        my $e = $elements->[$x];
        if ($e <= $elem) {
            return !0 if $e == $elem;
            $lx = $x + 1;
        }
        else {
            $hx = $x - 1;
        }
    }
    return !1;
}

# $cmp = $ds1->compare($ds2);
sub compare {
    my ($this, $that) = @_;
    my $order = $this->order;
    my $cmp   = $order <=> $that->order;
    return $cmp if $cmp;
    my $lx = $this->[_F_INDEX_MIN];
    my $rx = $that->[_F_INDEX_MIN];
    my $le = $this->[_F_ELEMENTS];
    my $re = $that->[_F_ELEMENTS];
    foreach my $i (0 .. $order) {
        $cmp = $le->[$lx] <=> $re->[$rx];
        return $cmp if $cmp;
        ++$lx <= $order or $lx = 0;
        ++$rx <= $order or $rx = 0;
    }
    return 0;
}

# $bool = $ds1->same_plane($ds2);
sub same_plane {
    my ($this, $that) = @_;
    my $order = $this->order;
    return !1 if $order != $that->order;
    my $le = $this->[_F_ELEMENTS];
    my $re = $that->[_F_ELEMENTS];
    my $delta0 = $re->[0] - $le->[0];
    if (!$delta0) {
        foreach my $x (1 .. $order) {
            return !1 if $re->[$x] != $le->[$x];
        }
        return !0;
    }
    my $modulus = $this->modulus;
    my $delta1  = $delta0 < 0? $delta0 + $modulus: $delta0 - $modulus;
    foreach my $x (1 .. $order) {
        my $delta = $re->[$x] - $le->[$x];
        return !1 if $delta != $delta0 && $delta != $delta1;
    }
    return !0;
}

# @e = $ds1->common_elements($ds2);
sub common_elements {
    my ($this, $that) = @_;
    my $order = $this->order;
    my @common = ();
    return @common if $order != $that->order;
    my $li = 0;
    my $ri = 0;
    my $lx = $this->[_F_INDEX_MIN];
    my $rx = $that->[_F_INDEX_MIN];
    my $le = $this->[_F_ELEMENTS];
    my $re = $that->[_F_ELEMENTS];
    my $lv = $le->[$lx];
    my $rv = $re->[$rx];
    while (1) {
        my $cmp = $lv <=> $rv;
        push @common, $lv if !$cmp;
        if ($cmp <= 0) {
            ++$lx <= $order or $lx = 0;
            last if ++$li > $order;
            $lv = $le->[$lx];
        }
        if ($cmp >= 0) {
            ++$rx <= $order or $rx = 0;
            last if ++$ri > $order;
            $rv = $re->[$rx];
        }
    }
    return @common;
}

# ($factor, $delta) = $ds1->find_linear_map($ds2);
sub find_linear_map {
    my ($this, $that) = @_;
    my $factor = _find_factor($this, $that);
    my $delta  = _delta_f($this, $factor, $that);
    return ($factor, $delta);
}

# @factor_delta_pairs = $ds1->find_all_linear_maps($ds2);
sub find_all_linear_maps {
    my ($this, $that) = @_;
    my $f1 = eval { _find_factor($this, $that) };
    return () if !defined $f1;
    my $modulus = $this->modulus;
    return
        sort { $a->[0] <=> $b->[0] }
        map {
            my $f = mulmod($f1, $_, $modulus);
            my $d = _delta_f($this, $f, $that);
            [$f, $d]
        } $this->multipliers;
}

1;
__END__