Beefy Boxes and Bandwidth Generously Provided by pair Networks
XP is just a number
 
PerlMonks  

Egyptian fractions

by jimt (Chaplain)
on Aug 24, 2006 at 17:51 UTC ( [id://569404]=obfuscated: print w/replies, xml ) Need Help??

Probably not my best obfuscation, but I'm proud of some parts of it. Takes a fraction as an argument and returns the egyptian fraction for that number.

For example:

clark:~/Desktop jim$ ./egypt.pl 19/20 19/20 = 1/2 + 1/3 + 1/9 + 1/180 clark:~/Desktop jim$ ./egypt.pl 18/20 18/20 = 1/2 + 1/3 + 1/15
sub egypt {my$f=2;do{return(--$f,egypt($_[0]/$f))unless$_[0]%$f++} while$_[0]>1}print"$ARGV[0] = ",join(' + ',pharaoh(hiero(split '/', $ARGV[0]))),"\n";sub hiero{my(%nf,%df);$nf{$_}++for egypt(shift);++ $df{$_}for egypt(shift);do{($df{$_},$nf{$_})=($df{$_}-$nf{$_},$nf{$_} -$df{$_})if$df{$_}}for keys%nf;return(eval join('*',1,(map{($_)x $nf{$_}}keys%nf)),eval join('*',1,(map{($_)x$df{$_}}keys%df)))}sub pharaoh {return$_[0]==1?"$_[0]/$_[1]":("1/".(int($_[1]/$_[0])+1), pharaoh(hiero($_[0]*(int($_[1]/$_[0])+1)-$_[1],$_[1]*(int($_[1]/ $_[0])+1))))}

Since I know people like to know how these things work, I'll explain it. First, as always, we'll re-format to pretty it up.

sub egypt { my $f = 2; do {return (--$f, egypt($_[0] / $f)) unless $_[0] % $f++} while $_[0] > 1; } print "$ARGV[0] = " , join( ' + ', pharaoh(hiero(split '/', $ARGV[0])) ), "\n"; sub hiero { my (%nf, %df); $nf{$_}++ for egypt(shift); $df{$_}++ for egypt(shift); do {($df{$_}, $nf{$_}) = ($df{$_} - $nf{$_}, $nf{$_} - $df{$_}) if $df{$_} } for keys %nf; return (eval join('*', 1, (map {($_) x $nf{$_}} keys %nf)), eval join('*', 1, (map {($_) x $df{$_}} keys %df))); } sub pharaoh { return $_[0] == 1 ? "$_[0]/$_[1]" : ( "1/" . (int($_[1] / $_[0]) + 1), pharaoh(hiero($_[0] * (int($_[1] / $_[0]) + 1) - $_[1], $_[1] * (int($_[1] / $_[0]) + 1)))) }

The first function is arbitrarily named 'egypt' to make it blend it. It returns the prime factorization of whatever number is passed into it. I'm fairly proud of this one.

sub egypt { my $f = 2; do {return (--$f, egypt($_[0] / $f)) unless $_[0] % $f++} while $_[0] > 1; }

We initialize an arbitrary variable ($f) to 2 (smallest possible divisor), then, as long as the number we've passed in is > 1, check to see if it's divisible by $f. If it is, then return a list starting with $f and appending the recursed prime factorization of the number / $f. Note that we post increment $f to let us jump up to the next number (if the modulus is nonzero), so we pre-decrement it before returning so we actually have the one that was the divisor.

print "$ARGV[0] = " , join( ' + ', pharaoh(hiero(split '/', $ARGV[0])) ), "\n";

The next block does the algorithm to find egyptian fractions, as defined up in that wikipedia link. It just calls the pharaoh function with the LCD fraction and displays the results. Simple.

sub hiero { my (%nf, %df); $nf{$_}++ for egypt(shift); $df{$_}++ for egypt(shift); do {($df{$_}, $nf{$_}) = ($df{$_} - $nf{$_}, $nf{$_} - $df{$_}) if $df{$_} } for keys %nf; return (eval join('*', 1, (map {($_) x $nf{$_}} keys %nf)), eval join('*', 1, (map {($_) x $df{$_}} keys %df))); }

The hiero function just reduces a numerator & denominator to lowest terms. It's easy. First, we set up some hashes. Next, we populate the prime factorization of the numerator and denominator into each hash, using the number of appearances of the number as the value. So the prime factorization of 8 is 2,2,2, and we end up with 2 => 3 in our hash.

Then simply loop over the keys in the numerator hash. If that key exists in the denominator hash, then delete from the opposite hash the number of times the value occurs for you. For example, if you passed in 4/8, you'd end up with 2 => 2 in the numerator hash and 2 => 3 in the denominator hash. 2 is in both, so we subtract 3 from the numerator hash's value and 2 from the denominator hash's value, ending up with 2 => -1 in the numerator hash and 2 => 0 in the denominator hash.

This ensures that any common factors are removed from the numerator and denominator. Next, all we need to do is map out the factors and repeat them the number of times they appear in the hash. So 2 => 3 turns into (2,2,2). We toss on a 1 to this list to avoid errors when there are no factors left. Then we just join those lists with asterisks and eval 'em, yielding the LCD.

sub pharaoh { return $_[0] == 1 ? "$_[0]/$_[1]" : ( "1/" . (int($_[1] / $_[0]) + 1), pharaoh(hiero($_[0] * (int($_[1] / $_[0]) + 1) - $_[1], $_[1] * (int($_[1] / $_[0]) + 1)))) }

pharaoh here is where the magic is. It takes a numerator and denominator as arguments, then recursively builds up the list of unit fractions to return. First of all, if the numerator is 1, then return numerator / denominator (we're done). Otherwise, we apply the algorithm and return a list consisting of int( 1 / (denominator / numerator) ) + 1, and the result of calling pharaoh on the LCD form of the fraction minus the value we just calculated. Further explanation is left as an exercise to the reader, but I'd recommend assigning this (int($_1 / $_[0]) + 1) to a variable and replacing in the function, it'll be clearer.

Replies are listed 'Best First'.
Re: Egyptian fractions (Golf Anyone?)
by Limbic~Region (Chancellor) on Aug 25, 2006 at 15:22 UTC
    jimt,
    Obfu is not my thing, but I do love interesting math problems. The following is more of a golf (156) than an obfu:
    use Math::Pari 'lcm';my($n,$d)=split/\//,shift;sub p{print"1/$_[0] "}w +hile(1){my $u=int$d/$n+1;p$u;my $l=lcm$d,$u;$n=$n*$l/$d-$l/$u;$d=$l;p($d),die$/if +$n==1}
    Golfing with strictures and warnings is a handicap so (146):
    ($n,$d)=split/\//,shift;sub p{print"1/$_[0] "}while(1){use Math::Pari lcm;$u=int$d/$n+1;p$u;$l=lcm$d,$u;$n=$n*$l/$d-$l/$u;$d=$l;p($d),die$/i +f$n==1}
    Since I am not a golfer either, anyone interested in the starting code can use it to improve things.
    #!/usr/bin/perl use strict; use warnings; use Math::Pari qw/gcd lcm/; my ($num, $den) = split m|/|, $ARGV[0]; my @result; while (1) { my $unit = int($den / $num + 1); push @result, $unit; my $lcm = lcm($den, $unit); $num = ($num * ($lcm / $den)) - ($lcm / $unit); $den = $lcm; push @result, $den and last if $num == 1; } print "1/$_ " for @result;

    Cheers - L~R

      Hi Limbic~Region,

      A couple of small changes will get you down to 142:

      ($n,$d)=split'/',shift;sub p{print"1/$_[0] "}{use Math::Pari lcm;$u=int$d/$n+1;p$u;$l=lcm$d,$u;$n=$n*$l/$d-$l/$u;$d=$l;p($d),die$/i +f$n==1;redo}

      s''(q.S:$/9=(T1';s;(..)(..);$..=substr+crypt($1,$2),2,3;eg;print$..$/
        Refinement of this particular track at 130:
        ($n,$d)=split'/',pop;sub p{print"1/@_ "}{use Math::Pari lcm;p$u=int$d/$n+1;$l=lcm$d,$u;$n=$n*$l/$d-($d=$l)/$u;$n-1?redo: p$d+die$/}
      81 chars:
      $ perl -e '($n,$d)=split"/",pop;{1while++$x<$d/$n;warn"1/$x\n";$n=$n*$ +x-$d;$d*=$x;redo if$n}' 18/20 1/2 1/3 1/15
      Update: changing the entire structure to a C-style for-loop saves 2 chars, so this gives me a 79-char solution:
      for(($n,$d)=split"/",pop;$n;){1while++$x<$d/$n;warn"1/$x\n";$n=$n*$x-$ +d;$d*=$x}
      This uses the greedy heuristic, so it doesn't find the "best" Egyptian fractions. $x keeps increasing until 1/$x is <= our current fraction. Then we print 1/$x and subtract from the current fraction. But there is no need to put this fraction in lowest terms, which saves a lot of strokes.

      If POSIX::ceil were available, 1while++$x<$d/$n is really just $x=ceil($d/$n) ..

      BTW, this 58-char solution would work if it weren't for floating-point error:

      perl -e '$f=eval pop;{1while++$x<1/$f;warn"1/$x\n";redo if$f-=1/$x}' 1 +9/20
      I don't have Math::Pari where I am, but it may be possible to use it in the above approach (it treats rationals with absolute precision).

      blokhead

      Ah, no, Limbic~Region,

      your's ain't golf. No modules, lest I would say

      use Junk;do

      - 11 chars.. but that could be golfed down to use J;d - 7 chars ;-)

      blokhead's solutions look neat, but they hog my cpu... (with 2355/12344 - still not finished in 3/4 hour ;-)

      #!/usr/bin/perl -l pop=~m|/|;($f,$g)=($`,$');sub d{int($_[1]/$_[0]+1)}sub g{($x,$y)=@_;($x,$y)=($y,$x% $y)while$y;$x} sub re{($p ,$e,$r,$l)=@_;($p,$l)=($p*$l-$e*$r,$e*$l);$g=g($p,$l); for($p,$l){$_/=$g};($p,$l)}while($f>1){push@o,"1/".d(# $f,$g);($f,$g)=re($f,$g,1,d($f,$g));}print join' + ',# @o,"$f/$g";# ungolfed and thus not for production use!

      way too long...

      <update> golfed down a bit... (198 chars, counting newlines).

      ($z,$n)=($_=pop)=~/(.+)\/(.+)/;$s='==';for(;;){$m=int($n/$z+1);$_ .=" $s 1/".($z==1?$n:$m);$z<=1&&last;($z,$n)=($z*$m-$n,$m*$n);($x ,$y)=($z,$n);($x,$y)=($y,$x%$y)while$y;$z/=$x;$n/=$x;$s='+'}print
      It computes egyptian fractions like this
      qwurx [shmem] ~> perl -l egy.pl 2355/12344 2355/12344 == 1/6 + 1/42 + 1/3282 + 1/15755059 + 1/744665636525384

      in no time... </update>

      --shmem

      _($_=" "x(1<<5)."?\n".q·/)Oo.  G°\        /
                                    /\_¯/(q    /
      ----------------------------  \__(m.====·.(_("always off the crowd"))."·
      ");sub _{s./.($e="'Itrs `mnsgdq Gdbj O`qkdq")=~y/"-y/#-z/;$e.e && print}
        A modification of blokhead's idea... what golf can go without regular expressions? 72 strokes:
        for($_=pop;/\//,$`;$_=$x*$`-$'.'/'.$'*$x){1while++$x<$'/$`;print"1/$x "}
        I counted the line break as one character, is that kosher?
        Update: I see also shmem had a similar thought, I had missed that before posting somehow.
        Update again: Changed the title of this node to not be dumb. Sorry all, thanks for the feedback those who messaged me!
        ~dewey
Re: Egyptian fractions
by dewey (Pilgrim) on Aug 25, 2006 at 02:30 UTC
    Nice work! A good explanation too, I will need to take this code apart before I really understand it thoroughly. For some reason hiero catches my fancy in particular...
    ~dewey

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: obfuscated [id://569404]
Approved by Limbic~Region
Front-paged by shmem
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others surveying the Monastery: (7)
As of 2024-03-28 08:26 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found