Math : Fraction, Transform 0.54 in 27/50

[Le Soldat Inconnu]

Hi,

I do algo to get fraction of number, like get 27/50 with value 0.54

Limit is for x/y, x and y have to be smallest than 10000 to don't have problem

If you arrive to find wrong result, said to me, i try it a lot of time to be sure it's work in all case. Thanks

Code: Select all

Procedure.s Fraction(Valeur.d)
	Protected NewList Reste.q()
	Protected Entier.q, Numerateur.q, Denominateur.q, Reel.d, Signe.q
	
	Reel.d = Abs(Valeur)
	If Reel <> Valeur
		Signe = -1
	Else
		Signe = 1
	EndIf
	Entier = IntQ(ValD(StrD(Reel, 4)))
	; Debug Reel
	; Debug Entier
	While ValD(StrD(Reel, 4)) <> Entier 
		AddElement(Reste())
		Reste() = Entier
		Reel = 1 / (Reel - Entier)
		Entier = IntQ(ValD(StrD(Reel, 4)))
		; Debug Reel
		; Debug Entier
		If ListSize(Reste()) > 32 Or Entier > 10000
			ProcedureReturn ""
		EndIf
	Wend
	
	If ListSize(Reste()) = 0
		ProcedureReturn ""
	EndIf
	
	; Debug ""
	
	Numerateur = Entier
	Denominateur = 1
	; Debug Str(Numerateur) + "/" + Str(Denominateur)
	
	Repeat 
		Swap Numerateur, Denominateur
		Numerateur = Reste() * Denominateur + Numerateur
		; Debug Str(Numerateur) + "/" + Str(Denominateur)
	Until PreviousElement(Reste()) = 0
	
	; Debug ""
	
	If ValD(StrD(Valeur - Signe * Numerateur / Denominateur, 12)) = 0 And Numerateur <= 10000 And Denominateur <= 10000 ; And Str(Denominateur) <> LSet("1", Len(Str(Denominateur)), "0")
		Debug Str(Signe * Numerateur) + "/" + Str(Denominateur)
		ProcedureReturn Str(Signe * Numerateur) + "/" + Str(Denominateur)
	Else
		ProcedureReturn ""
	EndIf
	
EndProcedure

Debug Fraction(0.154)

Speed test

Code: Select all

Procedure.s Fraction(Valeur.d)
	Protected NewList Reste.q()
	Protected Entier.q, Numerateur.q, Denominateur.q, Reel.d, Signe.q
	
	Reel.d = Abs(Valeur)
	If Reel <> Valeur
		Signe = -1
	Else
		Signe = 1
	EndIf
	Entier = IntQ(ValD(StrD(Reel, 4)))
	; Debug Reel
	; Debug Entier
	While ValD(StrD(Reel, 4)) <> Entier 
		AddElement(Reste())
		Reste() = Entier
		Reel = 1 / (Reel - Entier)
		Entier = IntQ(ValD(StrD(Reel, 4)))
		; Debug Reel
		; Debug Entier
		If ListSize(Reste()) > 32 Or Entier > 10000
			ProcedureReturn ""
		EndIf
	Wend
	
	If ListSize(Reste()) = 0
		ProcedureReturn ""
	EndIf
	
	; Debug ""
	
	Numerateur = Entier
	Denominateur = 1
	; Debug Str(Numerateur) + "/" + Str(Denominateur)
	
	Repeat 
		Swap Numerateur, Denominateur
		Numerateur = Reste() * Denominateur + Numerateur
		; Debug Str(Numerateur) + "/" + Str(Denominateur)
	Until PreviousElement(Reste()) = 0
	
	; Debug ""
	
	If ValD(StrD(Valeur - Signe * Numerateur / Denominateur, 12)) = 0 And Numerateur <= 10000 And Denominateur <= 10000 ; And Str(Denominateur) <> LSet("1", Len(Str(Denominateur)), "0")
		Debug Str(Signe * Numerateur) + "/" + Str(Denominateur)
		ProcedureReturn Str(Signe * Numerateur) + "/" + Str(Denominateur)
	Else
		ProcedureReturn ""
	EndIf
	
EndProcedure

Ok = 0
Max = 10000
Temps1 = ElapsedMilliseconds()
For n = 1 To Max
	Repeat
		Numerateur = Random(9999) + 1
		Denominateur = Random(9999) + 1
	Until Numerateur % Denominateur <> 0 ; Si la division n'est pas entiere
	Debug ""
	Debug Str(Numerateur) + "/" + Str(Denominateur)
	Reel.d = Numerateur / Denominateur
	Texte.s = Fraction(Reel)
	If Texte <> ""
		Ok + 1
	EndIf
