Reverse Greedy Issue

[Juggler_IN]

I need help with the following logic ---
Let N/QR be a fraction in lowest terms with Q > 1 and GCD(Q,R) = 1, and let x,y be the smallest positive solution of Nx - Qy = R. Then 1/Qx is a term in the expansion, with remainder y/Rx. Expand this remainder and iterate until reaching a unit fraction remainder.

Take the fraction 31/311. Here we have N = 31, Q = 311, R = 1. (Notice that the denominator QR = 311, and we require Q > 1 and also that Q,R must have no common factor. Therefore, we set Q = 311 and R = 1.) So for this first round, we need the smallest integer solution of 31x - 311y = 1. It’s easy to see by trial and error that the smallest solution is x = 301, y = 30. Therefore, the first term in our expansion is 1/Qx = 1/[(311)(301)] = 1/93611, with the remainder y/Rx = 30/[(1)(301)]. That is: 31/311 = 1/93611 + 30/301.

If we apply the same procedure to 30/301, we should get: 30/301 = 1/688 + 11/112.

With n=31, d=311, my code gives ReverseGreedy(31,311) = 1/93611 + 30/301 (Correct). But with n=30, d=301, my code gives ReverseGreedy(30,301) = 1/87591 + 29/291 (Incorrect); it should be 1/688 + 11/112.

How?

Function ReverseGreedy(n As Double, d As Double) As String
    
    Dim sUnit As String
    Dim x As Double, y As Double, q As Double
    Dim bFlag As Boolean

    sUnit = ""

    ' Smallest positive solution of Nx - Qy = 1.
    For x = 1 To d
        For y = 1 To n
            If (n * x - d * y) = 1 Then
                bFlag = True
                Exit For
            End If
        Next y
        If bFlag Then Exit For
    Next x

    q = d * x

    sUnit = "1/" & q & " + " & y & "/" & x

    ReverseGreedy = sUnit

End Function

 

[StephenCrump]

Your code only allows for R=1. The first case works because 311 is prime, i.e. R = 1.

In the 30/301 case, Q=43 and R=7 (not 1). Lowest integer solution for x is 16, so first term of expansion is 1/Qx = 1/(43 x 16) = 1/688.

 

[Juggler_IN]

@StephenCrump; How do I adjust the code ... attaching a function for Prime factors.

Function PrimeF(n As Long) As Variant

    Dim r As String
    Dim f As Long, i As Long
    Dim a As Variant

    i = 0
    ReDim a(0&)
    f = 2

    r = ""

    Do While n > 1
        If n Mod f = 0 Then
            r = r & f & ", "
            n = n / f
            a(i) = f
            i = i + 1
            ReDim Preserve a(0 To i)
        Else
            f = f + 1
        End If
    Loop

    If Len(r) > 0 Then r = Left(r, Len(r) - 2)

    ReDim Preserve a(0 To i - 1)

    PrimeF = a

End Function

 

[StephenCrump]

I hadn't encountered "greedy algorithms" being used in the context of fractions.

So some background first (which looks like where you're getting your examples): https://www.mathpages.com/home/kmath150/kmath150.htm

Here's an interesting algorithm for expanding any fraction into a sum of unit fractions. The method seems to give quite nice expansions that look rather "Egyptian". In fact, it reproduces many of the expansions listed in the Rhind Papyrus, etc. But, in addition, it easily generates a nice four-term expansion of 31/311 and in general seems to work well even on "difficult" fractions. Still, it's so simple that it could have been known to the ancient Egyptians.

It's fairly simple to edit your ReverseGreedy function to take a third parameter, r, and to use this in place of the fixed value 1. Rather than using nested loops, you could also use the Extended Euclidean Algorithm to get x and y more directly.

It would be useful then to make it recursive.

But first you have a decision on how to split the denominator into Q and r. With 301 it was straightforward, as 301 factorises neatly into 43 x 7. But what about 735, say, which is 7x7x5x3?

The Mathpages link notes: One interesting aspect of the algorithm described above is that, in order to yield unique results, we must include a rule for factoring the denominators. Interestingly, with Q = 1 it's equivalent to the Greedy Method, whereas with R = 1 it's equivalent to the Continued Fraction Method. Thus it appears that a wide variety of methods can be expressed in terms of the above algorithm together with some fixed rule for factoring the denominators.

From a historical standpoint, it seems (from what little I know of the surviving documentary evidence) that the ancient Egyptians didn't use any single rule for generating their unit fraction expansions. As a result, the small quantity of material must be sub-divided into even smaller sets of examples of any particular method, which makes it almost impossible to build conclusive arguments about the thought processes they may have used.

 

[Juggler_IN]

@StephenCrump; "But first you have a decision on how to split the denominator into Q and r. With 301 it was straightforward, as 301 factorises neatly into 43 x 7. But what about 735, say, which is 7x7x5x3?" ... an important point indeed ... Just to be able to manually take a fraction to its unit state (eg. 31/311 = 1/16 + 1/28 + 1/688 + 1/93611), we can pass both q & r ... like 49, 15 for 735 or 35, 21 ... and tabulate the unit fractions?

 

[StephenCrump]

Try this code. It's only lightly tested, but it replicates the several examples quoted in the Mathpages link above.

	A	B	C	D	E	F	G
1	Input	Output
2	163	1	1	1	1	1
3	431	118,525	55	350	42	3
4	Check
5	0.3782	0.0000	0.0182	0.0029	0.0238	0.3333
6	0.3782
Sheet1

Cell Formulas
Range	Formula
B5:G5	B5	=IFERROR(B2/B3,"")
A5	A5	=A2/A3
A6	A6	=SUM(B5:K5)

Dim MaxTerms As Long
Dim Fractions() As Variant
Sub Test()

    Dim N As Long, Q As Long, r As Long, f As Long
    
    MaxTerms = 10 'say
    ReDim Fractions(1 To 2, 1 To MaxTerms)
    
    N = Range("A2").Value
    f = Range("A3").Value
    Q = GetFactor(f)
    r = f / Q
    
    Call ReverseGreedy(N, Q, r, 1)
    
    Range("B2").Resize(2, MaxTerms) = Fractions
    
End Sub
Sub ReverseGreedy(N As Long, Q As Long, r As Long, k As Long)

    Dim x As Long, y As Long, f As Long, g As Long
    
    x = Application.RoundUp((Q + r) / N, 0)
    
    Do Until (N * x - r) Mod Q = 0
        x = x + 1
    Loop

    y = (N * x - r) / Q
    Fractions(1, k) = 1
    Fractions(2, k) = Q * x
    g = GCD(y, r * x)
    Fractions(1, k + 1) = y / g
    Fractions(2, k + 1) = r * x / g
    
    f = GetFactor(r * x / g)
    If k < MaxTerms - 1 And y / g > 1 Then Call ReverseGreedy(y / g, f, r * x / g / f, k + 1)

End Sub
Function GetFactor(ByVal N As Long) As Long

    Dim i As Long, iMax As Long, count() As Long
    
    ReDim count(1 To N)
    i = 2
    
    Do Until i > N
        If N Mod i = 0 Then
            N = N / i
            count(i) = count(i) + 1
            If i > iMax Then iMax = i
        Else
            i = i + 1
        End If
    Loop

    GetFactor = iMax ^ count(iMax)

End Function
Function GCD(ByVal A As Long, ByVal B As Long) As Long

    'We know B will be larger in these examples
    Dim r As Long
    
    Do Until A = 0
        r = B Mod A
        B = A
        A = r
    Loop
    
    GCD = B

End Function

 

[Juggler_IN]

@StephenCrump; Bingo! I appreciate your help.