Skip to content

Commit

Permalink
- Added the Number egypt_greedy(p/q) method.
Browse files Browse the repository at this point in the history
Returns the array of denominators, such that `egypt_greedy(p/q).sum{|d| 1/d }` equals `p/q`.

Example:

	say egypt_greedy(9/10)      #=> [2, 3, 15]
	say egypt_greedy(5/121)     #=> [25, 757, 763309, 873960180913, 1527612795642093418846225]
  • Loading branch information
trizen committed Oct 2, 2024
1 parent ec8e01d commit 262b1a5
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 0 deletions.
38 changes: 38 additions & 0 deletions lib/Sidef/Types/Number/Number.pm
Original file line number Diff line number Diff line change
Expand Up @@ -8587,6 +8587,44 @@ package Sidef::Types::Number::Number {
return ((bless \$q1), (bless \$q2));
}

sub egypt_greedy {
my ($pq) = @_;

$pq = _any2mpq($$pq) // return undef;

my $n = Math::GMPz::Rmpz_init();
my $d = Math::GMPz::Rmpz_init();

Math::GMPq::Rmpq_get_num($n, $pq);
Math::GMPq::Rmpq_get_den($d, $pq);

my @list;

state $g = Math::GMPz::Rmpz_init_nobless();

while (Math::GMPz::Rmpz_cmp_ui($n, 1) != 0) {

my $x = Math::GMPz::Rmpz_init();
Math::GMPz::Rmpz_div($x, $d, $n);
Math::GMPz::Rmpz_add_ui($x, $x, 1);

push @list, $x;

Math::GMPz::Rmpz_mul($n, $n, $x);
Math::GMPz::Rmpz_sub($n, $n, $d);
Math::GMPz::Rmpz_mul($d, $d, $x);
Math::GMPz::Rmpz_gcd($g, $n, $d);

if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
Math::GMPz::Rmpz_divexact($n, $n, $g);
Math::GMPz::Rmpz_divexact($d, $d, $g);
}
}

push @list, $d;
Sidef::Types::Array::Array->new([map { bless \$_ } @list]);
}

sub dump {
my ($x) = @_;

Expand Down
13 changes: 13 additions & 0 deletions lib/Sidef/Types/Number/Number.pod
Original file line number Diff line number Diff line change
Expand Up @@ -2960,6 +2960,19 @@ Aliases: I<exponential_divisors>

=cut

=head2 egypt_greedy

egypt_greedy(p/q)

Greedy algorithm for Egyptian fraction expansion (also called the Fibonacci-Sylvester algorithm): at each step, extract the largest unit fraction less than the target and replace the target with the remainder.

say egypt_greedy(9/10) #=> [2, 3, 15]
say egypt_greedy(5/121) #=> [25, 757, 763309, 873960180913, 1527612795642093418846225]

Returns the array of denominators, such that C<egypt_greedy(p/q).sum{|d| 1/d }> equals C<p/q>.

=cut

=head2 ei

ei(x)
Expand Down
6 changes: 6 additions & 0 deletions scripts/Tests/number_methods.sf
Original file line number Diff line number Diff line change
Expand Up @@ -1526,4 +1526,10 @@ assert_eq([farey_neighbors(100, 89/99)], [80/89, 9/10])
assert_eq([farey_neighbors(100, 50/79)], [31/49, 19/30])
assert_eq([farey_neighbors(100, 19/30)], [50/79, 45/71])

assert_eq(egypt_greedy(9/10), %n[2, 3, 15])
assert_eq(egypt_greedy(5/121), %n[25, 757, 763309, 873960180913, 1527612795642093418846225])

assert_eq(egypt_greedy(43/97).sum{.inv}, 43/97)
assert_eq(egypt_greedy(97/43).sum{.inv}, 97/43)

say "** Tests passed!"

0 comments on commit 262b1a5

Please sign in to comment.