Next
Temps2 = ElapsedMilliseconds()

MessageRequester("Fraction", "Result OK : " + Str(Ok) + "/" + Str(Max) + Chr(10) + "Time to convert " + Str(Max) + " values = " + Str(Temps2 - Temps1) + "ms" + Chr(10) + "Time to do 1 converion = " + StrD((Temps2 - Temps1)/Max, 3) + "ms")

[KJ67]

Thank you, this is useful.

[c4s]

So this is a conversion from a float number to a fraction. I don*'t think all those Str() and Val() are very performant. Look what I found: http://stackoverflow.com/questions/9572 ... -fractions
The first posted code looks really fast: http://www.ics.uci.edu/~eppstein/numth/frap.c

Maybe someone is able to translate this c source into PureBasic?

[drgolf]

Hello, here is the port of the c code :
http://www.ics.uci.edu/~eppstein/numth/frap.c

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

Dim m.i(2,2)
Define x.d,startx.d,maxden.i,ai.i,t.i
;
startx=0.54 ; float number to approximate
x=startx
maxden=50 ; maximum denominator
;
m(0,0)=1: m(1,1) = 1
m(0,1)=0: m(1,0) = 0
;
While (m(1,0) *  (Int(x) + m(1,1)) <= maxden) 
  ai=Int(x)
  
 	t = m(0,0) * ai + m(0,1)
 	m(0,1) = m(0,0)
 	m(0,0) = t
 	t = m(1,0) * ai + m(1,1)
 	m(1,1) = m(1,0)
 	m(1,0) = t
 	If x-ai=0
 	  Break ;     // AF: division by zero
 	  EndIf
 	  x = 1/(x - ai)
 	If(x>$7FFFFFFF) 
 	  Break;  // AF: representation failure
 	  EndIf
Wend 
;
;
Debug Str(m(0,0))+"/"+Str(m(1,0))+"  error : "+Str(startx - (m(0,0) / m(1,0)))
 	   
; 
;     /* now try other possibility */
     ai = (maxden - m(1,1)) / m(1,0)
     m(0,0) = m(0,0) * ai + m(0,1)
     m(1,0) = m(1,0) * ai + m(1,1)

Debug Str(m(0,0))+"/"+Str(m(1,0))+"  error : "+Str(startx - (m(0,0) / m(1,0)))



; /*
; ** find rational approximation To given real number
; ** David Eppstein / UC Irvine / 8 Aug 1993
; **
; ** With corrections from Arno Formella, May 2008
; **
; ** usage: a.out r d
; **   r is real number To approx
; **   d is the maximum denominator allowed
; **
; ** based on the theory of continued fractions
; ** If x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...)))
; ** then best approximation is found by truncating this series
; ** (With some adjustments in the last term).
; **
; ** Note the fraction can be recovered As the first column of the matrix
; **  ( a1 1 ) ( a2 1 ) ( a3 1 ) ...
; **  ( 1  0 ) ( 1  0 ) ( 1  0 )
; ** Instead of keeping the sequence of continued fraction terms,
; ** we just keep the last partial product of these matrices.
; */
; 
; #include <stdio.h>
; 
; main(ac, av)
; int ac;
; char ** av;
; {
;     double atof();
;     int atoi();
;     void exit();
; 
;     long m[2][2];
;     double x, startx;
;     long maxden;
;     long ai;
; 
;     /* Read command line arguments */
;     If (ac != 3) {
; 	fprintf(stderr, "usage: %s r d\n",av[0]);  // AF: argument missing
; 	exit(1);
;     }
;     startx = x = atof(av[1]);
;     maxden = atoi(av[2]);
; 
;     /* initialize matrix */
;     m[0][0] = m[1][1] = 1;
;     m[0][1] = m[1][0] = 0;
; 
;     /* loop finding terms Until denom gets too big */
;     While (m[1][0] *  ( ai = (long)x ) + m[1][1] <= maxden) {
; 	long t;
; 	t = m[0][0] * ai + m[0][1];
; 	m[0][1] = m[0][0];
; 	m[0][0] = t;
; 	t = m[1][0] * ai + m[1][1];
; 	m[1][1] = m[1][0];
; 	m[1][0] = t;
;         If(x==(double)ai) Break;     // AF: division by zero
; 	x = 1/(x - (double) ai);
;         If(x>(double)0x7FFFFFFF) Break;  // AF: representation failure
;     } 
; 
;     /* now remaining x is between 0 And 1/ai */
;     /* approx As either 0 Or 1/m where m is max that will fit in maxden */
;     /* first try zero */
;     printf("%ld/%ld, error = %e\n", m[0][0], m[1][0],
; 	   startx - ((double) m[0][0] / (double) m[1][0]));
; 
;     /* now try other possibility */
;     ai = (maxden - m[1][1]) / m[1][0];
;     m[0][0] = m[0][0] * ai + m[0][1];
;     m[1][0] = m[1][0] * ai + m[1][1];
;     printf("%ld/%ld, error = %e\n", m[0][0], m[1][0],
; 	   startx - ((double) m[0][0] / (double) m[1][0]));
; }
; 

