File Coverage

blib/lib/Math/DifferenceSet/Planar.pm

Criterion	Covered	Total	%
statement	835	835	100.0
branch	288	296	97.3
condition	60	69	86.9
subroutine	133	133	100.0
pod	79	79	100.0
total	1395	1412	98.8

package Math::DifferenceSet::Planar;

use strict;
use warnings;
use Carp qw(croak);
use Math::DifferenceSet::Planar::Data;
use Math::Prime::Util qw(
    is_prime is_power is_prime_power euler_phi factor_exp gcd
    mulmod addmod invmod powmod divmod logint
);

# 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_ZETA        =>  4;  # "zeta" value
use constant _F_ETA         =>  5;  # "eta" value, initially undef
use constant _F_THETA       =>  6;  # translation amount from canonical set
use constant _F_PRINC_ELEMS =>  7;  # principal elements arrayref
use constant _F_SUPPL_ELEMS =>  8;  # supplemental elements arrayref
use constant _F_LAMBDA      =>  9;  # plane logarithm value or undef
use constant _F_ELEMENTS    => 10;  # elements arrayref, initially undef
use constant _F_X_START     => 11;  # index of start element in elements
use constant _F_X_GAP       => 12;  # index of max gap element in elements
use constant _NFIELDS       => 13;

# usable native integer bits, typically 63 or 31
use constant _NATIVE_BITS     => logint(~0, 2);
# max order safe to use with native integer arithmetic
use constant _MAX_SMALL_ORDER => int( sqrt(2)*((1<<(_NATIVE_BITS>>1))-0.5) );

*canonize = \&lex_canonize;

our $VERSION = '1.003';

our $_MAX_ENUM_COUNT  = 32768;          # limit for stored rotator set size
our $_MAX_MEMO_COUNT  = 4096;           # limit for memoized values
our $_USE_SPACES_DB   = 1;              # enable looking up rotators

my $current_data = undef;               # current M::D::P::Data object

my %memo_n_planes  = ();                # memoized n_planes values
my @memo_np_orders = ();                # memoized orders FIFO
my %memo_rotators  = ();                # memoized rotators arrayrefs
my @memo_ro_orders = ();                # memoized orders FIFO

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

# 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  |
#             +-------+-------+--   --+-------+
#
# Otherwise, return undef.
#
sub _rotators {
    my ($class, $order, $base, $exponent, $modulus) = @_;
    return $memo_rotators{$order} if exists $memo_rotators{$order};
    my $rotators = undef;
    my $space = $_USE_SPACES_DB && $class->_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 = [$radices, $depths, $inverses];
    }
    elsif (_n_planes($order, $exponent, $modulus) <= $_MAX_ENUM_COUNT) {
        my $mult = 3 * $exponent;
        my $sieve = '1' x $modulus;
        substr($sieve, 1, 1) = '0';
        my $e = 1;
        for (2 .. $mult) {
            $e = mulmod($e, $base, $modulus);
            substr($sieve, $e, 1) = '0';
        }
        my @rot = (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, $x, 1) = '0';
                my $e = $x;
                for (2 .. $mult) {
                    $e = mulmod($e, $base, $modulus);
                    substr($sieve, $e, 1) = '0';
                }
                push @rot, $x;
                last if $mult == @rot;
            }
        }
        $rotators = \@rot;
    }
    $memo_rotators{$order} = $rotators;
    push @memo_ro_orders, $order;
    delete $memo_rotators{shift @memo_ro_orders}
        while $_MAX_MEMO_COUNT < @memo_ro_orders;
    return $rotators;
}

# iterative rotator base generator, slow, but memory efficient
sub _sequential_rotators {
    my ($order, $base, $exponent, $modulus) = @_;
    my $n_planes = _n_planes($order, $exponent, $modulus);
    my $mult = 3 * $exponent;
    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);
            }
            my $e = $x;
            for (2 .. $mult) {
                $e = mulmod($e, $base, $modulus);
                next ELEMENT if $e < $x;
            }
            ++$mx;
            return $x;
        }
    };
}

# structured rotator base iterator, time and space efficient
sub _structured_rotators {
    my ($modulus, $rotators) = @_;
    my ($radices, $depths, $inverses) = @{$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;
    };
}

# format a planar difference set space description, like '7^3 [2^6 5^2]'
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 0 if logint($base, 2) * $exponent >= _NATIVE_BITS;
    my $power = 1;
    while ($exponent) {
        $power *= $base if 1 & $exponent;
        $exponent >>= 1 and $base *= $base;
    }
    return $power;
}

