File Coverage

blib/lib/Math/NumSeq/PythagoreanHypots.pm

Criterion	Covered	Total	%
statement	84	86	97.6
branch	20	22	90.9
condition	8	9	88.8
subroutine	17	17	100.0
pod	4	4	100.0
total	133	138	96.3

# Copyright 2010, 2011, 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 .

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

use vars '$VERSION', '@ISA';
$VERSION = 72;

use Math::NumSeq 7; # v.7 for _is_infinite()
@ISA = ('Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;

use Math::NumSeq::Primes;
use Math::Prime::XS 0.23 'is_prime'; # version 0.23 fix for 1928099

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


# use constant name => Math::NumSeq::__('Pyathagorean Hypotenuses');
use constant description => Math::NumSeq::__('The hypotenuses of Pythagorean triples, ie. integers C for which there\'s some A>=1,B>=1 satisfying A^2+B^2=C^2.  Primitive hypotenuses are where A,B have no common factor.');
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;
use constant values_min => 5;
use constant i_start => 1;

use constant parameter_info_array =>
  [ {
     name        => 'pythagorean_type',
     type        => 'enum',
     default     => 'all',
     choices     => ['all','primitive'],
     choices_display => [Math::NumSeq::__('All'),
                         Math::NumSeq::__('Primitive'),
                        ],
    },
  ];

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

# cf A002144 - primes 4n+1, the primitive elements of hypots x!=y
#              -1 is a quadratic residue ...
#    A002365 - the "y" of prime "c" ??
#    A002366 - the "x" of prime "c" ??
#    A046083 - the "a" smaller number, ordered by "c"
#    A046084 - the "b" second number, ordered by "c"
#
#    A008846 - primitives, x,y no common factor
#    A004613 - all prime factors are 4n+1, is 1 then primitive hypots
#
#    A009000 - hypots with repetitions
#    A009012 - "b" second number, ordered by "b", with repetitions
#
my %oeis_anum = (all       => 'A009003',  # distinct a!=b and a,b>0
                 primitive => 'A008846',
                );
# OEIS-Catalogue: A009003
# OEIS-Catalogue: A008846 pythagorean_type=primitive
sub oeis_anum {
  my ($self) = @_;
  return $oeis_anum{$self->{'pythagorean_type'}};
}

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

sub rewind {
  my ($self) = @_;
  $self->{'array'} = [];
  $self->{'i'} = $self->i_start;
  $self->{'hi'} = 1;
}

sub next {
  my ($self) = @_;
  my $array = $self->{'array'};
  for (;;) {
    if (defined (my $value = shift @$array)) {
      return ($self->{'i'}++,
              $value);
    }
    my $lo = $self->{'hi'} + 1;
    $self->{'hi'} = my $hi = $lo + 1000;
    @$array = _hypots_block ($self, $lo, $hi);
  }
}

sub _hypots_block {
  my ($self, $lo, $hi) = @_;

  if ($self->{'pythagorean_type'} eq 'primitive') {
    return grep {$self->pred($_)} $lo .. $hi;

  } else {
    my %hypots;
    foreach my $p (grep {($_ & 3)==1}
                   Math::NumSeq::Primes::_primes_list(2, $hi)) {
      @hypots{map {$_*$p}
                (int(($lo+$p-1)/$p) .. int($hi/$p))
              } = ();
    }
    return sort {$a<=>$b} keys %hypots;
  }
}

sub pred {
  my ($self, $value) = @_;
  ### pred: $value

  if ($value < 5
      || _is_infinite($value)
      || $value != int($value)) {
    return 0;
  }

  my $pythagorean_type = $self->{'pythagorean_type'};
  ### $pythagorean_type

  unless ($value % 2) {
    ### even ...
    if ($pythagorean_type ne 'all') {
      ### primitive and prime never even ...
      return 0;
    }
    do {
      $value /= 2;
    } until ($value % 2);
  }
  unless ($value <= 0xFFFF_FFFF) {
    return undef;
  }
  $value = "$value"; # numize Math::BigInt for speed

  my $limit = int(sqrt($value));

  for (my $p = 3; $p <= $limit; $p += 2) {
    if (($value % $p) == 0) {
      if (($p & 3) == 1) {
        ### found 4k+1 prime: $p
        if ($pythagorean_type eq 'all') {
          return 1;
        }
      } else {
        ### found 4k+3 prime: $p
        if ($pythagorean_type eq 'primitive') {
          return 0;
        }
      }

      do {
        $value /= $p;
      } while (($value % $p) == 0);

      $limit = int(sqrt($value));  # new smaller limit
      ### divide out prime: "$p new limit $limit"
    }
  }

  # $value now 1 or prime
  if ($pythagorean_type eq 'primitive') {
    # all, check this isn't a 4k+3 prime
    return ($value % 4) != 3;
  } else {
    # all, last chance to see a 4k+1 prime
    return ($value > 1
            && ($value % 4) == 1);
  }
}

1;
__END__