This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.23
in repository libmath-prime-util-perl.

commit 9393ea90eaf6313e4e8e15d4cb4e46d190fdba51
Author: Dana Jacobsen <d...@acm.org>
Date:   Wed Feb 27 19:34:18 2013 -0800

    Add consec int lcm, documentation changes
---
 Changes                |   4 ++
 TODO                   |   1 +
 lib/Math/Prime/Util.pm | 159 ++++++++++++++++++++++++++++++++++++-------------
 3 files changed, 122 insertions(+), 42 deletions(-)

diff --git a/Changes b/Changes
index c2f25eb..60d17ca 100644
--- a/Changes
+++ b/Changes
@@ -11,6 +11,10 @@ Revision history for Perl extension Math::Prime::Util.
     - PP Zeta code replaced (for no-MPFR, non-bignums) with new series.  The
       new code is much more accurate for small values, and *much* faster.
 
+    - Add consecutive_integer_lcm function, just like MPU::GMP's.
+
+    - Add examples for first and second Chebyshev functions.
+
 0.22 26 February 2013
 
     - Move main factor loop out of xs and into factor.c.
diff --git a/TODO b/TODO
index 0acb07b..31fa191 100644
--- a/TODO
+++ b/TODO
@@ -38,6 +38,7 @@
 - More efficient Mertens.  The current version has poor growth.
 
 - Add first and second Chebyshev functions.  See Planat and Solé (2011).
