Enhance precision (without extra lib)

Just starting out? Need help? Post your questions and find answers here.
User avatar
Michael Vogel
Addict
Addict
Posts: 2680
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Enhance precision (without extra lib)

Post by Michael Vogel »

When calculating some triangle, I got wrong values under some circumstances. Does anyone know if certain math functions in Purebasic have a higher or lower "quality"?

Problems in hab, hac and so on...

Code: Select all

Define.d a,b,c,s
Define.d area
Define.d ha,hb,hc
Define.d wa,wb,wc

a=223.606797749979
b=200
c=100

s=(a+b+c)/2
area=Sqr(s*(s-a)*(s-b)*(s-c))
ha=area*2/a
hb=area*2/b
hc=area*2/c

wa=Degree(ACos((a*a-b*b-c*c)/(-2*b*c)))
wb=Degree(ACos((b*b-c*c-a*a)/(-2*c*a)))
wc=Degree(ACos((c*c-a*a-b*b)/(-2*a*b)))

hac=Sqr(b*b-ha*ha)
hab=Sqr(c*c-ha*ha)
hba=Sqr(c*c-hb*hb)
hbc=Sqr(a*a-hb*hb)
hcb=Sqr(a*a-hc*hc)
hca=Sqr(b*b-hc*hc)

ShowVariableViewer()

