bigfloat.pm

视频监控网络部分的协议ddns,的模块的实现代码,请大家大胆指正.

  if (defined $params[0])
    {
    $x->bround($params[0],$params[2]);		# then round accordingly
    }
  else
    {
    $x->bfround($params[1],$params[2]);		# then round accordingly
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;

  $x;
  }

sub _len_to_steps
  {
  # Given D (digits in decimal), compute N so that N! (N factorial) is
  # at least D digits long. D should be at least 50.
  my $d = shift;

  # two constants for the Ramanujan estimate of ln(N!)
  my $lg2 = log(2 * 3.14159265) / 2;
  my $lg10 = log(10);

  # D = 50 => N => 42, so L = 40 and R = 50
  my $l = 40; my $r = $d;

  # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
  $l = $l->numify if ref($l);
  $r = $r->numify if ref($r);
  $lg2 = $lg2->numify if ref($lg2);
  $lg10 = $lg10->numify if ref($lg10);

  # binary search for the right value (could this be written as the reverse of lg(n!)?)
  while ($r - $l > 1)
    {
    my $n = int(($r - $l) / 2) + $l;
    my $ramanujan = 
      int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
    $ramanujan > $d ? $r = $n : $l = $n;
    }
  $l;
  }

sub bnok
  {
  # Calculate n over k (binomial coefficient or "choose" function) as integer.
  # set up parameters
  my ($self,$x,$y,@r) = (ref($_[0]),@_);

  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($self,$x,$y,@r) = objectify(2,@_);
    }

  return $x if $x->modify('bnok');

  return $x->bnan() if $x->is_nan() || $y->is_nan();
  return $x->binf() if $x->is_inf();

  my $u = $x->as_int();
  $u->bnok($y->as_int());

  $x->{_m} = $u->{value};
  $x->{_e} = $MBI->_zero();
  $x->{_es} = '+';
  $x->{sign} = '+';
  $x->bnorm(@r);
  }

