File Coverage

lib/Math/Big.pm

Criterion	Covered	Total	%
statement	359	412	87.1
branch	93	150	62.0
condition	13	29	44.8
subroutine	25	26	96.1
pod	19	19	100.0
total	509	636	80.0

#############################################################################
# Math/Big.pm -- usefull routines with Big numbers (BigInt/BigFloat)

package Math::Big;

require 5.006002;	# anything lower is simple untested

use strict;
use warnings;

use Math::BigInt '1.97';
use Math::BigFloat;
use Exporter;

our $VERSION   = '1.14';
our @ISA       = qw( Exporter );
our @EXPORT_OK = qw( primes fibonacci base hailstone factorial
                     euler bernoulli pi log
                     tan cos sin cosh sinh arctan arctanh arcsin arcsinh
                  );
our @F;						# for fibonacci()

# some often used constants:
my $four    = Math::BigFloat->new(4);
my $sixteen = Math::BigFloat->new(16);
my $fone    = Math::BigFloat->bone();		# pi
my $one     = Math::BigInt->bone();		# hailstone, sin, cos etc
my $two     = Math::BigInt->new(2);		# hailstone, sin, cos etc
my $three   = Math::BigInt->new(3);		# hailstone

my $five         = Math::BigFloat->new(5);	# for pi
my $twothreenine = Math::BigFloat->new(239);	# for pi

# In scalar context this returns the prime count (# of primes <= N).
# In array context it returns a list of primes from 2 to N.
sub primes {
  my $end = shift;
  return unless defined $end;
  $end = $end->numify() if ref($end) =~ /^Math::Big/;
  if ($end < 2) { return !wantarray ? Math::BigInt->bzero() : (); }
  if ($end < 3) { return !wantarray ? $one->copy : ($two->copy); }
  if ($end < 5) { return !wantarray ? $two->copy : ($two->copy, $three->copy); }

  $end-- unless ($end & 1);
  my $s_end = $end >> 1;
  my $whole = int( ($end>>1) / 15);
  # Be conservative.  This would result in terabytes of array output.
  die "Cannot return $end primes!" if $whole > 1_145_324_612;  # ~32 GB string
  my $sieve = "100010010010110" . "011010010010110" x $whole;
  substr($sieve, $s_end+1) = ''; # Clip to the right number of entries
  my ($n, $limit) = ( 7, int(sqrt($end)) );
  while ( $n <= $limit ) {
    for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
      substr($sieve, $s, 1) = '1';
    }
    do { $n += 2 } while substr($sieve, $n>>1, 1);
  }

  return Math::BigInt->new(1 + $sieve =~ tr/0//) if !wantarray;

  my @primes = (2, 3, 5);
  $n = 7-2;
  foreach my $s (split("0", substr($sieve, 3), -1)) {
    $n += 2 + 2 * length($s);
    push @primes, $n if $n <= $end;
  }
  return map { Math::BigInt->new($_) } @primes;
}

sub fibonacci
  {
  my $n = shift; return unless defined $n;

  if (ref($n))
    {
    return if $n->sign() ne '+';		# < 0, NaN, inf
    }
  else
    {
    return if $n < 0;
    }

  #####################
  # list context
  if (wantarray)
    {
    # make a scalar (if $n doesn't fit into a scalar, the resulting array
    # will be too big, anyway)
    $n = $n->numify() if ref($n);

    my @fib = (Math::BigInt -> bzero(),Math::BigInt -> bone(),Math::BigInt -> bone);
    my $i = 3;							# no BigInt
    while ($i <= $n)
      {
      $fib[$i] = $fib[$i-1]+$fib[$i-2]; $i++;
      }
    return @fib;
    }
  #####################
  # scalar context

  _fibonacci_fast($n);
  }

my $F;

