CSI T 2 0 1 3 Rational Arithm etic w ith Floating Point Vaclav - PowerPoint PPT Presentation

CSI T 2 0 1 3 Rational Arithm etic w ith Floating Point Vaclav Skala University of West Bohemia, Plzen, Czech Republic VSB Technical University, Ostrava, Czech Republic http://www.VaclavSkala.eu Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 1

CSI T 2 0 1 3 Plzen ( Pilsen) City Plzen is an old city [first records of Plzen castle 976] city of culture, industry, and brewery. City, where today’s beer fermentation process was invented that is why today’s beers are called Pilsner - world wide Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 2

CSI T 2 0 1 3 University of W est Bohem ia 1 7 5 3 0 students + 9 8 7 PhD students Com puter Science and Engineering Mathem atics (+ Geomatics) Physics Cybernetics Mechanics (Computational) • Over 5 0 % of income from research and application projects • NTIS project (investment of 64 mil. EUR) • 2 nd in the ranking of Czech technical / informatics faculties 2009, 2012 Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 3

CSI T 2 0 1 3 “Real science” in the XXI century Courtesy of Czech Film, Barrandov Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 4

CSI T 2 0 1 3 Num erical system s • Binary system is used nearly exclusively • Octal & hexadecimal representation is used • If we would be direct descendants of tetrapods – we would have a great advantage – “simple counting in hexadecimal system” Nam e Base Digits E m in E m ax BI NARY B 16 Half 2 10+1 − 14 15 B 32 Single 2 23+1 − 126 127 B 64 Double 2 52+1 − 1022 1023 B 128 Quad 2 112+1 − 16382 16383 DECI MAL D 32 10 7 − 95 96 D 64 10 16 − 383 384 D 128 10 34 − 6143 6144 IEEE 758-2008 standard Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 5

CSI T 2 0 1 3 Mathem atically perfect algorithm s fail due to instability Main issues • stability, robustness of algorithms • acceptable speed • linear speedup – results depends on HW, CPU …. parameters ! Num erical stability • limited precision of float / double • tests A ? B with floats if A = B then ….. else ….. ; if A = 0 then ….. else …. should be forbidden in programming languages • division operation should be removed or postponed to the last moment if possible - “blue screens”, system resets Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 6

CSI T 2 0 1 1 3 T Typical exam pl es of in stability y s in E 3 • inter rsection o of 2 lines E 2 or a E 3 • point t lies on a line in plane in Ax + + By + C C = 0 o or Ax + B By + Cz + D = 0 • detec ction if a a line inte ersects a a polygon n, touche es a vert tex or pa asses throu ugh Typical T problem m double x x = -1; d double p p = ….; d w hile ( x x < +1) w { { if (x == p) C Console.O Out.Write eLine(” * *** ”) x += = p; } } /* if p / p = 0.1 t then no output, if p = 0 0.25 the en expec cted outp put */ A Amman, Jorda an 2013 Vacla v Skala http://www.V VaclavSkala.e eu 7

CSI T 2 0 1 1 3 Delauna D ay triang gulation n & Voro onoi dia agram P Point ins ide of a circle giv ven by th hree poin nts – pro oblems w with mesh hing p points in regular rectangu ular grid . I It can be e seen th hat the D DT & VD is ver V ry sensit tive to a point p position change S ?? ? ?? ROBU USTNES A Amman, Jorda an 2013 Vacla v Skala http://www.V VaclavSkala.e eu 8

CSI T 2 0 1 3 Floating point 4 • Not all numbers are represented � � 1 � correctly 1 � 2 � 3 � 5 � 3 � • Logarithmic arithmetic … • Continuous fractions � � �3; 7,15,1,292,1,1,1,2,1,3,1 … � • Interval arithmetic Generally NOT valid x + y = [a + c, b + d] x = [ a , b ] identities due to limited x - y = [a - d, b - c] y = [ c , d ] precision x × y = [min(ac, ad, bc, bd), max(ac, ad, bc, bd)] • ��� � � � ��� � � � 1 x / y = [min(a/c, a/d, b/c, b/d), [ � � � � � ] max(a/c, a/d, b/c, b/d)] if y ≠ 0 • x � � y � � �� � ���� � �� Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 9