+  MPU::GMP has primorial and consecutive_integer_lcm which are what we need.
 
 - Keep speeding up factoring for very large values.
 
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 7e76610..4c494e9 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -25,7 +25,7 @@ our @EXPORT_OK = qw(
                      nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
                      random_prime random_ndigit_prime random_nbit_prime
                      random_strong_prime random_maurer_prime
-                     primorial pn_primorial
+                     primorial pn_primorial consecutive_integer_lcm
                      factor all_factors
                      moebius mertens euler_phi jordan_totient exp_mangoldt
                      divisor_sum
@@ -435,23 +435,25 @@ sub primes {
 # It will generate primes with more bits, but it slows down a lot.  The
 # time variation becomes quite extreme once bit sizes get over 6000 or so.
 #
-# Random timings for 1M calls:
-#   0.032   system rand
-#   0.20    Math::Random::MT::Auto
-#   0.96    Math::Random::Secure  irand  (with Math::Random::ISAAC::XS)
-#   1.71    Math::Random::Secure         (with Math::Random::ISAAC::XS)
-#   1.90   *Bytes::Random::Secure irand  (with Math::Random::ISAAC::XS)
-#   2.40    Math::Random::Secure  irand
-#   3.37    Math::Random::Secure
-#   2.21   *Bytes::Random::Secure        (with Math::Random::ISAAC::XS)
-#   3.32   *Bytes::Random::Secure irand
-#   3.62   *Bytes::Random::Secure
-# 123.8    *Crypt::Random         irand
-# BRS: sub irand { return unpack("L", random_bytes(4)); }
-#      sub rand  {return ($_[0]||1)*(unpack("L",random_bytes(4))/4294967296.0)}
-# CR:  makerandom(Size=>32, Uniform=>1, Strength=>0)
-#      haveged daemon running to stop /dev/random blocking
-# Both BRS and CR have more features that this isn't measuring.
+# Random timings for 10M calls:
+#    1.92    system rand
+#    2.62    Math::Random::MT::Auto
+#   12.0     Math::Random::Secure           w/ISAAC::XS
+#   12.6     Bytes::Random::Secure OO       w/ISAAC::XS
+#   31.1     Bytes::Random::Secure OO
+#   44.5     Bytes::Random::Secure function w/ISAAC::XS
+#   44.8     Math::Random::Secure
+#   71.5     Bytes::Random::Secure function
+# 1840       Crypt::Random
+#
+# time perl -E 'sub irand {int(rand(4294967296));} irand() for 1..10000000;'
+# time perl -E 'use Math::Random::MT::Auto qw/irand/; irand() for 1..10000000;'
+# time perl -E 'use Math::Random::Secure qw/irand/; irand() for 1..10000000;'
+# time perl -E 'use Bytes::Random::Secure qw/random_bytes/; sub irand {return 
unpack("L",random_bytes(4));} irand() for 1..10000000;'
+# time perl -E 'use Bytes::Random::Secure; my $rng = 
Bytes::Random::Secure->new(); sub irand {return $rng->irand;} irand() for 
1..10000000;'
+# time perl -E 'use Crypt::Random qw/makerandom/; sub irand 
{makerandom(Size=>32, Uniform=>1, Strength=>0)} irand() for 1..100_000;'
+# > haveged daemon running to stop /dev/random blocking
+# > Both BRS and CR have more features that this isn't measuring.
 #
 # To verify distribution:
 #   perl -Iblib/lib -Iblib/arch -MMath::Prime::Util=:all -E 'my %freq; 
$n=1000000; $freq{random_nbit_prime(6)}++ for (1..$n); printf("%4d %6.3f%%\n", 
$_, 100.0*$freq{$_}/$n) for sort {$a<=>$b} keys %freq;'
@@ -1033,6 +1035,44 @@ sub pn_primorial {
   return primorial(nth_prime($n));
 }
 
+sub consecutive_integer_lcm {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return 1 if $n <= 1;
+
+  my $pn = 1;
+  if ($n >= (($_Config{'maxbits'} == 32) ? 22 : 46)) {
+    if (!defined $Math::BigInt::VERSION) {
+      eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
+      or do { croak "Cannot load Math::BigInt"; };
+    }
+    $pn = Math::BigInt->bone();
+  }
+  # Ensure we use their type
+  $pn = $_[0]->copy->bone() if ref($_[0]) eq 'Math::BigInt';
+
+  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::consecutive_integer_lcm) {
+    my $clcm = Math::Prime::Util::GMP::consecutive_integer_lcm($n);
+    return int($clcm) unless ref($pn) eq 'Math::BigInt';
+    return $pn->copy->badd("$clcm");
+  }
+
+  my @primes = @{primes(2,$n)};
+  while (@primes) {
+    my $p = $primes[0];
+    my $pmin = int($n/$p);
+    last if $p > $pmin;
+    my $p_power = $p*$p;
+    $p_power *= $p while $p_power <= $pmin;
+    $pn *= $p_power;
+    shift @primes;
+  }
+  foreach my $p (@primes) {
+    $pn *= $p;
+  }
+  return $pn;
+}
+
 
 sub all_factors {
   my $n = shift;
@@ -1994,15 +2034,16 @@ Version 0.22
   $sigma  = divisor_sum( $n );
   $sigma2 = divisor_sum( $n, sub { $_[0]*$_[0] } );
 
-  # The primorial n# (product of all primes <= n)
-  say "15# (2*3*5*7*11*13) is ", primorial(15);
-  # The primorial p(n)# (product of first n primes)
-  say "P(9)# (2*3*5*7*11*13*17*19*23) is ", pn_primorial(9);
+  # primorial n#, primorial p(n)#, and lcm
+  say "The product of primes below 47 is ",     primorial(47);
+  say "The product of the first 47 primes is ", pn_primorial(47);
+  say "lcm(1..1000) is ", consecutive_integer_lcm(1000);
 
   # Ei, li, and Riemann R functions
-  my $ei = ExponentialIntegral($x);    # $x a real: $x != 0
-  my $li = LogarithmicIntegral($x);    # $x a real: $x >= 0
-  my $R  = RiemannR($x)                # $x a real: $x > 0
+  my $ei   = ExponentialIntegral($x);   # $x a real: $x != 0
+  my $li   = LogarithmicIntegral($x);   # $x a real: $x >= 0
+  my $R    = RiemannR($x)               # $x a real: $x > 0
+  my $Zeta = RiemannZeta($x)            # $x a real: $x >= 0
 
 
   # Precalculate a sieve, possibly speeding up later work.
@@ -2063,7 +2104,7 @@ L<Math::Prime::Util::GMP> if you plan to use bigints with 
this module, as
 it will make it run much faster.
 
 
-Some of the functions, notably:
+Some of the functions, including:
 
   factor
   is_prime
@@ -2072,6 +2113,10 @@ Some of the functions, notably:
   next_prime
   prev_prime
   nth_prime
+  moebius
+  mertens
+  euler_phi
+  exp_mangoldt
 
 work very fast (under 1 microsecond) on small inputs, but the wrappers for
 input validation and bigint support take more time than the function itself.
@@ -2569,6 +2614,16 @@ The result will be a L<Math::BigInt> object if it is 
larger than the native
 bit size.
 
 
+=head2 consecutive_integer_lcm
+
+  $lcm = consecutive_integer_lcm($n);
+
+Given an unsigned integer argument, returns the least common multiple of all
+integers from 1 to C<n>.  This can be done by manipulation of the primes up
+to C<n>, resulting in much faster and memory-friendly results than using
+a factorial.
+
+
 =head2 random_prime
 
   my $small_prime = random_prime(1000);      # random prime <= limit
@@ -2605,15 +2660,19 @@ a non-blocking seed.  This will create good quality 
random numbers, so there
 should be little reason to change unless one is generating long-term keys,
 where using the blocking random source may be preferred.
 
-Examples of irand functions:
+Examples of various ways to set your own irand function:
 
   # Math::Random::Secure.  Uses ISAAC and strong seed methods.
   use Math::Random::Secure;
   prime_set_config(irand => \&Math::Random::Secure::irand);
 
-  # Bytes::Random::Secure.  Also uses ISAAC and strong seed methods.
-  use Bytes::Random::Secure qw/random_bytes/;
-  prime_set_config(irand => sub { return unpack("L", random_bytes(4)); });
+  # Bytes::Random::Secure (OO interface with full control of options):
+  use Bytes::Random::Secure ();
+  BEGIN {
+    my $rng = Bytes::Random::Secure->new( Bits => 512 );
+    sub irand { return $rng->irand; }
+  }
+  prime_set_config(irand => \&irand);
 
   # Crypt::Random.  Uses Pari and /dev/random.  Very slow.
   use Crypt::Random qw/makerandom/;
@@ -3008,7 +3067,7 @@ the input is less than 1, results will be calculated as 
C<Ei(ln x)>.
 
   my $z = RiemannZeta($s);
 
-Given a floating point input C<s> where C<s E<gt> 0>, returns the floating
+Given a floating point input C<s> where C<s E<gt>= 0>, returns the floating
 point value of ζ(s)-1, where ζ(s) is the Riemann zeta function.  One is
 subtracted to ensure maximum precision for large values of C<s>.  The zeta
 function is the sum from k=1 to infinity of C<1 / k^s>.  This function only
@@ -3019,7 +3078,8 @@ were Math::BigFloat objects.
 
 For non-BigInt/BigFloat objects, the result should be accurate to at least 14
 digits.  The XS code uses a rational Chebyshev approximation between 0.5 and 5,
-and a series for larger values.
+and a series for other values.  The PP code uses an identical series for all
+values.
 
 For BigInt / BigFloat objects, we first check to see if the Math::MPFR module
 is installed.  If so, then it is used, as it will return results much faster
@@ -3068,6 +3128,19 @@ Print some primes above 64-bit range:
     # Similar code using Pari:
     # perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI 
"100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while 
($p <= $end) { say $p; $p = nextprime($p+1); }'
 
+The first Chebyshev function θ(x):
+
+  use List::Util qw/sum/;
+  sub chebyshev1 { 0+sum( map { log($_) } @{primes(shift)} ) }
+  #  or using primorial and BigFloat:
+  sub chebyshev1a { Math::BigFloat->new(primorial(shift))->blog }
+
+The second Chebyshev function ψ(x):
+
+  use List::Util qw/sum/;
+  sub chebyshev2 { 0+sum( map { log(exp_mangoldt($_)) } 1 .. shift ) }
+  #  or using cilcm and BigFloat:
+  sub chebyshev2a { Math::BigFloat->new(consecutive_integer_lcm(shift))->blog }
 
 Project Euler, problem 3 (Largest prime factor):
 
@@ -3151,17 +3224,19 @@ Here's the best way for PE187.  Under 30 milliseconds 
from the command line:
 I have not completed testing all the functions near the word size limit
 (e.g. C<2^32> for 32-bit machines).  Please report any problems you find.
 
-Perl versions earlier than 5.8.0 have issues with 64-bit that show up in the
-factoring tests.  The test suite will try to determine if your Perl is broken.
-If you use later versions of Perl, or Perl 5.6.2 32-bit, or Perl 5.6.2 64-bit
-and keep numbers below C<~ 2^52>, then everything works.  The best solution is
-to update to a more recent Perl.
+Perl versions earlier than 5.8.0 have a rather broken 64-bit implementation,
+in that the values are actually stored as doubles.  Hence any value larger
+than C<~ 2^49> will start losing bottom bits.  This causes numerous functions
+to not work properly.  The test suite will try to determine if your Perl is
+broken (this only applies to really old versions of Perl compiled for 64-bit
+when using numbers larger than C<~ 2^49>).  The best solution is updating to
+a more recent Perl.
 
 The module is thread-safe and should allow good concurrency on all platforms
-that support Perl threads except Win32 (Cygwin works).  With Win32, either
-don't use threads or make sure C<prime_precalc> is called before using
-C<primes>, C<prime_count>, or C<nth_prime> with large inputs.  This is B<only>
-an issue if you use non-Cygwin Win32 and call these routines from within
+that support Perl threads except Win32.  With Win32, either don't use threads
+or make sure C<prime_precalc> is called before using C<primes>,
+C<prime_count>, or C<nth_prime> with large inputs.  This is B<only>
+an issue if you use non-Cygwin Win32 B<and> call these routines from within
 Perl threads.
 
 
@@ -3203,7 +3278,7 @@ Perl modules, counting the primes to C<800_000_000> (800 
million):
        0.47  Math::Prime::Util::PP       0.14     Perl (Lehmer's method)
        0.9   Math::Prime::Util           0.01     mod-30 sieve
        2.9   Math::Prime::FastSieve      0.12     decent odd-number sieve
-      11.7   Math::Prime::XS             0.29     "" but needs a count API
+      11.7   Math::Prime::XS             0.26     needs some optimization
       15.0   Bit::Vector                 7.2
       48.9   Math::Prime::Util::PP       0.14     Perl (fastest I know of)
      170.0   Faster Perl sieve (net)     2012-01  array of odds

-- 
Alioth's /usr/local/bin/git-commit-notice on 
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to