# check whether order is small enough and calculate modulus
sub _modulus {
    my ($order) = @_;
    if ($order > _MAX_SMALL_ORDER) {
        croak "order $order too large for this platform\n";
    }
    return +($order + 1) * $order + 1;
}

# calculate minimal equivalent of a factor
sub _min_factor {
    my ($factor, $base, $modulus) = @_;
    my $f0 = $factor;
    my $f  = mulmod($f0, $base, $modulus);
    while ($f != $f0) {
        $factor = $f if $f < $factor;
        $f = mulmod($f, $base, $modulus);
    }
    return $factor;
}

# boolean whether an ordered strictly increasing list contains an element
sub _ol_contains {
    my ($haystack, $needle) = @_;
    my $lx = 0;
    my $hx = $#{$haystack};
    while ($lx <= $hx) {
        my $x = ($lx + $hx) >> 1;
        my $v = $haystack->[$x];
        my $cmp = $needle <=> $v;
        return !0 if !$cmp;
        if ($cmp < 0) {
            $hx = $x - 1;
        }
        else {
            $lx = $x + 1;
        }
    }
    return !1;
}

# calculate number of planes with memoization
sub _n_planes {
    my ($order, $exponent, $modulus) = @_;
    return $memo_n_planes{$order} if exists $memo_n_planes{$order};
    my $n_planes = $memo_n_planes{$order} =
        euler_phi($modulus) / (3 *$exponent);
    push @memo_np_orders, $order;
    delete $memo_n_planes{shift @memo_np_orders}
        while $_MAX_MEMO_COUNT < @memo_np_orders;
    return $n_planes;
}

# arrange elements of a planar difference set in numerical order
# return arrayref of sorted list, start index, gap index
sub _sort_elements {
    my $modulus  = shift;
    my @elements = sort { $a <=> $b } @_;
    my $xs = my $xl = $#elements;
    my $xg = my $xh = 0;
    my $lo = $elements[$xl] - $modulus;
    my $hi = $elements[0];
    my $mg = my $ml = my $d = $hi - $lo;
    while ($xh < $#elements) {
        $xl = $xh++;
        ($lo, $hi) = ($hi, $elements[$xh]);
        $d = $hi - $lo;
        $ml = $d, $xs = $xl if $d < $ml;
        $mg = $d, $xg = $xh if $d > $mg;
    }
    croak "duplicate element: $elements[$xs]" if !$ml;
    croak 'delta 1 elements missing'          if $ml > 1;
    return (\@elements, $xs, $xg);
}

# generate fill elements (0, 1 or 2 values)
sub _fill_elements {
    my ($order, $modulus) = @_;
    my $type = $order % 3;
    return ()  if $type == 2;
    return (0) if $type == 0;
    my $m3 = $modulus / 3;
    return ($m3, $m3 << 1);
}

# recall elements or generate them, updating object attributes
sub _elements {
    my ($this) = @_;
    my $elements = $this->[_F_ELEMENTS];
    return $elements if $elements;
    my (           $order,   $base,   $exponent,   $modulus,   $theta) =
        @{$this}[_F_ORDER, _F_BASE, _F_EXPONENT, _F_MODULUS, _F_THETA];
    my $mult = 3 * $exponent;
    my @elem = ();
    foreach my $e0 (@{$this->[_F_PRINC_ELEMS]}) {
        my $e = $e0;
        push @elem, addmod($e0, $theta, $modulus);
        for (2 .. $mult) {
            $e = mulmod($e, $base, $modulus);
            push @elem, addmod($e, $theta, $modulus);
        }
    }
    foreach my $e0 (@{$this->[_F_SUPPL_ELEMS]}) {
        push @elem, addmod($e0, $theta, $modulus);
        my $e = mulmod($e0, $base, $modulus);
        while ($e != $e0) {
           push @elem, addmod($e, $theta, $modulus);
           $e = mulmod($e, $base, $modulus);
        }
    }
    foreach my $e0 (_fill_elements($order, $modulus)) {
        push @elem, addmod($e0, $theta, $modulus);
    }
    ($elements) = @{$this}[_F_ELEMENTS, _F_X_START, _F_X_GAP] =
        _sort_elements($modulus, @elem);
    return $elements;
}

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