OpenWindow(0,0,0,0,0,"",#PB_Window_Invisible)
Repeat : Until WaitWindowEvent()=#PB_Event_CloseWindow
superadnim
Enthusiast
Enthusiast
Posts: 480
Joined: Thu Jul 27, 2006 4:06 am

Re: Enhance precision (without extra lib)

Post by superadnim »

hi you should try using EnableExplicit
the compiler is defining hac, hab, hba, hbc, hcb and hca as integer by default because no type has been defined

for higher resolution computation ive always been told bigint is the key so integer math is more accurate but for most uses or purposes double floating point are enough.

the variable viewer has type icons if you look at the variables that are not working for you they show I in red instead of D in green.

now for some reason i get NaN so i copied a sqrt function off the forum but it is much slower

Code: Select all

EnableExplicit

Procedure.d Sqrt(n.d)
  ProcedureReturn Pow(10,Log10(n)/2.0)
EndProcedure

Define.d a,b,c,s
Define.d area
Define.d ha,hb,hc
Define.d wa,wb,wc

Define.d hac, hab, hba, hbc, hcb, hca

a=223.606797749979
b=200
c=100

s=(a+b+c)/2
area=Sqrt(s*(s-a)*(s-b)*(s-c))
ha=area*2/a
hb=area*2/b
hc=area*2/c

wa=Degree(ACos((a*a-b*b-c*c)/(-2*b*c)))
wb=Degree(ACos((b*b-c*c-a*a)/(-2*c*a)))
wc=Degree(ACos((c*c-a*a-b*b)/(-2*a*b)))

hac=Sqrt(b*b-ha*ha)
hab=Sqrt(c*c-ha*ha)
hba=Sqrt(c*c-hb*hb)
hbc=Sqrt(a*a-hb*hb)
hcb=Sqrt(a*a-hc*hc)
hca=Sqrt(b*b-hc*hc)

ShowVariableViewer()

OpenWindow(0,0,0,0,0,"",#PB_Window_Invisible)
Repeat : Until WaitWindowEvent()=#PB_Event_CloseWindow
sorry if i cant help more than this

:lol: should I bash the keyboard and give up?
:?
User avatar
Michael Vogel
Addict
Addict
Posts: 2680
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Re: Enhance precision (without extra lib)

Post by Michael Vogel »

What a shame, the missing definition has been because of my 'optimization' of the code snippet :oops:

The NaN error appears because the square root of a negative number (because the triangle height gets a value higher than the side)...
...so sqrt is doing it''s job in that situation.

Reduced the snippet to just it's first two calculation steps...

Code: Select all

Define.d a,b,c,s
Define.d area
Define.d ha,hb,hc
Define.d hab,hba,hca
Define.d hac,hbc,hcb
Define.d wa,wb,wc
Define.d z

Procedure.d Sqrt(n.d)
	ProcedureReturn Pow(10,Log10(n)/2.0)
EndProcedure
Procedure.d Frac(innumber.d)
	ProcedureReturn innumber - Int (innumber)
EndProcedure
Procedure.d Mul(a.d,b.d)
	ProcedureReturn Int(a)*Int(b)+frac(a)*frac(b)+Int(a)*frac(b)+frac(a)*Int(b)
EndProcedure
Procedure Check(d.d)
	
	#Show=15
	
	Debug StrD(d,#show)+"  sqr:"+StrD(Sqr(d),#show)+"  sqrt:"+StrD(sqrt(d),#show)
	
EndProcedure

a=223.606797749979
b=200
c=100

s=(a+b+c)/2
check(s*(s-a)*(s-b)*(s-c))
Check((s*s*s*s-s*a*s*s-s*s*b*s+s*a*b*s)-(s*s*s*c-s*a*s*c-s*s*b*c+s*a*b*c))
Check(Mul(Mul(s,(s-a)),mul((s-b),(s-c))))
End
User avatar
Michael Vogel
Addict
Addict
Posts: 2680
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Re: Enhance precision (without extra lib)

Post by Michael Vogel »

Hm, thought I get rid of some problems when slpitting the integer and fraction part so I would be able do all calculations with integers...
...but the precision got already lost when setting the fraction part of a variable...

Code: Select all

#Precision=10000000000000000
#Precision_=100000000

Structure qdot
	number.q
	fraction.q
EndStructure

Procedure qSet(*qdot.qdot,value.d)

	*qdot\number=Int(value)
	*qdot\fraction=(value-*qdot\number)*#precision

EndProcedure
Procedure.d qGet(*qdot.qdot)

	ProcedureReturn *qdot\number+(*qdot\fraction/#Precision)

EndProcedure

Procedure qAdd(*result.qdot,*qa.qdot,*qb.qdot)
	
	*result\fraction=*qa\fraction+*qb\fraction
	*result\number=*qa\number+*qb\number+*result\fraction/#Precision
	*result\fraction%#Precision
	; ProcedureReturn qGet(*result)

EndProcedure
Procedure qSub(*result.qdot,*qa.qdot,*qb.qdot)
	
	*result\fraction=*qa\fraction-*qb\fraction
	*result\number=*qa\number-*qb\number+*result\fraction/#Precision
	*result\fraction%#Precision
	; ProcedureReturn qGet(*result)

EndProcedure
Procedure qMul(*result.qdot,*qa.qdot,*qb.qdot)
	
	Protected.q mul,high,low
	
	*result\number=*qa\number * *qb\number
	*result\fraction=0
	
	Debug "a "+Str(*qa\number)+" . "+RSet(Str(*qa\fraction),16,"0")
	Debug "a "+Str(*qb\number)+" . "+RSet(Str(*qb\fraction),16,"0")
	
	mul=(*qa\fraction/#Precision_)*(*qb\fraction/#Precision_)
	*result\number+mul/#Precision
	*result\fraction+mul%#Precision
	Debug "+ "+Str(mul%#Precision)
	Debug "= "+Str(*result\fraction)
	
	mul=*qa\number * *qb\fraction/#Precision_
	*result\number+mul/#Precision_
	*result\fraction+(mul%#Precision_)*#Precision_
	Debug "+ "+Str(mul%#Precision_)
	Debug "= "+Str(*result\fraction)
	
	mul=*qb\number * *qa\fraction/#Precision_
	*result\number+mul/#Precision_
	*result\fraction+(mul%#Precision_)*#Precision_
	Debug "+ "+Str(mul%#Precision_)
	Debug "= "+Str(*result\fraction)
	
	
EndProcedure

q.qdot
qa.qdot
qb.qdot
qc.qdot
qSet(qa,223.606797749979)
qSet(qb,200)
qSet(qc,100)
s.d
qAdd(q,qa,qb)
qAdd(q,q,qc)

a.d=3.6
b.d=4.8
qSet(qa,a)
qSet(qb,b)
qMul(q,qa,qb)
Debug a*b
Debug qGet(q)

Maybe I should include GMP, does anyone has actual lib/dll files which could be included to a Purebasic program to have a create a single exe file?
User avatar
STARGÅTE
Addict
Addict
Posts: 2090
Joined: Thu Jan 10, 2008 1:30 pm
Location: Germany, Glienicke
Contact:

Re: Enhance precision (without extra lib)

Post by STARGÅTE »

What exactly is your idea?
PB 6.01 ― Win 10, 21H2 ― Ryzen 9 3900X, 32 GB ― NVIDIA GeForce RTX 3080 ― Vivaldi 6.0 ― www.unionbytes.de
Lizard - Script language for symbolic calculations and moreTypeface - Sprite-based font include/module
User avatar
Michael Vogel
Addict
Addict
Posts: 2680
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Re: Enhance precision (without extra lib)

Post by Michael Vogel »

I've written some programs based on floating point arithmetics (GPS tools and also some for astronomy and math). Actually I needed a quick solution to be able to check distance measurement data - so I wrote the small program triangle this week...

Screenshot

...which works fine until a value will have been entered where some trivial operations lead to a unprecise (and illegal) result. If you'd like to have a look my program, you only have to press F1 to change alpha to 90° and this will lead to the result as seen in the first (better third) post here...

...no problem to "repair" this special case in that program, but it's a workaround and I am unsure if there are more "bad" numbers to be aware of and if my other tools are also at risk :?
User avatar
STARGÅTE
Addict
Addict
Posts: 2090
Joined: Thu Jan 10, 2008 1:30 pm
Location: Germany, Glienicke
Contact:

Re: Enhance precision (without extra lib)

Post by STARGÅTE »

I can not start the programm, it crashes.

What do I have to see in the Screenshot?
And what is the message of the code in post 3?
PB 6.01 ― Win 10, 21H2 ― Ryzen 9 3900X, 32 GB ― NVIDIA GeForce RTX 3080 ― Vivaldi 6.0 ― www.unionbytes.de
Lizard - Script language for symbolic calculations and moreTypeface - Sprite-based font include/module
User avatar
Michael Vogel
Addict
Addict
Posts: 2680
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Re: Enhance precision (without extra lib)

Post by Michael Vogel »

STARGÅTE wrote:I can not start the programm, it crashes.
Absolutely weird, there's nothing special inside and it works fine here - double checked it that moment (Win 8.1 64Bit)
STARGÅTE wrote:What do I have to see in the Screenshot?
That's how it look like with "working" numbers...
Of course, wrong results won't look like that, here's the screenshot when using an 90° angle: Image
STARGÅTE wrote:And what is the message of the code in post 3?
It should demonstrate that the trivial calculation of 's' can be transformed and gives different answers, but all of them show wrong digits.
User avatar
STARGÅTE
Addict
Addict
Posts: 2090
Joined: Thu Jan 10, 2008 1:30 pm
Location: Germany, Glienicke
Contact:

Re: Enhance precision (without extra lib)

Post by STARGÅTE »

>> Absolutely weird, there's nothing special inside and it works fine here - double checked it that moment (Win 8.1 64Bit)
Now it works. Don't know what happens.

>> It should demonstrate that the trivial calculation of 's' can be transformed and gives different answers, but all of them show wrong digits.
Of cause. Even with GMP you would end up with such observation. The digit position depends on your set precision.
It is a limitation of arbitrary but still limited precision of real numbers.
Especially for irrational numbers like many root-numbers.

The problem I see is, that you use in the first post a function like ACos() which is only valid between -1 and 1.
Since multiply floating point arithmetics with valid a, b and c (for a triangle) can also lead to numbers like 1.00001
Therefore I use stable functions like:

Code: Select all

Procedure.d ArcSin(Value.d)
	If Value < -1.0
		ProcedureReturn -#PI/2
	ElseIf Value > 1.0
		ProcedureReturn #PI/2
	Else
		ProcedureReturn ASin(Value)
	EndIf
EndProcedure

Procedure.d ArcCos(Value.d)
	If Value < -1.0
		ProcedureReturn #PI/2
	ElseIf Value > 1.0
		ProcedureReturn 0
	Else
		ProcedureReturn ACos(Value)
	EndIf
EndProcedure

Procedure.d Sqrt(Value.d)
	If Value < 0.0
		ProcedureReturn 0.0
	Else
		ProcedureReturn Sqr(Value)
	EndIf
EndProcedure

; Because 0.85 != (1.0-0.15) in double precision.

Debug Degree(ASin(0.85/(1.0-0.15)))
Debug Degree(ArcSin(0.85/(1.0-0.15)))

Debug Degree(ACos(0.85/(1.0-0.15)))
Debug Degree(ArcCos(0.85/(1.0-0.15)))

Debug Sqr((1.0-0.15)-0.85)
Debug Sqrt((1.0-0.15)-0.85)
Btw: The GMP Lib has no build-in functions like cos() or sin().
PB 6.01 ― Win 10, 21H2 ― Ryzen 9 3900X, 32 GB ― NVIDIA GeForce RTX 3080 ― Vivaldi 6.0 ― www.unionbytes.de
Lizard - Script language for symbolic calculations and moreTypeface - Sprite-based font include/module
User avatar
Michael Vogel
Addict
Addict
Posts: 2680
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Re: Enhance precision (without extra lib)

Post by Michael Vogel »

Thanks for investing your time in this story and you're right, that innacurate values are leading to illegal results.

As far as I remember - had math at the university some hundred years ago - there shouldn't be a legal triangle producing a value <1 or >1 for the term u=(a*a-b*b-c*c)/(-2*b*c) within Acos(u) as it could be written as (a/b*a/c-b/c-c/b)/2...

However, I am happy to have another "repair kit" to "clean" illegal values. I am doing so with lines like hac=Sqr(Max(b*b-ha*ha,0)) because the heights get (slightly) higher values as their corresponding sides, which is impossible of course.

Played around a while with different approachs to create "better" math functions for add/sub/mul/div and sqrt but in fact the original routines are not that bad, it's a kind of bad luck that numerous errors never get wiped out (i fact they are increasing with each calculation).

Code: Select all


Procedure error(value.d,info.s)
	Protected s.s
	Protected v.d
	Protected e.q
	v=value
	If v<0
		v=-v
	EndIf
	If v>1
		v-Int(v)
	EndIf
	s=Mid(StrD(v,20),5,18)
	e=1000000000000000000-Val(s)
	Debug ".."+Mid(StrD(v,2),2,3)+" : "+info+RSet(Str(Abs(e)),10)+" delta"	
EndProcedure

a.d=4.6
b.d=3.8
error(a,"a")
error(b,"b")
error(a-b,"-")
error(b-a,"-")
error(a+b,"+")
error((b+a)*(a+b),"*")
User avatar
Michael Vogel
Addict
Addict
Posts: 2680
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Re: Enhance precision (without extra lib)

Post by Michael Vogel »

Just if anyone who can't remember how a triangle looks like, I've finished this windows only program...
...have fun.
Post Reply