| Subject: | Incorrect values, bigints, etc. |
1. goedel(10000) should be 2. It instead comes out as 11. Fix:
Use "sort {$a <=> $b}" in goedel sub.
2. goedel(4999) should be 24821251455656250000. It gets turned into a double and no longer is usable. This is a pretty significant issue.
There are various ways to fix this. Here is a way to do it, which vastly simplifies the function as well as being faster. Note it changes the dependencies:
use Math::Prime::Util qw/nth_prime/;
sub goedel {
my $n = shift;
my %opts = (
q/offset/ => 0,
q/reverse/ => 0,
@_ );
croak "n should be a non-negative integer"
if $n eq '' || $n =~ tr/0123456789//c;
my $offset = $opts{'offset'};
croak "offset should be a non-negative integer"
if $offset eq '' || $offset =~ tr/0123456789//c;
my $one = $n - $n + 1;
my $g = $one;
$n = reverse $n if $opts{'reverse'};
foreach my $i (1 .. length($n)) {
$g *= ($one * nth_prime($i)) ** (substr($n, $i-1, 1)+$offset);
}
$g;
}
No need for reduce, max, pairwise, pow cache, _next_prime, etc. If the input is a BigInt, then the result is also a BigInt (with the same back end). There are other solutions from Math::Prime::Util other than using nth_prime (e.g. forprimes or prime_iterator) but this works fine unless we're sending in billion digit numbers (the exponentiation and product will kill us long before we hti that point).
$one will be either an integer 1 or a BigInt 1 based on the type of the input. If it is a bigint, then multiplying it times the nth prime turns that into a bigint so the exponentiation is done as a bigint. If one were really concerned about performance it would be easy enough to split, e.g.:
if (ref($n) eq 'Math::BigInt') {
# ... it's a bigint, so do appropriate stuff
} else {
# ... it's not a bigint, so don't.
}
Alternately, we could try to do the right thing based on the input:
use Math::BigInt try => "GMP,Pari";
and then do bigint processing entirely based on the size of the input and whether we are 32-bit or 64-bit.