# identify a plane by its rotator value with respect to a given plane
# (setting its lambda value if possible)
sub _log {
    my ($this, $ref) = @_;
    my $modulus    =  $this->[_F_MODULUS];
    my $delta      = -$this->[_F_THETA];
    my $ref_lambda =   $ref->[_F_LAMBDA];
    my $factor = 0;
    my %this_e = ();
    my $this_elements = $this->_elements;
    foreach my $e (@{$this_elements}) {
        $this_e{$delta? addmod($e, $delta, $modulus): $e} = 1;
    }
    if (@{$ref->[_F_PRINC_ELEMS]}) {
        my $inv_r = invmod($ref->[_F_PRINC_ELEMS]->[0], $modulus);
        ELEM:
        foreach my $o (@{$this->[_F_PRINC_ELEMS]}) {
            my $ro = mulmod($inv_r, $o, $modulus);
            foreach my $e (@{$ref->[_F_PRINC_ELEMS]}) {
                next ELEM if !exists $this_e{ mulmod($e, $ro, $modulus) };
            }
            $factor = $ro;
            last;
        }
    }
    else {
        my $ri = $this->iterate_rotators;
        ROT:
        while (my $ro = $ri->()) {
            foreach my $e (@{$ref->[_F_SUPPL_ELEMS]}) {
                next ROT if !exists $this_e{ mulmod($e, $ro, $modulus) };
            }
            $factor = $ro;
            last;
        }
    }
    croak 'unaligned sets' if !$factor;
    my $base = $this->[_F_BASE];
    if ($ref_lambda) {
        $this->[_F_LAMBDA] =
            _min_factor(
                mulmod($ref_lambda, $factor, $modulus), $base, $modulus
            );
    }
    return $factor;
}

# $factor = _find_factor($ds1, $ds2);
sub _find_factor {
    my ($this, $that) = @_;
    return 1 if $this == $that;
    my $order = $this->order;
    croak 'sets of same size expected' if $order != $that->order;
    my $log_this = $this->[_F_LAMBDA];
    my $log_that = $that->[_F_LAMBDA];
    if (!$log_that) {
        $log_that = _log($that, $this);
    }
    elsif (!$log_this) {
        $log_this = _log($this, $that);
    }
    return $log_this? divmod($log_that, $log_this, $this->modulus): $log_that;
}

# 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 $elements = $that->_elements;
    my $s        = $elements->[$that->[_F_X_START]];
    return addmod($s, -mulmod($x, $factor, $modulus), $modulus);
}

# $bool = _is_mult($factor, $base, $mult, $modulus);
sub _is_mult {
    my ($factor, $base, $mult, $modulus) = @_;
    return !0 if $factor == $base;
    my $p = $base;
    for (2 .. $mult-1) {
        $p = mulmod($p, $base, $modulus);
        return !0 if $factor == $p;
    }
    return !1;
}

# ($order, $base, $exponent, $key) = _order_from_params(@_);
sub _order_from_params {
    my ($order, $exponent) = @_;
    my $base = undef;
    if (defined $exponent) {
        $base  = $order;
        croak "order base $base is not a prime" if !is_prime($base);
        $order = _pow($base, $exponent);
        croak "order $base ** $exponent too large for this platform"
            if !$order || $order > _MAX_SMALL_ORDER;
    }
    else {
        croak "order $order too large for this platform"
            if $order > _MAX_SMALL_ORDER;
        $exponent = is_prime_power($order, \$base);
        croak "order $order is not a prime power" if !$exponent;
    }
    return ($order, $base, $exponent);
}

# ($this, $order, $base, $exponent, $modulus) = _full_params(@_);
sub _full_params {
    my $this = shift;
    my ($order, $base, $exponent, $modulus);
    if (@_) {
        ($order, $base, $exponent) = _order_from_params(@_);
        $modulus = _modulus($order);
    }
    else {
        croak 'parameters expected if called as a class method' if !ref $this;
        (              $order,   $base,   $exponent,   $modulus) =
            @{$this}[_F_ORDER, _F_BASE, _F_EXPONENT, _F_MODULUS];
    }
    return ($this, $order, $base, $exponent, $modulus);
}

# $bool = $class->_known_ref('ref_std', $base, $exponent);
sub _known_ref {
    my ($class, $attribute, $base, $exponent) = @_;
    my $order = defined($exponent)? _pow($base, $exponent): $base;
    return !1 if !$order || $order > $class->_data->max_order;
    my $pds   = $class->_data->get($order, 'base', $attribute);
    return
        $pds && (!defined($exponent) || $base == $pds->base) &&
        $pds->$attribute != 0;
}

# identity transformation
sub _no_change { $_[0] }

