Egyptian Fractions

Question

Overview:

From Wikipedia: An Egyptian fraction is the sum of distinct unit fractions. That is, each fraction in the expression has a numerator equal to 1 and a denominator that is a positive integer, and all the denominators differ from each other. The value of an expression of this type is a positive rational number a/b. Every positive rational number can be represented by an Egyptian fraction.

Challenge:

Write the shortest function that will return the values of all the denominators for the smallest set of unit fractions that add up to a given fraction.

Rules/Constraints:

• Input will be two positive integer values.
  • This can be on STDIN, argv, comma separated, space delimited, or any other method you prefer.
• The first input value shall be the numerator and the second input value the denominator.
• The first input value shall be less than the second.
• The output may include a value(s) that exceeds the memory limitations of your system/language (RAM, MAX_INT, or whatever other code/system constraints exist). If this happens, simply truncate the result at the highest possible value and note that somehow (i.e. ...).
• The output should be able to handle a denominator value up to at least 2,147,483,647 (2³¹-1, signed 32-bit int).
  • A higher value (long, etc.) is perfectly acceptable.
• The output shall be a listing of all values of denominators of the smallest set of unit fractions found (or the fractions themselves, i.e. 1/2).
• The output shall be ordered ascending according to the value of the denominator (descending by the value of the fraction).
• The output can be delimited any way you wish, but there must be some character between so as to differentiate one value from the next.
• This is code golf, so the shortest solution wins.

Exmaples:

• Input 1:

  43, 48

• Output 1:

  2, 3, 16

• Input 2:

  8/11

• Output 2:

  1/2 1/6 1/22 1/66

• Input 3:

  5 121

• Output 3:

  33 121 363

Answer 1

Python 2, 169 167 chars

x,y=input()
def R(n,a,b):
 if n<2:return[b/a][b%a:]
 for m in range((b+a-1)/a,b*n/a):
  L=R(n-1,a*m-b,m*b)
  if L:return[m]+L
n=L=0
while not L:n+=1;L=R(n,x,y)
print L

Takes comma-separated args on stdin and prints a python list on stdout.

$ echo 8,11 | ./egypt.py 
[2, 5, 37, 4070]

Answer 2

Common Lisp, 137 chars

