view release on metacpan or  search on metacpan

lib/Math/NumSeq/ReRound.pm  view on Meta::CPAN

  ### ReRound pred(): $value

  my $extra_multiples = $self->{'extra_multiples'};

  if (_is_infinite($value)) {
    return undef;
  }
  if ($value <= 1 || $value != int($value)) {
    return ($value == 1);
  }

  # special case m=1 stepping down to an even number
  if (($value -= $extra_multiples) % 2) {
    return 0;
  }

  my $m = 2;
  while ($value > $m) {
    ### at: "value=$value  m=$m"

    if (($value -= $extra_multiples*$m) <= 0) {
      ### no, negative: $value
      return 0;
    }
    ### subtract to: "value=$value"

    ### rem: "modulus=".($m+1)." rem ".($value%($m+1))
    my $rem;
    if (($rem = ($value % ($m+1))) == $m) {
      ### no, remainder: "rem=$rem  modulus=".($m+1)
      return 0;
    }

    $value -= $rem;
    $m++;
  }

  ### final ...
  ### $value
  ### $m

  return ($value == $m);
}

# # v1.02 for leading underscore
use constant 1.02 _PI => 4*atan2(1,1); # similar to Math::Complex pi()

# value = m*1+m*2+m*3+...+m*i
# value = m*i*(i+1)/2
# i*i + i - 2v/m = 0
# i = (-1 + sqrt(1 + 4*2v/m)) / 2
#   = -1/2 + sqrt(1 + 4*2v/m)/2
#   = -1/2 + sqrt(4 + 2v/m)
# i -> sqrt(2v/m)
#
# as extra_multiples dominates the estimate changes from
#     extra_multiples==0         est = sqrt(pi*value)
# up to
#     extra_multiples==m large   est = sqrt(2*value/m)
#
# What formula that progressively morphs pi to 2/m ?
#
sub value_to_i_estimate {
  my ($self, $value) = @_;
  if ($value < 0) { return 0; }

  if ($self->{'extra_multiples'} == 0) {
    # use sqrt(pi)~=296/167 to preserve Math::BigInt
    return int( (sqrt(int($value)) * 296) / 167);
  } else {
    return int(sqrt(int (2 * $value / $self->{'extra_multiples'})));

    # return int(sqrt(int ($value * _PI/($self->{'extra_multiples'}+1)**2))
    #            + ($self->{'extra_multiples'} > 0
    #               ? -0.5 + sqrt(4 + 2*$value/$self->{'extra_multiples'})
    #               : 0));
  }



  # # large extra_multiples
  # return int(sqrt(int (2 * $value
  #                      / ($self->{'extra_multiples'}+1))));
  #
  # return int(sqrt(int (($value * 355)
  #                      / (113 * ($self->{'extra_multiples'}+1)))));
  #
  #
  # # extra_multiples==14
  # return int(sqrt(int (_PI * $value)) * 429/2048); #        429=13*11*3
  #
  # # extra_multiples==12
  # return int(sqrt(int (_PI * $value)) * 462/2048); # 231/1024) # 11*7*3=231
  #
  # # extra_multiples==10
  # return int(sqrt(int (_PI * $value))) * 126/512; # 63/256;  # 7*3*2=126
  #
  # # extra_multiples==8
  # return int(sqrt(int (_PI * $value))) * 35/128;             # 7*5=35
  #
  # # extra_multiples==6
  # return int(sqrt(int (_PI * $value))) * 10/32; # 5/16       # 5*2=10
  #
  # # extra_multiples==4
  # return int(sqrt(int (_PI * $value))) * 3/8;
  #
  # # extra_multiples==2
  # return int(sqrt(int (_PI * $value))) * 1/2;
  #
  # # extra_multiples==0
  # return int(sqrt(int (_PI * $value))) * 1/1;
  #
  #  2 1
  #  4 11
  #  6 1010
  #  8 100011
  # 10 1111110
  # 12 111001110
  # 14 110101101

  #

lib/Math/NumSeq/ReRound.pm  view on Meta::CPAN

=head2 Predicate

The rounding procedure can be reversed to test for a ReRound value.

    for i=2,3,4,etc
      value -= extra_multiples*i
      if value < 0 then not a ReRound

      remainder = value mod i
      if remainder==i-1 then not a ReRound

      value -= remainder    # round down to multiple of i

    stop when value <= i
    is a ReRound if value==i, and i is its index

For example to test 28, it's a multiple of 2, so ok for the final rounding
step.  It's predecessor in the rounding steps was a multiple of 3, so round
down to a multiple of 3 which is 27.  The predecessor of 27 was a multiple
of 4 so round down to 24.  But at that point there's a contradiction because
if 24 was the value then it's already a multiple of 3 and so wouldn't have
gone up to 27.  This case where a round-down gives a multiple of both i and
i-1 is identified by the remainder = value % i == i-1.  If value is already
a multiple of i-1 then subtracting an i-1 would leave it still so.

=head2 Value to i Estimate

For the default sequence as noted above the values grow roughly as

    value ~= i^2 / pi

so can be estimated as

    i ~= sqrt(value) * sqrt(pi)

There's no need for high accuracy in pi.  The current code uses an
approximation sqrt(pi)~=296/167 for faster calculation if the value is a
C<Math::BigInt> or C<Math::BigRat>.

    i ~= sqrt(value) * 296/167

C<extra_multiples =E<gt> m> adds a fixed amount i*m at each step, for a
total

    value_EM = m + 2*m + 3*m + 4*m + ... + i*m
             = m * i*(i+1)/2

As m becomes large this part of the value dominates, so

    value =~ m * i*(i+1)/2     for large m

    i =~ sqrt(value*2/m)

As m increases the factor sqrt(pi) progressively morphs towards sqrt(2/m).
For m even it might be sqrt(pi)/2, 3/8*sqrt(pi), 5/16*sqrt(pi),
35/128*sqrt(pi), etc, and for m odd it might be 2*sqrt(1/pi),
4/3*sqrt(1/pi), 16/15*sqrt(1/pi), 32/35*sqrt(1/pi), 256/315*sqrt(1/pi), etc.
What's the pattern?

The current code uses the "large m" formula for any mE<gt>0, which is no
more than roughly a factor 1.25 too big.

=head1 SEE ALSO

L<Math::NumSeq>,
L<Math::NumSeq::ReReplace>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 2011, 2012, 2013, 2014, 2016, 2019, 2020 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