[Le Soldat Inconnu]

Test of algo in C :
0.3155 = 11/35
hum, sorry, it's not good 11/35 = 0.3142857142857143

with mine
0.3155 = 631/2000

[c4s]

You have to set a higher value for "maxden". For example 10000 like in your code.
Then the result is 59/187 = 0,315508... Well yours is still better but not as fast.

http://www.webmath.com/dec2fract.html
Maybe this approach is the best one?

Btw: Thanks drgolf!

[Arctic Fox]

Here is another method posted by Xombie
http://www.purebasic.fr/english/viewtopic.php?t=18845

[Le Soldat Inconnu]

I do optimisation of my code (3 time more fast)

I give good result for all fraction [1 to 10000]/[1 to 10000] (My function always checks fraction found with original value)
To get fraction, it's take 0.050 ms on my computer so it's really fast

@C4S : I find nothing more fast than ValD(StrD(Value, NbDecimals)) to round value with X decimals, if somebody have idea to do this fastest

@Artic fox : i try but result don't look good with this method too

Bye

[Michael Vogel]

Hi, drgolf

your code looks fine and is quite fast, but which of the results should be taken? In the example below, the first line showing the famous 355/133 approximation for pi is the better result compared to the second line. Maybe an epsilon value could be added to break as soon a value is close enough to the given value...

Michael

Hello, here is the port of the c code :
http://www.ics.uci.edu/~eppstein/numth/frap.c

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

Dim m.i(2,2)
Define x.d,startx.d,maxden.i,ai.i,t.i
;
startx=#PI
x=startx
maxden=5000;

m(0,0)=1: m(1,1) = 1
m(0,1)=0: m(1,0) = 0
;
While (m(1,0) *  (Int(x) + m(1,1)) <= maxden) 
  ai=Int(x)
  
 	t = m(0,0) * ai + m(0,1)
 	m(0,1) = m(0,0)
 	m(0,0) = t
 	t = m(1,0) * ai + m(1,1)
 	m(1,1) = m(1,0)
 	m(1,0) = t
 	If x-ai=0
 	  Break ;     // AF: division by zero
 	  EndIf
 	  x = 1/(x - ai)
 	If(x>$7FFFFFFF) 
 	  Break;  // AF: representation failure
 	  EndIf
Wend 
;
;
Debug Str(m(0,0))+"/"+Str(m(1,0))+"  error : "+StrD(startx - (m(0,0) / m(1,0)))
 	   
; 
;     /* now try other possibility */
     ai = (maxden - m(1,1)) / m(1,0)
     m(0,0) = m(0,0) * ai + m(0,1)
     m(1,0) = m(1,0) * ai + m(1,1)

Debug Str(m(0,0))+"/"+Str(m(1,0))+"  error : "+StrD(startx - (m(0,0) / m(1,0)))

[Trond]

I can't get the code from the original post to work properly. With some numbers it returns no result at all. Try #pi, 1, 2, 3, 1.45345, they all give an empty result.

This is my try. It's really inaccurate, but also small and fast:

Code: Select all