BEGIN
  {
  # Sequence of the first few fibonacci numbers:

  #     0,1,2,3,4,5,6,7, 8, 9, 10,11,12, 13, 14, 15, 16, 17,  18,  19,  20,
  @F = (0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765,
  #     21,
        10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229,
  #     30,
        832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352,
  #     37,                                                 42,
        24157817, 39088169, 63245986, 102334155, 165580141, 267914296,
  );
#433494437
#701408733
#1134903170
#1836311903
#2971215073
#4807526976
#7778742049
#12586269025
  for (my $i = 0; $i < @F; $i++)
    {
    $F[$i] = Math::BigInt->new($F[$i]);
    }
  }

sub _fibonacci_fast
  {
  my $x = shift; $x = $x->numify() if ref($x);
  return $F[$x] if $x < @F;

  # Knuth, TAOCR Vol 1, Third Edition, p. 81
  # F(n+m) = Fm * Fn+1 + Fm-1 * Fn

  # if m is set to n+1, we get:
  # F(n+n+1) = F(n+1) * Fn+1 + Fn * Fn
  # F(n*2+1) = F(n+1) ^ 2 + Fn ^ 2

  # so to know Fx, we must know F((x-1)/2), which only works for odd x
  # Fortunately:
  # Fx+1 = F(x) + F(x-1)
  # when x is even, then are x+1 and x-1 odd and can be calculated by the
  # same means, and from this we get Fx.

  # starting with level 0 at Fn we fill a hash with the different n we need
  # to calculate all Fn of the previous level. Here is an example for F1000:

  # To calculate F1000, we need F999 and F1001 (since 1000 is even)
  # To calculate F999, we need F((999-1)/2) and F((999+1)/2), this are 499
  # and 500. Sine for F1001 we need 500 and 501, we need F(500), F(501) and
  # F(499) to calculate F(1001) and F(999), and from these we can get F(1000).
  # For 500, we need 499 and 501, both are already needed.
  # For 501, we need 250 and 251. An so on and on until all values at a level
  # are under 17.
  # For the deepest level we use a table-lookup. The other levels are then
  # calulated backwards, until we arive at the top and the result is then in
  # level 0.

  # level
  #   0        1         2           3    and so on
  # 1000 ->   999   ->  499 <-  ->  249
  #    |	|---->  500  |
  #    |-->  1001   ->  501 <-  ->  250
  #                       |------>  251

  my @fibo;
  $fibo[0]->{$x} = undef;		# mark our final result as needed
  # if $x is even we need these two, too
  if ($x % 1 == 0)
    {
    $fibo[0]->{$x-1} = undef; $fibo[0]->{$x+1} = undef;
    }
  # XXX
  # for statistics
  my $steps = 0; my $sum = 0; my $add = 0; my $mul = 0;
  my $level = 0;
  my $high = 1;				# keep going?
  my ($t,$t1,$f);			# helper variables
  while ($high > 0)
    {
    $level++;				# next level
    $high = 0;				# count of results > @F
#      print "at level $level (high=$high)\n";
    foreach $f (keys %{$fibo[$level-1]})
      {
      $steps ++;
      if (($f & 1) == 0)		# odd/even?
        {
        # if it is even, add $f-1 and $f+1 to last level
        # if not existing in last level, we must add
        # ($f-1-1)/2 & ($f-1-1/2)+1 to the next level, too
	$t = $f-1;
        if (!exists $fibo[$level-1]->{$t})
          {
          $fibo[$level-1]->{$t} = undef; $t--; $t /= 2;			# $t is odd
          $fibo[$level]->{$t} = undef; $fibo[$level]->{$t+1} = undef;
          $high = 1 if $t+1 >= @F;	# any value not in table?
          }
	$t = $f+1;
        if (!exists $fibo[$level-1]->{$t})
          {
          $fibo[$level-1]->{$t} = undef; $t--; $t /= 2;			# $t is odd
          $fibo[$level]->{$t} = undef; $fibo[$level]->{$t+1} = undef;
          $high = 1 if $t+1 >= @F;	# any value not in table?
          }
#        print "$f even: ",$f-1," ",$f+1," in level ",$level-1,"\n";
        }
      else
        {
        # else add ($_-1)/2and ($_-1)/2 + 1 to this level
        $t = $f-1; $t /= 2;
        $fibo[$level]->{$t} = undef; $fibo[$level]->{$t+1} = undef;
        $high = 1 if $t+1 >= @F;	# any value not in table?
#       print "$_ odd: $t ",$t+1," in level $level (high = $high)\n";
        }
      }
    }
  # now we must fill our structure backwards with the results, combining them.
  # numbers in the last level can be looked up:
  foreach $f (keys %{$fibo[$level]})
    {
    # this inserts Math::BigInt objects, making the math below work:
    $fibo[$level]->{$f} = $F[$f];
    }
#  use Data::Dumper; print Dumper(\@fibo);
  my $l = $level;				# for statistics
  while ($level > 0)
    {
    $level--;
    #$sum += scalar keys %{$fibo[$level]};	# for statistics
    # first do the odd ones
    my $fibo_level = $fibo[$level];
    foreach $f (keys %$fibo_level)
      {
      next if ($f & 1) == 0;
      $t = $f-1; $t /= 2; my $t1 = $t+1;
      $t = $fibo[$level+1]->{$t};
      $t1 = $fibo[$level+1]->{$t1};
      #$fibo_level->{$f} = $t*$t+$t1*$t1;
      $fibo_level->{$f} = $t->copy->bmul($t)->badd($t1->copy()->bmul($t1));
      #$mul += 2; $add ++;			# for statistics
      }
    # now the even ones
    foreach $f (keys %{$fibo[$level]})
      {
      next if ($f & 1) != 0;
      $fibo[$level]->{$f} = $fibo[$level]->{$f+1} - $fibo[$level]->{$f-1};
      #$add ++;					# for statistics
      }
    }
#  print "sum $sum level $l => ",$sum/$l," steps $steps adds $add muls $mul\n";
  $fibo[0]->{$x};
  }

