#!/usr/bin/perl -w

# Copyright 2012, 2013 Kevin Ryde

# This file is part of Math-PlanePath.
#
# Math-PlanePath 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-PlanePath 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-PlanePath.  If not, see <http://www.gnu.org/licenses/>.


use 5.010;
use strict;
use List::Util 'min', 'max';
use Math::PlanePath::FactorRationals;

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

{
  foreach my $n (1 .. 20) {
    print Math::PlanePath::FactorRationals::_pos_to_pn__negabinary($n),",";
  }
  exit 0;
}

{
  # different pos=49 numbers got=69 want=88, and more diff
  # N=50 = 5*5*2
  my $path = Math::PlanePath::FactorRationals->new;
  foreach my $x (1 .. 50) {
    my $n = $path->xy_to_n(1,$x);
    print "$x   $n\n";
  }
  exit 0;
}


# Return ($good, $prime,$exp, $prime,$exp,...).
# $good is true if a full factorization is found.
# $good is false if cannot factorize because $n is too big or infinite.
#
# If $n==0 or $n==1 then there are no prime factors and the return is
# $good=1 and an empty list of primes.
#
sub INPROGRESS_prime_factors_and_exps {
  my ($n) = @_;
  ### _prime_factors(): $n

  unless ($n >= 0) {
    return 0;
  }
  if (_is_infinite($n)) {
    return 0;
  }

  # if ($n <= 0xFFFF_FFFF) {
  #   return (1, prime_factors($n));
  # }

  my @ret;
  unless ($n % 2) {
    my $count = 0;
    do {
      $count++;
      $n /= 2;
    } until ($n % 2);
    push @ret, 2, $count;
  }

  # Stop at when prime $p reaches $limit and when no prime factor has been
  # found for the last 20 attempted $p.  Stopping only after a run of no
  # factors found allows big primorials 2*3*5*7*13*... to be divided out.
  # If the divisions are making progress reducing $i then continue.
  #
  # Would like $p and $gap to count primes, not just odd numbers.  Perhaps
  # a table of small primes.  The first gap of 36 odds between primes
  # occurs at prime=31469.  cf A000230 smallest prime p for gap 2n.

  my $limit = 10_000 / (_blog2_estimate($n) || 1);
  my $gap = 0;
  for (my $p = 3; $gap < 36 || $p <= $limit ; $p += 2) {
    if ($n % $p) {
      $gap++;
    } else {
      do  {
        ### prime: $p
        $n /= $p;
        push @ret, $p;
      } until ($n % $p);

      if ($n <= 1) {
        ### all factors found ...
        return (1, @ret);
      }
      # if ($n < 0xFFFF_FFFF) {
      #   ### remaining factors by XS ...
      #   return (1, @ret, prime_factors($n));
      # }
      $gap = 0;
    }
  }
  return 0;  # factors too big
}