# $it = $class->_iterate_refs('ref_std', 'zeta_canonize', 10, 20);
sub _iterate_refs {
    my ($class, $attribute, $transform, @minmax) = @_;
    my $dit = $class->_data->iterate_refs($attribute, @minmax);
    return sub {
        my $pds = $dit->();
        return undef if !$pds;
        my $this = eval { $class->_from_pds($pds->order, $pds) };
        return $this && $this->multiply($pds->$attribute)->$transform;
    };
}

# $ds = Math::DifferenceSet::Planar->_from_pds($order, $pds);
sub _from_pds {
    my ($class, $order, $pds) = @_;
    my $modulus  = _modulus($order);
    my $base     = $pds->base;
    my $exponent = $base == $order? 1: logint($order, $base);
    my $main     = $pds->main_elements;
    my (@princ, @suppl) = ();
    foreach my $e (@{$main}) {
        if (gcd($modulus, $e) == 1) {
            push @princ, $e;
        }
        else {
            push @suppl, $e;
        }
    }
    my $lambda = undef;
    if (my $log = $pds->ref_std) {
        $lambda =
            $log == 1? 1: _min_factor(invmod($log, $modulus), $base, $modulus);
    }

    return bless [
        $order,
        $base,
        $exponent,
        $modulus,
        0,                      # zeta
        0,                      # eta
        0,                      # theta
        \@princ,
        \@suppl,
        $lambda,
        undef,                  # elements
        undef,                  # index_start
        undef,                  # index_gap
    ], $class;
}

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

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

sub set_database {
    my $class = shift;
    $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 !1 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 = shift;
    my ($order, $base, $exponent) = _order_from_params(@_);
    my $pds = $class->_data->get($order);
    if (!$pds) {
        my $key = join q[, ], @_;
        croak "PDS($key) not available";
    }
    return $class->_from_pds($order, $pds);
}

sub lex_reference {
    my $class = shift;
    my ($order, $base, $exponent) = _order_from_params(@_);
    my $pds = $class->_data->get($order);
    if ($pds && (my $lambda = $pds->ref_lex)) {
        return $class->_from_pds($order, $pds)->multiply($lambda)->lex_canonize;
    }
    return undef;
}

sub gap_reference {
    my $class = shift;
    my ($order, $base, $exponent) = _order_from_params(@_);
    my $pds = $class->_data->get($order);
    if ($pds && (my $lambda = $pds->ref_gap)) {
        return
            $class->_from_pds($order, $pds)->multiply($lambda)->gap_canonize;
    }
    return undef;
}

sub std_reference {
    my $class = shift;
    my ($order, $base, $exponent) = _order_from_params(@_);
    my $pds = $class->_data->get($order);
    if ($pds && (my $lambda = $pds->ref_std)) {
        return
            $class->_from_pds($order, $pds)->zeta_canonize->multiply($lambda);
    }
    return undef;
}

# $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  = _modulus($order);
    if (grep { $_ < 0 || $modulus <= $_ } @_) {
        my $max = $modulus - 1;
        croak "element values inside range 0..$max expected";
    }
    my ($elements, $index_start, $index_gap) = _sort_elements($modulus, @_);
    my $n_mult = 3 * $exponent;

    # find zeta and theta
    my ($lx, $ux, $c) = (0, 0, 0);
    my ($e0, $e2, $e3) = @{$elements}[$index_start, 0, 0];
    my $is_m3 = $order % 3 == 1;
    my $de    = $is_m3? $modulus/3: $order + 1;
    my $bogus = 0;
    while ($c != $de) {
        if ($c < $de) {
            $ux = 0 if ++$ux > $order;
            $e3 = $elements->[$ux];
            $bogus = 1, last if $ux == $lx;
        }
        else {
            $bogus = 1, last if ++$lx > $order;
            $e2 = $elements->[$lx];
        }
        $c = $e3 < $e2? $modulus + $e3 - $e2: $e3 - $e2;
    }
    croak "delta $de elements missing\n" if $bogus;
    my ($zeta, $theta);
    if ($is_m3) {
        $theta = addmod($e3, $de, $modulus);
        $zeta  = mulmod($theta, $order - 1, $modulus);
    }
    else {
        $zeta  = addmod(mulmod($e3, $order, $modulus), -$e0, $modulus);
        $theta = $zeta && divmod($zeta, $order - 1, $modulus);
    }

    my @elems =
        sort { $a <=> $b } map { addmod($_, -$theta, $modulus) } @{$elements};
    my @princ = ();
    my @suppl = ();
    my %todo  = map {($_ => 1)} @elems;
    foreach my $start (@elems) {
        next if !exists $todo{$start};
        delete $todo{$start};
        my $this  = mulmod($start, $base, $modulus);
        my $count = 1;
        while ($this != $start) {
            if (!defined delete $todo{$this}) {
                croak
                    "bogus set: prime divisor $base of order $order " .
                    'is not a multiplier';
            }
            ++$count;
            $this = mulmod($this, $base, $modulus);
        }
        if ($count == $n_mult) {
            if (gcd($start, $modulus) == 1) {
                push @princ, $start;
            }
            else {
                push @suppl, $start;
            }
        }
        elsif ($count >= 3) {
            push @suppl, $start;
        }
    }
    my $eta = $exponent == 1? $zeta: undef;

    return bless [
        $order,
        $base,
        $exponent,
        $modulus,
        $zeta,
        $eta,
        $theta,
        \@princ,
        \@suppl,
        undef,                  # lambda
        $elements,
        $index_start,
        $index_gap,
    ], $class;
}

