Math : Fraction, Transform 0.54 in 27/50

[Michael Vogel]

Simple, fast (?) and exact

But too exact for me, I prefer solutions allowing a (small) delta

Code: Select all

Procedure.q gcd(m.q,n.q)

       Protected k.q = 0

       If n < 0 : n * (-1) : EndIf
       If m < 0 : m * (-1) : EndIf

       While (m&1=0 And n&1=0)
               m>>1
               n>>1
               k+1
       Wend
       If n&1=0
               Swap m, n
       EndIf

       While m <> 0
               While (m&1=0)
                       m >> 1
               Wend

               If m < n
                       Swap m, n
               EndIf
               m - n
       Wend

       ProcedureReturn n * (1 << k)

EndProcedure
Procedure DoFract(d.d)

       Protected z.q
       Protected n.q=1
       Protected g.q=0

       Protected dummy.s=StrD(d)

       Protected i=FindString(dummy,".",1)
       If i
               z=Val(ReplaceString(dummy,".",""))
               n=Val(LSet("1",Len(dummy)-i+1,"0"))
       EndIf

       g=gcd(z,n)

       z/g
       n/g

       Debug StrD(d)+" := "+Str(z)+" / "+Str(n)

EndProcedure

For i=1 To 10
       DoFract(ValD(Left(StrD(#PI),i)))
Next i

[Le Soldat Inconnu]

@Michael Vogel : I just try with this
DoFract(57/67)
and it's give back
0.8507462687 := 8507462687 / 10000000000
Not good too



I do Calculator in this moment, here a screen

Like you see, i arrive to find fraction with Pi, Sqr(2) and Sqr(3)
It's simple, if my algo (see first post) doesn't give result for a value, i try to have result this value/Pi. and if it"s give back result, i know i have to multiply my fraction by Pi

[Michael Vogel]

@Michael Vogel : I just try with this
DoFract(57/67)
and it's give back
0.8507462687 := 8507462687 / 10000000000
Not good too
[...]
Like you see, i arrive to find fraction with Pi, Sqr(2) and Sqr(3)
It's simple, if my algo (see first post) doesn't give result for a value, i try to have result this value/Pi. and if it"s give back result, i know i have to multiply my fraction by Pi

Your calculator looks brilliant

The 57/67 point just show that some (periodical) values have no exact representation in a floating point number (btw I have not known, that StrD(value) truncates the value and eliminates some ciphers as well ) -- I know, I should have put some quotes around the word exact in my posting

Michael

PS found a calculator program which does exactly what I would do, it's here: http://www.hpmuseum.org/software/41/dec2fr.htm
PSS I've stolen the code of drgolf and modified it to have a function which returns the closer value....

Code: Select all

; ** find rational approximation To given real number
; ** David Eppstein / UC Irvine / 8 Aug 1993
; **
; ** With corrections from Arno Formella, May 2008
; **
; ** PureBasic portage by DrGolf, Mai 2010

Procedure.s Fraction(start.d,flag.l=#True,maxden.q=10000)

	Protected x.d
	Protected ai.q
	Protected t.q
	Protected p0.q
	Protected p1.q
	Protected q0.q
	Protected q1.q

	p1=1: q1 = 1
	p0=0: q0 = 0

	x=start

	While (q0 *  (Int(x) + q1) <= maxden)

		ai=Int(x)

		t = p1*ai+p0
		p0 = p1
		p1 = t

		t = q0*ai+q1
		q1 = q0
		q0 = t

		If x-ai=0
			Debug "Zero Exit..."
			Break ;     // AF: division by zero
		EndIf

		x = 1/(x - ai)
		If (x>$7FFFFFFF)
			Debug "Exit..."
			Break;  // AF: representation failure
		EndIf

	Wend

	; Solution 2...
	ai = (maxden - q1) / q0
	p0 = p1 * ai + p0
	q1 = q0 * ai + q1

	If Abs(start - (p1 / q0)) > Abs(start - (p0 / q1))
		;Debug "Solution 2 prefered..."
		p1=p0
		q0=q1
	EndIf

	;Debug Str(p1)+"/"+Str(q0)+"  error : "+StrD(start - (p1 / q0))
	If flag And p1>q0
		ProcedureReturn Str(Int(p1/q0))+" + "+Str(p1%q0)+" / "+Str(q0)
	Else
		ProcedureReturn Str(p1)+" / "+Str(q0)
	EndIf

EndProcedure

Debug Fraction(#PI,0)
Debug Fraction(#PI,1,100000)

[Al_the_dutch]

Thanks, I was looking for this and it is very useful. But.. I found a tiny error in it... Brackets... Here is my version.

Code: Select all

; https://www.purebasic.fr/english/viewtopic.php?t=42091&start=15

; ** find rational approximation To given real number
; ** David Eppstein / UC Irvine / 8 Aug 1993
; ** With corrections from Arno Formella, May 2008
; ** PureBasic portage by DrGolf, Mai 2010
; ** Al_the_dutch July 2024

EnableExplicit

Procedure.s Fraction(Start.d, MaxDen.q=10000)

	Protected x.d
	Protected.q x_int, t, NomAlt= 0, Nom = 1, DeNom = 0, DeNomAlt = 1

	x = Start
	x_int = Int(x)
	
	;While DeNom * x_int + DeNomAlt <= MaxDen; OK!!
	While DeNom * (x_int + DeNomAlt) <= MaxDen; Originally 	While (q0 *  (Int(x) + q1) <= maxden); not correct
	   ;                  =                          = Remove brackets

		t = Nom * x_int + NomAlt
		NomAlt = Nom
		Nom = t

		t = DeNom * x_int + DeNomAlt
		DeNomAlt = DeNom
		DeNom = t

		If x - x_int = 0
			Debug "Zero Exit..."
			Break ;     // AF: division by zero
		EndIf

		x = 1/(x - x_int)
		If (x>$7FFFFFFF)
			Debug "Exit..."
			Break;  // AF: representation failure
		EndIf
		
		x_int = Int(x)
		
	Wend

	; Solution 2...
	x_int = (MaxDen - DeNomAlt) / DeNom
	NomAlt = Nom * x_int + NomAlt
	DeNomAlt = DeNom * x_int + DeNomAlt

	If Abs(Start - (Nom / DeNom)) > Abs(Start - (NomAlt / DeNomAlt))
		;Debug "Solution 2 prefered..."
		Nom = NomAlt
		DeNom = DeNomAlt
	EndIf

	;Debug Str(Nom)+"/"+Str(DeNom)+"  error : "+StrD(Start - (Nom / DeNom))
	ProcedureReturn Str(Nom)+" / "+Str(DeNom) + " " + StrD(Start - Nom/DeNom, 15)
	
EndProcedure

Procedure.s Iif(Test.b, ok$, nok$)
   
   If Test
      ProcedureReturn ok$
   EndIf
   ProcedureReturn nok$
   
EndProcedure

Define q.d = #PI; https://math.stackexchange.com/questions/3506435/best-possible-rational-approximation-of-pi
; https://mathworld.wolfram.com/PiApproximations.html
; Convergents of the pi continued fractions are the simplest approximants To pi. 
; The first few are given by 3, 22/7, 333/106, 355/113, 103993/33102, 104348/33215, ... (OEIS A002485 And A002486), 
; which are good To 0, 2, 4, 6, 9, 9, 9, 10, 11, 11, 12, 13, ... (OEIS A114526) decimal digits, respectively.

Define k.a, l.a, MaxDen.q, Diff0.d = 1.0, Diff1.d, Fract$

;Testing - You only expect Smaller or Equal.
For k = 1 To 6
   For l = 1 To 9
      MaxDen = l* Pow(10, k)
      Fract$ = Fraction(q, MaxDen)
      Diff1 = ValD(Right(Fract$, 18))
      Debug "MaxDen : " + MaxDen + "  Frac: "+ Fract$ + Iif(Bool(Abs(Diff1) <= Abs(Diff0)), " Smaller Or Equal", " NO GOOD!")
      Diff0 = Diff1
   Next l   
Next k