sub base
  {
  my ($number,$base) = @_;

  $number = Math::BigInt->new($number) unless ref $number;
  $base = Math::BigInt->new($base) unless ref $base;

  return if $number < $base;
  my $n = Math::BigInt->new(0);
  my $trial = $base;
  # 9 = 2**3 + 1
  while ($trial < $number)
    {
    $trial *= $base; $n++;
    }
  $trial /= $base; $a = $number - $trial;
  ($n,$a);
  }

sub to_base
  {
  # after an idea by Tilghman Lesher
  my ($x, $base, $alphabet) = @_;

  $x = Math::BigInt->new($x) unless ref $x;

  return '0' if $x->is_zero();

  # setup defaults:
  $base = 2 unless defined $base;
  my @digits = $alphabet ? split //, $alphabet : ('0' .. '9', 'A' .. 'Z');

  if ($base > scalar(@digits))
    {
    require Carp;
    Carp::carp("Base $base higher base than number of digits (" . scalar @digits . ") in alphabet");
    }

  if (!$x->is_pos())
    {
    require Carp;
    Carp::carp("to_base() needs a positive number");
    }

  my $o = $x->copy();
  my $r;

  my $result = '';
  while (!$o->is_zero)
    {
    ($o, $r) = $o->bdiv($base);
    $result = $digits[$r] . $result;
    }

  $result;
  }