# $ds = Math::DifferenceSet::Planar->from_elements(
#   0, 1, 3, 9, 27, 49, 56, 61, 77, 81
# );
sub from_elements {
    my $class = shift;
    my $this  = $class->from_elements_fast(@_);
    if(my $ref = eval { $class->new(@_ - 1) }) {
        eval {
            my ($factor, $delta) = $ref->find_linear_map($this);
            !$ref->multiply($factor)->translate($delta)->compare($this)
        } or
        croak 'apparently not a planar difference set';
    }
    return $this;
}

# $ds = Math::DifferenceSet::Planar->from_lambda($order, $lambda);
# $ds = Math::DifferenceSet::Planar->from_lambda($order, $lambda, $theta);
sub from_lambda {
    my ($class, $order, $lambda, $theta) = @_;
    my ($base, $exponent);
    $exponent = is_prime_power($order, \$base);
    croak "this implementation cannot handle order $order" if !$exponent;
    my $modulus = _modulus($order);
    croak "impossible lambda value $lambda" if gcd($modulus, $lambda) != 1;
    my $l = mulmod($lambda, $base, $modulus);
    while ($l != $lambda) {
        croak "non-canonical lambda value $lambda" if $l < $lambda;
        $l = mulmod($l, $base, $modulus);
    }
    croak "non-canonical theta value $theta"
        if $theta && ($theta < 0 || $modulus <= $theta);
    my $ref = $class->std_reference($order);
    croak "reference set of order $order not available" if !$ref;
    my $this = $ref->multiply($lambda);
    return $theta? $this->translate($theta): $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 = _modulus($order);
    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//;
}

# $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;
        my $this = $class->_from_pds($pds->order, $pds);
        return $this;
    };
}

sub available_min_order { $_[0]->_data->min_order }
sub available_max_order { $_[0]->_data->max_order }
sub available_count     { $_[0]->_data->count     }

sub known_std_ref {
    my $class = shift;
    return $class->_known_ref('ref_std', @_);
}

sub known_lex_ref {
    my $class = shift;
    return $class->_known_ref('ref_lex', @_);
}

sub known_gap_ref {
    my $class = shift;
    return $class->_known_ref('ref_gap', @_);
}

sub iterate_known_std_refs {
    my ($class, @minmax) = @_;
    return $class->_iterate_refs('ref_std', '_no_change', @minmax);
}

sub iterate_known_lex_refs {
    my ($class, @minmax) = @_;
    return $class->_iterate_refs('ref_lex', 'lex_canonize', @minmax);
}

sub iterate_known_gap_refs {
    my ($class, @minmax) = @_;
    return $class->_iterate_refs('ref_gap', 'gap_canonize', @minmax);
}

sub known_std_ref_min_order { $_[0]->_data->ref_min_order('ref_std') }
sub known_std_ref_max_order { $_[0]->_data->ref_max_order('ref_std') }
sub known_std_ref_count     { $_[0]->_data->ref_count(    'ref_std') }
sub known_lex_ref_min_order { $_[0]->_data->ref_min_order('ref_lex') }
sub known_lex_ref_max_order { $_[0]->_data->ref_max_order('ref_lex') }
sub known_lex_ref_count     { $_[0]->_data->ref_count(    'ref_lex') }
sub known_gap_ref_min_order { $_[0]->_data->ref_min_order('ref_gap') }
sub known_gap_ref_max_order { $_[0]->_data->ref_max_order('ref_gap') }
sub known_gap_ref_count     { $_[0]->_data->ref_count(    'ref_gap') }