CSI T 2 0 1 3 Statem ents like if <float> = <float> then …. or if <float> ≠ <float> then …. should not be allow ed in languages Quadratic equation ���√� � ���� �� � � �� � � � 0 usually solved as � �,� � �� If � � � 4�� then � � ��� � ��������� � � 4�� �/2 � � � � � � � � � � � ⁄ to get m ore reliable results. Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 10

CSI T 2 0 1 3 at � � 77617 , � � 33096 Function value com putation ���, �� � 333.75� � � � � �11� � � � � � � � 121� � � 2� � 5.5� � � �/�2�� � � 6.33835 10 �� single precision � � 1,1726039400532 double precision � � 1,1726039400531786318588349045201838 extended precision The correct result is “somewhere” in the interval of [ �0,82739605994682136814116509547981629����, �0,82739605994682136814116509547981629����� Exact solution ���, �� � �2 � � 2� � � 54767 66192 Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 11

CSI T 2 0 1 3 � � � �� � �� � � Num erical com putations � Hilbert’s Matrix ����� � �� � � �� � ��1� ��� �� � � � 1� �� � � � 1 � �� � � � 1 � �� � � � 2 � �� � � 1 � � � � � � � 1.0E+19 1.0E+17 1.0E+15 1.0E+13 1.0E+11 1.0E+09 1.0E+07 Error 1.0E+05 1.0E+03 1.0E+01 1.0E-01 1.0E-03 ε 1.0E-05 1.0E-07 ε p 1.0E-09 ξ 1.0E-11 1.0E-13 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Order of the Hilbert matrix Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 12

CSI T 2 0 1 3 Projective Space X = [ X, Y ] T X ∈ E 2 w c P 2 D(P ) 2 p D( ρ x = [ x, y: w ] T x ∈ P 2 2 E 2 D(E ) w=1 c=1 x D(p) ρ Conversion : Y A X B X = x / w Y = y / w & w ≠ 0 x y a b (a) (b) If w = 0 then x represents “an ideal point” - a point in infinity, i.e. it is a directional vector. The Euclidean space E 2 is represented as a plane w = 1 . Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 13

CSI T 2 0 1 3 Points and vectors • Vectors are “freely m ovable” – not having a fixed position � � � �� � , � � : 0� � • Points are not “freely m ovable” – they are fixed to an origin of the current coordinate system � � � �� � , � � : � � � � and � � � �� � , � � : � � � � usually in textbooks � � � � � � 1 A vector � � � � � � � in the Euclidean coordinate system - CORRECT Horrible “construction” DO NOT USE I T – I T I S TOTALLY W RONG � � � � � � � � �� � � � � , � � � � � : � � � � � � � � �� � � � � , � � � � � : 1 � 1� � � �� � � � � , � � � � � : 0� � This was presented as “How a vector” is constructed in the projective space � � in a textbook!! W RONG, W RONG, W RONG This construction has been found in SW as well !! Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 14

CSI T 2 0 1 3 Points and vectors A vector given by two points in the projective space � � � � � � � � �� � � � � � � � � , � � � � � � � � � : � � � � � � This is the CORRECT SOLUTI ON , but what is the interpretation? A “difference” of coordinates of tw o points is a vector in the m athem atical m eaning and � � � � is a “scaling” factor actually In the projective representation (if the vector length matters) � � � � � � � � �� � � � � � � � � , � � � � � � � � � : � � � � � � � �� � � � � � � � � , � � � � � � � � � � : 0� � � � � � � � � We have to strictly distinguish if we are working with points, i.e. vector as a data structure represents the coordinates, or with a vector in the mathematical meaning stored in a vector data structure. VECTORS x FRAMES Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 15

CSI T 2 0 1 3 Duality For simplicity, let us w c 2 2 P D(P ) consider a line p defined p ρ D( ) as: 2 E 2 D(E ) w=1 c=1 aX + bY + c = 0 x D(p) ρ We can multiply it by Y A X B w ≠ 0 and we get: ax + by + cw = 0 x y a b p T x = 0 i.e. (a) (b) p = [ a, b: c ] T x = [ x, y: w ] T =[ wX, wY: w ] T A line p ∈ E 2 is actually a plane in the projective space P 2 (point [ 0,0:0 ] T excluded) Amman, Jordan 2013 Vaclav Skala http://www.VaclavSkala.eu 16