I am a supporter of
St. Joseph's hospice.
 If you find this site useful or if it helped you, consider a small donation to
St. Joseph's, please.

Information on
St. Joseph's

JustGiving - Sponsor me now!



"If you don't know where you are going, you might wind up someplace else." [Yogy Berra]

Which rational number is a good proxy of π (3.1415926...)? Enter in cell A1 ‘=pi()', in cell B1 your maximal denominator (for example 10), and in cells C1:D1 ‘=NRN(A1,B1)' as array formula (with CTRL + SHIFT + ENTER). You will get in C1:D1 22 and 7. That means: 22/7 is the nearest rational number to π with a denominator not higher than 10. For 1000 in B1 you would get 355/113.

This algorithm does NOT necessarily find the nearest rational number to a given floating point number with a given maximal denominator and a maximal absolute error. The good message is, though, that it would then return a #NUM! error. In this case please try an individual maximal absolute error.

The author's (Oliver Aberth) original intention was to support exact computation with rational numbers, for example solving a set of linear equations with rational coefficients.

Please read my disclaimer.

Option Explicit

Function sbNRN(dFloat As Double, lMaxDen As Long, _
    Optional dMaxErr As Double = -1#) As Variant
'Computes nearest rational number to dFloat with a maximal denominator
'lMaxDen and a maximal absolute error dMaxErr and returns result as a
'variant Nominator / Denominator.
'Select two adjacent cells in spreadsheet and enter as array formula.
'See: Oliver Aberth, A method for exact computation with rational numbers,
'     JCAM, vol 4, no. 4, 1978
'Reverse(moc.liborplus.www) V1.0 10-Nov-2019

Dim lSgn As Long
Dim dB As Double
Dim lA As Long
Dim lP1 As Long, lP2 As Long, lP3 As Long
Dim lQ1 As Long, lQ2 As Long, lQ3 As Long
Dim vR As Variant

If TypeName(Application.Caller) <> "Range" Then
   sbNRN = CVErr(xlErrRef)
   Exit Function
End If

With Application.Caller
ReDim vR(1 To .Rows.Count, 1 To .Columns.Count)

If dMaxErr = -1# Then dMaxErr = 1# / (2# * CDbl(lMaxDen) ^ 2)
lSgn = Sgn(dFloat)
lP1 = 0
lP2 = 1
lQ1 = 1
lQ2 = 0
dB = Abs(dFloat)

Do While lMaxDen > lQ2
    lA = Int(dB)
    lP3 = lA * lP2 + lP1
    lQ3 = lA * lQ2 + lQ1
    If Abs(dB - CDbl(lA)) < 0.00000000000001 Then
        Exit Do
    End If
    dB = 1# / (dB - CDbl(lA))
    lP1 = lP2
    lP2 = lP3
    lQ1 = lQ2
    lQ2 = lQ3

If lQ3 > lMaxDen Then
    lQ3 = lQ2
    lP3 = lP2
    If lQ2 > lMaxDen Then
        lQ3 = lQ1
        lP3 = lP1
    End If
End If

'If absolute error exceeds 1/2Q^2 then Aberth's lemma p. 286 might not apply.
'But the user can override this and check the result himself.
If Abs(dFloat - lSgn * lP3 / lQ3) > dMaxErr Then
    sbNRN = CVErr(xlErrNum)
    Exit Function
End If

vR(1, 1) = lSgn * lP3

If .Rows.Count > 1 Then
    vR(2, 1) = lQ3
    vR(1, 2) = lQ3
End If

sbNRN = vR
End With

End Function

Sulprobil   Get it done   Contact   Disclaimer   Impressum   Download