Math-AnyNum
view release on metacpan or search on metacpan
[ADDITIONS]
- acmp(x, y): absolute comparison of `x` and `y`.
- polygonal(n, k): returns the nth k-gonal number.
- polygonal_root(n, k): returns the k-gonal root of `n`.
- polygonal_root2(n, k): returns the second k-gonal root of `n`.
- ipolygonal_root(n, k): returns the integer k-gonal root of `n`.
- ipolygonal_root2(n, k): returns the second integer k-gonal root of `n`.
- is_polygonal(n, k): returns 1 when `n` is a k-gonal number.
- is_polygonal2(n, k): returns 1 when `n` is a second k-gonal number.
- faulhaber_sum(n, p): computes 1^p + 2^p + 3^p + ... + n^p, using Faulhaber's formula.
0.12 2017-09-18
[ADDITIONS]
- Added the `rat_approx(n)` function, which returns the smallest rational approximation for a given real number `n`.
[IMPROVEMENTS]
- The newly added functions in Math::MPFR-3.36, Rmpfr_q_div() and Rmpfr_z_div(), are now used by Math::AnyNum.
Changes
examples/agm_pi.pl
examples/arithmetic_coding.pl
examples/bernoulli_numbers_from_primes.pl
examples/bernoulli_numbers_recursive.pl
examples/bernoulli_seidel.pl
examples/binary_arithmetic_coding.pl
examples/binradix_arithmetic_coding.pl
examples/BPSW_primality_test.pl
examples/computing_pi.pl
examples/faulhaber_s_formula.pl
examples/fibonacci.pl
examples/fibonacci_validation.pl
examples/halley_s_method.pl
examples/inverse_of_factorial.pl
examples/inverse_of_fibonacci.pl
examples/is_power.pl
examples/krzysztof_reformulated_zeta_function.pl
examples/lambert_W.pl
examples/mandelbrot_set.pl
examples/miller_rabin_primality_test.pl
examples/newton_s_method.pl
examples/partial_sums_of_sigma_function.pl
examples/pi_machin.pl
examples/PSW_primality_test.pl
examples/rsa_algorithm.pl
examples/solve_pell_equation.pl
examples/tac-compressor.pl
SHA256 cad84790401acdebfda5284a92eb4d9bc90713580a99157fbf4ed9fd4b2cf117 examples/BPSW_primality_test.pl
SHA256 d5d7de9e034fc39bd68bedba6613bc471b1f8261cad65462dd9aa700b25d2281 examples/PSW_primality_test.pl
SHA256 fbec24c8676d4bc209671250e728a2af88167bd728666acc5949d3960669b25d examples/agm_pi.pl
SHA256 3aabd02bb5cc7433dc2d5fe54ff9be1a988101c0716c44cb1f6ad114f65d5cf4 examples/arithmetic_coding.pl
SHA256 c62077e52b27ae92221cabf2b2ec87a080f570a4b0620cac4c236c6ac9f30725 examples/bernoulli_numbers_from_primes.pl
SHA256 24520cd2f3aaa72f9a08b5d8c0d61873bfdd2f8483555aabf39d07de4d16712b examples/bernoulli_numbers_recursive.pl
SHA256 b3e1e01e348f5b3da1e4c99215d665945712fc7279078695ae58a275412a88b9 examples/bernoulli_seidel.pl
SHA256 5c952dcb0dd038f40e5e7a856e0a7e9b61ee0ff2060fe7a6900f8aa38afa2762 examples/binary_arithmetic_coding.pl
SHA256 6eb2cabd317574ed1cf5dba030d12e2502f5efe6bdb329de72df8f78206cc0c7 examples/binradix_arithmetic_coding.pl
SHA256 e0f2b1392386cea7d33d3228ab0a297e07ab52d73c6d5b19badaa5258092fd1f examples/computing_pi.pl
SHA256 f45437d321a5ebdabe509d6b527b7792290895f00a832666acea82c222c1136f examples/faulhaber_s_formula.pl
SHA256 111044be67e44022c41ffcc2db10e44fc7a0913ad3259f4edd09c41180d56437 examples/fibonacci.pl
SHA256 594ed1e376be79787f538cae08e794e711b27a61a490e25611691f660609c955 examples/fibonacci_validation.pl
SHA256 3628a9c81a7b80d3d155deeca96bccf804da471a9c61cd87a9b16429565b8ee8 examples/halley_s_method.pl
SHA256 5d0a51ac61f12e00019de06bda13999f98052ebf53f30e20881eba608aff771e examples/inverse_of_factorial.pl
SHA256 6a5799e4e70e4976b93d85e7f2f4761e3f543057d642fff803270d4c392ddd47 examples/inverse_of_fibonacci.pl
SHA256 d7929cef15a9ce122b5ea663f6ef61307d0e75f739e658221e368ff00ee06f0e examples/is_power.pl
SHA256 a5527a9482d2396b36fda3b397cbff46285a0e7e94d521c4a19fce1497907908 examples/krzysztof_reformulated_zeta_function.pl
SHA256 36e8e8e4b15ece8f380566f6786153752728ff54a9b814132cd0078d22e61563 examples/lambert_W.pl
SHA256 55c45179689f1297a3faaf9a64697b29be7042140e5110316bd0c8109ff75e8c examples/mandelbrot_set.pl
SHA256 ef57b7f63ab12f4fc430e53de87249fa135ca21491c075455c9cd4804ecaeb8a examples/miller_rabin_primality_test.pl
SHA256 928280c185647c6a0590a1860ee8ab75c5a7ae8d6dc65a6fad835d2e010ab95e examples/newton_s_method.pl
SHA256 8834273ee7ec15ce567284aea912b1d1e4148541a18491db47ab6ed3809a9364 examples/partial_sums_of_sigma_function.pl
SHA256 ecce3ea0f64463bc6e7e727fddf0a6d09c94b8cfbdb5c4c347453717b2f030b3 examples/pi_machin.pl
SHA256 55a514b7d591f48f3d523a2c672b03ef579bac5b69a55933d8bffa32562d9c15 examples/rsa_algorithm.pl
SHA256 fe055c47ffaea5aeac9ebe14709cea2082a9b00db5e495bbec22c95b18c0db71 examples/solve_pell_equation.pl
SHA256 c41413eaa426f5e0d3b688d83311caf1767921ef4faca58211bba49b160cf867 examples/tac-compressor.pl
SHA256 c8d9cfdbc321af860c95e0c49118bb1bea278f1c0557a1fa312bd673335947f4 examples/tribonacci.pl
examples/faulhaber_s_formula.pl view on Meta::CPAN
#!/usr/bin/perl
# The formula for calculating the sum of consecutive
# numbers raised to a given power, such as:
# 1^p + 2^p + 3^p + ... + n^p
# where p is a positive integer.
# See also:
# https://en.wikipedia.org/wiki/Faulhaber%27s_formula
use 5.020;
use strict;
use warnings;
use lib qw(../lib);
use experimental qw(signatures);
use Math::AnyNum qw(:overload binomial factorial faulhaber_sum);
# This function returns the n-th Bernoulli number
examples/faulhaber_s_formula.pl view on Meta::CPAN
foreach my $i (1 .. $n) {
foreach my $k (0 .. $i - 1) {
$B[$i] //= 0;
$B[$i] -= $B[$k] / factorial($i - $k + 1);
}
}
return $B[-1] * factorial($#B);
}
# The Faulhaber's formula
# See: https://en.wikipedia.org/wiki/Faulhaber%27s_formula
sub faulhaber_s_formula ($n, $p) {
my $sum = 0;
for my $k (0 .. $p) {
$sum += (-1)**$k * binomial($p + 1, $k) * bernoulli_number($k) * $n**($p - $k + 1);
}
return $sum / ($p + 1);
}
# Alternate expression using Bernoulli polynomials
# See: https://en.wikipedia.org/wiki/Faulhaber%27s_formula#Alternate_expressions
sub bernoulli_polynomials ($n, $x) {
my $sum = 0;
for my $k (0 .. $n) {
$sum += binomial($n, $k) * bernoulli_number($n - $k) * $x**$k;
}
return $sum;
}
sub faulhaber_s_formula_2 ($n, $p) {
(bernoulli_polynomials($p + 1, $n + 1) - bernoulli_number($p + 1)) / ($p + 1);
}
foreach my $m (1 .. 15) {
my $n = int(rand(100));
my $t1 = faulhaber_s_formula($n, $m);
my $t2 = faulhaber_s_formula_2($n, $m);
my $t3 = faulhaber_sum($n, $m);
say "Sum_{k=1..$n} k^$m = $t1";
die "error: $t1 != $t2" if ($t1 != $t2);
die "error: $t1 != $t3" if ($t1 != $t3);
}
__END__
Sum_{k=1..79} k^1 = 3160
examples/partial_sums_of_sigma_function.pl view on Meta::CPAN
#!/usr/bin/perl
# Author: Daniel "Trizen" Șuteu
# Date: 10 Novermber 2018
# https://github.com/trizen
# A generalized formula with O(sqrt(n)) complexity for computing the partial-sum of `k^m * sigma_j(k)`, for `1 <= k <= n`:
#
# Sum_{k=1..n} k^m * sigma_j(k)
#
# for any fixed integers m >= 0 and j >= 0.
# Formula:
# Sum_{k=1..n} k^m * sigma_j(k) = Sum_{k=1..floor(sqrt(n))} F(m, k) * (F(m+j, floor(n/k)) - F(m+j, floor(n/(k+1))))
# + Sum_{k=1..floor(n/(floor(sqrt(n))+1))} k^(m+j) * F(m, floor(n/k))
#
# where F(n,x) is Faulhaber's formula for `Sum_{k=1..x} k^n`, defined in terms of Bernoulli polynomials as:
#
# F(n, x) = (Bernoulli(n+1, x+1) - Bernoulli(n+1)) / (n+1)
#
# where Bernoulli(n,x) are the Bernoulli polynomials and Bernoulli(n) is the n-th Bernoulli number.
# See also:
# https://en.wikipedia.org/wiki/Divisor_function
# https://en.wikipedia.org/wiki/Faulhaber%27s_formula
# https://en.wikipedia.org/wiki/Bernoulli_polynomials
# https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html
use 5.020;
use strict;
use warnings;
use lib qw(../lib);
use experimental qw(signatures);
use Math::AnyNum qw(isqrt idiv ipow bernoulli faulhaber_sum dirichlet_sum);
sub faulhaber_partial_sum_of_sigma($n, $m, $j) { # using Faulhaber's formula
my $total = 0;
my $s = isqrt($n);
my $u = int($n / ($s + 1));
for my $k (1 .. $s) {
$total += faulhaber_sum($k, $m) * (
faulhaber_sum(int($n/$k), $m+$j)
- faulhaber_sum(int($n/($k+1)), $m+$j)
examples/zeta_2n.pl view on Meta::CPAN
#!/usr/bin/perl
# Calculate zeta(2n) using a closed-form formula.
# See also:
# https://en.wikipedia.org/wiki/Riemann_zeta_function
use 5.010;
use strict;
use warnings;
use lib qw(../lib);
use Memoize qw(memoize);
examples/zeta_2n_fast.pl view on Meta::CPAN
#!/usr/bin/perl
# Calculate zeta(2n) using a closed-form formula.
# See also:
# https://en.wikipedia.org/wiki/Riemann_zeta_function
use 5.010;
use strict;
use warnings;
use lib qw(../lib);
use experimental qw(signatures);
lib/Math/AnyNum.pm view on Meta::CPAN
}
Math::GMPz::Rmpz_fits_ulong_p($k)
? Math::GMPz::Rmpz_powm_ui($r, $n, Math::GMPz::Rmpz_get_ui($k), $m)
: Math::GMPz::Rmpz_powm($r, $n, $k, $m);
bless \$r;
}
#
## Geometric summation formula
## https://en.wikipedia.org/wiki/Geometric_series
#
sub geometric_sum ($$) {
my ($n, $r) = @_;
$n = ref($n) eq __PACKAGE__ ? $$n : _star2obj($n);
$r = ref($r) eq __PACKAGE__ ? $$r : _star2obj($r);
bless \__div__(__sub__(__pow__($r, __add__($n, 1)), 1), __sub__($r, 1));
}
#
## Faulhaber's summation formula
## https://en.wikipedia.org/wiki/Faulhaber%27s_formula
#
sub faulhaber_sum ($$) {
my ($n, $p) = @_;
$n = _star2mpz($n) // goto &nan;
$p = _star2ui($p) // goto &nan;
if ($p == 0) {
return bless \$n;
lib/Math/AnyNum.pod view on Meta::CPAN
=head2 kronecker
kronecker(n, m) #=> Scalar
Returns the Kronecker symbol I<(n|m)>, which is a generalization of the Jacobi symbol for all integers I<m>.
=head2 faulhaber_sum
faulhaber_sum(n, k) #=> Int | NaN
Computes the power sum C<1^k + 2^k + 3^k +...+ n^k>, using Faulhaber's formula.
The value for C<k> must be a non-negative integer. Returns NaN otherwise.
Example:
faulhaber_sum(5, 2) = 1^2 + 2^2 + 3^2 + 4^2 + 5^2 = 55
=head2 geometric_sum
geometric_sum(n, r) #=> Any
Computes the geometric sum C<1 + r + r^2 + r^3 + ... + r^n>, using the following formula:
geometric_sum(n, r) = (r^(n+1) - 1) / (r - 1)
Example:
geometric_sum(5, 8) = 8^0 + 8^1 + 8^2 + 8^3 + 8^4 + 8^5 = 37449
=head2 dirichlet_sum
dirichlet_sum(n, \&f, \&g, \&F, \&G) #=> Int | NaN
lib/Math/AnyNum.pod view on Meta::CPAN
sub { faulhaber_sum($_[0], 2) }, # G(n) = Sum_{k=1..n} g(k)
)
=head2 harmonic | harmfrac
harmonic(n) #=> Rat | NaN
Returns the n-th Harmonic number C<H_n>. The harmonic numbers are the sum of
reciprocals of the first C<n> natural numbers: C<1 + 1/2 + 1/3 + ... + 1/n>.
For values greater than 7000, binary splitting (Fredrik Johansson's elegant formulation) is used.
=head2 secant_number
secant_number(n) #=> Int
Returns the n-th secant number (A000364), starting with C<secant_number(0) = 1>.
=head2 tangent_number
tangent_number(n) #=> Int
( run in 0.616 second using v1.01-cache-2.11-cpan-26ccb49234f )