sub hailstone
  {
  # return in list context the hailstone sequence, in scalar context the
  # number of steps to reach 1
  my ($n) = @_;

  $n = Math::BigInt->new($n) unless ref $n;

  return if $n->is_nan() || $n->is_negative();

  # Use the Math::BigInt lib directly for more speed, since all numbers
  # involved are positive integers.

  my $lib = Math::BigInt->config()->{lib};
  $n = $n->{value};
  my $three_ = $three->{value};
  my $two_ = $two->{value};

  if (wantarray)
    {
    my @seq;
    while (! $lib->_is_one($n))
      {
      # push @seq, Math::BigInt->new( $lib->_str($n) );
      push @seq, bless { value => $lib->_copy($n), sign => '+' }, "Math::BigInt";

      # was: ($n->is_odd()) ? ($n = $n * 3 + 1) : ($n = $n / 2);
      if ($lib->_is_odd($n))
        {
        $n = $lib->_mul ($n, $three_); $n = $lib->_inc ($n);

        # We now know that $n is at least 10 ( (3 * 3) + 1 ) because $n > 1
        # before we entered, and since $n was odd, it must have been at least
        # 3. So the next step is $n /= 2:
        push @seq, bless { value => $lib->_copy($n), sign => '+' }, "Math::BigInt";
        # this is better, but slower:
        #push @seq, Math::BigInt->new( $lib->_str($n) );
        # next step is $n /= 2 as usual (we save the else {} block, too)
        }
      $n = $lib->_div($n, $two_);
      }
    push @seq, Math::BigInt->bone();
    return @seq;
    }

  my $i = 1;
  while (! $lib->_is_one($n))
    {
    $i++;
    # was: ($n->is_odd()) ? ($n = $n * 3 + 1) : ($n = $n / 2);
    if ($lib->_is_odd($n))
      {
      $n = $lib->_mul ($n, $three_); $n = $lib->_inc ($n);

      # We now know that $n is at least 10 ( (3 * 3) + 1 ) because $n > 1
      # before we entered, and since $n was odd, it must have been at least 3.
      # So the next step is $n /= 2 as usual (we save the else {} block, too).
      $i++;			# one more (we know that $n cannot be 1)
      }
    $n = $lib->_div($n, $two_);
    }
  Math::BigInt->new($i);
  }

sub factorial
  {
  # calculate n! - use Math::BigInt bfac() for speed
  my ($n) = shift;

  if (ref($n))
    {
    $n->copy()->bfac();
    }
  else
    {
    Math::BigInt->new($n)->bfac();
    }
  }

sub bernoulli
  {
  # returns the nth Bernoulli number. In scalar context as Math::BigFloat
  # fraction, in list context as two Math:BigFloats, which, if divided, give
  # the same result. The series runs this:
  # 1/6, 1/30, 1/42, 1/30, 5/66, 691/2730, etc

  # Since I do not have yet a way to compute this, I have a table of the
  # first 40. So bernoulli(41) will fail for now.

  my $n = shift;

  return if $n < 0;
  my @table_1 = ( 1,1, -1,2 );					# 0, 1
  my @table = (
                1,6, -1,30, 1,42, -1,30, 5,66, -691,2730,	# 2, 4,
                7,6, -3617,510, 43867,798,
		-174611,330,
                854513,138,
		'-236364091',2730,
		'8553103',6,
                '-23749461029',870,
                '8615841276005',14322,
		'-7709321041217',510,
		'2577687858367',6,
		'-26315271553053477373',1919190,
		'2929993913841559',6,
		'-261082718496449122051',13530,			# 40
              );
  my ($a,$b);
  if ($n < 2)
    {
    $a = Math::BigFloat->new($table_1[$n*2]);
    $b = Math::BigFloat->new($table_1[$n*2+1]);
    }
  # n is odd:
  elsif (($n & 1) == 1)
    {
    $a = Math::BigFloat->bzero();
    $b = Math::BigFloat->bone();
    }
  elsif ($n <= 40)
    {
    $n -= 2;
    $a = Math::BigFloat->new($table[$n]);
    $b = Math::BigFloat->new($table[$n+1]);
    }
  else
    {
    die 'Bernoulli numbers over 40 not yet implemented.' if $n > 40;
    }
  wantarray ? ($a,$b): $a/$b;
  }

sub euler
  {
  # Calculate Euler's number.
  # first argument is x, so that result is e ** x
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = $_[0];
  $x = Math::BigFloat->new($x) if !ref($x) || (!$x->isa('Math::BigFloat'));

  $x->bexp($_[1]);
  }