Procedure.s ToFraction(V.d, Precision = 1000)
  Protected T.q, N.q, I.q
  T = Int(V)*Precision + Round((V-Int(V))*Precision, #PB_Round_Nearest)
  N = Precision
  
  I = 1
  Repeat
    I + 1
    If T % I = 0 And N % I = 0
      T / I
      N / I
      I = 1
    EndIf
  Until I >= T Or I >= N
  
  ProcedureReturn Str(T) + "/" + Str(N)
EndProcedure

[Le Soldat Inconnu]

@Trond :

Explication
#Pi no result because no fraction
1, 2, 3, 4, etc... no result because there no fraction for this number
1.45345, the fraction is 29069/20000.I increase my system to get value to [1 to 100000]/[1 to 100000]. It's don't give result because fraction have to big value and i put limit

For your code, i try "ToFraction(647/845, 10000)" and it's give back 3871/5000
We have to be carreful with Round( because it don't use double but only float

here the code with increase limit to 100000

Code: Select all

Procedure.s Fraction(Valeur.d)
	Protected NewList Reste.q()
	Protected Entier.q, Numerateur.q, Denominateur.q, Reel.d, Signe.q
	
	Reel.d = Abs(Valeur)
	If Reel <> Valeur
		Signe = -1
	Else
		Signe = 1
	EndIf
	Entier = IntQ(ValD(StrD(Reel, 5)))
	; Debug Reel
	; Debug Entier
	While ValD(StrD(Reel, 5)) <> Entier 
		AddElement(Reste())
		Reste() = Entier
		Reel = 1 / (Reel - Entier)
		Entier = IntQ(ValD(StrD(Reel, 5)))
		; Debug Reel
		; Debug Entier
		If ListSize(Reste()) > 32 Or Entier > 100000
			ProcedureReturn ""
		EndIf
	Wend
	
	If ListSize(Reste()) = 0
		ProcedureReturn ""
	EndIf
	
	; Debug ""
	
	Numerateur = Entier
	Denominateur = 1
	; Debug Str(Numerateur) + "/" + Str(Denominateur)
	
	Repeat 
		Swap Numerateur, Denominateur
		Numerateur = Reste() * Denominateur + Numerateur
		; Debug Str(Numerateur) + "/" + Str(Denominateur)
	Until PreviousElement(Reste()) = 0
	
	; Debug ""
	
	If ValD(StrD(Valeur - Signe * Numerateur / Denominateur, 12)) = 0 And Numerateur <= 100000 And Denominateur <= 100000 ; And Str(Denominateur) <> LSet("1", Len(Str(Denominateur)), "0")
		Debug Str(Signe * Numerateur) + "/" + Str(Denominateur)
		ProcedureReturn Str(Signe * Numerateur) + "/" + Str(Denominateur)
	Else
		ProcedureReturn ""
	EndIf
	
EndProcedure

Debug Fraction(1.45345)

[Arctic Fox]

1, 2, 3, 4, etc... no result because there no fraction for this number

1/1, 2/1, 3/1, 4/1 ... right?

[Michael Vogel]

1, 2, 3, 4, etc... no result because there no fraction for this number

1/1, 2/1, 3/1, 4/1 ... right?

pi (\/¯2 etc.) cannot be represented by a fraction, but the rounded value(s) for #PI (and Sqr(2) etc.) can, of course!

Converting the floating point values of a computer to a fraction means you should beware of approximations because you never handle with exact values (1/10 have also be rounded in the binary system), so that's why such a routine should also find a nice result for #PI as well

Dozens of years ago I have done such a routine for my HP41 (a calculator ) maybe I can find the program to convert it to PB (a totally new programming technique ) – or soething like that can be also found on sites, like hpcalc.org etc.

Michael

[Trond]

1, 2, 3, 4, etc... no result because there no fraction for this number

1/1, 2/1, 3/1, 4/1 ... right?

Yes.

[Trond]

More accurate, but not very fast:

Code: Select all

Procedure.s ToFraction2(V.d, MaxDenominator = 10000)
  Protected MinN, MaxN, BestN, Add
  Protected MinDiff.d = 100
  Protected MaxDiff.d = 0
  
  MaxDenominator - 1
  For N = 1 To MaxDenominator
    Diff.d = Abs(V*N-Int(V*N))
    If Diff < MinDiff
      MinDiff = Diff
      MinN = N
      If Diff < 0.000001
        Break
      EndIf
    EndIf
    If Diff > MaxDiff
      MaxDiff = Diff
      MaxN = N
      If Diff > 0.99999
        Break
      EndIf
    EndIf
  Next
  If MinDiff < 1-MaxDiff
    BestN = MinN
  Else
    BestN = MaxN
    Add = 1
  EndIf
  
  Debug ""
  Debug StrD(V)
  Debug StrD(Int(V*BestN+Add)/BestN)
  
  ProcedureReturn Str(Int(V*BestN+Add)) + "/" + Str(BestN)
EndProcedure


Debug ToFraction2(1)
Debug ToFraction2(2)
Debug ToFraction2(2.5)
Debug ToFraction2(2.75)
Debug ToFraction2(#PI, 30)
Debug ToFraction2(#PI, 200)
Debug ToFraction2(1.22)
Debug ToFraction2(1.45345)