# $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;
# $z = $ds->zeta;
# $t = $ds->theta;
# $l = $ds->lambda;
sub order          { $_[0]->[_F_ORDER   ] }
sub order_base     { $_[0]->[_F_BASE    ] }
sub order_exponent { $_[0]->[_F_EXPONENT] }
sub modulus        { $_[0]->[_F_MODULUS ] }
sub zeta           { $_[0]->[_F_ZETA    ] }
sub theta          { $_[0]->[_F_THETA   ] }
sub lambda         { $_[0]->[_F_LAMBDA  ] }

sub min_element {
    my ($this) = @_;
    my $elements = $this->_elements;
    return $elements->[0];
}

sub max_element {
    my ($this) = @_;
    my $elements = $this->_elements;
    return $elements->[-1];
}

sub n_planes {
    my ($this, $order, $base, $exponent, $modulus) = _full_params(@_);
    return _n_planes($order, $exponent, $modulus)
}

# @e  = $ds->elements;
sub elements {
    my ($this) = @_;
    my $elements = $this->_elements;
    return 0+@{$elements} if !wantarray;
    my $x_start = $this->[_F_X_START];
    return @{$elements}[$x_start .. $#{$elements}, 0 .. $x_start-1];
}

# $e0 = $ds->element(0);
sub element {
    my ($this, $index) = @_;
    my $n_elems  = $this->[_F_ORDER] + 1;
    return undef if $index < -$n_elems || $n_elems <= $index;
    my $elements = $this->_elements;
    my $x_eff    = $this->[_F_X_START] + $index;
    $x_eff -= $n_elems if $x_eff >= $n_elems;
    return $elements->[$x_eff];
}

# @e  = $ds->elements_sorted;
sub elements_sorted {
    my ($this) = @_;
    my $elements = $this->_elements;
    return @{$elements};
}

# $e0 = $ds->element_sorted(0);
sub element_sorted {
    my ($this, $index) = @_;
    my $elements = $this->_elements;
    return $elements->[$index];
}

# $ds1 = $ds->translate(1);
sub translate {
    my ($this, $delta) = @_;
    my $modulus = $this->[_F_MODULUS];
    $delta %= $modulus;
    return $this if !$delta;
    my $that = bless [@{$this}], ref $this;
    my $elements = $this->[_F_ELEMENTS];
    if ($elements) {
        my $lim = $modulus - $delta;
        my @elems = my @wrap = ();
        foreach my $e (@{$elements}) {
            if ($e < $lim) {
                push @wrap, $e + $delta;
            }
            else {
                push @elems, $e - $lim;
            }
        }
        my $dx = @elems;
        push @elems, @wrap;
        my $x_start = addmod($this->[_F_X_START], $dx, 0+@elems);
        my $x_gap   = addmod($this->[_F_X_GAP],   $dx, 0+@elems);
        $that->[_F_ELEMENTS] = \@elems;
        $that->[_F_X_START]  = $x_start;
        $that->[_F_X_GAP]    = $x_gap;
    }
    my ($order, $zeta, $theta) = @{$this}[_F_ORDER, _F_ZETA, _F_THETA];
    $that->[_F_ZETA] =
        addmod($zeta, mulmod($delta, $order-1, $modulus), $modulus);
    $that->[_F_THETA] = addmod($theta, $delta, $modulus);
    if (defined (my $e = $that->[_F_ETA])) {
        my $base = $this->[_F_BASE];
        $that->[_F_ETA] =
            addmod($e, mulmod($delta, $base-1, $modulus), $modulus);
    }
    return $that;
}

# $ds2 = $ds->canonize;
# $ds2 = $ds->lex_canonize;
sub lex_canonize {
    my ($this) = @_;
    my $elements = $this->_elements;
    return $this->translate(- $elements->[$this->[_F_X_START]]);
}

# $ds2 = $ds->gap_canonize;
sub gap_canonize {
    my ($this) = @_;
    my $elements = $this->_elements;
    return $this->translate(- $elements->[$this->[_F_X_GAP]]);
}

# $ds2 = $ds->zeta_canonize;
sub zeta_canonize {
    my ($this) = @_;
    return $this->translate(- $this->[_F_THETA]);
}

# $it  = $ds->iterate_rotators;
# while (my $m = $it->()) {
#   ...
# }
sub iterate_rotators {
    my ($this, $order, $base, $exponent, $modulus) = _full_params(@_);
    my $rotators = $this->_rotators($order, $base, $exponent, $modulus);
    return _sequential_rotators($order, $base, $exponent, $modulus)
        if !$rotators;
    return _structured_rotators($modulus, $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)->lex_canonize: undef;
    };
}