sub sin
  {
  # calculate sinus
  # first argument is x, so that result is sin(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:      x^3   x^5   x^7   x^9
  #    sin = x - --- + --- - --- + --- ...
  # 		  3!    5!    7!    9!

  # difference for each term is thus x^2 and 1,2

  my $sin = $x->copy(); my $last = 0;
  my $sign = 1;				# start with -=
  my $x2 = $x * $x; 			# X ^ 2, difference between terms
  my $over = $x2 * $x; 			# X ^ 3
  my $below = Math::BigFloat->new(6); my $factorial = Math::BigFloat->new(4);
  while ($sin->bcmp($last) != 0) # no $x-$last > $diff because bdiv() limit on accuracy
    {
    $last = $sin->copy();
    if ($sign == 0)
      {
      $sin += $over->copy()->bdiv($below,$d);
      }
    else
      {
      $sin -= $over->copy()->bdiv($below,$d);
      }
    $sign = 1-$sign;					# alternate
    $over *= $x2;					# $x*$x
    $below *= $factorial; $factorial++;			# n*(n+1)
    $below *= $factorial; $factorial++;
    }
  $sin->bround($d-1);
  }

sub cos
  {
  # calculate cosinus
  # first argument is x, so that result is cos(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:      x^2   x^4   x^6   x^8
  #    cos = 1 - --- + --- - --- + --- ...
  # 		  2!    4!    6!    8!

  # difference for each term is thus x^2 and 1,2

  my $cos = Math::BigFloat->bone(); my $last = 0;
  my $over = $x * $x;			# X ^ 2
  my $x2 = $over->copy();		# X ^ 2; difference between terms
  my $sign = 1;				# start with -=
  my $below = Math::BigFloat->new(2); my $factorial = Math::BigFloat->new(3);
  while ($cos->bcmp($last) != 0) # no $x-$last > $diff because bdiv() limit on accuracy
    {
    $last = $cos->copy();
    if ($sign == 0)
      {
      $cos += $over->copy()->bdiv($below,$d);
      }
    else
      {
      $cos -= $over->copy()->bdiv($below,$d);
      }
    $sign = 1-$sign;					# alternate
    $over *= $x2;					# $x*$x
    $below *= $factorial; $factorial++;			# n*(n+1)
    $below *= $factorial; $factorial++;
    }
  $cos->round($d-1);
  }

sub tan
  {
  # calculate tangens
  # first argument is x, so that result is tan(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:  1         2            3            4           5

  #		      x^3          x^5          x^7          x^9
  #    tan = x + 1 * -----  + 2 * ----- + 17 * ----- + 62 * ----- ...
  # 		       3           15           315         2835
  #
  #  2^2n * ( 2^2n - 1) * Bn * x^(2n-1)          256*255 * 1 * x^7   17
  #  ---------------------------------- : n=4:  ----------------- = --- * x^7
  #               (2n)!                            40320 * 30       315
  #
  # 8! = 40320, B4 (Bernoully number 4) = 1/30

  # for each term we need: 2^2n, but if we have 2^2(n-1) we use n = (n-1)*2
  # 2 copy, 7 bmul, 2 bdiv, 3 badd, 1 bernoulli

  my $tan = $x->copy(); my $last = 0;
  my $x2 = $x*$x;
  my $over = $x2*$x;
  my $below = Math::BigFloat->new(24);	 	# (1*2*3*4) (2n)!
  my $factorial = Math::BigFloat->new(5);	# for next (2n)!
  my $two_n = Math::BigFloat->new(16);	 	# 2^2n
  my $two_factor = Math::BigFloat->new(4); 	# 2^2(n+1) = $two_n * $two_factor
  my ($b,$b1,$b2); $b = 4;
  while ($tan->bcmp($last) != 0) # no $x-$last > $diff because bdiv() limit on accuracy
    {
    $last = $tan->copy();
    ($b1,$b2) = bernoulli($b);
    $tan += $over->copy()->bmul($two_n)->bmul($two_n - $fone)->bmul($b1->babs())->bdiv($below,$d)->bdiv($b2,$d);
    $over *= $x2;				# x^3, x^5 etc
    $below *= $factorial; $factorial++;		# n*(n+1)
    $below *= $factorial; $factorial++;
    $two_n *= $two_factor;			# 2^2(n+1) = 2^2n * 4
    $b += 2;					# next bernoulli index
    last if $b > 40;				# safeguard
    }
  $tan->round($d-1);
  }

