# Copyright 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-NumSeq.
#
# Math-NumSeq is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-NumSeq is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-NumSeq.  If not, see <http://www.gnu.org/licenses/>.

package Math::NumSeq::LongFractionPrimes;
use 5.004;
use strict;

use vars '$VERSION', '@ISA';
$VERSION = 71;
use Math::NumSeq;
@ISA = ('Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;

use Math::NumSeq::PrimeFactorCount;;
*_prime_factors = \&Math::NumSeq::PrimeFactorCount::_prime_factors;

use Math::NumSeq::Primes;
use Math::NumSeq::Squares;

# uncomment this to run the ### lines
#use Smart::Comments;


# use constant name => Math::NumSeq::__('...');
use constant description => Math::NumSeq::__('Primes where the fraction 1/p has digits repeating in full-length period p-1.  This is when the radix is a primitive root modulo p.');
use constant default_i_start => 1;
use constant characteristic_integer => 1;

use Math::NumSeq::Base::Digits
  'parameter_info_array';   # radix parameter

sub values_min {
  my ($self) = @_;
  unless (exists $self->{'values_min'}) {
    for (my $value = 2; ; $value++) {
      if ($self->pred($value)) {
        $self->{'values_min'} = $value;
        last;
      }
    }
  }
  return $self->{'values_min'};
}
sub values_max {
  my ($self) = @_;
  unless (exists $self->{'values_max'}) {
    $self->{'values_max'} = undef;
    if ($self->{'perfect_square'}) {
      for (my $value = $self->{'radix'}**2; $value >= 2; $value--) {
        if ($self->pred($value)) {
          $self->{'values_max'} = $value;
          last;
        }
      }
    }
  }
  return $self->{'values_max'};
}

#------------------------------------------------------------------------------

my @oeis_anum;
$oeis_anum[10] = 'A006883';
# OEIS-Catalogue: A006883

sub oeis_anum {
  my ($self) = @_;
  ### oeis_anum: $self
  return $oeis_anum[$self->{'radix'}];
}

#------------------------------------------------------------------------------

sub rewind {
  my ($self) = @_;
  $self->{'primes'} = Math::NumSeq::Primes->new;
  $self->{'i'} = $self->i_start;
  $self->{'perfect_square'} = Math::NumSeq::Squares->pred($self->{'radix'});
}

sub next {
  my ($self) = @_;
  ### LongFractionPrimes next() ...
  my $radix = $self->{'radix'};
  for (;;) {
    (undef, my $prime) = $self->{'primes'}->next;
    # FIXME: 3 ?
    if (_is_primitive_root ($self->{'radix'}, $prime) && $prime != 3) {
      return ($self->{'i'}++, $prime);
    }
    if ($self->{'perfect_square'} && $prime > $radix) {
      return;
    }
  }
}

sub pred {
  my ($self, $value) = @_;
  ### LongFractionPrimes pred(): "$value"
  return (Math::NumSeq::Primes->pred($value)
          && _is_primitive_root ($self->{'radix'}, $value));
}

sub _is_primitive_root {
  my ($base, $modulus) = @_;

  if (_is_infinite($modulus)) {
    return undef;
  }
  if ($modulus < 2) {
    return 0;
  }

  my $exponent = $modulus - 1;
  my ($good, @primes) = _prime_factors($exponent);
  return undef unless $good;

  my $prev_p = 0;
  while (defined (my $p = shift @primes)) {
    next if $p == $prev_p;
    $prev_p = $p;

    ### $p
    ### div: $exponent/$p
    ### powmod: _powmod($base, $exponent/$p, $modulus)

    if (_powmod($base, $exponent/$p, $modulus) <= 1) {
      return 0;
    }
  }

  # my $power = $base;
  # foreach (1 .. $value-2) {
  #   $power %= $value;
  #   if ($power == 1) {
  #     ### no, at: $_
  #     return 0;
  #   }
  #   $power *= $base;
  # }

  return 1;
}

sub _powmod {
  my ($base, $exponent, $modulus) = @_;

  my @exponent = _bit_split_hightolow($exponent)
    or return 1;

  my $power = $base;
  shift @exponent; # high 1 bit

  while (defined (my $bit = shift @exponent)) {  # high to low
    $power *= $power;
    $power %= $modulus;
    if ($bit) {
      $power *= $base;
      $power %= $modulus;
    }
  }
  return $power;
}

sub _bit_split_hightolow {
  my ($n) = @_;
  # ### _bit_split_hightolow(): "$n"

  if (ref $n) {
    if ($n->isa('Math::BigInt')
        && $n->can('as_bin')) {
      ### BigInt: $n->as_bin
      return split //, substr($n->as_bin,2);
    }
  }
  my @bits;
  while ($n) {
    push @bits, $n % 2;
    $n = int($n/2);
  }
  return reverse @bits;
}

1;
__END__

=for stopwords Ryde Math-NumSeq

=head1 NAME

Math::NumSeq::LongFractionPrimes -- primes for which fraction 1/p has a long period

=head1 SYNOPSIS

 use Math::NumSeq::LongFractionPrimes;
 my $seq = Math::NumSeq::LongFractionPrimes->new;
 my ($i, $value) = $seq->next;

=head1 DESCRIPTION

I<In progress ...>

This is the primes for which fraction 1/p written out in decimal has digits
repeating in period p-1.

    2, 3, 7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, ...
    starting i=1 for prime=2

For example 1/7=0.142857142857142857... is runs of 7-1=6 repeating digits
"142857", so 7 is in the sequence.  On the other hand 1/11=0.09090909... is
only 2 repeating digits, so is not in the sequence.

A prime p has full period p-1 digits when the base 10 is a primitive root
modulo p, meaning that 10 mod p, 100 mod p, 1000 mod p, ..., 10^(p-2) mod p,
are all != 1.

=head2 Radix

An optional C<radix> parameter selects a base other than decimal.

If C<radix> is a square, 4,9,16,etc then only primes 2 and 3 have period p-1
and the sequence stops after them.

=head1 FUNCTIONS

See L<Math::NumSeq/FUNCTIONS> for behaviour common to all sequence classes.

=over 4

=item C<$seq = Math::NumSeq::LongFractionPrimes-E<gt>new ()>

Create and return a new sequence object.

=item C<$value = $seq-E<gt>pred($value)>

Return true if C<$value> is a prime and fraction 1/p has digit period p-1.

=item C<$i = $seq-E<gt>i_start ()>

Return 1, the first term in the sequence being at i=1.

=back

=head1 SEE ALSO

L<Math::NumSeq>,
L<Math::NumSeq::Primes>

=head1 HOME PAGE

L<http://user42.tuxfamily.org/math-numseq/index.html>

=head1 LICENSE

Copyright 2012, 2013, 2014 Kevin Ryde

Math-NumSeq is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation; either version 3, or (at your option) any later
version.

Math-NumSeq is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
more details.

You should have received a copy of the GNU General Public License along with
Math-NumSeq.  If not, see <http://www.gnu.org/licenses/>.

=cut