sub iterate_planes_zc {
    my ($this) = @_;
    my $ref  = $this->zeta_canonize;
    my $r_it = $ref->iterate_rotators;
    return sub {
        my $r = $r_it->();
        return $r? $ref->multiply($r): undef;
    };
}

sub iterate_principal_planes {
    my ($this) = @_;
    my $ref = $this->zeta_canonize;
    my (      $princ_elems,   $modulus) =
    @{$ref}[_F_PRINC_ELEMS, _F_MODULUS];
    my $pex = 0;
    return sub {
        return undef if $pex >= @{$princ_elems};
        my $r = invmod($princ_elems->[$pex++], $modulus);
        return $ref->multiply($r)->lex_canonize;
    };
}

sub iterate_principal_planes_zc {
    my ($this) = @_;
    my $ref = $this->zeta_canonize;
    my (      $princ_elems,   $modulus) =
    @{$ref}[_F_PRINC_ELEMS, _F_MODULUS];
    my $pex = 0;
    return sub {
        return undef if $pex >= @{$princ_elems};
        my $r = invmod($princ_elems->[$pex++], $modulus);
        return $ref->multiply($r);
    };
}

# @pm = $ds->multipliers;
sub multipliers {
    my ($this, $order, $base, $exponent, $modulus) = _full_params(@_);
    my $n_mult = 3 * $exponent;
    return $n_mult if !wantarray;
    my @mult = (1, $base);
    my $x = $base;
    while (@mult < $n_mult) {
        $x = mulmod($x, $base, $modulus);
        push @mult, $x;
    }
    return @mult;
}

# $ds3 = $ds->multiply($m);
sub multiply {
    my ($this, $factor) = @_;
    my (           $order,   $base,   $exponent,   $modulus) =
        @{$this}[_F_ORDER, _F_BASE, _F_EXPONENT, _F_MODULUS];
    $factor %= $modulus;
    return $this if 1 == $factor;
    croak "factor $factor is not coprime to modulus"
        if gcd($modulus, $factor) != 1;
    my $theta1 = $this->[_F_THETA];
    my $theta  = $theta1 && mulmod($theta1, $factor, $modulus);
    my $mult = 3 * $exponent;
    return $this->translate($theta - $theta1) if
        _is_mult($factor, $base, $mult, $modulus);
    my (           $zeta,   $eta,   $lambda) =
        @{$this}[_F_ZETA, _F_ETA, _F_LAMBDA];
    $zeta &&= mulmod($zeta, $factor, $modulus);
    $eta  &&= mulmod($eta,  $factor, $modulus);
    my @princ = ();
    my @suppl = ();
    foreach my $e (@{$this->[_F_PRINC_ELEMS]}) {
        my $p0 = my $p = mulmod($e, $factor, $modulus);
        for (2 .. $mult) {
            $p = mulmod($p, $base, $modulus);
            $p0 = $p if $p0 > $p;
        }
        push @princ, $p0;
    }
    @princ = sort { $a <=> $b } @princ;
    foreach my $e (@{$this->[_F_SUPPL_ELEMS]}) {
        push @suppl,
            _min_factor(mulmod($e, $factor, $modulus), $base, $modulus);
    }
    @suppl = sort { $a <=> $b } @suppl;
    if ($lambda) {
        $lambda =
            _min_factor(mulmod($lambda, $factor, $modulus), $base, $modulus);
    }
    return bless [
        $order,
        $base,
        $exponent,
        $modulus,
        $zeta,
        $eta,
        $theta,
        \@princ,
        \@suppl,
        $lambda,
        undef,                          # elements
        undef,                          # index_start
        undef,                          # index_gap
    ], ref $this;
}

# ($e1, $e2) = $ds->find_delta($delta);
sub find_delta {
    my ($this, $delta) = @_;
    my $order    = $this->order;
    my $modulus  = $this->modulus;
    my $elements = $this->_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);
}

# $e0 = $ds->start_element;
sub start_element {
    my ($this) = @_;
    my $elements = $this->_elements;
    return $elements->[$this->[_F_X_START]];
}

# ($e1, $e2) = $ds->peak_elements
sub peak_elements {
    my ($this) = @_;
    return $this->find_delta($this->modulus >> 1);
}

