File Coverage

blib/lib/Math/PlanePath/FactorRationals.pm

Criterion	Covered	Total	%
statement	131	203	64.5
branch	33	62	53.2
condition	5	11	45.4
subroutine	23	28	82.1
pod	4	4	100.0
total	196	308	63.6

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 .


# Multiples of prime make grid.

# [13] L. S. Johnston, Denumerability of the rational number system, Amer. Math. Monthly, 55 (Feb.
#      1948), no. 2, 65-70.
# www.jstor.org/stable/2305738


# prime factors q1,..qk of n
# f(m/n) = m^2*n^2/ (q1q2...qk)

# Kevin McCrimmon, 1960
#
# integer prod p[i]^a[i] -> rational prod p[i]^b[i]
# b[i] = a[2i-1] if a[2i-1]!=0
#    b[1]=a[1], b[2]=a[3], b[3]=a[5]
# b[i] = -a[2k] if a[2i-1]=0 and is kth such
#
# b[i] = f(a[i]) where f(n) = (-1)^(n+1) * floor((n+1)/2)
#   f(0) =  0
#   f(1) =  1
#   f(2) = -1
#   f(3) =  2
#   f(4) = -2

# Gerald Freilich, 1965
#
# f(n) = n/2      if n even n>=0
#      = -(n+1)/2 if n odd n>0
# f(0)=0/2      =  0
# f(1)=-(1+1)/2 = -1
# f(2)=2/2      =  1
# f(3)=-(3+1)/2 = -2
# f(4)=4/2      =  2

# Yoram Sagher, "Counting the rationals", American Math Monthly, Nov 1989,
# page 823.  http://www.jstor.org/stable/2324846
#
# m = p1^e1.p2^e2...
# n = q1^f1.q2^f2...
# f(m/n) = p1^2e1.p2^2e2... . q1^(2f1-1).q2^(2f2-1)...
# so     0 -> 0              0 ->  0
#    num 1 -> 2              1 -> -1
#        2 -> 4              2 ->  1
# den -1 1 -> 2*1-1 = 1      3 -> -2
#     -2 2 -> 2*2-1 = 3      4 ->  2

# Umberto Cerruti, "Ordinare i razionali Gli alberi di Keplero e di
# Calkin-Wilf", following T.J. Heard
#   B(2k)=-k   even=negative and zero
#   B(2k-1)=k  odd=positive
#   which is Y/X invert
# B(0 =2*0)   =  0
# B(1 =2*1-1) =  1
# B(2 =2*1)   = -1
# B(3 =2*2-1) =  2
# B(4 =2*2)   = -2


package Math::PlanePath::FactorRationals;
use 5.004;
use strict;
use Carp 'croak';
use List::Util 'min';
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

use vars '$VERSION', '@ISA';
$VERSION = 128;
use Math::PlanePath;
@ISA = ('Math::PlanePath');

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'digit_join_lowtohigh';

use Math::PlanePath::CoprimeColumns;
*_coprime = \&Math::PlanePath::CoprimeColumns::_coprime;

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


# Not yet.
use constant parameter_info_array =>
  [ { name      => 'factor_coding',
      display   => 'Sign Encoding',
      type      => 'enum',
      default   => 'even/odd',
      choices         => ['even/odd','odd/even',
                          'negabinary','revbinary',
                         ],
      choices_display => ['Even/Odd','Odd/Even',
                          'Negabinary','Revbinary',
                          ],
    },
  ];

use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant x_minimum => 1;
use constant y_minimum => 1;
use constant gcdxy_maximum => 1;  # no common factor
use constant absdy_minimum => 1;

# factor_coding=even/odd
# factor_coding=odd/even
#   dir_minimum_dxdy() suspect dir approaches 0.
#   Eg. N=5324   = 2^2.11^3     dx=3,dy=92   0.97925
#       N=642735 = 3^5.23^2     dX=45 dY=4    Dir4=0.05644
#         642736 = 2^4.17^2.139
#   dir_maximum_dxdy() suspect approaches 360 degrees
#   use constant dir_maximum_dxdy => (0,0);  # the default
#
# factor_coding=negabinary
#   dir_minimum_dxdy() = East 1,0 at N=1
#   dir_maximum_dxdy() believe approaches 360 degrees
#   Eg. N=40=2^3.5 X=5, Y=2
#       N=41=41    X=41, Y=1
#   N=multiple 8 and solitary primes, followed by N+1=prime is dX=big, dY=-1
#
# factor_coding=revbinary
#   dir_maximum_dxdy() approaches 360 degrees  dY=-1, dX=big
#   Eg. N=7208=2^3*17*53 X=17*53  Y=2
#       N=7209=3^4*89    X=3^4*89 Y=1
#                       dX=6308  dY=-1