{
  my @primes = (2,3,5,7);
  sub _extend_primes {
    for (my $p = $primes[-1] + 2; ; $p += 2) {
      if (_is_prime($p)) {
        push @primes, $p;
        return;
      }
    }
  }
  sub _is_prime {
    my ($n) = @_;
    my $limit = int(sqrt($n));
    for (my $i = 0; ; $i++) {
      if ($i > $#primes) { _extend_primes(); }
      my $prime = $primes[$i];
      if ($n % $prime == 0) { return 0; }
      if ($prime > $limit) { return 1; }
    }
  }

  # $aref is an arrayref of prime exponents, [a,b,c,...]
  # Return their product 2**a * 3**b * 5**c * ...
  #
  sub _factors_join {
    my ($aref, $zero) = @_;
    ### _factors_join(): $aref
    my $n = $zero + 1;
    for (my $i = 0; $i <= $#$aref; $i++) {
      if ($i > $#primes) { _extend_primes(); }
      $n *= ($primes[$i] + $zero) ** $aref->[$i];
    }
    ### join: $n
    return $n;
  }

  # Return an arrayref of prime exponents of $n.
  # Eg. [a,b,c,...] for $n == 2**a * 3**b * 5**c * ...
  sub _factors_split {
    my ($n) = @_;
    ### _factors_split(): $n
    my @ret;
    for (my $i = 0; $n > 1; $i++) {
      if ($i > 6541) {
        ### stop, primes too big ...
        return;
      }
      if ($i > $#primes) { _extend_primes(); }

      my $count = 0;
      while ($n % $primes[$i] == 0) {
        $n /= $primes[$i];
        $count++;
      }
      push @ret, $count;
    }
    return \@ret;
  }

  # ### f: 2*3*3*5*19
  # ### f: _factors_split(2*3*3*5*19)
  # ### f: _factors_join(_factors_split(2*3*3*5*19),0)


  # factor_coding => 'spread'

  # "spread"
  # if ($self->{'factor_coding'} eq 'spread') {
  #   # N = 2^e1 * 3^e2 * 5^e3 * 7^e4 * 11^e5 * 13^e6 * 17^e7
  #   # X = 2^e1 * 3^e3 * 5^e5 * 7^e7,  Y = 1
  #   #
  #   # X = 2^e1        * 5^e5          e3=0,e7=0
  #   # Y =        3^e2        * 7^e4
  #   #
  #   # X=1,0,1
  #   # Y=0,0,0
  #   # 22 = 1,0,0,0,1
  #   # num = 1,0,1 = 2*5 = 10
  #   #
  #   my $xexps = _factors_split($x)
  #     or return undef;  # overflow
  #   my $yexps = _factors_split($y)
  #     or return undef;  # overflow
  #   ### $xexps
  #   ### $yexps
  #
  #   my @nexps;
  #   my $denpos = -1; # to store first at $nexps[1]
  #   while (@$xexps || @$yexps) {
  #     my $xexp = shift @$xexps || 0;
  #     my $yexp = shift @$yexps || 0;
  #     ### @nexps
  #     ### $xexp
  #     ### $yexp
  #     push @nexps, $xexp, 0;
  #     if ($xexp) {
  #       if ($yexp) {
  #         ### X,Y common factor ...
  #         return undef;
  #       }
  #     } else {
  #       ### den store to: "denpos=".($denpos+2)."  yexp=$yexp"
  #       $nexps[$denpos+=2] = $yexp;
  #     }
  #   }
  #   ### @nexps
  #   return (_factors_join(\@nexps, $x*0*$y));
  #
  # } els

  # if ($self->{'factor_coding'} eq 'spread') {
  #   # N = 2^e1 * 3^e2 * 5^e3 * 7^e4 * 11^e5 * 13^e6 * 17^e7
  #   # X = 2^e1 * 3^e3 * 5^e5 * 7^e7,  Y = 1
  #   #
  #   # X = 2^e1        * 5^e5          e3=0,e7=0
  #   # Y =        3^e2        * 7^e4
  #   #
  #   # 22 = 1,0,0,0,1
  #   # num = 1,0,1 = 2*5 = 10
  #   # den = 0
  #   #
  #   my $nexps = _factors_split($n)
  #     or return;  # too big
  #   ### $nexps
  #   my @dens;
  #   my (@xexps, @yexps);
  #   while (@$nexps || @dens) {
  #     my $exp = shift @$nexps;
  #     if (@$nexps)  {
  #       push @dens, shift @$nexps;
  #     }
  #
  #     if ($exp) {
  #       ### to num: $exp
  #       push @xexps, $exp;
  #       push @yexps, 0;
  #     } else {
  #       ### zero take den: $dens[0]
  #       push @xexps, 0;
  #       push @yexps, shift @dens;
  #     }
  #   }
  #   ### @xexps
  #   ### @yexps
  #   return (_factors_join(\@xexps,$zero),
  #           _factors_join(\@yexps,$zero));
  #
  # } else
}

{
  # reversing binary, max factor=3
  # 0 0  0  fac=0
  # 1 1  1  fac=1
  # 2 2  2  fac=1
  # 3 -1  3  fac=3
  # 4 4  4  fac=
  # 5 -3  5  fac=
  # 6 -2  6  fac=3
  # 7 3  7  fac=
  # 8 8  8  fac=
  # 9 -7  9  fac=
  # 10 -6  10  fac=
  # 11 7  11  fac=
  # 12 -4  12  fac=3
  # 13 5  13  fac=
  # 14 6  14  fac=
  # 15 -5  15  fac=3
  # 16 16  16  fac=

  my $max_fac = 0;
  foreach my $n (0 .. 2**20) {
    my $pn = Math::PlanePath::FactorRationals::_pos_to_pn__revbinary($n);
    my $ninv = Math::PlanePath::FactorRationals::_pn_to_pos__revbinary($pn);

    my $fac = $n / abs($pn||1);
    if ($fac >= $max_fac) {
      $max_fac = $fac;
    } else {
      $fac = '';
    }
    print "$n $pn  $ninv  fac=$fac\n";

    die unless $ninv == $n;
  }
  print "\n";
  exit 0;
}

{
  # negabinary, max factor approach 5
  my %rev;
  my $max_fac = 0;
  foreach my $n (0 .. 2**20) {
    my $power = 1;
    my $nega = 0;
    for (my $bit = 1; $bit <= $n; $bit <<= 1) {
      if ($n & $bit) {
        $nega += $power;
      }
      $power *= -2;
    }
    my $fnega = Math::PlanePath::FactorRationals::_pos_to_pn__negabinary($n);
    my $ninv = Math::PlanePath::FactorRationals::_pn_to_pos__negabinary($nega);

    my $fac = -$n / ($nega||1);
    if ($fac > $max_fac) {
      $max_fac = $fac;
      print "$n $nega   $fnega $ninv  fac=$fac\n";
    } else {
      $fac = '';
    }
    $rev{$nega} = $n;
  }
  print "\n";
  exit 0;
  foreach my $nega (sort {$a<=>$b} keys %rev) {
    my $n = $rev{$nega};
    print "$nega $n\n";
  }
  exit 0;
}