# ($e1, $e2, $delta) = $ds->largest_gap;
sub largest_gap {
    my ($this) = @_;
    my $elements = $this->_elements;
    my $x2 = $this->[_F_X_GAP];
    my ($e1, $e2) = @{$elements}[$x2 - 1, $x2];
    my $delta = $x2? $e2 - $e1: $this->modulus + $e2 - $e1;
    return ($e1, $e2, $delta);
}

# $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 ($x) = $this->find_delta( invmod($p, $m) );
        my $s   = $this->[_F_ELEMENTS]->[$this->[_F_X_START]];
        $eta = $this->[_F_ETA] = addmod(mulmod($x, $p, $m), -$s, $m);
    }
    return $eta;
}

sub plane_principal_elements    { @{$_[0]->[_F_PRINC_ELEMS]} }
sub plane_supplemental_elements { @{$_[0]->[_F_SUPPL_ELEMS]} }

sub plane_fill_elements {
    my ($this) = @_;
    my @fe = _fill_elements(@{$this}[_F_ORDER, _F_MODULUS]);
    return @fe;
}

sub plane_derived_elements_of {
    my ($this, @elem) = @_;
    my ($base, $modulus) = @{$this}[_F_BASE, _F_MODULUS];
    my @de = ();
    foreach my $e0 (@elem) {
        my $e = mulmod($e0, $base, $modulus);
        while ($e != $e0) {
            push @de, $e;
            $e = mulmod($e, $base, $modulus);
        }
    }
    return @de;
}

sub plane_unit_elements {
    my ($this) = @_;
    if (!wantarray) {
        return 3 * $this->[_F_EXPONENT] * @{$this->[_F_PRINC_ELEMS]};
    }
    return
        map { ($_, $this->plane_derived_elements_of($_)) }
        @{$this->[_F_PRINC_ELEMS]};
}

sub plane_nonunit_elements {
    my ($this) = @_;
    if (!wantarray) {
        return $this->[_F_ORDER] + 1 - $this->plane_unit_elements;
    }
    return
        (
            (
                map { @{$_} }
                sort { @{$b} <=> @{$a} || $a->[0] <=> $b->[0] }
                map { [$_, $this->plane_derived_elements_of($_)] }
                @{$this->[_F_SUPPL_ELEMS]}
            ),
            _fill_elements(@{$this}[_F_ORDER, _F_MODULUS])
        );
}

# $bool = $ds->contains($e)
sub contains {
    my ($this, $elem) = @_;
    my $elements  = $this->_elements;
    return _ol_contains($elements, $elem);
}

# $cmp = $ds1->compare($ds2);
sub compare {
    my ($this, $that) = @_;
    my $order = $this->order;
    my $cmp   = $order <=> $that->order;
    return $cmp if $cmp;
    my $le = $this->_elements;
    my $re = $that->_elements;
    foreach my $x (0 .. $order) {
        $cmp = $le->[$x] <=> $re->[$x];
        return $cmp if $cmp;
    }
    return 0;
}

# $cmp = $ds1->compare_topdown($ds2);
sub compare_topdown {
    my ($this, $that) = @_;
    my $order = $this->order;
    my $cmp   = $order <=> $that->order;
    return $cmp if $cmp;
    my $le = $this->_elements;
    my $re = $that->_elements;
    my $x = $order;
    while ($x >= 0) {
        $cmp = $le->[$x] <=> $re->[$x];
        return $cmp if $cmp;
        --$x;
    }
    return 0;
}

# $bool = $ds1->same_plane($ds2);
sub same_plane {
    my ($this, $that) = @_;
    my $order = $this->order;
    return !1 if $order != $that->order;
    my $l1 = $this->[_F_LAMBDA];
    my $l2 = $that->[_F_LAMBDA];
    return $l1 == $l2 if $l1 && $l2;
    my $le = $this->[_F_PRINC_ELEMS];
    my $re;
    if (@{$le}) {
        $re = $that->[_F_PRINC_ELEMS];
    }
    else {
        $le = $this->[_F_SUPPL_ELEMS];
        $re = $that->[_F_SUPPL_ELEMS];
    }
    foreach my $x (0 .. $#{$le}) {
        return !1 if $re->[$x] != $le->[$x];
    }
    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 $le = $this->_elements;
    my $re = $that->_elements;
    my $lv = $le->[0];
    my $rv = $re->[0];
    while (1) {
        my $cmp = $lv <=> $rv;
        push @common, $lv if !$cmp;
        if ($cmp <= 0) {
            last if ++$li > $order;
            $lv = $le->[$li];
        }
        if ($cmp >= 0) {
            last if ++$ri > $order;
            $rv = $re->[$ri];
        }
    }
    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__