sub bexp
  {
  # Calculate e ** X (Euler's number to the power of X)
  my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('bexp');

  return $x->binf() if $x->{sign} eq '+inf';
  return $x->bzero() if $x->{sign} eq '-inf';

  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my ($scale,@params);
  ($x,@params) = $x->_find_round_parameters($a,$p,$r);

  # also takes care of the "error in _find_round_parameters?" case
  return $x if $x->{sign} eq 'NaN';

  # no rounding at all, so must use fallback
  if (scalar @params == 0)
    {
    # simulate old behaviour
    $params[0] = $self->div_scale();	# and round to it as accuracy
    $params[1] = undef;			# P = undef
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r;			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    }
  else
    {
    # the 4 below is empirical, and there might be cases where it's not enough...
    $scale = abs($params[0] || $params[1]) + 4;	# take whatever is defined
    }

  return $x->bone(@params) if $x->is_zero();

  if (!$x->isa('Math::BigFloat'))
    {
    $x = Math::BigFloat->new($x);
    $self = ref($x);
    }
  
  # when user set globals, they would interfere with our calculation, so
  # disable them and later re-enable them
  no strict 'refs';
  my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
  my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
  # we also need to disable any set A or P on $x (_find_round_parameters took
  # them already into account), since these would interfere, too
  delete $x->{_a}; delete $x->{_p};
  # need to disable $upgrade in BigInt, to avoid deep recursion
  local $Math::BigInt::upgrade = undef;
  local $Math::BigFloat::downgrade = undef;

  my $x_org = $x->copy();

  # We use the following Taylor series:

  #           x    x^2   x^3   x^4
  #  e = 1 + --- + --- + --- + --- ...
  #           1!    2!    3!    4!

  # The difference for each term is X and N, which would result in:
  # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term

  # But it is faster to compute exp(1) and then raising it to the
  # given power, esp. if $x is really big and an integer because:

  #  * The numerator is always 1, making the computation faster
  #  * the series converges faster in the case of x == 1
  #  * We can also easily check when we have reached our limit: when the
  #    term to be added is smaller than "1E$scale", we can stop - f.i.
  #    scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
  #  * we can compute the *exact* result by simulating bigrat math:

  #  1   1    gcd(3,4) = 1    1*24 + 1*6    5
  #  - + -                  = ---------- =  --                 
  #  6   24                      6*24       24

  # We do not compute the gcd() here, but simple do:
  #  1   1    1*24 + 1*6   30
  #  - + -  = --------- =  --                 
  #  6   24       6*24     144

  # In general:
  #  a   c    a*d + c*b 	and note that c is always 1 and d = (b*f)
  #  - + -  = ---------
  #  b   d       b*d

  # This leads to:         which can be reduced by b to:
  #  a   1     a*b*f + b    a*f + 1
  #  - + -   = --------- =  -------
  #  b   b*f     b*b*f        b*f

  # The first terms in the series are:

  # 1     1    1    1    1    1     1     1     13700
  # -- + -- + -- + -- + -- + --- + --- + ---- = -----
  # 1     1    2    6   24   120   720   5040   5040

  # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!

  if ($scale <= 75)
    {
    # set $x directly from a cached string form
    $x->{_m} = $MBI->_new(
    "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
    $x->{sign} = '+';
    $x->{_es} = '-';
    $x->{_e} = $MBI->_new(79);
    }
  else
    {
    # compute A and B so that e = A / B.
 
    # After some terms we end up with this, so we use it as a starting point:
    my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
    my $F = $MBI->_new(42); my $step = 42;

    # Compute how many steps we need to take to get $A and $B sufficiently big
    my $steps = _len_to_steps($scale - 4);
#    print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
    while ($step++ <= $steps)
      {
      # calculate $a * $f + 1
      $A = $MBI->_mul($A, $F);
      $A = $MBI->_inc($A);
      # increment f
      $F = $MBI->_inc($F);
      }
    # compute $B as factorial of $steps (this is faster than doing it manually)
    my $B = $MBI->_fac($MBI->_new($steps));
    
#  print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";

    # compute A/B with $scale digits in the result (truncate, not round)
    $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
    $A = $MBI->_div( $A, $B );

    $x->{_m} = $A;
    $x->{sign} = '+';
    $x->{_es} = '-';
    $x->{_e} = $MBI->_new($scale);
    }

  # $x contains now an estimate of e, with some surplus digits, so we can round
  if (!$x_org->is_one())
    {
    # raise $x to the wanted power and round it in one step:
    $x->bpow($x_org, @params);
    }
  else
    {
    # else just round the already computed result
    delete $x->{_a}; delete $x->{_p};
    # shortcut to not run through _find_round_parameters again
    if (defined $params[0])
      {
      $x->bround($params[0],$params[2]);		# then round accordingly
      }
    else
      {
      $x->bfround($params[1],$params[2]);		# then round accordingly
      }
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;

  $x;						# return modified $x
  }

sub _log
  {
  # internal log function to calculate ln() based on Taylor series.
  # Modifies $x in place.
  my ($self,$x,$scale) = @_;

  # in case of $x == 1, result is 0
  return $x->bzero() if $x->is_one();

  # XXX TODO: rewrite this in a similiar manner to bexp()

  # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log

  # u = x-1, v = x+1
  #              _                               _
  # Taylor:     |    u    1   u^3   1   u^5       |
  # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
  #             |_   v    3   v^3   5   v^5      _|

  # This takes much more steps to calculate the result and is thus not used
  # u = x-1
  #              _                               _
  # Taylor:     |    u    1   u^2   1   u^3       |
  # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
  #             |_   x    2   x^2   3   x^3      _|

  my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);

  $v = $x->copy(); $v->binc();		# v = x+1
  $x->bdec(); $u = $x->copy();		# u = x-1; x = x-1
  $x->bdiv($v,$scale);			# first term: u/v
  $below = $v->copy();
  $over = $u->copy();
  $u *= $u; $v *= $v;				# u^2, v^2
  $below->bmul($v);				# u^3, v^3
  $over->bmul($u);
  $factor = $self->new(3); $f = $self->new(2);

  my $steps = 0 if DEBUG;  
  $limit = $self->new("1E-". ($scale-1));
  while (3 < 5)
    {
    # we calculate the next term, and add it to the last
    # when the next term is below our limit, it won't affect the outcome
    # anymore, so we stop

    # calculating the next term simple from over/below will result in quite
    # a time hog if the input has many digits, since over and below will
    # accumulate more and more digits, and the result will also have many
    # digits, but in the end it is rounded to $scale digits anyway. So if we
    # round $over and $below first, we save a lot of time for the division
    # (not with log(1.2345), but try log (123**123) to see what I mean. This
    # can introduce a rounding error if the division result would be f.i.
    # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
    # if we truncated $over and $below we might get 0.12345. Does this matter
    # for the end result? So we give $over and $below 4 more digits to be
    # on the safe side (unscientific error handling as usual... :+D

    $next = $over->copy->bround($scale+4)->bdiv(
      $below->copy->bmul($factor)->bround($scale+4), 
      $scale);

## old version:    
##    $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);

    last if $next->bacmp($limit) <= 0;

    delete $next->{_a}; delete $next->{_p};
    $x->badd($next);
    # calculate things for the next term
    $over *= $u; $below *= $v; $factor->badd($f);
    if (DEBUG)
      {
      $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
      }
    }
  print "took $steps steps\n" if DEBUG;
  $x->bmul($f);					# $x *= 2
  }

sub _log_10
  {
  # Internal log function based on reducing input to the range of 0.1 .. 9.99
  # and then "correcting" the result to the proper one. Modifies $x in place.
  my ($self,$x,$scale) = @_;

  # Taking blog() from numbers greater than 10 takes a *very long* time, so we
  # break the computation down into parts based on the observation that:
  #  blog(X*Y) = blog(X) + blog(Y)
  # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
  # $x is the faster it gets. Since 2*$x takes about 10 times as
  # long, we make it faster by about a factor of 100 by dividing $x by 10.

  # The same observation is valid for numbers smaller than 0.1, e.g. computing
  # log(1) is fastest, and the further away we get from 1, the longer it takes.
  # So we also 'break' this down by multiplying $x with 10 and subtract the
  # log(10) afterwards to get the correct result.

  # To get $x even closer to 1, we also divide by 2 and then use log(2) to
  # correct for this. For instance if $x is 2.4, we use the formula:
  #  blog(2.4 * 2) == blog (1.2) + blog(2)
  # and thus calculate only blog(1.2) and blog(2), which is faster in total
  # than calculating blog(2.4).

  # In addition, the values for blog(2) and blog(10) are cached.

  # Calculate nr of digits before dot:
  my $dbd = $MBI->_num($x->{_e});
  $dbd = -$dbd if $x->{_es} eq '-';
  $dbd += $MBI->_len($x->{_m});

  # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
  # infinite recursion

  my $calc = 1;					# do some calculation?

  # disable the shortcut for 10, since we need log(10) and this would recurse
  # infinitely deep
  if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
    {
    $dbd = 0;					# disable shortcut
    # we can use the cached value in these cases
    if ($scale <= $LOG_10_A)
      {
      $x->bzero(); $x->badd($LOG_10);		# modify $x in place
      $calc = 0; 				# no need to calc, but round
      }
    # if we can't use the shortcut, we continue normally
    }
  else
    {
    # disable the shortcut for 2, since we maybe have it cached
    if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
      {
      $dbd = 0;					# disable shortcut
      # we can use the cached value in these cases
      if ($scale <= $LOG_2_A)
        {
        $x->bzero(); $x->badd($LOG_2);		# modify $x in place
        $calc = 0; 				# no need to calc, but round
        }
      # if we can't use the shortcut, we continue normally
      }
    }

  # if $x = 0.1, we know the result must be 0-log(10)
  if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
      $MBI->_is_one($x->{_m}))
    {
    $dbd = 0;					# disable shortcut
    # we can use the cached value in these cases
    if ($scale <= $LOG_10_A)
      {
      $x->bzero(); $x->bsub($LOG_10);
      $calc = 0; 				# no need to calc, but round
      }
    }

  return if $calc == 0;				# already have the result

  # default: these correction factors are undef and thus not used
  my $l_10;				# value of ln(10) to A of $scale
  my $l_2;				# value of ln(2) to A of $scale

  my $two = $self->new(2);

  # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
  # so don't do this shortcut for 1 or 0
  if (($dbd > 1) || ($dbd < 0))
    {
    # convert our cached value to an object if not already (avoid doing this
    # at import() time, since not everybody needs this)
    $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;

    #print "x = $x, dbd = $dbd, calc = $calc\n";
    # got more than one digit before the dot, or more than one zero after the
    # dot, so do:
    #  log(123)    == log(1.23) + log(10) * 2
    #  log(0.0123) == log(1.23) - log(10) * 2
  
    if ($scale <= $LOG_10_A)
      {
      # use cached value
      $l_10 = $LOG_10->copy();		# copy for mul
      }
    else
      {
      # else: slower, compute and cache result
      # also disable downgrade for this code path