sub sinh
  {
  # calculate sinus hyperbolicus
  # first argument is x, so that result is sinh(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:       x^3   x^5   x^7
  #    sinh = x + --- + --- + --- ...
  # 	           3!    5!    7!

  # difference for each term is thus x^2 and 1,2

  my $sinh = $x->copy(); my $last = 0;
  my $x2 = $x*$x;
  my $over = $x2 * $x; my $below = Math::BigFloat->new(6); my $factorial = Math::BigFloat->new(4);
  while ($sinh->bcmp($last)) # no $x-$last > $diff because bdiv() limit on accuracy
    {
    $last = $sinh->copy();
    $sinh += $over->copy()->bdiv($below,$d);
    $over *= $x2;					# $x*$x
    $below *= $factorial; $factorial++;			# n*(n+1)
    $below *= $factorial; $factorial++;
    }
  $sinh->bround($d-1);
  }

sub cosh
  {
  # calculate cosinus hyperbolicus
  # first argument is x, so that result is cosh(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:       x^2   x^4   x^6
  #    cosh = x + --- + --- + --- ...
  # 	           2!    4!    6!

  # difference for each term is thus x^2 and 1,2

  my $cosh = Math::BigFloat->bone(); my $last = 0;
  my $x2 = $x*$x;
  my $over = $x2; my $below = Math::BigFloat->new(); my $factorial = Math::BigFloat->new(3);
  while ($cosh->bcmp($last)) # no $x-$last > $diff because bdiv() limit on accuracy
    {
    $last = $cosh->copy();
    $cosh += $over->copy()->bdiv($below,$d);
    $over *= $x2;					# $x*$x
    $below *= $factorial; $factorial++;			# n*(n+1)
    $below *= $factorial; $factorial++;
    }
  $cosh->bround($d-1);
  }

sub arctan
  {
  # calculate arcus tangens
  # first argument is x, so that result is arctan(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:      x^3   x^5   x^7   x^9
  # arctan = x - --- + --- - --- + --- ...
  # 		  3     5    7      9

  # difference for each term is thus x^2 and 2:
  # 2 copy, 1 bmul, 1 badd, 1 bdiv

  my $arctan = $x->copy(); my $last = 0;
  my $x2 = $x*$x;
  my $over = $x2*$x; my $below = Math::BigFloat->new(3); my $add = Math::BigFloat->new(2);
  my $sign = 1;
  while ($arctan->bcmp($last)) # no $x-$last > $diff because bdiv() limit on A
    {
    $last = $arctan->copy();
    if ($sign == 0)
      {
      $arctan += $over->copy()->bdiv($below,$d);
      }
    else
      {
      $arctan -= $over->copy()->bdiv($below,$d);
      }
    $sign = 1-$sign;					# alternate
    $over *= $x2;					# $x*$x
    $below += $add;
    }
  $arctan->bround($d-1);
  }

sub arctanh
  {
  # calculate arcus tangens hyperbolicus
  # first argument is x, so that result is arctanh(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:       x^3   x^5   x^7   x^9
  # arctanh = x + --- + --- + --- + --- + ...
  # 	 	   3     5    7      9

  # difference for each term is thus x^2 and 2:
  # 2 copy, 1 bmul, 1 badd, 1 bdiv

  my $arctanh = $x->copy(); my $last = 0;
  my $x2 = $x*$x;
  my $over = $x2*$x; my $below = Math::BigFloat->new(3); my $add = Math::BigFloat->new(2);
  while ($arctanh->bcmp($last)) # no $x-$last > $diff because bdiv() limit on A
    {
    $last = $arctanh->copy();
    $arctanh += $over->copy()->bdiv($below,$d);
    $over *= $x2;					# $x*$x
    $below += $add;
    }
  $arctanh->bround($d-1);
  }