#------------------------------------------------------------------------------
# even/odd

# $n>=0, return a positive if even or negative if odd
#   $n==0  return  0
#   $n==1  return -1
#   $n==2  return +1
#   $n==3  return -2
#   $n==4  return +2
sub _pos_to_pn__even_odd {
  my ($n) = @_;
  return ($n % 2 ? -1-$n : $n) / 2;
}

# # $n is positive or negative, return even for positive or odd for negative.
# #   $n==0   return 0
# #   $n==-1  return 1
# #   $n==+1  return 2
# #   $n==-2  return 3
# #   $n==+2  return 4
# sub _pn_to_pos__even_odd {
#   my ($n) = @_;
#   return ($n >= 0 ? 2*$n : -1-2*$n);
# }

#------------------------------------------------------------------------------
# odd/even

# $n>=0, return a positive if even or negative if odd
#   $n==0  return  0
#   $n==1  return +1
#   $n==2  return -1
#   $n==3  return +2
#   $n==4  return -2
sub _pos_to_pn__odd_even {
  my ($n) = @_;
  return ($n % 2 ? $n+1 : -$n) / 2;
}

# # $n is positive or negative, return odd for positive or even for negative.
# #   $n==0   return 0
# #   $n==+1  return 1
# #   $n==-1  return 2
# #   $n==+2  return 3
# #   $n==-2  return 4
# sub _pn_to_pos__odd_even {
#   my ($n) = @_;
#   return ($n <= 0 ? -2*$n : 2*$n-1);
# }

#------------------------------------------------------------------------------
# negabinary

sub _pn_to_pos__negabinary {
  my ($n) = @_;
  my @bits;
  while ($n) {
    my $bit = ($n % 2);
    push @bits, $bit;
    $n -= $bit;
    $n /= 2;
    $n = -$n;
  }
  return digit_join_lowtohigh(\@bits, 2,
                              $n); # zero
}
sub _pos_to_pn__negabinary {
  my ($n) = @_;
  return (($n & 0x55555555) - ($n & 0xAAAAAAAA));
}

#------------------------------------------------------------------------------
# revbinary
# A065620  pos -> pn
# A065621  pn(+ve) -> pos
# A048724  pn(-ve) -> pos        n XOR 2n
# A048725  A048724 twice
# cf
# A073122  minimizing by taking +/- powers  cf A072219 A072339

# rev = 2^e0 - 2^e1 + 2^e2 - 2^e3 + ... + (-1)^t*2^et
#       0 <= e0 < e1 < e2 ...
sub _pos_to_pn__revbinary {
  my ($n) = @_;
  my $sign = 1;
  my $ret = 0;
  for (my $bit = 1; $bit <= $n; $bit *= 2) {
    if ($n & $bit) {
      $ret += $bit*$sign;
      $sign = -$sign;
    }
  }
  return $ret;
}
sub _pn_to_pos__revbinary {
  my ($n) = @_;
  my @bits;
  while ($n) {
    my $bit = ($n % 2);
    push @bits, $bit;
    $n -= $bit;
    $n /= 2;
    if ($bit) {
      $n = -$n;
    }
  }
  return digit_join_lowtohigh(\@bits, 2,
                              $n); # zero
}

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

my %factor_coding__pos_to_pn = ('even/odd' => \&_pos_to_pn__even_odd,
                                'odd/even' => \&_pos_to_pn__odd_even,
                                negabinary => \&_pos_to_pn__negabinary,
                                revbinary  => \&_pos_to_pn__revbinary,
                               );
my %factor_coding__pn_to_pos = (# 'even/odd' => \&_pn_to_pos__even_odd,
                                # 'odd/even' => \&_pn_to_pos__odd_even,
                                negabinary => \&_pn_to_pos__negabinary,
                                revbinary  => \&_pn_to_pos__revbinary,
                               );