(defun z(n)(labels((i(n s r)(cond((= n 0)r)((< n(/ 1 s))(i n(ceiling(/ 1 n))r))(t(i(- n(/ 1 s))(1+ s)(cons s r))))))(reverse(i n 2'()))))

(z 43/48) -> (2 3 16)

(z 8/11) -> (2 5 37 4070)

(z 5/121) -> (25 757 763309 873960180913 1527612795642093418846225)

No need to worry about huge numbers, or handling fractional notation!

Answer 3

PHP 82 bytes

<?for(fscanf(STDIN,"%d%d",$a,$b);$a;)++$i<$b/$a||printf("$i ",$a=$a*$i-$b,$b*=$i);

This could be made shorter, but the current numerator and denominator need to be keep as whole numbers to avoid floating point rounding error (instead of keeping the current fraction).

Sample usage:

$ echo 43 48 | php egyptian-fraction.php
2 3 16
$ echo 8 11 | php egyptian-fraction.php
2 5 37 4070

Answer 4

C, 163 177 chars

6/6: At last, the program now correctly handles truncation in all cases. It took a lot more chars than I was hoping for, but it was worth it. The program should 100% adhere to the problem requirements now.

d[99],c,z;
r(p,q,n,i){for(c=n+q%p<2,i=q/p;c?d[c++]=i,0:++i<n*q/p;)q>~0U/2/i?c=2:r(i*p-q,i*q,n-1);}
main(a,b){for(scanf("%d%d",&a,&b);!c;r(a,b,++z));while(--c)printf("%d\n",d[c]);}

The program takes the numerator and denominator on standard input. The denominators are printed to standard output, one per line. Truncated output is indicated by printing a zero denominator at the end of the list:

$ ./a.out
2020 2064
2
3
7
402
242004

$ ./a.out
6745 7604
2
3
19
937
1007747
0

The denominators in the second example sum to 95485142815 / 107645519046, which differs from 6745 / 7604 by roughly 1e-14.

Answer 5

Python, 61 chars

Input from STDIN, comma separated.
Output to STDOUT, newline separated.
Doesn't always return the shortest representation (e.g. for 5/121).

a,b=input()
while a:
    i=(b+a-1)/a
    print"1/%d"%i
    a,b=a*i-b,i*b

Characters counted without unneeded newlines (i.e. joining all lines within the while using ;).
The fraction is a/b.
i is b/a rounded up, so I know 1/i <= a/b.
After printing 1/i, I replace a/b with a/b - 1/i, which is (a*i-b)/(i*b).

Answer 6

C, 94 bytes

n,d,i;main(){scanf("%i%i",&n,&d);for(i=1;n>0&++i>0;){if(n*i>=d)printf("%i ",i),n=n*i-d,d*=i;}}

Try It Online

edit: A shorter version of the code was posted in the comments so I replaced it. Thanks!

Answer 7

Stax, 18 bytes

é├WüsOΩ↨÷╬6H╒σw[▐â

Run and debug it

At each step, it tries to minimize the subsequent numerator. It seems to work, but I can't prove it.

Answer 8

APL (Dyalog Extended), 7 bytes

⌂efract

Try it online!

From some reason, Dyalog APL has a builtin for this.

Takes the numerator and denominator on the left and the right

Answer 9

Husk, 8 bytes

ḟo=⁰ΣṖİ\

Try it online!

Larger denominators take very long, as we are calculating the powerset of an infinite list.

Takes input as a fraction a/b.

Explanation

ḟo=⁰ΣṖİ\
      İ\ Infinite list [1,1/2,1/3...]
     Ṗ   Powerset
ḟo       first sublist which satisfies:
    Σ    sum
  =⁰     equals input?

Answer 10

AXIOM, 753 bytes

L==>List FRAC INT
macro  M(q)==if c<m or(c=m and m<999 and reduce(max,map(denom,q))<xv)then(m:=c;a:=q;xv:=reduce(max,map(denom,a)))
f(x,n)==(y:=x;a:L:=[];c:=0;q:=denom x;q:=q^4;for i in n.. repeat((c:=c+1)>50=>(a:=[];break);1/i>y=>1;member?(1/i,a)=>1;a:=concat(a,1/i);(y:=y-1/i)=0=>break;numer(y)=1 and ~member?(y,a)=>(a:=concat(a,y);break);(i:=floor(1/y))>q=>(a:=[];break));a)
h(x:FRAC INT):L==(a:L:=[];x>1=>a;numer(x)=1=>[x];n:=max(2,floor(1/x));xv:=m:=999;d:=denom x;zd:=divisors d;z:=copy zd;for i in 2..30 repeat z:=concat(z,i*zd);d:=min(10*d,n+9*m);for i in n..d repeat((c:=maxIndex(b:=f(x,i)))=0=>1;c>m+1=>1;M(b);v:=reduce(+,delete(b,1));for j in z repeat((c:=1+maxIndex(q:=f(v,j)))=1=>1;member?(b.1,q)=>1;q:=concat(b.1,q);M(q)));reverse(sort a))

The idea would be apply the "Greedy Algorithm" with different initial points, and save the list that has minimum length. But not always it would find the min solution with less difined: "array A will be less than array B if and only if A has few elements of B, or if the number of elements of A is the same of number of elements of B, than A it is less than B if the more little element of A is bigger as number, than the more little element of B". Ungolfed and test

-- this would be the "Greedy Algorithm"
fracR(x,n)==
   y:=x;a:L:=[];c:=0;q:=denom x;q:=q^4
   for i in n.. repeat
      (c:=c+1)>50   =>(a:=[];break)
      1/i>y         =>1
      member?(1/i,a)=>1
      a:=concat(a,1/i)
      (y:=y-1/i)=0  =>break
      numer(y)=1 and ~member?(y,a)=>(a:=concat(a,y);break)
      (i:=floor(1/y))>q           =>(a:=[];break)
   a

-- Return one List a=[1/x1,...,1/xn] with xn PI and x=r/s=reduce(+,a) or return [] for fail
Frazione2SommaReciproci(x:FRAC INT):L==
    a:L:=[]
    x>1       =>a
    numer(x)=1=>[x]
    n:=max(2,floor(1/x));xv:=m:=999;d:=denom x;zd:=divisors d;z:=copy zd
    for i in 2..30 repeat z:=concat(z,i*zd)
    d:=min(10*d,n+9*m) 
    for i in n..d repeat
        (c:=maxIndex(b:=fracR(x,i)))=0=>1 
        c>m+1                         =>1
        M(b)
        v:=reduce(+,delete(b,1))
        for j in z repeat
              (c:=1+maxIndex(q:=fracR(v,j)))=1=>1
              member?(b.1,q)                  =>1
              q:=concat(b.1,q)
              M(q) 
    reverse(sort a)

(7) -> [[i,h(i)] for i in [1/23,2/23,43/48,8/11,5/121,2020/2064,6745/7604,77/79,732/733]]
   (7)
      1   1      2   1  1      43  1 1  1      8  1 1  1  1
   [[--,[--]], [--,[--,---]], [--,[-,-,--]], [--,[-,-,--,--]],
     23  23     23  12 276     48  2 3 16     11  2 6 22 66
      5    1  1   1      505  1 1 1  1    1
    [---,[--,---,---]], [---,[-,-,-,---,----]],
     121  33 121 363     516  2 3 7 602 1204
     6745  1 1  1  1    1      1       77  1 1 1  1  1   1
    [----,[-,-,--,---,-----,------]], [--,[-,-,-,--,---,---]],
     7604  2 3 19 950 72238 570300     79  2 3 8 79 474 632
     732  1 1 1  1   1    1     1
    [---,[-,-,-,--,----,-----,-----]]]
     733  2 3 7 45 7330 20524 26388
                                                      Type: List List Any
       Time: 0.07 (IN) + 200.50 (EV) + 0.03 (OT) + 9.28 (GC) = 209.88 sec
(8) -> h(124547787/123456789456123456)
   (8)
        1             1                         1
   [---------, ---------------, ---------------------------------,
    991247326  140441667310032  613970685539400439432280360548704
                                     1
    -------------------------------------------------------------------]
    3855153765004125533560441957890277453240310786542602992016409976384
                                              Type: List Fraction Integer
                     Time: 17.73 (EV) + 0.02 (OT) + 1.08 (GC) = 18.83 sec
(9) -> h(27538/27539)
         1 1 1  1  1    1      1        1
   (9)  [-,-,-,--,---,-----,------,----------]
         2 3 7 52 225 10332 826170 1100871525
                                              Type: List Fraction Integer
                     Time: 0.02 (IN) + 28.08 (EV) + 1.28 (GC) = 29.38 sec

reference and numbers from: http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fractions/egyptian.html

for add something, this below would be the one optimized for find min length fraction that has the max denominator less (and not optimized for lengh)

L==>List FRAC INT

-- this would be the "Greedy Algorithm"
fracR(x,n)==
   y:=x;a:L:=[];c:=0;q:=denom x;q:=q^20
   for i in n.. repeat
      (c:=c+1)>1000  =>(a:=[];break)
      1/i>y          =>1
      member?(1/i,a) =>1
      a:=concat(a,1/i)
      (y:=y-1/i)=0  =>break
      numer(y)=1 and ~member?(y,a)=>(a:=concat(a,y);break)
      (i:=floor(1/y))>q           =>(a:=[];break)
   a

-- Return one List a=[1/x1,...,1/xn] with xn PI and x=r/s=reduce(+,a) or return [] for fail
Frazione2SommaReciproci(x:FRAC INT):L==
    a:L:=[]
    x>1       =>a
    numer(x)=1=>[x]
    n:=max(2,floor(1/x));xv:=m:=999;d:=denom x;zd:=divisors d;z:=copy zd; 
    w1:= if d>1.e10 then 1000 else 300; w2:= if d>1.e10 then 1000 else if d>1.e7 then 600 else if d>1.e5 then 500 else if d>1.e3 then 400 else 100;
    for i in 2..w1 repeat(mt:=(i*zd)::List PI;mv:=[yy for yy in mt|yy>=n];z:=sort(removeDuplicates(concat(z,mv)));#z>w2=>break)
    for i in z repeat
        (c:=maxIndex(b:=fracR(x,i)))=0=>1 
        c>m+1                         =>1
        if c<m or(c=m and m<999 and reduce(max,map(denom,b))<xv)then(m:=c;a:=b;xv:=reduce(max,map(denom,a)))
        v:=reduce(+,delete(b,1))
        for j in z repeat
              (c:=1+maxIndex(q:=fracR(v,j)))=1=>1
              member?(b.1,q)                  =>1
              q:=concat(b.1,q)
              if c<m or(c=m and m<999 and reduce(max,map(denom,q))<xv)then(m:=c;a:=q;xv:=reduce(max,map(denom,a)))
    reverse(sort a)

the results:

(5) -> [[i,Frazione2SommaReciproci(i)] for i in [1/23,2/23,43/48,8/11,5/121,2020/2064,6745/7604,77/79,732/733]]
   (5)
      1   1      2   1  1      43  1 1  1      8  1 1  1  1
   [[--,[--]], [--,[--,---]], [--,[-,-,--]], [--,[-,-,--,--]],
     23  23     23  12 276     48  2 3 16     11  2 6 22 66
      5    1  1   1      505  1 1 1  1    1
    [---,[--,---,---]], [---,[-,-,-,---,----]],
     121  33 121 363     516  2 3 7 602 1204
     6745  1 1  1  1    1      1       77  1 1 1  1  1   1
    [----,[-,-,--,---,-----,------]], [--,[-,-,-,--,---,---]],
     7604  2 3 19 950 72238 570300     79  2 3 8 79 474 632
     732  1 1 1  1   1    1     1
    [---,[-,-,-,--,----,-----,-----]]]
     733  2 3 7 45 7330 20524 26388
                                                      Type: List List Any
                     Time: 0.08 (IN) + 53.45 (EV) + 3.03 (GC) = 56.57 sec
(6) -> Frazione2SommaReciproci(124547787/123456789456123456)
   (6)
        1            1               1                  1
   [---------, ------------, ----------------, -------------------,
    994074172  347757767307  2764751529594496  1142210063701888512
                      1
    -------------------------------------]
    2531144929865351036156388364636113408
                                              Type: List Fraction Integer
         Time: 0.15 (IN) + 78.30 (EV) + 0.02 (OT) + 5.28 (GC) = 83.75 sec
(7) -> Frazione2SommaReciproci(27538/27539)
         1 1 1  1   1     1       1       1
   (7)  [-,-,-,--,----,-------,-------,-------]
         2 3 7 43 1935 3717765 5204871 7105062
                                              Type: List Fraction Integer
                     Time: 0.05 (IN) + 45.43 (EV) + 2.42 (GC) = 47.90 sec

It seems many good denominators have as factor divisors of the input fraction denominator.

Answer 11

APL(NARS), 2502 bytes

fdn←{1∧÷⍵}⋄fnm←{1∧⍵}⋄ffl←{m←⎕ct⋄⎕ct←0⋄r←⌊⍵⋄⎕ct←m⋄r}⋄divisori←{a[⍋a←{∪×/¨{0=≢⍵:⊂⍬⋄s,(⊂1⌷⍵),¨s←∇1↓⍵}π⍵}⍵]}

r←frRF w;x;y;c;q;i;j
(x i)←w⋄i-←1⋄y←x⋄r←⍬⋄c←0⋄q←fdn x⋄q←q*20
i+←1
→4×⍳∼1000<c+←1⋄→6
j←÷i⋄→2×⍳j>y⋄→2×⍳(⊂j)∊r⋄r←r,(⊂j)⋄y←y-j⋄→0×⍳y=0⋄→5×⍳1≠fnm y⋄→5×⍳(⊂y)∊r⋄r←r,⊂y⋄→0
→2×⍳∼q<i←ffl ÷y
r←⍬

r←fr2SumF x;n;xv;m;d;zd;z;i;b;c;t;v;j;k;q;w1;w2;t;b1
z←r←⍬⋄→0×⍳1≤ffl x
:if 1=fnm x⋄r←,⊂x⋄→0⋄:endif
n←2⌈ffl÷x⋄xv←m←999⋄d←fdn x⋄zd←divisori d
w1←1000⋄w2←50⋄:if d>1.e10⋄w2←700⋄:elseif d>1.e7⋄w2←600⋄:elseif d>1.e5⋄w2←500⋄:elseif d>1.e3⋄w2←400⋄:elseif d>1.e2⋄w2←100⋄:endif
:for i :in ⍳w1⋄z←∪z∪k/⍨{⍵≥n}¨k←i×zd⋄:if w2<≢z⋄:leave⋄:endif⋄:endfor
z←∪z∪zd⋄z←z[⍋z]
:for i :in z
    :if 0=c←≢b←frRF x i ⋄:continue⋄:endif
    :if      c>m+1      ⋄:continue⋄:endif
    :if      c<m        ⋄m←c⋄r←b⋄xv←⌈/fdn¨b
    :elseif (c=m)∧(m<999)
         :if xv>t←⌈/fdn¨b⋄m←c⋄r←b⋄xv←t⋄:endif
    :endif
    :if c≤2⋄:continue⋄:endif
    v←↑+/1↓b⋄b1←(⊂↑b)
    :for j :in z
       :if 1=c←1+≢q←frRF v j⋄:continue⋄:endif
       :if        b1∊q      ⋄:continue⋄:endif
       q←b1,q
       :if  c<m⋄m←c⋄r←q     ⋄xv←⌈/fdn¨q
       :elseif (c=m)∧(m<999)
           :if xv>t←⌈/fdn¨q⋄m←c⋄r←q⋄xv←t⋄:endif
       :endif
    :endfor
:endfor
→0×⍳1≥≢r⋄r←r[⍋fdn¨r]

Traslation from AXIOM code for this problem, to APL, using for the first time (for me) the fraction type (that is bignum...).

103r233 means the fraction 103/233. Test:

  ⎕fmt fr2SumF 1r23
┌1────┐
│ 1r23│
└~────┘
  ⎕fmt fr2SumF 2r23
┌2──────────┐
│ 1r12 1r276│
└~──────────┘
  ⎕fmt fr2SumF 43r48
┌3────────────┐
│ 1r2 1r3 1r16│
└~────────────┘
  fr2SumF 8r11
1r2 1r6 1r22 1r66 
  fr2SumF 5r121
1r33 1r121 1r363 
  fr2SumF 2020r2064
1r2 1r3 1r7 1r602 1r1204 
  fr2SumF 6745r7604
1r2 1r3 1r19 1r950 1r72238 1r570300 
  fr2SumF 77r79
1r2 1r3 1r8 1r79 1r474 1r632 
  fr2SumF 732r733
1r2 1r3 1r7 1r45 1r7330 1r20524 1r26388 
  fr2SumF 27538r27539
1r2 1r3 1r7 1r43 1r1935 1r3717765 1r5204871 1r7105062 
  fr2SumF 124547787r123456789456123456
1r994074172 1r347757767307 1r2764751529594496 1r1142210063701888512 
  1r2531144929865351036156388364636113408 

Answer 12

GAP, 340 bytes

f:=function(b,l)local x,R;if l=1 then if IsInt(b)then
return[[b]];fi;return[];fi;R:=[];for x in[Int(b)+1..Int(l*b)]do
UniteSet(R,List(Filtered(f(b*x/(x-b),l-1),d->not x in d),s->UnionSet([x],s)));od;return
R;end;e:=function(n,d)local l,R;l:=0;repeat l:=l+1;R:=f(d/n,l);until R<>[];return
Filtered(R,s->s[l]=Minimum(List(R,s->s[l])))[1];end;

Try it online!

The entry function is e, which takes two parameters, the numerator and denominator in that order, as required. 2 characters could be golfed (replace n,d with q and d/n with 1/q) if I may make e take the fraction n/d direct, and 2 more if I may make e take its reciprocal d/n. But this goes against the question.

The helper function f returns a list of the ways to represent 1/b in l distinct terms.