sub arcsin
  {
  # calculate arcus sinus
  # first argument is x, so that result is arcsin(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:      1 * x^3   1 * 3 * x^5   1 * 3 * 5 * x^7
  # arcsin = x + ------- + ----------- + --------------- + ...
  # 		 2 *  3    2 * 4 *  5    2 * 4 * 6 *   7

  # difference for each term is thus x^2 and two muls (fac1, fac2):
  # 3 copy, 3 bmul, 1 bdiv, 3 badd

  my $arcsin = $x->copy(); my $last = 0;
  my $x2 = $x*$x;
  my $over = $x2*$x; my $below = Math::BigFloat->new(6);
  my $fac1 = Math::BigFloat->new(1);
  my $fac2 = Math::BigFloat->new(2);
  my $two = Math::BigFloat->new(2);
  while ($arcsin->bcmp($last)) # no $x-$last > $diff because bdiv() limit on A
    {
    $last = $arcsin->copy();
    $arcsin += $over->copy()->bmul($fac1)->bdiv($below->copy->bmul($fac2),$d);
    $over *= $x2;					# $x*$x
    $below += $one;
    $fac1 += $two;
    $fac2 += $two;
    }
  $arcsin->bround($d-1);
  }

sub arcsinh
  {
  # calculate arcus sinus hyperbolicus
  # first argument is x, so that result is arcsinh(x)
  # Second argument is accuracy (number of significant digits), it
  # stops when at least so much plus one digits are 'stable' and then
  # rounds it. Default is 42.
  my $x = shift; $x = 0 if !defined $x;
  my $d = abs(shift || 42); $d = abs($d)+1;

  $x = Math::BigFloat->new($x) if ref($x) ne 'Math::BigFloat';

  # taylor:      1 * x^3   1 * 3 * x^5   1 * 3 * 5 * x^7
  # arcsin = x - ------- + ----------- - --------------- + ...
  # 		 2 *  3    2 * 4 *  5    2 * 4 * 6 *   7

  # difference for each term is thus x^2 and two muls (fac1, fac2):
  # 3 copy, 3 bmul, 1 bdiv, 3 badd

  my $arcsinh = $x->copy(); my $last = 0;
  my $x2 = $x*$x; my $sign = 0;
  my $over = $x2*$x; my $below = 6;
  my $fac1 = Math::BigInt->new(1);
  my $fac2 = Math::BigInt->new(2);
  while ($arcsinh ne $last) # no $x-$last > $diff because bdiv() limit on A
    {
    $last = $arcsinh->copy();
    if ($sign == 0)
      {
      $arcsinh += $over->copy()->bmul(
        $fac1)->bdiv($below->copy->bmul($fac2),$d);
      }
    else
      {
      $arcsinh -= $over->copy()->bmul(
        $fac1)->bdiv($below->copy->bmul($fac2),$d);
      }
    $over *= $x2;					# $x*$x
    $below += $one;
    $fac1 += $two;
    $fac2 += $two;
    }
  $arcsinh->round($d-1);
  }

sub log
  {
  my ($x,$base,$d) = @_;

  my $y;
  if (!ref($x) || !$x->isa('Math::BigFloat'))
    {
    $y = Math::BigFloat->new($x);
    }
  else
    {
    $y = $x->copy();
    }
  $y->blog($base,$d);
  $y;
  }

sub pi
  {
  # calculate PI (as suggested by Robert Creager)
  my $digits = abs(shift || 1024);

  my $d = $digits+5;

  my $pi =  $sixteen * arctan( scalar $fone->copy()->bdiv($five,$d), $d )
             - $four * arctan( scalar $fone->copy()->bdiv($twothreenine,$d), $d);
  $pi->bround($digits+1);	# +1 for the "3."
  }

1;

__END__