sub new {
  my $self = shift->SUPER::new(@_);

  my $factor_coding = ($self->{'factor_coding'} ||= 'even/odd');
  $factor_coding__pos_to_pn{$factor_coding}
    or croak "Unrecognised factor_coding: ",$factor_coding;

  return $self;
}

sub n_to_xy {
  my ($self, $n) = @_;
  ### FactorRationals n_to_xy(): "$n"

  if ($n < 1) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  # what to do for fractional $n?
  {
    my $int = int($n);
    if ($n != $int) {
      ### frac ...
      my $frac = $n - $int;  # inherit possible BigFloat/BigRat
      my ($x1,$y1) = $self->n_to_xy($int);
      my ($x2,$y2) = $self->n_to_xy($int+1);
      my $dx = $x2-$x1;
      my $dy = $y2-$y1;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
    $n = $int;
  }

  my $zero = $n * 0;

  my $pos_to_pn = $factor_coding__pos_to_pn{$self->{'factor_coding'}};
  my $x = my $y = ($n * 0) + 1;  # inherit bignum 1
  my ($limit,$overflow) = _limit($n);
  ### $limit
  my $divisor = 2;
  my $dstep = 1;
  while ($divisor <= $limit) {
    if (($n % $divisor) == 0) {
      my $count = 0;
      for (;;) {
        $count++;
        $n /= $divisor;
        if ($n % $divisor) {
          my $pn = &$pos_to_pn($count);
          ### $count
          ### $pn
          my $pow = ($divisor+$zero) ** abs($pn);
          if ($pn >= 0) {
            $x *= $pow;
          } else {
            $y *= $pow;
          }
          last;
        }
      }
      ($limit,$overflow) = _limit($n);
      ### $limit
    }
    $divisor += $dstep;
    $dstep = 2;
  }
  if ($overflow) {
    ### n too big ...
    return;
  }

  ### remaining $n is prime, count=1: "n=$n"
  my $pn = &$pos_to_pn(1);
  ### $pn
  my $pow = $n ** abs($pn);
  if ($pn >= 0) {
    $x *= $pow;
  } else {
    $y *= $pow;
  }

  ### result: "$x, $y"
  return ($x, $y);
}

sub xy_to_n {
  my ($self, $x, $y) = @_;

  $x = round_nearest ($x);
  $y = round_nearest ($y);
  ### FactorRationals xy_to_n(): "x=$x y=$y"

  if ($x < 1 || $y < 1) {
    return undef;  # negatives and -infinity
  }
  if (is_infinite($x)) { return $x; } # +infinity or nan
  if (is_infinite($y)) { return $y; } # +infinity or nan

  if ($self->{'factor_coding'} eq 'negabinary'
      || $self->{'factor_coding'} eq 'revbinary') {
    ### negabinary or revbinary ...
    my $pn_to_pos = $factor_coding__pn_to_pos{$self->{'factor_coding'}};
    my $n = 1;
    my $zero = $x * 0 * $y;

    # Factorize both $x and $y and apply their pn_to_pos encoded powers to
    # make $n.  A common factor between $x and $y is noticed if $divisor
    # divides both.

    my ($limit,$overflow) = _limit(max($x,$y));
    my $dstep = 1;
    for (my $divisor = 2; $divisor <= $limit; $divisor += $dstep, $dstep=2) {
      my $count = 0;
      if ($x % $divisor == 0) {
        if ($y % $divisor == 0) {
          return undef;  # common factor
        }
        while ($x % $divisor == 0) {
          $count++;
          $x /= $divisor;  # mutate loop variable
        }
      } elsif ($y % $divisor == 0) {
        while ($y % $divisor == 0) {
          $count--;
          $y /= $divisor;  # mutate loop variable
        }
      } else {
        next;
      }

      # Here $count > 0 if from $x or $count < 0 if from $y.
      ### $count
      ### pn: &$pn_to_pos($count)

      $count = &$pn_to_pos($count);
      $n *= ($divisor+$zero) ** $count;

      # new search limit, perhaps smaller than before
      ($limit,$overflow) = _limit(max($x,$y));
    }

    if ($overflow) {
      ### x,y too big to find all primes ...
      return undef;
    }

    # Here $x and $y are primes.
    if ($x > 1 && $x == $y) {
      ### common factor final remaining prime x,y ...
      return undef;
    }

    # $x is power p^1 which is negabinary=1 or revbinary=1 so multiply into
    # $n.  $y is power p^-1 and -1 is negabinary=3 so cube and multiply into
    # $n.
    $n *= $x;
    $n *= $y*$y*$y;

    return $n;

  } else {
    ### assert: $self->{'factor_coding'} eq 'even/odd' || $self->{'factor_coding'} eq 'odd/even'
    if ($self->{'factor_coding'} eq 'odd/even') {
      ($x,$y) = ($y,$x);
    }

    # Factorize $y so as to make an odd power of its primes.  Only need to
    # divide out one copy of each prime, but by dividing out them all the
    # $limit to search up to is reduced, usually by a lot.
    #
    # $ymult is $y with one copy of each prime factor divided out.
    # $ychop is $y with all primes divided out as they're found.
    # $y itself is unchanged.
    #
    my $ychop = my $ymult = $y;

    my ($limit,$overflow) = _limit($ychop);
    my $dstep = 1;
    for (my $divisor = 2; $divisor <= $limit; $divisor += $dstep, $dstep=2) {
      next if $ychop % $divisor;

      if ($x % $divisor == 0) {
        ### common factor with X ...
        return undef;
      }
      $ymult /= $divisor;           # one of $divisor divided out
      do {
        $ychop /= $divisor;         # all of $divisor divided out
      } until ($ychop % $divisor);
      ($limit,$overflow) = _limit($ychop);  # new lower $limit, perhaps
    }

    if ($overflow) {
      return undef; # Y too big to find all primes
    }

    # remaining $ychop is a prime, or $ychop==1
    if ($ychop > 1) {
      if ($x % $ychop == 0) {
        ### common factor with X ...
        return undef;
      }
      $ymult /= $ychop;
    }

    return $x*$x * $y*$ymult;
  }
}

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

# all rationals X,Y >= 1 with no common factor
use Math::PlanePath::DiagonalRationals;
*xy_is_visited = Math::PlanePath::DiagonalRationals->can('xy_is_visited');

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

# even/odd
#   X=2^10 -> N=2^20 is X^2
#   Y=3 -> N=3
#   Y=3^2 -> N=3^3
#   Y=3^3 -> N=3^5
#   Y=3^4 -> N=3^7
#   Y*Y / distinct prime factors
#
# negabinary
#   X=prime^2 -> N=prime^6       is X^3
#   X=prime^6 -> N=prime^26      is X^4.33
#   maximum 101010...10110 -> 1101010...10 approaches factor 5
#   same for negatives
#
# revbinary
#   X=prime^k -> N=prime^(3k)    ix X^3

# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### rect_to_n_range()

  $x1 = round_nearest ($x1);
  $y1 = round_nearest ($y1);
  $x2 = round_nearest ($x2);
  $y2 = round_nearest ($y2);

  my $n = max($x1,$x2) * max($y1,$y2);
  my $n_squared = $n * $n;
  return (1,
          ($self->{'factor_coding'} eq 'negabinary'
           ? ($n_squared*$n_squared) * $n     # X^5*Y^5
           : $self->{'factor_coding'} eq 'revbinary'
           ? $n_squared * $n                  # X^3*Y^3
           # even/odd, odd/even
           : $n_squared));                    # X^2*Y^2
}


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

# _limit() returns ($limit,$overflow).
#
# $limit is the biggest divisor to attempt trial division of $n.  If $n <
# 2^32 then $limit=sqrt($n) and that will find all primes.  If $n >= 2^32
# then $limit is smaller than sqrt($n), being calculated from the length of
# $n so as to make a roughly constant amount of time doing divisions.  But
# $limit is always at least 50 so as to divide by primes up to 50.
#
# $overflow is a boolean, true if $n is too big to search for all primes and
# $limit is something smaller than sqrt($n).  $overflow is false if $limit
# has not been capped and is enough to find all primes.
#
sub _limit {
  my ($n) = @_;
  my $limit = int(sqrt($n));
  my $cap = max (int(65536 * 10 / length($n)),
                 50);
  if ($limit > $cap) {
    return ($cap, 1);
  } else {
    return ($limit, 0);
  }
}

1;
__END__