This document is an edited reprint of the paper What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg, publishe

Transcription

1 모든컴퓨터과학자가 알아야할 부동소수점의모든것 What Every Computer Scientist Should Know About Floating-Point Arithmetic Jaebum Lee J B i

2 This document is an edited reprint of the paper What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg, published in the March 1991 issue of Comptuer Surveys. Copyright 1991, Association for Computing Machinery, Inc., reprinted by permission. This document is translated by Jaebum Lee, with many other contributors, in a part of open-computer-book distribution project. 이문서의모든내용은 ncg_goldberg.html 에공개되어있는모든컴퓨터과학자가알아야하는부동소수점연산의모든것 (What Every Computer Scientist Should Know About Floating-Point Arithmetic) 을번역한것입니다. 양질의무료컴퓨터관련문서를한국어로공급하기위한프로젝트의일환으로, 이문서는 tistory.com 에서무료로배포되고있습니다. ( 오픈북프로젝트에관해자세한내용은 을참조하시기바랍니다 ) 문서전체번역은이재범 (kev0960@gmail.com) 외수많은사람들이수정에도움을주셨습니다. 원문서에서의주석은모두각주 (footnote) 로달았고, 특별한설명이더필요하다고생각된부분은역자주를책여백 (marginnote) 에달았습니다. 또한, 마지막의참고문헌목록은생략하였으니, 혹시참고문헌을찾고싶은분들은원문을참조하시기바랍니다. 부동소수점연산 (Floating-point arithmetic) 은소수의사람들만이다루는주제로여겨져왔습니다. 이는많은컴퓨터시스템상에서부동소수점연산이사용된다는사실을생각해보았을때매우놀랄만한일인데요, 거의모든컴퓨터언어에서는부동소수점데이터타입이있고, PC 를비롯한슈퍼컴퓨터들에까지거의모든장치에는부동소수점가속기 (floating-point accelerator) 들이있습니다. 많은컴파일러들은부동소수점알고리즘을컴파일하는데사용될것이며, 사실상모든운영체제에서는오버플로우 (overflow) 와같은부동소수점예외 (exception) 를처리할수있어야만합니다. 이문서에서는부동소수점의여러특징중에서컴퓨터시스템의설계에가장직접적으로영향을주었던부분에대해다룰것입니다. 일단, 그시작으로, 부동소수점표기법의설립배경과, 반올림오차, 그리고 IEEE 부동소수점표준에대해다룰것이며, 끝부분에서는어떻게컴퓨터제작자들이부동소수점을더잘지원할수있는지에대해설명하도록할것입니다. 중요키워드 : 비정규화수 (Denormalized number), 예외 (Exception), 부동소수점 (Floating-point), 부동소수점표준 (Floating-point standard), 점진적언더플로우 (Gradual underflow), 보호숫자 (Guard digit), NaN, 오버플로우 (Overflow), 상대오차 (Relative error), 반올림오차 (Rounding error), 반올림모드 (Rounding mode), ulp, 언더플로우 (Underflow) ii

3 차례 차례 들어가며 iii v 제 1 장 반올림오차 (Rounding Error) 부동소수점형식 (Floating-point Formats) 상대오차와 Ulp 보호숫자 (Guard Digits) 지워짐 (Cancellation) 정확한반올림연산 제 2 장 IEEE 표준 형식과연산들 밑수 (base) 정밀도 (Precision) 지수 연산들 (operations) 특별한값들 NaN 들 무한대 부호있는영 (signed zero) 정규화되지않는수 예외, 플래그, 예외처리기 (Exceptions, Flags, Trap handlers) 예외처리기 (Trap handler) 반올림방식들 (Rounding Modes) 플래그 (Flag) 제 3 장 시스템적측면 명령어집합 (instruction sets) 언어와컴파일러 (Languages and Compilers) IEEE 표준 iii

4 3.2 최적화기 처리 (Handling) 제 4 장자세한증명들 반올림오차 이진수를십진수로변환하기 덧셈에서의오차 제 5 장끝마치며 감사의말 iv

5 들어가며 컴퓨터시스템제작자들은종종부동소수점연산에대한지식을요구로합니다. 하지만, 놀랍게도부동소수점에대해자세하게설명한것들은극히드물지요. 심지어부동소수점에대해다루는몇안되는책들중하나인 Pat Sterbenz 의 Floating-Point Computation 은절판된지오래입니다. 이문서는부동소수점연산의여러특징중에서컴퓨터시스템설계에직접적으로영향을주는것들에대해다루고자합니다. 이를위해크게 3 개의부분으로문서를나누었는데, 첫번째반올림오차에선부동소수점덧셈, 뺄셈, 곱셈, 나눗셈시어떠한방식으로반올림을하느냐에따라달라지는양상을살펴볼것입니다. 또한, 반올림오차, ulps, 그리고상대오차를계산하는두가지방식을다룰것입니다. 두번째 IEEE 표준에선 IEEE 부동소수점표준에대해이야기하는데, 이는상업하드웨어업체들에게빠르게채택되고있는방식중하나입니다. 세번째시스템적측면에서는부동소수점을컴퓨터시스템설계측면에서살펴볼것인데, 예를들어명령어세트 (instruction set) 설계, 컴파일러최적화, 예외처리등에대해다룰것입니다. 참고로, 이문서에서저는간단한산수로쉽게보일수있는것들에대해서는증명하지않고넘어가도록하였습니다. 특히모든 ( 단순한 ) 문장에대해시시콜콜하게집고넘어가는것은불필요하다고생각되기때문이지요. 글의논지에서벗어나는자질구레한증명들은뒤에자세한것들을참고하시기바랍니다. 참고로이장에서등장하는여러정리들의증명뒤에는 로표시되어있는데, 증명을생략했을경우그냥정리만써놓았습니다. v

6

7 제 1 장 반올림오차 (Rounding Error) 무한히많은실수들을유한개의비트로표현하기위해서는근사적으로표현해야합니다. 비록정수개수도무한하다해도, 대부분의프로그램에서사용하는스케일의정수들은 32 비트이내로저장할수있습니다. 반면에, 실수간의연산에서는 32 비트이내로계산결과를정확하게표현할수있는경우가드뭅니다. 따라서, 부동소수점연산의결과는반드시반올림을통해유한비트표현방식으로나타내져야만합니다. 그렇기때문에반올림오차는부동소수점연산에서의어쩔수없는특징이기도하지요. 상대오차와 Ulp 에서는어떻게반올림오차를측정하는지에대해설명할것입니다. 어차피대부분의부동소수점연산에서반올림오차가발생하므로, 산술연산작업에서필요이상보다살짝더발생하는오차에대해관심을가져야할까요? 보호숫자 (Guard Digits) 에서는두인접한수를뺄때발생하는오차를줄이기위한방법으로도입된보호숫자라는개념에대해설명할것입니다. 보호숫자라는개념은 IBM 에서중요하게생각해온개념으로 1968년에 IBM 의 System/360 아키텍쳐의배정밀도 (double precision) 형식에추가된이후로모든장치들에게도입되었습니다. 이보호숫자의효용성에대해설명하기위해 2 가지예제를살펴볼것입니다. IEEE 표준에서는단순히보호숫자를필수적으로사용해야한다는것보다더많은것을요구합니다. IEEE 표준에서는덧셈, 뺄셈, 곱셈, 나눗셈, 그리고제곱근연산에대한알고리즘을제공하며, 이표준에따라만들어진장치에서는연산시제공된알고리즘을사용한결과와같은결과가나오도록요구합니다. 따라서, 어떤프로그램이다른장치로옮겨졌을때, 만일그장치가 IEEE 표준을지원하는장치라면, 기본연산결과는장치에관계없이모두같을것입니다. 따라서 IEEE 표준을지원하는장치끼리매우쉽게프로그램을포팅 (porting) 할수있게되지요. 더자세한사항은정확히반올림되는연산들섹션에서살펴볼것입니다. 1.1 부동소수점형식 (Floating-point Formats) 그동안실수를표현하는방식으로꽤많은종류의형태가제안되었지만, 현재까지 가장많이사용하는방식은부동소수점표현입니다. 1) 부동소수점방식은 ( 항상 부동소수점이란, 말그대로소수점의위치가고정된것이아닌, 단어그래로둥둥떠다니며위치가바뀐다는의미입니다 1

8 짝수라고가정되는 ) 밑수 (base) β 와정밀도 (precision) p 를인자로가지며표현하게됩니다. 예를들어서 β = 10 이고 p = 3 이라면, 0.1 을표현했을때 이렇게표기할수있게됩니다. 만일 β = 2 이고 p = 24 라면, 0.1 은정확하게표현할수없고, 단순히근사적으로 로표시할수밖에없습니다. 일반적으로, 부동소수점수는 ±d.dd...d β e 로표현되며, d.dd...d 부분의경우가수 (significand) 라고불리며 p 자리수입니다. 즉, 정확히말하면 ±d 0.d 1 d 2...d p 1 β e 로나타내겠지요. ± (d 0 + d 1 β d p 1 β (p 1) )β e, (0 d i β) (1.1) 여기서부동소수점수라는단어는, 식 1.1 의형태로정확히표현되는수를지칭하는말로사용하도록하겠습니다. 부동소수점의또다른중요한인자로지수의크기범위를나타내는값들이있는데, 이를각각 θ max 와 θ min 이라고지칭합니다. 따라서 β p 개의가능한가수들과, θ max θ min + 1 개의가능한지수들을표현하기하기위해서는총 [log 2 (e max e min + 1)] + [log 2 (β p )] + 1 비트가필요로하며, 끝에붙은 +1 은부호비트를위해필요한추가적인비트를나타냅니다. 모든실수가부동소수점형식으로표현될수없는이유는크게두가지가있는데, 가장흔한이유로, 0.1 을예로들을수있습니다. 비록 0.1 은십진수로유한개의자리수로표현할수있지만, 이진법으로표현시에는무한개의자리수가필요하게됩니다. 따라서 β = 2 일경우, 0.1 은두개의부동소수점수사이에위치하게되지만, 그들둘다 0.1 은아닙니다. 다른경우로는부동소수점으로표현가능한범위를벗어난경우인데, 만일절대값이 β β emax 보다크거나, 1.0 β e min 보다작을경우에부동소수점으로표현할수없게됩니다. 이문서에서는후자보다는전자로인해표현될수없는실수들에대해다룰것이지만, 후자의경우에대해서는무한대와정규화되지않는수을참고하시기바랍니다. 부동소수점표현은반드시유일할필요는없습니다. 예를들어서 과 은모두 0.1 을표현하지요. 만일맨앞자리가 0 이아니라면 ( 즉, 식 1.1 에서 d 0 0 ) 부동소수점표현이정규화 (normalized) 되었다고말합니다. 따라서, 는정규화되있는것이고, 는정규화되어있지않은것입니다. β = 2, p = 3, θ min = 1, θ max = 2 인경우, 그림 1.1 에서나타난것처럼 16개의정규화된부동소수점이가능합니다. 하지만, 이렇게정규화된부동소수점들만을사용하기로제한하게된다면, 0 을표현하지못한다는문제점이생기게됩니다. 이러한문제를해결하기위한가장자연스러운방법은, 1.0 β e min 로표현하는것인데, 이러한 1) 그외의제안된방식으로는 floating slash 와 signed logarithm 등의표현방식이있습니다. [Matula and Kornerup 1985; Swartzlander and Alexopoulos 1975]. 2

9 그림 1.1: β = 2, p = 3, e min = 1, e max = 2 일때의정규화된수들 방법을사용할경우만일지수를 k 비트로표현한다면총 2 k 1 만큼의수만을지수로사용할수있게됩니다 ( 왜냐하면나머지하나는 0 을위해남겨놓았으므로 ) 참고로부동소수점에서의 는부동소수점표기중일부지, 실제로부동소수점수의곱셈연산을의미하는것이아닙니다. 예를들어서 ( ) ( ) 에서는오직한번의부동소수점곱셈이발생하게됩니다. 1.2 상대오차와 Ulp 부동소수점연산시에반올림오차는필연적으로발생하는것이기때문에그오차를측정하는것은매우중요한일입니다. 일단여기서기본적으로 β = 10, p = 3 인부동소수점들을고려하도록합시다. 만일부동소수점연산의결과가 이고, 정확한결과는.0314 였다면, 마지막자리에서 2 의오차가있었음을알수있습니다. 마찬가지로만일정확한값은 인데, 부동소수점표현으로 로나타났다면, 마지막자리에서.159 의오차가발생하였음을알수있습니다. 따라서일반적으로, 만일실수 z 를부동소수점수 d.d...d β e 로근사시켰다면, 마지막자리에서발생하는오차는 d.d...d (z/β e ) β p 1 임을알수있습니다. 2) 3) Ulps 는 Units in the last place 의약자로마지막자리의단위를의미합니다. 부동소수점연산의결과가참값과최대한으로가깝다하더라도, 그오차는최대.5ulp 까지날수있습니다. 부동소수점수와참값과의차이를측정하는다른방식으로상대오차 (relative error) 가있는데, 이는단순히두수의차이를참값으로나눈것입니다. 예를들어, 를 으로근사하였다면, 상대오차의크기는.00159/ 가됩니다..5ulp 에해당하는오차의크기를계산하기위해서, 만일어떤실수가부동소수점수 d.dd...dd β e 로표현되었다고하면,.5ulp 의크기는 β β e 이됩니다. 여기서 β 은 β/2 입니다. 부동소수점수의가수가 p 자리라면, 그오차에는 p 개의 0 이오게됩니다. 따라서,.5ulp 에해당하는오차는 ((β/2)β p ) β e 가됩니다. d.dd...dd β e 꼴의모든숫자들이같은크기의절대오차 (absolute error) 을가지게되지만, 상대오차의경우 ((β/2)β p ) β e /β e 와 ((β/2)β p ) β e /β e+1 의범위를가지게됩니다. 2) 여기서 z 가 β e max + 1 보다크거나 β e min 보다작은경우는고려하지않습니다. 물론, 추후에따로명시하기전까지이러한범위에있는수들은고려하지않습니다. 3) z 을 z 를부동소수점수로근사한것이라가정하면, d.d...d (z/β e ) β p 1 는 z z /uls(z ) 이라한것과같습니다. 참고로오차를계산하는좀더정확한공식은 z z /uls(z) 입니다 - Ed. 3

10 워블 (Wobble) 의사전적 정의은불안정함, 떨림, 동요등입니다. 이말은, 1 2 β p 1 2 ulp β 2 β p (1.2) 따라서.5ulp 에대응하는상대오차의경우, β 의배수로값이바뀔수있으며, 이 인자를워블 (wobble) 이라부릅니다. 식에서의상한을 ϵ = (β/2)β p 로정의하면, 언제나어떤실수를가장가까운부동소수점수로반올림하였을때상대오차는항상 ϵ 이하며, 이를머신입실론 (machine epsilon) 이라부릅니다. 위예제에서상대오차는.00159/ 였는데, 이와같이상대오차의 크기는매우작으므로, 보통 ϵ 의배수로많이나타냅니다. 이경우 ϵ = (β/2)β p = 5(10) 3 =.005 이므로, 상대오차의경우 (.00159/ )/.005)ϵ 0.1ϵ 이됩니다. 상대오차와 ulps 의차이를설명하기위해실수 x = 를생각합시다. 이는 x = 로근사되며, 오차는 0.5ulps 이고, 상대오차는 0.8ϵ 입니다. 다음에, 만일 8 x 을계산한다면, 참값은 8x = 98.8 이지만, 계산된값은 8 x = 입니다. 이경우, 오차는 4.0ulps 가되었지만, 상대오차는 0.8ϵ 으로그대로입니다. ulps 로 측정된오차의크기는 8 배나커졌지만, 상대오차는 0.8ϵ 로그대로이지요. 일반적으로, 밑수가 β 일때, 상대오차를일정하게유지시킨다면 ulps 로측정된오차는 β 의배수로 늘어날수있지만, 역으로식에서도나타나듯이, ulps 오차를고정시킨다면상대오차는 β 의배수로늘어날수있습니다. 반올림으로발생한오차를측정하는가장자연스러운방법은 ulps 로측정하는 것입니다. 예를들어, 가장가까운부동소수점수로반올림할때발생하는오차는.5ulp 보다작거나같겠지요. 하지만, 여러가지연산들이포함되어있는식에서의반올림 오차를분석할때, 상대오차를측정하는것이좀더나은방법입니다. 반올림오차 에서, 여러식들을연산할때의상대오차를계산하는방법을살펴볼것입니다. ϵ 이 가장가까운부동소수점수로반올림하는효과를 β 의배수로뻥튀기할수있기 때문에작은 β 를사용하는것이오차를측정하는데좀더정확합니다. 만일반올림오차의차수 (order) 에만관심을가진다면 ulps 와 ϵ 은단순히 β 배 차이이기때문에서로바꿔가며사용할수있습니다. 예를들어서만일부동소수점 수가 n ulps 만큼의오차가있다면, 이말은 logβn 만큼의자리를신뢰할수없다는 뜻입니다. 만일, 상대오차가 nϵ 이라면, contaminated digits log βn (1.3) 1.3 보호숫자 (Guard Digits) 두부동소수점수의차이를계산하는한가지방법으로, 단순히두수의차이를계산한뒤에가장가까운부동소수점수로반올림하는방법이있습니다. 이방법은두수의차이가매우크다면, 비효율적인방법입니다. 예를들어 p = 3, 이라면, 계산결과 4

11 x = y = x y = 로반올림하면, 다시 가됩니다. 이모든자리를다계산하기보다는, 부동소수점하드웨어는보통고정된개수의자리수에서만연산을하게됩니다. 만일, p 개의자리수만보관한다고생각하고, 차이계산시더작은수를오른쪽으로쉬프트시킨다고생각하면 ( 쉬프트시뒤자리수는버려지게된다 ), 는 x = y = x y = 가됩니다. 이결과는차이를계산한뒤반올림한것과정확히같지요. 이번에는다른예제를살펴봅시다 의경우 x = y = x y = 정확한값은.17 이므로, 계산된값과그차이는무려 30ulps 에달합니다. Theorem 인자가 β 와 p 인부동소수점형식을사용할때, p 자리를사용해서두부동소수점수의차이를계산한다면, 그결과의상대오차는 β 1 까지커질수있다. Proof. x y 를계산할때상대오차가 β 1 인경우는 x = 이고, y =.ρρ...ρ 일때, ρ = β 1 이면된다. 여기서 y 는 p 자리수이고, 정확한차이는 x y = β p 이다. 하지만, 오직 p 자리만큼만사용한다고하였으므로, y 의가장오른쪽자리수는떨어져나가게된다. 따라서, 계산된차이는 β p+1 이고, 오차는 β p β p+1 = β p (β 1) 이고, 상대오차는 β p (β 1)/β p = β 1 이되므로성립한다. 나머지모든뺄셈에서의가능한상대오차가 β 1 보다작다는사실은자명하게보일수있습니다. β = 2 일대, 상대오차는결과값과동일하게될수있으며, β = 10 이면, 상대오차는결과값에 9 배나커질수있습니다. 다르게말하면, β = 2 일때, 식 1.3 에따르면, 신뢰할수없는자리수의경우 log 2 (1/ϵ) = log 2 (2 p ) = p 라는말입니다. 다시말해, 모든 p 자리수가틀렸다는이야기이지요. 만일추가적으로 1 개자리수가이와같은상황을방지하기위해더해졌다고가정해봅시다. ( 보호숫자 ) 이렇게되면뺄셈시에, 더작은수는 p + 1 개자리수로분할되며, 뺄셈의결과는 p 자리수로반올림됩니다. 보호숫자를도입하게될경우, 위예제는 5

12 x = y = x y = 가되고결과는정확합니다. 1 개의보호숫자를통해, 상대오차는 ϵ 보다커질수있으며, 에서나타나는것처럼, x = y = x y = 이결과는 102 로반올림되며, 정확한결과인 에비교해볼때상대오차는.006 이고이는 ϵ =.005 보다크게됩니다. 대부분의경우, 계산결과의상대오차는 ϵ 보다살짝커질수있으며정확히말하자면, Theorem 만일인자 β, p 를사용한부동소수점형식에서, x 와 y 사이의뺄셈을 p + 1 자리의수 (1 개의보호숫자 ) 로수행하였다면, 그결과의상대오차는 2ϵ 보다작다. 이에대한증명은반올림오차에서할것입니다. 참고로위증명에서 x 와 y 가양수인지음수인지상관이없기때문에뺄셈에대해증명한다면자동으로덧셈에대해서증명이됩니다. 1.4 지워짐 (Cancellation) 참고로숫자가 의미없 게 되었다라는말은, 영 어에서오염되었다 (contaminated) 라고합니다. 이차방정식의근의공식 b± b 2 4ac 2a 앞선내용에따르면보호숫자없이는매우가까운두수의뺄셈은엄청나게큰상대오차를발생시킨다는사실을알수있었습니다. 다시말해, 뺄셈이포함된식은매우큰상대오차를발생시키게되는데, 심지어는너무나상대오차가켜져서, 모든숫자들이의미없게만들어버릴수도있었습니다. ( 정리 1.3.1). 두개의가까운수를뺄때, 앞자리숫자 (most significant digits) 들은서로지워지게됩니다. 이때두가지방식의지워짐 (cancellation) 이가능한데, 이는각각나쁜 (catastrophic) 과착한 (benign) 지워짐이있습니다. 먼저나쁜지워짐 (catastrophic cancellation) 은, 두피연산자들에이미반올림오차가적용되어있을때발생하게됩니다. 예를들어 2차방정식의근의공식을계산할때, b 2 4ac 를계산하게됩니다. 이때, b 2 와 4ac 에는이미부동소수점곱셈을수행하였기때문에반올림오차가반영된상태입니다. 만일이들이가장가까운부동소수점으로반올림된다고할때 ( 따라서오차가.5ulp 이내 ), 이들사이에뺄셈을한다면, 두피연산자차이가적기때문에 ( 이렇다고앞에서가정함 ) 매우정확한자리수들이모두지워지고, 뒤에반올림으로인해신뢰할수없는자리수들만남겨놓게됩니다. 따라서, 그오차는기존의.5ulps 에서수 ulps 로뛰어오르게됩니다. 예를들어서 b = 3.34, a = 1.22, c = 2.28 이라고가정해봅시다. b 2 4ac 의정확한값은

13 입니다. 그런데, b 2 은 11.2 로반올림되고, 4ac 는 11.1 로반올림되서계산결과가.1 이나왔다면, 비록뺄셈에서 이정확히.1 이되었다고해도이는무려 700 ulps 의오차가됩니다. 착한지워짐 (benign cancellation) 은두개의정확한수의뺄셈시에발생하게됩니다. 만일 x, y 모두에반올림오차가없다면, 정리 에의해보호숫자를이용해서뺄셈수행시 x y 는매우작은상대오차 (2ϵ 이하 ) 를가지게됩니다. 나쁜지워짐의가능성을내포하고있는공식을때로는이를포함하고있지않도록변형할수있습니다. 예를들어서아래와같은이차방정식의근의공식을살펴봅시다 = 이고, 의마지막자리가 이므로 700 ulps 가된다. r 1 = b + b 2 4ac, r 2 = b b 2 4ac 2a 2a (1.4) 만일 b 2 >> ac 라서 b 2 4ac 에서지워짐이발생하지않는다면, 단순히 b 2 4ac b 라고볼수있습니다. 하지만이경우덧셈 ( 뺄셈 ) 에서나쁜지워짐이발생하게 되는데, 이를피하기위해서는식 1.4 의분자분모에 b b 2 4ac 를곱해서, r 1 = 2c b b 2 4ac, r 2 = 2c b + b 2 4ac (1.5) 를요도할수있습니다. 만일 b 2 >> ac 이고 b > 0 이라면, r 1 을식을이용해서계산하는것은나쁜지워짐을발생시키게됩니다. 하지만반대로식를이용해서 r 1 을계산하고식을이용해서 r 2 를계산하면이를막을수있습니다. 만일 b < 0 의경우반대로하면되겠지요. x 2 y 2 또한나쁜지워짐을야기시킬수있는식입니다. 각각을제곱해서곱하는것보다는 (x y)(x + y) 를계산하는것이더정확합니다. 이차방정식근의공식과는달리, 뺄셈에있음에도불구하고그뺄셈이착한뺄셈이기때문에큰반올림오차없이처리할수있습니다. 정리 에의해 x y 의상대오차는최대 2ϵ 이기때문이죠. x + y 도마찬가지입니다. 두개의작은상대오차를가지는부동소수점수를곱한그결과역시작은상대오차 ( 섹션반올림오차참고 ) 를내므로효율적이라볼수있습니다. 계산된결과와실제정확한결과를구분하기위해서다음과같은기호를사용하도록하겠습니다. 예를들어서 x y 를실제로정확한차이라고한다면, x y 는컴퓨터로계산된차이가됩니다. ( 즉, 반올림오차를포함한 ) 비슷하게도,,, 등이각각덧셈, 곱셈, 나눗셈을의미합니다. 또한, 대문자로쓰여진함수들, 예를들어 LN(x), SQRT (x) 는컴퓨터로계산된값을의미하고, 원래쓰이던표현 (ln x, x) 등은정확히계산된값을의미합니다. 비록 (x y) (x y) 는 x 2 y 2 의훌륭한근사라고해도, 만일 x, y 자체가실제값 ˆx, ŷ 의근사값일수도있습니다. 예를들어 ˆx, ŷ 가 10진수로정확히표현되지만이진수로는정확히표현되지않는값일수도있습니다. 이경우비록 x y 는 x y 의 7

14 좋은근사일수있지만 ˆx ŷ 에비해큰상대오차를가질수있기때문에 (x+y)(x y) 의정확도는 x 2 y 2 에비해크게효과적이지않을수있습니다. 사실 (x + y)(x y) 나 x 2 y 2 의경우계산하는데걸리는시간이거의비슷하기때문에, 그나마좀더정확한전자가선호되겠지만, 일반적인경우나쁜지워짐을착한지워짐으로바꾼다고해서근본적인문제 ( 입력값자체가근사값이라는사실 ) 을해결할수없습니다. 하지만, 지워짐자체를아예발생하게하지않는다면, 입력되는데이터가근사값이라고해도매우의미있을것입니다. x 2 y 2 은 (x + y)(x y) 로씀으로써나쁜지워짐을착한지워짐으로바꿀수있었기때문에정확도를향상시킬수있었습니다. 다음예제로, 나쁜지워짐을가지고있지만모두착한지워짐만으로바꿀수있는흥미로운예를하나살펴볼것입니다. 삼각형의넓이는아래식 1.6 에따라각변의길이 a, b, c 를통해표현할수있습니다. A = s(s a)(s b)(s c), where s = (a + b + c)/2 (1.6) 예를들어삼각형이매우납작하다고가정합시다 (a b + c). 즉, s a 가됩니다. 따라서, 식 1.6 에서의 (s a) 항은근접한두값을빼게되것이고, 심지어그들중하나는반올림오차를포함하고있을수도있습니다. a = 9.0, b = c = 4.53 인경우 s 의정확한값은 9.03 이고 A 는 입니다. s 의계산된값은 9.05 로오직 2 ulps 의오차가있지만 A 의계산된결과는 3.04 로, 무려 70 ulps 의오차가발생하게됩니다. 아래는식 1.6 을재구성해서납작한삼각형에경우에도정확히계산할수있도록수정한것입니다. [Kahan 1986] A = (a + (b + c))(c (a b))(c + (a b))(a + (b c)), a b c (1.7) 4 만일 a, b, c 가 a b c 를만족하지않더라도, 식 1.7 자체는순환적이기때문에그냥알파벳순서만살짝바꿔주면됩니다. 위로수정된결과를이용하면계산된값은 2.35 가나오며이는오직 1 ulps 의오차입니다. 위예제만을보더라도알겠지만일반적인경우에도식 1.6 보다는식 1.7 가훨씬나은모습을보여줍니다. Theorem 식 1.7 를사용하였을때발생하는반올림오차는최대 11ϵ 이다. ( 단, 뺄셈은보호숫자를이용하였고, e.005, 그리고제곱근의계산오차는 0.5ulp 이내이다.) ϵ <.005 라는조건은대부분의부동소수점시스템에서만족하고있습니다. 예를들어서 β = 2, p 8 나 β = 10, p 3 만되도, ϵ <.005 를만족합니다. 정리 에서의식에서상대오차를고려할때, 위식자체가부동소수점연산을통해처리된다는것을이해하셔야합니다. 따라서우리가실제로계산하는식은 8

15 SQRT ((a (b c)) (c (a b)) (c (a b)) (a (b c))) 4 (1.8) 보통오차상한선은매우크게잡히는것이일반적인데, 예를들어정리 에서는넉넉하게 11ϵ 정도로잡았는데, 위예제에서계산되었던 2.35는실제값 에대해상대오차가 0.7ϵ 정도밖에안됩니다. 이렇게오차상한선을크게잡는이유는나중에수치적으로증명할때조금더편리하기때문입니다. 마지막예로 x << 1 일때의 (1 + x) n 가착한지움을살펴보도록수정하는것을살펴봅시다. 참고로이식은금융계산에서자주쓰이게되는데, 예컨대매일 100 달러를은행계좌에집어넣은다고가정하고, 연이율이 6% 라고생각합시다. 만일 n = 365 이고, i =.06 이라면 1년이지난후에계좌에누적된돈의양은 100 (1+i/n) n 1 1/n 달러가됩니다. β = 2, p = 24 일때컴퓨터로계산해본다면그결과 $ 가나오는데, 실제로정확한결과는 $ 로 $1.40 의오차가발생하게됩니다. 이오차가발생하는이유는매우간단한데, 1 + i/n 은 1 에 을더하는것이므로, 너무나값이작기때문에 i/n 의뒷자리비트들이잘려나가게됩니다. 특히이반올림오차는 1 + i/n 이 n 승이될때더확대됩니다. 이문제많은 (1+i/n) n 은 e n ln(1+i/n) 로다시쓸수있는데, 이번에문제는 ln(1+x) 를작은 x 에대해연산하는것이문제가됩니다. 이근사식을이용하는한가지방법으로 ln(1 + x) x 임을생각할수있는데, 이경우, 계산된식이 $ 이되어오히려기존의오차였던 $1.40 보다증가한 $3.21 이됩니다. 하지만다행이도 ln(1 + x) 를정확하게계산할수있는방법이있는데이는정리 에나타나있습니다. [Hewlett-Packard 1982] 이공식을이용하면 $ 로오직 2 센트오차만발생하게되죠. Theorem 만일 ln(1 + x) 를아래공식을이용해계산하면 x if 1 x = 1 ln(1 + x) = x ln(1+x) (1+x) 1 if 1 x 1 일때상대오차는최대 5ϵ 이다. ( 단, 0 x < 3/4 이고, 뺄셈시보호숫자를사용하고, e < 0.1, ln 을 0.5ulp 이내로계산한다.) 사실이공식은위를만족하는어떤 x 값에대해서도적용가능하지만, 실제로우리가관심을가지는것은그대로 ln(1 + x) 를계산시에나쁜지워짐이발생하는 x << 1 일경우입니다. 위식자체가매우이상하게보일지라도, 왜이식이작동하는지는매우 9

16 간단하게알수있습니다. ln(1 + x) 를 x( ln(1 + x) ) = xµ(x) x 여기서문제가되는부분은 µ(x) 를정확하게계산할수있느냐인데, 왜냐하면 µ(x) = ln(1+x) x 를계산할때 1 + x 를연산시큰반올림오차가생기기때문입니다. 하지만여기서중요한사실은 µ 의값은 ln(1 + x) x 이기에거의일정하다는사실 입니다. 따라서 x 를살짝바꾸는것은큰오차를야기하지않습니다. 다시말해, 만일 x x 라면, xµ( x) 를계산하는것은 xµ(x) = ln(1 + x) 의좋은근사가될것이라는 것이지요. 그런데, x 와 x + 1 모두정확하게계산할수있는 x 를찾을수있나요? 바로 x = (1 x) 1 을사용하는것입니다. 왜냐하면 x + 1 은정확히 1 x 가될것이기 때문이죠. 이지워짐 (Cancellation) 에서다루어왔던내용들을정리해보자면, 차이가거의 없지만정확하게알려진두수의뺄셈시에 ( 착한지움 ) 보호숫자를도입한다면일정 정확도가보장된다는것입니다. 또한, 때로는공식자체가부정확한결과를낼수있지 만, 착한지움을이용한다면같은공식이라도다르게연산함으로써높은정확도를낼 수있게됩니다. 하지만이역시뺄셈을수행시보호숫자를도입한연산이보장될때 이야기입니다. 사실보호숫자를도입하는데필요한비용은크지않습니다. 왜냐하면 이는단순히 1 개비트정도를추가하는것뿐이기때문이죠. 예를들어서 54 비트 배정밀도누산기 (adder) 의경우 1 개비트를추가하는데전체비용에 2% 정도밖에 되지않습니다. 이정도비용으로는여러분은공식 1.6 과같은삼각형넓이를구하는 알고리즘이라던지, ln(1 + x) 을계산하는알고리즘을사용할수있게되는것입니다. 대부분의현대컴퓨터들이보호숫자를사용하는가운데, 일부시스템 (Cray System 등 ) 에서는사용하지않고있습니다. 1.5 정확한반올림연산 보호숫자를이용한부동소수점연산은정확하게계산한다음에가장가까운부동소수점수로반올림하는것보다는정확하지않습니다. 이러한방식으로계산하는것을정확히반올림되었다 (exactly rounded) 라고부릅니다 4). 정리 이전에바로나온예제가보여주듯이, 보호숫자를이용한연산은언제나정확한반올림된결과를보여주지는않습니다. 이전장에서보호숫자를이용하여여러알고리즘을처리하는것을살펴보았다면이번장에서는정확한반올림을요구로하는알고리즘들을살펴볼것입니다. 아직까지반올림에대한정확한정의를내리지않았습니다. 사실반올림자체는정확히반을가루는예외적인경우를제외하고는매우직관적입니다. 더가까운쪽으로올리거나버리면되는것인데, 문제는그중간에수가위치할경우입니다 는 12 4) 혹은 correctly rounded 라고불리기도합니다 10

17 로내림이되어야할까요아니면 13 으로올림이되어야할까요. 어떤학교에서는아마 0,1,2,3,4 는내림을해야하고 5,6,7,8,9 는올림을해야한다고가르쳐야할것입니다. 이와같은방식은 Digital Equipment Corporation 의 VAX 컴퓨터들이채택하고있는방식이기도합니다. 반올림을수행하는다른방식으로는마지막자리가 5 인경우, 올림을하던내림을하던두가지가능성이있기때문에, 50% 의경우에올림을하고나머지 50% 의경우의내림을수행하도록하면된다는것입니다. 즉이방식에서는항상반올림된수가짝수가되도록반올림을하게결정하면올림과내림모두공평하게처리할수있게됩니다.. 즉, 12.5 는 13 이되는것이아니라 12 가되고, 13.5 는 13 이아니라 14 가됩니다. 둘중어느방식이효과적인지는 Reiser and Knuth [1975] 에서다루고있으며, 그결과항상짝수로반올림하는것이효과적이라결정났습니다. Theorem x, y 를부동소수점수라하고, x 0 = x, x 1 = (x 0 y) y,..., x n = (x n 1 y) y 라정의합시다. 만일 와 가모두정확한짝수반올림연산을수행된다면, 모든 n 에대해 x n = x 이거나, n 1 에대해 x n = x 1 을만족한다. ( 단, n 1) 위정리 의내용을명확히하자면, 예를들어 β = 10, p = 3 이고 x = 1.00, y =.555 라고가정해봅시다. 반올림시에, 위순열은 x 0 y = 1.56, x 1 = = 1.01, x 1 y = = 1.57,... 가되고, 각순열마다 x n 의크기는 x n = 9.45(n 845) 일때까지.01 씩늘어나게됩니다. 반면에짝수반올림시에, x n 은항상 1.00 이됩니다. 이예제를통해알수있는점은, 그냥반올림을사용한다면계산결과는조금씩증가하게되는데짝수반올림을사용하는경우그러한현상이나타나지않는다는점입니다. 따라서이제부터이문서전체에서항상짝수반올림만을사용할것이며, 그냥반올림이라칭할시에, 짝수반올림을의미하는것입니다. 정확한반올림을사용하는한가지경우는바로다정밀도 (multiple precision) 연산입니다. 다정밀도연산은보통두가지방식으로구현하는데, 크기가큰가수를가지고있는부동소수점수를사용한다던지 ( 예를들어 WORD 의배열로보관한다는등 ), 다정밀도부동소수점수를보통의부동소수점수의배열의합으로써나타내면됩니다. 여기서다룰방식은후자로, 보통의부동소수점들을배열로구현하는것은고급언어에서쉽게코딩할수있기때문입니다. 다만, 이를위해서는정확한반올림연산을필요로합니다. 이러한시스템에서, 곱셈연산은다음과같이구현됩니다. 예를들어서다정밀도부동소수점수 x, y 에대해서각각이 x = x h + x l 이고 y = y h + y l 이라한다면, 둘의곱은 여기서 WORD 는컴퓨터자료형의한종류로보통 32비트 (4바이트) 정수를의미합니다. xy = x h y h + x h y l + x l y h + x l y l 로주어집니다. 만일 x, y 의가수가 p 비트라면, 우변의각각의항들역시가수의비트가 p 비트일것이고, x l, x h, y l, y h 각각은 p/2 비트로나타낼수있게됩니다. 만일 p 가짝수라면, 간단히 x 0.x 1...x p/2 1 과 x p/2...x p 1 의합으로나타낼수 11

18 있겠지요. 하지만, p 가홀수라면전자와같이간단히두개로쪼갤수없게되는데, 하지만음수를사용함으로써비트하나를늘려서짝수비트로만들수있습니다. 예를들어서만일 β = 2, p = 5, x = 일때 x 는 x h =.11 와 x l = 로쪼갤수있습니다. 수를쪼개는방식은이방식말고도여러방법이있는데, Dekker [1971] 덕분에매우간단히계산할수있게되었습니다. 다만, 이방법으로는 1 개이상의보호숫자가필요합니다. Theorem p 를부동소수점정밀도 ( 단 β > 2 일때, p 는항상짝수 ), 그리고부동소수점연산시에정확히반올림된다고가정합시다. 만일 k = [p/2] 라정의하고, m = β k + 1 이라정의하면, x 는 x = x h + x l 로분할할수있으며, 이때 x h = (m x) (m x x), x l = x x h 이고각각의 x i 는 [p/2] 비트의정밀도로표현가능하다. 이기법은후의증명에서많이사용됩니다. 증명시 크기를조정하여 라는말은모두 β 의멱수를곱해서가수는유지시키고, 지수만바꾼다는의미입니다. 많은경우지수는크게관련이없기때문에이기법을사용하여원하는범위내로주어진수를위치시킬수있습니다 이정리 가어떻게작동하는지확인하기위해 β = 10, p = 4, b = 3.476, a = 3.463, c = 라해봅시다. b 2 ac 를가장가까운부동소수점수로반올림하면 이되는데, 이때 b b = 12.08, a c = 이므로따라서계산된 b 2 ac 는.03 입니다. 이는무려오차가 480ulps 나되지요. 정리 을이용하면, b = , a = , c = 로나타낼수잇습니다. 이변형된형태를이용하여계산을수행해보면, b 2 = 가되고 ( 이시점에서는아직그값이계산된상태는아님 ), ac = ( ) = 가됩니다. 이제, b 2 ac 를각각의항에대한뺄셈을수행한다면그결과 = 이되며이는정확히반올림된결과와같게됩니다. 정리 이정확한반올림을필요로한다는것을보여주기위해서한가지예로 p = 3, β = 2, x = 7 일때를생각해봅시다. m = 5, mx = 35, m x = 32 가되는데, 만일한개의보호숫자로뺄셈을수행하게된다면 (m x) x = 28 이됩니다. 따라서 x h = 4, x l = 3 이되어 x l 은 [p/2] = 1 비트로표현할수없게되죠. 정확한반올림의마지막예를위해, m 을 10 으로나누는경우를생각해봅시다. 그결과는실제로 m/10 한결과와부동소수점수와일반적으로다를것입니다. 하지만놀랍게도 β = 2 인경우에, m/10 한결과에 10 을곱한다면다시 m 으로돌아오게됩니다. 사실더일반적인경우에도이는성립하는데 (Kahan 이이를멋지게보였습니다 ) 증명은기발하지만궁금하지않은독자들은 IEEE 표준으로넘어가시기바랍니다. Theorem β = 2 일때, m, n 이모두정수이고, m < 2 p 1 이고 n = 2 i + 2 j 의특별한꼴일때, 연산시정확한반올림을한다면 (m n) n = m 을만족한다. Proof. 2의멱수를곱하는것은아무런문제가되지않는데왜냐하면이는가수를바꾸는것이아니라단순히지수만바꾸는것이기때문입니다. 따라서만일 q = m/n 이라면, 2 의멱수를곱해서 2 p 1 n < 2 p 가되고, m 역시마찬가지로해서 1/2 < q < 1 가 12

19 되도록조정할수있게됩니다. 따라서 2 p 2 < m < 2 p 가되죠. 이때, m 이 p 개의비트를가지게되므로, 소수점오른쪽으로최대한개의비트만가질수있게됩니다. m 의부호를바꾸는것은증명에영향을주지않으므로 q > 0 이라가정합시다. 만일 ˆq = m n 이라정의한다면, 위정리를증명하기위해서는 nˆq m 1 4 (1.9) 임을보이면충분합니다. 왜냐하면 m 은소수점뒤로최대 1 개의비트만있으므로 nˆq 는자동으로 m 으로반올림될것이기때문입니다. 예외적으로 nˆq m = 1/4 가되는경우를살펴보자면, 앞서조건에서조정하기전 m 은 m < 2 p 1 이였기때문에, 조정된이후의 m 의마지막비트역시 0 이됩니다. 따라서절반으로나뉜경우역시 m 으로반올림됩니다. ( 항상짝수로반올림됨 ) 이제 q =.q 1 q 2..., 이라정의하고, ˆq =.q 1 q 2...q p 1 이라정의합시다. 이제 nˆq m 을계산하기위해서는먼저 ˆq q = N/2 p+1 m/n 을계산하도록합시다.. 이때, N 은홀수인정수이다. n = 2 i +2 j 이고 2 p 1 n < 2 p 이므로어떤 k p 2 에대해반드시 n = 2 p k 가되어야만합니다. 따라서 따라서, ˆq q = nn 2p+1 m n2 p+1 = (2p 1 k + 1)N 2 p+1 k m n2 p+1 k 여기서분자는당연히정수이고, N 이홀수이기때문에, 분자는홀수가됩니다. ˆq q 1 n2 p+1 k 가성립합니다. 이때 q < ˆq 라가정하면 ( 그반대도역시비슷하다 ), nˆq < m 이 성립하고, 따라서 m nˆq = m nˆq = n(q ˆq) n(q (ˆq 2 p 1 )) 따라서정리가증명됩니다. n(2 p 1 1 n2 p+1 k ) = (2 p k )2 p 1 2 p 1+k = 1 4 위정리는 2 가아닌임의의밑수 β 에대해서도 2 i + 2 j 가아닌 β i + β j 에대해서 성립하게됩니다. 다만, β 가커지게되면, 그꼴의좀더드물어질뿐입니다. 우리는이제다음질문에답할수있게되었습니다. 과연, 기초연산에서필요한 것보다살짝조금더반올림오차가발생하는것에대해신경을써야할까? 그에대한 답은, 신경을써야한다이고, 왜냐하면기초연산을정확하게하는일이야말로여러공 13

20 식들이작은상대오차를야기하도록사용할수있기때문입니다. 지워짐 (Cancellation) 에서우리는정확한결과를내기위해보호숫자를필요로하는여러가지알고리즘에대해논의하였습니다. 만일그공식들에게주어지는입력값자체에불확실성이있다면정리 과 의내용들은아마쓸모가없어질것입니다. 착한지워짐이었던 x y 가, x, y 자체가근사값이되다면나쁜지워짐이될수있기때문이죠. 하지만, 정확한연산자체는부정확한입력값에대해서도매우유용한데, 왜냐하면정리 과 에서나타났듯이정확한관계를구축할수있기때문입니다. 이들은비록부동소수점변수가참값의근사치라도매우유용하게사용할수있습니다. 14

21 제 2 장 IEEE 표준 IEEE 에는두개의부동소수점연산표준이있습니다. IEEE 754 는이진표준으로단정밀도 (single precision) 에는 β = 2, p = 24 를, 배정밀도 (double precision) 에는 p = 53 을요구합니다 [IEEE 1987]. 또한이표준에는단정밀도와배정밀도에서어떠한식으로비트를배열할지도명시하고있습니다. IEEE 854 는 754 와는달리 β = 2 와 β = 10 을허용하고, 어떠한식으로부동소수점수가비트들로인코딩되야하는지명시하지않았습니다 [Cody et al. 1984]. 또한명확히 p 로어떠한값을써야하는지지정하지는않았지만, 단, 배정밀도에서가능한 p 값들을제한해놓았습니다. 앞으로 IEEE 표준이란말은두개의표준을모두지칭하는데사용하도록하겠습니다. 이장에서는 IEEE 표준에대해간단하게짚고넘어가도록하겠습니다. 각각의섹션에서는 IEEE 표준의한측면에대해살펴보도록하며, 왜이러한것을정했는지도이야기하겠습니다. 이문서는 IEEE 표준이가장좋은부동소수점표준이라고주장하기위해만들어진것이아니라, 단지어떻게하면 IEEE 표준을잘이해하고사용할수있는지설명해주는것입니다. 더자세하게알고싶은분들은 IEEE 표준문서를직접참고하시기바랍니다. [IEEE 1987; Cody et al. 1984] 2.1 형식과연산들 밑수 (base) 왜 IEEE 854 에서 β = 10 을허용하는지는명백합니다. 우리인간이십진수를사용해서계산을하기때문입니다. 따라서 β = 10 은특히계산기에서사용되기에매우적합한데, 왜냐하면계산기에서각각의계산결과는모두십진수로표시되기때문입니다. 하지만 IEEE 854 에서밑이 10 이아니라면반드시 2 를사용하도록정해놓은이유는있습니다. 상대오차와 Ulp 에서한가지이유를이야기했었었는데. 왜냐하면.5ulps 의반올림오차에해당하는상대오차가 β 의배수로증가되기때문입니다. 뿐만아니라, 밑수가커지게된다면상대적으로유효정밀도가떨어진다는점이있는데, 예를들어 β = 16, p = 1 인부동소수점과 β = 2, p = 4 인경우를각각살펴보도록 15

22 합시다. 이둘은모두동일하게 4 비트의가수를가지게됩니다. 이를이용해서 15/8 을 계산한다고생각해봅시다. β = 2 인경우에 15 는 으로표현되고, 15/8 은 이되서정확하게계산됩니다. 하지만 β = 16 일때 15 는 F 16 0 이고 (F 는 16 진수로 15) 15/8 은 로오직 1 개비트만정확하게나옵니다. 일반적으로 16 진수는최대 3 개의비트를손실하게되므로, p 16 진수자리수는유효정밀도가 이진수의경우 4p 였겠지만 4p 3 비트까지떨어지게됩니다. 따라서 β 가커지게되면 이러한문제들이발생하게되는데요, 그럼왜 IBM 은그들의 system/370 시스템에 β = 16 을사용하였을까요. 사실정확한이유는 IBM 만이알고있겠지만, 아마 2 개의이유를들수있겠습니다. 첫번째이유는더넓은범위의지수를이용할수있는 점인데요, IBM system/370 의단정밀도는 β = 16, p = 6 을사용하였습니다. 따라서 가수부분은 24 비트이고, 7 비트를지수부분이차지하게됬는데, 이를통해서 에서 = 2 28 까지의넓은범위의수를사용할수있게되었습니다. β = 2 일때 동일한지수범위를사용하려면 9 비트를지수로사용해야하고, 가수로오직 22 비트 밖에사용할수없게죕니다. 하지만, 앞에서도말했지만 β = 16 일때유효정밀도가 최대 21 비트까지떨어지므로오히러 β = 2 를사용했을안정적으로 23 비트의정밀도 (22 비트의가수지만뒤에서설명하겠지만추가적인비트로정밀도를한비트더얻을 수있음 ) 를얻을수있게되어서더안정적입니다. β = 16 을고른이유로또가능한것은바로쉬프트연산때문입니다. 두개의 부동소수점수를더할때, 만일그들의지수들이다르다면, 둘중하나의가수들이 쉬프트되어서소수점을맞춰야하므로연산속도가느려지게됩니다. β = 16, p = 1 시스템에서는 1 부터 15 까지의숫자들이모두같은지수를가지게되므로, ( 15 2 ) = 105 가지의쌍에대해덧셈시에쉬프트연산을할필요가없게됩니다. 반면에 β = 2, p = 4 시스템에서는 1 부터 15 까지의숫자들의지수들이 0 부터 3 까지에걸쳐있으므로 105 개의쌍중에서 70 개의쌍사이에서쉬프트연산이필요하게됩니다. 하지만대부분의현대하드웨어에서는쉬프트연산으로이한작업저하정도가 매우미미하기때문에워블을작게줄여주는 β = 2 를이용하는경우가많습니다. 특히 β = 2 를이용하게된다면 1 개의추가적인비트를가수에더사용할수있기때문에 1) 정밀도를더높일수있습니다. 모든부동소수점수가언제나정규화되기때문에, 가수의가장최고비트 ( 맨왼쪽 ) 는언제나 1 이되므로굳이이비트를저장하기위해 비트를낭비할필요가없게됩니다. 이러한트릭을이용하는것을, 숨겨진비트 (hidden bit) 를사용한다고말합니다. 하지만이트릭을이용하기위해서는부동소수점형식 (Floating-point Formats) 에서지적하였듯이, 0 을위한특별한처리가필요하게됩니다. 따라서 0 을표현하기위해지수가 e min 1 이고, 가수가모두 0 이라면, e min 1 이아니라, 0 을나타내도록정하였습니다. IEEE 754 단정밀도는 32 비트로구성되며, 1 비트는부호비트, 8 비트는지수비트, 그리고나머지 23 비트는가수비트를위해사용됩니다. 하지만숨겨진비트를이용하여, 1) 이사실은비록 Knuth ([1981], 211 쪽 ) 가 Konard Zuse 의공으로돌렸지만사실 Goldberd [1967] 이최초로발표하였습니다. 16

23 형식 인자 단정밀도 확장단정밀도 배정밀도 확장배정밀도 p e max > e min 지수비트수 형식전체비트수 표 2.1: 표 IEEE 형식인자들 비록 23 비트로인코드되지만, 실질적인가수부분은 24 비트 (p = 24) 가됩니다 정밀도 (Precision) IEEE 표준에서는 4 개의서로다른정밀도를지정합니다. 단정밀도, 배정밀도, 확장단정밀도 (single-extended), 확장배정밀도 (double-extended). IEEE 754 에서는대부분의부동소수점하드웨어에서제공하는단, 배정밀도와대응됩니다. 단정밀도는 32 비트워드를사용하고, 배정밀도는두개의인접한 32 비트워드를사용하게됩니다. 확장정밀도형식은약간더늘어난정밀도와지수범위를지원합니다. ( 표 참고 ) 확장정밀도의경우, IEEE 표준은얼마나많은비트들을더추가해야하는지에대한하한선만을지정하였습니다. 확장배정밀도의경우최소한 80 비트이상을사용하도록정하였는데, 표 에나온것은 79 비트만으로구현해도충분합니다. 그럼 1 비트를더추가적으로필요로하는이유는, 대부분의확장정밀도를지원하는하드웨어들이숨겨진비트를사용하지않기때문에 79 비트만으로충분하지만최소 80 비트이상을요구하고있습니다. 표준안은확장정밀도에좀더많은강조를하였는데 ( 반면에배정밀도지원은특별히강조하지않았지만 ), 표준안구현시에확장형식을꼭지원하도록강력히추천한다 (Implementations should support the extended format corresponding to the widest basic format supported..) 라며사용을권장하였습니다. 이렇게확장정밀도사용을도입한한가지중요한계기는바로계산기때문입니다. 계산기들은보통 10 개의숫자들만표시하지만, 내부적으로 13 개의숫자를사용해서처리하기때문입니다. 계산기의경우, exp 나 log 와같은함수들을 10 자리정확도로계산하기위해서는, 몇개의자리수가더필요하게됩니다. 따라서내부적으로 13 개의자리수를이용해서계산하되, 보여줄때는정확한 10 개비트만보여주게됩니다. IEEE 표준의확장정밀도역시이와비슷한역할을합니다. 이를이용해서라이브러리로하여금.5ulp 오차이내로단정밀도 ( 혹은배정밀도 ) 계산을수행할수있게됩니다. 하지만확장정밀도사용시중요한것이있는데, 바로이를사용시에반드시사용자에게투명하게공개해야된다는것입니다. 예를들어계산기의경우내부에서처리된값이화면에보여지는값과동일한정밀도로반올림되지않는다면, 그값을이용해서계산을계속하였을때, 내부에숨겨진숫자들때문에유저들의입장에서는불일치하게나타날수있기때문입니다. 17

24 확장정밀도를더설명하기위해, IEEE 754 단정밀도와 10진수사이의변환문제에대해생각해보도록합시다. 이상적으로, 단정밀도수들은표시될때충분한십진수숫자로표시되므로, 다시십진수로읽어들일때원래의단정밀도로복구할수있게됩니다. 사실 9 자리의십진수만이용해도유일한단정밀도이진수로복구할수있다는것이알려져있습니다 ( 이진수를십진수로변환하기참조 ) 어떤십진수를다시유일한이진표현으로변경할때, 1ulp 정도로작은반올림오차조차치명적일수있습니다. 왜냐하면이는아예다른답이되기때문이지요. 따라서이때효율적인알고리즘을위해서확장정밀도를사용하는것입니다. 만일확장단정밀도가사용가능하다면십진수를이진수로바꾸는매우직관적인방법이있습니다. 먼저 9 개의십진숫자들을 N 이라고합시다. ( 소수점을무시하고읽는다 ) 그럼표 에서 p 32 이고, 10 9 < 이므로 N 은정확히확장단정밀도로표현될수있습니다. 이제, 원래숫자의소수점위치를맞추기위해적절한 10 p 을찾습니다. 이제 10 p 을계산하는데, 만일 P 13 이라면, 이역시정확하게변환될수있습니다. 왜냐하면 = 이고 5 13 < 2 32 이기때문이죠. 이제마지막으로, N 과 10 P 를곱하면됩니다. 만일마지막작업이정확히수행되었다면, 가장가까운이진수로복구할수있게됩니다. 이진수를십진수로변환하기을보면어떻게마지막곱셈을수행하는지설명되어있습니다. 따라서 P 13 일때, 확장단정밀도형식을사용하는것은 9 자리십진수를정확히가장가까운이진수로변환하는것이가능합니다 ( 단, 정확히반올림된경우에 ). 만일 P > 13 이라면확장단정밀도로는항상정확히반올림된이진수로변환하는것이위알고리즘으로는충분하지않습니다. 하지만 Coonen [1984] 은이진수를십진수로변환하고다시원래이진수로복구하는데에는충분하다는사실을보였습니다. 만일배정밀도가지원된다면, 위알고리즘은확장단정밀도보다는배정밀도에서수행하여도됩니다. 하지만배정밀도를 17 자리십진수로변환하고다시역변환하는작업은확장배정밀도를필요로합니다 지수 bias 의사전적의미는편향, 치우침으로, 결과적으로지수값에 127 을뺌으로써그값이치우쳐져있다라는맥락에서이름붙인것이라생각하시면됩니다. 지수자체가양수혹은음수둘다될수있기때문에부호를나타내기위해특별한방법을필요로합니다. 보통부호있는수를나타내는방법으로두가지방법이있는데하나는부호와그절대값을따로보관 (sign/magnitude) 하는방법과 2 의보수표현 (2 s complement) 방식이있습니다. 부호비트를따로보관하는방법은 IEEE 형식에서가수부분을처리할대사용하는것입니다. 즉 1 개비트를부호를나타내는데사용하고, 나머지비트를그수의크기 ( 절대값 ) 을나타내는데사용하는것입니다. 한편 2 의보수표현은보통정수연산에서자주쓰이게됩니다. 이표현방식을사용하면 [ 2 p 1, 2 p 1 1] 범위의수들을모듈로 2 p 로나오는가장작은음이아는정수로표현하게됩니다. IEEE 이진표준은지수를처리하기위해이두방식모두사용하지않고, 그대신에바이어스 (biased) 표현을사용합니다. 예를들어단정밀도에서, 지수부분은 8 비트로 18

25 저장되게되는데, 이때바이어스는 127 입니다 ( 배정밀도는 1023). 이말은, 만일 k 가 부호없는정수로표현된지수비트라면, 실제로부동소수점수의지수값은 k 127 이라는뜻입니다. 이실제지수값은바이어스된지수인 k 와구분하기위해바이어스 되지않은지수 (unbiased exponent) 라고부릅니다. 표 에딷르면단정밀도는 e max = 127 이고 e min = 126 이었습니다. 왜 e min < e max 가되도록처리하였다면, 가장작은수의역수가 ( 1 2 e min ) 오버플로우 (overflow) 나지않게하기위함이었습니다. 비록가장큰수의역수가언더플로우 (underflow) 날수있지만, 많은경우언더플로우가오버플로우보다는덜심각한문제 이기때문입니다. 밑수 (base) 에서 e min 1 이 0 을표현하기위해사용된다고하였고, 특별한수들섹션에서 e max + 1 의사용예들을볼것이빈다. IEEE 단정밀도에서는, 바이어스된지수들의범위가 e min 1 = 127 에서 e max + 1 = 128 이라는것이고, 다시말해바이어스되지않는지수들의경우그범위가 0 에서 255 까지라는의미 입니다. 이는정확히부호없는 8 비트이진수가표현할수있는수들이지요 연산들 (operations) IEEE 표준에서는덧셈, 뺄셈, 곱셈, 그리고나눗셈의결과가정확히반올림되기를 요구하고있습니다. 즉, 연산이먼저정확히계산이된다음, 가장가까운부동소수점 수로반올림 ( 짝수반올림방식을사용해서 ) 되어야한다는말입니다. 이전의보호숫자 (Guard Digits) 에서, 두개의부동소수점들의지수차이가매우클때두수의합이나 차이를계산하는일은매우어려운일이라는것을말했었습니다. 따라서보호숫자를 도입해서작은상대오차를유지하면서효과적으로계산하였다고하였습니다. 하지만, 한개의보호숫자를이용해서계산하는것은정확히계산한다음에반올림하는것과 언제나같은결과를낼것이라보장할수없습니다. 하지만두번째보호숫자와, 세번째 스티키비트 (sticky bit) 를도입해총 3 개의보호숫자를이용해한개의보호숫자를 이용해서계산하는것보다살짝더많은비용으로정확히계산한후반올림하는효과를 낼수있습니다. [Goldberg 1990] 따라서이를통해표준안이효율적으로구현될수 있지요. 산술연산결과를명확히지정하는이유는바로프로그램들의이식성 (portability) 을위해서입니다. 같은 IEEE 산술표준을사용하는두개의기계사이에서작동하 는프로그램은반드시같은결과를내도록표준안이정해져있으며, 만일그결과가 다르다면산술연산에서의차이가아니라소프트웨어의버그라고할수있게됩니다. 세세하게표준안을정해놓은또다른이유로, IEEE 산술에서산술적으로올바른값을 낸다고증명된것들을 IEEE 표준을지원하는기계에서정확히작동할것이라보장할 수있기때문입니다. Brown [1981] 은모든부동소수점하드웨어에사용가능한공리 (axiom) 들을제안 하였습니다. 하지만이공리들은지워짐 (Cancellation) 이나정확한반올림연산에서 다루어왔던유용한알고리즘들을증명할수없었는데, 왜냐하면이들은하드웨어에 공통적으로들어있지않는특별한기능들 ( 보호숫자연산등 ) 을필요로하기때문이 오버플로우는 (overflow) 는단어그대로, 수의크기가너무커서보관할수있는범위를초과해버린다는의미입니다. 예를들어서 4비트정수형을보관하는자료형이있다면, 최대크기는 1111 인데, 만일여기에 1 을더하게된다면, 오버플로우가발생하며, 0000 이되버립니다. 이는엄청난계산상의문제가되겠지요. 언더플로우 (underflow) 는이정반대의경우입니다. 19

26 었습니다. 게다가 Brown 의공리들은연산이정확히수행되고반올림되는지단순히정의하는것보다훨씬복잡하였습니다. 따라서 Brown 의공리들로부터정리들을증명하는일은, 단순히작업들이정확히반올림된다고가정하고증명하는일보다훨씬어렵습니다. 부동소수점표준에서어느연산까지다루어야할지에대해서는아직의견이일치가되지않았습니다. 게다가, IEEE 표준에서는기본적인사칙연산외에도, 제곱근, 나머지연산 (modular), 그리고정확히반올림되는정수와부동소수점수와의변환등에대해명시하고있습니다. 뿐만아니라표준에서는내부형식과십진수사이의변환이정확히반올림되도록요구하였습니다. ( 아주큰수들은제외 ) Kulish 와 Mirnaker [1986] 은표준안에내적연산을포함시켜서그연산과정을정확히명시해달라고제안하였습니다. 그들은 IEEE 산술연산으로이용해서내적연사을계산했을때, 최종답안이꽤틀릴수있기때문이라고주장하였습니다. IEEE 표준안은표제작자의딜레마 (table maker s dilemma) 때문에초월함수들의값들이정확히반올림되도록요구하지않습니다. 예를들어여러분이지수함수의표를 4 자리정밀도로만든다고합시다. 그렇다면 exp(1.626) = 일때이를반올림을 으로해야할까요, 아니면 로해야할까요. 사실 exp(1.626) 을좀더적확히계산한다면, 이됩니다. 더정확히하면 더정확히하면 사실 exp 함수가초월함수이기대문에 exp(1.626) 값을정확히반올림하기위해서는 ddd 가될지 ddd 가될지모르기때문에무한히많은연산이필요하질도모릅니다. 따라서, IEEE는초월함수들의정밀도를정확하게구하는것은무의미하다고판단하였습니다. 이문제를해결하기위해제안된방법으로표준안에서이러한초월함수들의값을구하는알고리즘들을명확히지정하는방법이있는데, 사실모든하드웨어아키텍쳐위에서동일하게잘작동하는단일알고리즘은없기대문에결과적으로는표준안에서초월함수에대한어떠한필수적인요건은갖추지않게되었습니다. 유리근사 (Rational approximation), CORDIC 2), 그리고거대한표를이용하는것이현대에컴퓨터들에서초월함수값을계산하는기본적인방법이되었습니다. 각각의방법들은하드웨어에의존적이며, 모든하드웨어에대해잘작동하는알고리즘은아직까지없습니다. 2.2 특별한값들 몇몇의부동소수점하드웨어에서는모든비트패턴들이하나의유효한부동소수점수를나타냅니다. 예를들면 IBM 의 System/370 을들을수있죠. 하지만 VAM tm 의경우, 일부비트패턴들을특별한수를나타내기위해서예약해놓았는데, 이들을예약된피연산자 (reverved operand) 라고부릅니다. 이아이디어는 CDC 6600 으로거슬러올라가는데, 이시스템에서미정와무한대를나타내기위해특별한값들 - 2) CORDIC 은 Coordinate Rotation Digial Comptuer 의약자로, 쉬프트와덧셈연산을많이사용하는 ( 그리고극소수의덧셈과나눗셈 ) 초월함수계산방법입니다 [Walther 1971]. 20

27 지수 가수 나타내는값 e = e min 1 f = 0 ±0 e = e min 1 f 0 0.f 2 e min e min e e max - 1.f 2 e e = e max + 1 f = 0 ± e = e max + 1 f 0 NaN 표 2.2: 표 2.2 IEEE 754 의특별한값들 INDEFINITE, INFINITY 를사용하였습니다. IEEE 표준에서는이전통을계승하여서 NaN (Not a Number - 수가아님 ) 과무한대를도입하였습니다. 이러한특별한값들을도입하지않는한, 음수의제곱근을계산하는것과같은예외상황을연산을중지시키지않는한잘처리할수없게됩니다. IBM System/370 FORTRAN 의경우, 음수의제곱근과같은명령을처리할때기본값으로수행하는명령은오류메세지를표현하는것입니다. 게다가이시스템에서는모든비트패턴이하나의유효한수이기때문에, 예를들어 4 제곱근역시어떤부동소수점값을리턴하게됩니다. 이 System/370 FORTRAN 의경우 4 = 2 를리턴하였습니다. IEEE 산술에서는이러한상황에서 NaN 이리턴됩니다. IEEE 표준안에서는다음과같은몇개의특별한값들을정해놓았습니다. ( 표 2.2 참조 ) ±0, 비정규화된수, ±, 그리고 NaN 들입니다 ( 다음장에서설명하겠지만 NaN 은한개가아닙니다 ) 이러한특별한값들의지수부분은 e max + 1 이거나 e min 1 을가지게됩니다 NaN 들전통적으로 0/0 이나 1 과같은계산들은복구할수없는오류로, 계산을중단시키게만들었습니다. 하지만, 이러한값들이나타나도계산을계속해야만하는몇가지상황들이있는데요, 예를들어서함수의해 (zero) 를찾는루틴을만들었다고합시다. 어떤함수 f 의해를찾는 zero(f) 라는함수가있다면, 전통적으로이함수는 f 가정의되어있는구간 [a, b] 안에서해를찾도록구간에대한정보를따로요구를할것입니다. 따라서, 해찾는함수를정확히말하자면 zero(f, a, b) 라고할수있겠습니다. 하지만유용한해찾는함수라면, 구간 [a, b] 에대한정보와같은추가적인정보를필요로하지는않겠지요. 우리가해를찾는함수를만들었을때그함수가정의된구간까지따로입력해준다면매우귀찮은일이아닐수없습니다. 하지만, 이렇게 zero 힘수가해를찾기위해탐색을할때, 만일 f 의정의역을벗어나는값을넣게된다면, 아마도함수값을계산하는도중에 0/0 이나 1 과같은것들을계산할것이며, 이때전체계산과정을중단시키고불필요하게도해찾는과정역시중단되게되겠지요. 이렇나문제는 NaN 이라불리는특별한값을도입함으로서해결할수있습니다. 0/0 이나 1 과같은것을연산하였을때계산을중단시키기보다는 NaN 을리턴하도록말이지요. NaN 이발생할수있는사례들은에잘나타나있습니다. 따라서 zero(f) 가 f 의정의역밖을탐색하더라도, 이함수는 NaN 을리턴하여서, 함수가계산을 21

28 연산 NaN 이만들어지는경우 + + ( ) 0 / 0/0, / REM x REM 0, REM y x, x < 0 표 2.3: 표 NaN 을발생시키는연산들 중지시키지않고계속찾아나갈수있게된다는것이지요. 이러한사실을염두에둘때, 우리는 NaN 과보통의정상적인부동소수점수와의연산을했을때어떠한결과가나올지생각할수있습니다. 예를들어 f 가 return(-b + sqrt(d))/(2*a) 와같이정해졌다고해봅시다. 만일 d < 0 이라면 f 는 NaN 을리턴하게됩니다. 왜냐하면 d < 0 이면 sqrt(d) 가 NaN 이되므로 b + sqrt(d) 가 NaN 되도록해야하기때문입니다. 비슷하게도, 나눗셈시만일하나의피연산자가 NaN 이라면그나눈결과역시 NaN 이됩니다. 일반적으로 NaN 이어떤부동소수점연산을하게되면그결과는 NaN 이됩니다. 사용자가정의역범위를지정하지않아도안전하게해를찾는함수를설계하는또하나의방법은바로시그널 (signal) 을사용하는것입니다. 해찾는함수안에부동소수점예외를처리하기위한시그널핸들러 (signal handler) 을도입한다음에, 만일 f 를계산하는데정의역을벗어난값을계산해서부동소수점예외가발생하였다면, 시그널처리기에서이를받아처리한후다시계산을지속할수있습니다. 다만이방식을도입하였을때문제점은, 모든언어에서각기다른방법으로시그널을처리하는방법이있기때문에 ( 혹은아예없는경우도있다 ), 높은이식성을보여주지는못합니다. IEEE 754 에서 NaN 은주로지수가 e max + 1 이고 0 이아닌가수를가지는값으로구현됩니다. 정확히가수에어떠한값이들어가야하는문제는시스템에따라다르게설계되므로, 한개의유일한 NaN 이존재하는것은아니고, NaN 계열들이존재하게됩니다. 다만 NaN 과어떤부동소수점수와연산하였을때그결과역시 NaN 이되어야만합니다. 따라서만일어떠한긴연산의결과가 NaN 이였을때, 그시스템독립적인정보를포함하고있는가수를분석하여서언제가장첫번째 NaN 이발생하였는지를알수있게됩니다. 사실이마지막문장에는약간의어폐가있는데, 만일두피연산자모두 NaN 이라면그연산결과는둘중하나의 NaN 이므로가장먼저만들어진 NaN 을추적할수는없습니다 무한대 NaN 이 0/0 이나 1 과같은것들을처리하는데사용되었다면, 무한대는오버플로우가발생하였을때처리하는데사용됩니다. 무한대를따로도입하는것이단순히표현가능한최대의수를사용하는것보다더안전한방법입니다. 예를들어서 x 2 + y 2 를계산한다고해봅시다. 만일 β = 10, p = 3, e max = 98 이라고하고 x = 이고, y = 이라고한다면 x 2 는오버플로우가발생해서

29 이될것입니다. y 역시마찬가지로 x 2 + y 2 계산시각각은모두오버플로우가나서 전체결과역시오버플로우가발생한 이됩니다. 따라서최종연산결과는 = 이되고이는참값인 과전혀다르게나옵니다. 반면에 IEEE 산술을하게된다면 x 2 은 이고, 마찬가지로 y 2 역시 라서, 최종 결과는 가됩니다. 차라리이렇게처리하는것이, 참값과전혀다른부동소수점수를 리턴하는것보다더안전하다고볼수있습니다. 0 을 0 으로나누는것은 NaN 을발생시킵니다. 하지만 0 이아닌수를 0 으로 나누는것은무한대를발생시킵니다. 즉, 1/0 =, 1/0 = 가됩니다. 이두개의 상황을구분짓는이유는다음과같습니다. 만일 x 0 일때, f(x) 0, g(x) 0 이라고합시다. 그러면 f(x)/g(x) 는어떠한값도가질수있게됩니다. 예를들어서 f(x) = sin(x), g(x) = x 라고하면, x 0 일때 f(x)/g(x) 1 이됩니다. 반면에 f(x) = 1 cos(x) 라면, f(x)/g(x) 0 이되죠. 따라서 0/0 을생각해볼때, 이는그 어떤값도될수있습니다. 따라서 IEEE 표준에서는이를아예 NaN 이라고처리하였 습니다. 반면에 c > 0 이고, 임의의해석함수 f, g 에대해 f(x) c, g(x) 0 이라면 f(x)/g(x) ± 가됩니다. 따라서 IEEE 표준에서는 c 0 인이상 c/0 = ± 가 되도록정의하였습니다. 이때 의부호는 c 와 0 의부호에따라가도록하였습니다. 또한사용자들은플래그 (flag) 의상태를점검함으로써 의발생이오버플로우로인한 발생인지, 0 으로나눔으로인한발생인지확인할수있습니다. 전자의경우에서는 오버플로우플래그 (overflow flag) 가설정이되고, 후자의경우에서는제로플래그 (zero flag) 가설정이됩니다. 무한대를가지고있는연산에서그결과를결정하는것은꽤같단합니다. 무한대에 해당하는부분을 x 로바꾸고, x 를계산하면되지요. 따라서 3/ = 0 입니다. 왜냐하면 lim 3/x = 0 x 비슷하게도, 4 = 이고, = 가되빈다. 하지만극한값이존재하지 않는경우에그결과는 NaN 이되는데, 따라서 / 는 NaN 이됩니다. ( 표 에 여러예들이있습니다 ) 이러한논리는 0/0 이왜 NaN 이되었는지와일맥상통합니다. 어떠한것을연산을할때그중간에 NaN 이포함되어있다면전체연산결과또한 NaN 이됩니다. 하지만 ± 의경우, 1/ = 0 과같은규칙들때문에 가포함되어 있다고해도그연산결과가부동소수점수가될수있습니다. 예를들어 x/(x 2 + 1) 과같은함수를생각해봅시다. 사실이는안좋은식인데, 분모의 x 가 ββ e max/2 보 다커지게된다면오버플로우가발생해서 1/x 에가까운값이아니라 0 이라는값을 출력하게됩니다. 하지만위식을다시쓰면 1/(x + x 1 ) 로쓸수있는데, 이를통해 너무일찍오버플로우가발생하는일을막을수있게됩니다. 또한 x = 0 이여도 x = 0 : 1/( ) = 1/(0 + ) = 1/ = 0 이라는올바른결과를내게됩니다. 만일 무한대연산이없었다면이와같이식을변형하였을대따로 x = 0 인지를확인을해야 겠었지요. 무한대연산이있는덕택에여러특별한경우에따로체크를해줄필요가 없어졌습니다. 하지만그대신에 x/(x 2 +1) 처럼무한대에서이상한결과를내지않는지 23

30 확인해주어야만합니다 부호있는영 (signed zero) 0 은지수에 e min 1 을넣고가수부분이 0 인것으로표현됩니다. 부호비트에 2 개의서로다른값이들어갈수있기때문에 +0, 0 두가지의 0 이가능합니다. 만일 +0 과 0 가서로다른것으로취급된다면, 매우간단한 if(x == 0) 과같은작업들 조차 x 의부호에따라매우불확실한결과를내게됩니다. 따라서 IEEE 표준에서 는 0 < +0 이아니라 +0 = 0 이라고정해놓았습니다. 물론 0 의부호자체를 무시하는방법도있지만 IEEE 표준에서는그렇게하지않았습니다. 왜냐하면, 부호 있는 0 을포함한곱셈이나나눗셈시행시에, 보통의수에서부호를처리하듯이계산 결과가나타나기때문이지요. 따라서 3 ( + 0) = +0, +0/ 3 = 0 이됩니다. 만일 0 에부호가없었다면 1/(1/x) = x 와같은관계는 x = ± 일때성립하지않게 됩니다. 왜냐하면 0 에부호를두지않는다면 1/, 1/ + 모두 0 이되는데, 1/0 은 + 이기때문이지요. 따라서 0 에부호가없다면 1/(1/x) = x 와같은식을계산하는 과정에서오버플로우된정보의부호를읽어버리게되는결과를야기할수있습니다. 부호있는 0 을사용하는다른이유로언더플로우와 0 에서로그함수와같은불연속적인 함수를다룰때사용됩니다. IEEE 산술에서 log 0 = 이고 x < 0 일때 NaN 이라고 정의하는것은자연스럽습니다. 만일 x 가 0 으로언더플로우된매우작은음수라고 생각해봅시다. 그렇다면 0 에부호가있는덕택에 x 는음수가될것이므로, 로그함수는 NaN 을리턴하게되겠지요. 만일 0 에부호가없다면, 로그함수는이 0 이매우작은 음수가 0 으로언더플로우된건지양수가언더플로우된것인지구별할수없기대문에 그냥 를리턴하게됩니다. 다른예제로 0 에서불연속적인 signum 함수를생각할 수있는데요, 이함수는수의부호를리턴하는함수입니다. 부호있는 0 을사용하는가장흥미로운예제로복소수연산을들수있습니다. 간단한예를보자면, 1/z = 1/( z) 를생각해봅시다. 이식은 z 0 일때자명하게 성립합니다. 하지만 z = 1 일때에는, 좌변의경우 1/( 1) = 1 i 가되고, 우변 의경우 1/( 1) = 1/i = i 가되어서성립하지않게되죠. 따라서 1/z 1/( z) 가되버리는것입니다. 이러한문제가발생하게된근본적인이유는, 제곱근함수자체가 다중값함수이기 (multi-valued function) 때문에전체복소평면에서연속이되도록 값을선택할수없기때문입니다. 하지만, 모든음의실수를포함하는분기선 (branch) 을제외하고생각한다면 ( 이러한것을분기서법 (branch cut) 이라고합니다 ), 제곱근 함수를연속으로만들수있게됩니다. 그럼이제 x > 0 일때 x + i0 꼴의남은음의 실수들을어떻게처리해야할지문제가생기게됩니다. 여기서부호있는 0 은이문제를 해결하기위한완벽한방법이됩니다. x + i(+0) 꼴의수들은한개의부호 (i x) 를가지고분기선의반대쪽에있는 x + i( 0) 꼴의수는부호가 i x 가됩니다. 사실, 를계산하는자연스러운공식은이러한결과를내게되죠. 따라서 1/z = 1/( z) 로다시돌아가보면만일 z = 1 = 1 + i0 이라면, 24

31 1/z = 1/( 1+i0) = [( 1 i0)]/[( 1+i0)( 1 i0)] = ( 1 i0)/(( 1) ) = 1+i( 0) 따라서 1/z = 1 + i( 0) = i 이고, 1/( z) = 1/i = i 가되서등식이성립합니다. 따라서 IEEE 산술은모든 z 에대한항등식을성립시켜줍니다. 좀더복잡한예들은 Kahan [1987] 을참고하시기바랍니다. 비록 +0 과 0 을구분하는일은여러가지장점이있지만, 종종혼동을유발하기도합니다. 예를들어서부호있는영은 x = y 1/x = 1/y 관계식을만족시키지않기때문입니다. (x = +0, y = 0) 하지만 IEEE 위원회는부호있는 0 을사용하는데얻는이점이그단점을훨씬넘어섰다고판단하였습니다 정규화되지않는수예를들어서 β = 10, p = 3, e min = 98 일때 x = , y = 들은가능한최소의부동소수점수인 보다 10 배나큰지극히평범한부동소수점수들입니다. 그런데, 이들은이상한특징이있는데, x y 임에도불구하고 x y = 0 이된다는것이지요. 이러한문제가발생하는것은 x y = = 로정규화된수로표현하기에너무나작아서결국 0 이되버린다는것입니다. 하지만아래관계를유지하는일은매우중요합니다. x = y x y = 0 (2.1) 만일위관계가유지되지않는다면 if(x y then z = 1/(x-y) 와같은매우자명해보이는코드조차 0 으로나눗셈이발생하여오류가발생할수있씁니다. 이러한코드에서발생하는버그를찾아내는것은매우어려운일이지요. 만일식 2.1 가명확하다면매우신뢰할수있는부동소수점코드를작성할수있겠지만, 대부분의수에게만참이고일부수에서는거짓이라면, 아예사용할수없겠지요. IEEE 표준에서는식 2.1 을포함한많은유용한관계식을보장하기위해비정규화된수 (Denormalized numbers) 라는것을도입하였습니다. 사실이부분은표준에서도가장논쟁이많이되었던부분이였고, 실제로 754 가채택되기전까지많은지연이있었던부분입니다. IEEE 표준에호환이된다고주장하는많은고성능하드웨어들은비정규화수들을직접적으로지원하지않고, 비정규화가발생하였을대, IEEE 표준을지원하도록소프트웨어적인측면에서처리하도록하여고있습니다. 3) 비정규수라는아이디어를최초로제공한것은 Goldberd [1967] 로매우간단한아이디어입니다. 바로, 지수가 e min 일때, 가수는정규화될필요가없다는것입니다. 따라서 β = 10, p = 3, e min = 98 일때, 은더이상가장작은부동소수점수가아니게됩니다. 왜냐하면 3) 이러한사실은표준에서가장문제가되고있는부분이기도합니다. 언더플로우를자주발생시키는프로그램은예외를하드웨어보다소프트웨어상에서처리할때눈에띄는속도저하를느낄수있습니다. 25

32 그림 2.1: 점진적언더플로우를사용할때와사용하지않을때 역시부동소수점수가되기때문입니다. β = 2 이고숨겨진비트가사용되었을대, 약간의문제가있는데왜냐하면지수가 e min 여도언제나가수가 1.0 보다크게된다는것입니다 ( 암묵적으로맨앞의비트가 1 이니까 ). 그해결책은 0 을표현하는데사용했던것과동일한것으로표 2.2 에나와있습니다. 지수 e min 은비정규화수를나타내는데사용됩니다. 정확히말하면, 만일가수비트들이 b 1, b 2,..., b p 1 이고, 지수값이 e 라할때, 만일 e > e min 1 이면수는 1.b 1 b 2...b p 1 2 e 가되고, 만일 e = e min 1 이라면, 그수는 0.b 1 b 2...b p 1 2 e+1 로처리됩니다. 이때 +1 이필요한이유는비정규화수들의지수는 e min 이지 e min 1 이아니기때문입니다. 앞서 x = , y = 의예를다시생각해봅시다. 비정규화수를사용했다면 x y 는 0 이되지않고그대신에, 비정규화수인 이됩니다. 이러한동작은점진적언더플로우 (gradual underflow) 라고부릅니다. 그리고점진적언더플로우를이용하면식 2.1 을증명하기매우쉽습니다. 그림 2.1 는비정규화수들을보여주고있습니다. 먼저그림에서윗줄에는정규화된부동소수점수들을보여주고있는데, 0 과가장장근정규화된수간의차이인 1.0 β e min 을볼수있습니다. 따라서만일계산된결과가이사이안에떨어지게된다면바로 0 이되버리게된다는것입니다. 하지만아래줄을보면, 비정규화수를사용하였을때어떠현효과를볼수있는지알수있습니다. 그 간극 이메꿔져있기때문에, 계산결과크기가 1.0 β e min 이될지라도, 가장가까운비정규화수가되지 0 으로떨어지지않게됩니다. 수직선에비정규화수들이들어가게된다면, 두인접한수사이의간격은원래정규화수와동일하게작용합니다. 즉, 두수사이의간격은원래정규화수와같거나, β 의배수로달라지게되죠. 비정규화수가없다면, 그간격이 β p+1 β e min 에서 β e min 으로급격히차이가나게됩니다. 따라서많은알고리즘에서언더플로우문턱에서많은상대오차를보여주었던정규화된수들의경우, 점진적언더플로우를사용하게되면잘작동하게됩니다. 점진적언더플로우를사용하지않을경우정규화된입력들에대해간단한식인 x y 와같은것조차앞에서 x = , y = 의예에서보았듯이매우큰상대오차를보여주게됩니다. 아래예제 [Demmel 1984] 에서보여주듯이, 지워짐이없더라도큰상대오차는발생할수있습니다. 아래와같은두개의복소수 a + ib, c + id 26

33 의나눗셈을살펴봅시다. 자명하게도, 그식은 a + ib c + id ac + bd ad = c 2 + ibc + d2 c 2 + d 2 로주어집니다. 앞에서도이야기하였지만만일분모 c + id 의절대값이 ββ emax/2 보다커진다면실제계산값이범위안에있더라도오버플로우가나게됩니다. 따라서위 식을좀더제대로계산하는 Smith 의공식은 a + ib c + id = a+b(d/c) c+d(d/c) + i b a(d/c) c+d(d/c) b+a(c/d) d+c(c/d) i a+b(c/d) d+c(c/d) if ( d < c ) if ( d c ) (2.2) 이 Smith 의공식을적용하게되면 ( i10 98 )/( i( )) 이되어서점진적언더플로우인 0.5 를내게됩니다. 만일점진적언더플로우를사용하지않는다면 0.4 의계산결과는 0 으로내림이되어서오차가무려 100 ulps 에달할것입니다. 사실, 통상적인경우비정규화수들은오차의한계를 1.0 β e min 정도까지보장하게됩니다. 2.3 예외, 플래그, 예외처리기 (Exceptions, Flags, Trap handlers) 0 으로나누기나오버플로우와같은예외적인상황들이 IEEE 산술연산에서발생할경우, 디폴트동작은그냥결과값을전달하고연산을지속하는것입니다. 보통그결과값들은아마 0/0 이나 1 의경우 NaN 이겠고, 1/0 이나오버플로우의경우 가되겠지요. 이전섹션들에서이러한디폴트값들이어떻게사용될수있을지이야기하였습니다. 예외상황들이발생하였을때이러한값을리턴하는것뿐만이아니라, 상태플래그 (status flag) 역시설정이됩니다. IEEE 표준에서는유저들로하여금이렇나상태플래그의값을읽고쓸수있도록요구합니다. 이플래그들은한번설정이되면, 명시적으로클리어되기전까지값이유지됩니다. 플래그를검사하는것만이예컨대 가발생하였을때 1/0 때문인지, 아니면정말오버플로우로인한무한대인지구별할수있는유일한방법입니다. 때로는예외상황이발생하였을때작업을지속하는것이적절하지않을수도있습니다. 무한대에서 x/(x 2 + 1) 과같은예를들수있었는데요, 만일 x > ββ emax/2 가된다면분모가무한대가되서최종계산값이 0 이됩니다. 이전에도이야기했지만위식을 1/(x + x 1 ) 로다시쓰게된다면위문제를해결할수있겠지만, 언제나이렇게식을재구성하는것이항상답은아닙니다. IEEE 표준은예외처리기 (trap handler) 의도입을가능케하는것을강력하게추천하여는데요, 따라서예외가발생하였을때플래그를설정하는것이아니라예외핸들러가호출되도록하였습니다. 예외핸들러에서리턴된값을이용해서작업을적당히처리하게말이지요. 상태플래그를클리어하던지설정하던지결정하는것은예외핸들러의몫입니다. 27

34 예외 예외핸들러가사용되지않을때 예외핸들러에전달되는인자 오버플로우 ± 나 ±x max round(x2 α ) 언더플로우 0, ±2 e min, 비정규화수 round(x2 alpha ) 0 으로나눠짐 ± 피연산자 올바르지못한작업 NaN 피연산자 부정확함 round(x) round(x) 표 2.4: IEEE 754 의예외들. 참고로 x 는연산의정확한결과이고, 단정밀도의경우 α = 192, 배정밀도의경우 α = 1546 입니다. 그리고 x max = e max 입니다. 인터럽트 (interrupt) 란, 하드웨어상에서운영체제로보내는어떠한신호를뜻합니다. 주로이신호들은운영체제에서처리하는작업을중지하고이신호를먼저처리하기때문에인터럽트라불립니다. IEEE 표준은예외를 5 개의종류로분류하였습니다. 오버플로우, 언더플로우, 0 으로나눠짐, 올바르지못한작업 (invalid operation), 그리고부정확함 (inexact) 이죠. 각각의예외상황에따라독립적인상태플래그를가지게됩니다. 일단앞서말한 5개의상태플래그중 3 개는매우명백합니다. 올바르지못한작업은표 에서나타난상황들과임의의 NaN 들간의비교를의미합니다. 올바르지못한작업예외가발생하였을때디폴트로 NaN 이리턴되지만, 그역은참이아닙니다. 즉, 어떠한연산에서하나의피연산자가 NaN 이라면, 그결과는 NaN 이겠지만그연산이표 의조건을만족하지않는이상올바르지못한작업예외는발생하지않게됩니다. 부정확함예외 (inexact exception) 는부동소수점연산의결과가부정확할때발생합니다. 예를들어서 β = 10, p = 3 인시스템을생각해봅시다. 그럼, = 14.7 은정확하지만, = 15.0 은부정확하므로부정확함예외를발생시키게됩니다. ( = 15.05) 이진수를십진수로변환하기에서부정확함예외를사용하는알고리즘에대해설명할것입니다. 5 개예외의동작에대해서는표 2.3 에정리해놓았습니다. 사실부정확함예외가너무많이발생한다는구현상의문제가있습니다. 만일부동소수점하드웨어에플래그가있는대신에운영체제에인터럽트를보내는방식으로작동한다면, 너무나많은인터럽트가운영체제로전해져서문제가될수있기때문입니다. 이러한문제는소프트웨어적으로상태플래그를관리하는방식으로해결할수있습니다. 최초로예외가발생하였을때, 이에해당하는소프트웨어플래그를설정한뒤에, 부동소수점하드웨어로하여금그클래스의예외들을무시하도록설정하면되기때문입니다. 따라서그뒤에뒤따라나오는예외들은모두운영체제에인터럽틀르전달하지않고진행되며, 사용자가상태플래그를리셋하였을때다시예외를무시하지않도록설정하면됩니다 예외처리기 (Trap handler) 예외처리기를사용하는한가지명백한이유는기존의코드와호환성을위해서입니다. 옛날코들들의경우, 예외가발생하였을때작업이중지지되어서프로세스에예외처리기를설치할수있을것이라가정하고만들어져있습니다. 이는특히 do S until (x >= 100) 과같은코드를처리할때유용합니다. 왜냐하면 NaN 과수를 <,, >,, = ( 단 는아니다 ) 비교하는일은언제나 false 이므로위코드는 x 28

35 가 NaN 이된다면무한루프를돌게때문이죠. 또한예외처리기는오버플로우가발생할수있는 Π i=1 x i 와같은것을처리할때유용하게사용할수있습니다. 이에대한한가지해결책으로로그를사용하는것인데, exp log x i 를대신계산하는것입니다. 하지만이러한방식으로처리하였을때정밀도가더떨어지게되고, Π i=1 x i 를계산하는것보다비용이높아진다는문제가있습니다. 이문제를해결하는다른방법으로오버 / 언더플로우의개수를세는예외처리기를설치하여서해결할수있습니다 [Sterbenz 1974]. 그방법은다음과같습니다. 0 으로초기화된전역카운터가있고, 부분곱 p k = Π k i=1 x i 라정의한다음에, 어떤 k 에서오버플로우가발생한다면예외처리기에서카운터를 1 증가시키고, 지수를다시가능한범위안에들어오도록조정합니다. IEEE 754 단정밀도의경우 e max = 127 이므로, 만일 p k = 이된다면, 오버플로우가발생하게되서예외처리기가호출됩니다. 이때, 예외처리기는지수를단정밀도범위에포함되게조정하여서 p k 를 로조정합니다. ( 아래참고 ) 비슷하게도, p k 가언더플로우가나게되면카운터가하나줄어들며음의지수를단정밀도범위에포함되도록조정하게됩니다. 만일모든곱셈이수행된이후에, 카운터가 0 이라면그최종결과는그대로 p n 이될것입니다. 그런데, 카운터가양수라면그결과값은오버플로우된것이고음수라면언더플로우된것을알수있게됩니다. 만일부분곱에서한번이라도범위를초과하지않았더라면오류처리기는한번도호출되지않기때문에추가적인비용이들지않습니다. 만일그중간에오버 / 언더플로우가발생하게되었다고해도, 그계산결과는로그를사용하였을때보다더정확하게된다, 왜냐하면 p k 는 p k 1 로부터정확한곱셈을통해계산되었기때문입니다. IEEE 754 에서는예외처리기가호출될때, 인자로조정된결과 (wrapped around) 를전달하게됩니다. 여기서오버플로우의조정된결과라는의미는, 그결과값이마치무한대의정밀도를가지고있는것처럼계산된다음에 2 α 로나뉘어진후, 올바른정밀도로반올림된다는의미입니다. 언더플로우의경우그반대로 2 α 가곱해진다음에, 동일하게처리된다는것이지요. 여기서 α 는단정밀도의경우 192, 배정밀도의경우 1536 이됩니다. 이때문에앞선예에서 이 가된이유지요 반올림방식들 (Rounding Modes) IEEE 표준에서반올림은매계산시에항상발생하게됩니다. 기본값으로반올림의뜻은가장가까운값으로바꾼다는의미입니다. 표준에따르면 3 개의다른반올림방식을지원하도록하였는데요, 각각은 0 으로의반올림, + 를향한반올림, 그리고 를향한반올림입니다. 예를들어서, 정수로의변환연산을할경우, 를향한반올림은마치바닥함수 (floor function) 처럼작동하고, + 를향한반올림은마치천장함수 (ceiling function) 처럼작동하게됩니다. 반올림방식은오버플로우에도영향을주게되는데, 예를들어서 0 으로의반올림이나, 를향항반올림이작동중이라면, 양의오버플로우가발생하였을때, + 가아닌가장큰표현가능한수를리턴하도록바뀌게된다는것입니다. 마찬가지로음의오버플로우가발생하였을때, 만일 + 나 를향한반올림은쉽게말해내림 (floor) 이고, + 를향한반올림은올림과 (ceil) 동일합니다. 그리고 0 을향한반올림은말그대로정수부분만취하는것으로 (truncate), 어떤 x 의 sgn(x) y 와동일합니다. 29

36 구간산술연산이란, 수학적으로계산을할때, 반올림오차와수치적오차를고려하여어떠한구간안에포함시켜서신뢰할수있는계산결과를얻고자하는연산입니다. 예를들어서, 사람의키를어림을할때 2.0 미터로정확히세는것이아니라, 그키를 1.97 에서 2.03 미터사이라고구간을정한후이를가지고연산하는분야입니다. 0 으로의반올림이작동중이라면가장큰음의수를리턴하게됩니다. 이반올림방식의한가지응용예로구간산술연산 (interval arithmetic) 을들수있습니다. 예를들어서두수 x, y 의합은구간 [z, z] 안에있는데, 여기서 z 는 x y 가 로반올림된것이고, 반대로 z 는 x y 가 로반올림된것입니다. 만일반올림모드들을사용하지않는다면, 구간산술은보통 z = (x y)(1 ϵ) 과 z = (x y)(1 + ϵ) ( 단, ϵ 은머신입실론 ) 으로계산되며, 이는구간의크기를꽤크게잡는경향이있습니다. 구간연산의결과는항상구간으로표현되기때문에, 구간연산의입력값역시구간으로주어지는데요, 만일두개의구간 [x, x], [y, ȳ] 을더한다면, 그결과는 [z, z] 가되며, 이때 z 는 를향한반올림방식에서 x y 를구한것이고, z 는 + 를향한반올림방식에서 x ȳ 를구한것입니다. 부동소수점연산이구간산술연산으로수행된다면, 최종구간안에는반드시연산의정확한결과가포함됩니다. 구간의크기가매우커진다면 ( 사실대부분의경우가그렇습니다 ) 그구간연산의결과는별로유용하지못할것입니다. 왜냐하면참값이구넓은구간어디엔가있는것이기때문이죠 플래그 (Flag) IEEE 표준에는여러개의플래그들과모드들이있습니다. 앞에서도이야기하였지 만, 5 개의예외 ( 언더플로우, 오버플로우, 0 으로나눔, 올바르지못한연산, 부정확함 ) 각각에한개의상태플래그가할당되어있습니다. 또한모드들의경우 4 개로앞에서도 이야기하였듯이 0,, 를향한반올림과, 가장가까운수로의반올림 ( 짝수반올림 ) 말이지요. 이섹션에서는어덯게하면이모드들과플래그들을잘조합해서사용할수 있을지에대해다룰것입니다. 좀더복잡한예로 이진수를십진수로변환하기을 참고해보시기바랍니다. 예를들어서 x n 을계산하는루틴을짠다고생각해봅시다. 여기서 n 은양의정수 인데, 여러분은아래와같이간단히짤수있을것입니다. PositivePower ( x, n ) { while ( n i s even ) { x = x*x n = n/2 } u = x while ( true ) { n = n/2 i f ( n==0) return u x = x*x i f ( n i s odd ) u = u*x } 30

37 만일 n < 0 이라면, x n 을계산하기에좋은방법이 PositivePower(1/x, -n) 을 호출하기보다는 1/PositivePower(x, -n) 을계산하는것입니다. 왜냐하면전자의 경우이미반올림오차가생기는 1/x 를 n 번이나곱해야되기때문이죠. 반면에후자의 경우맨마지막연산에서한번의반올림오차만이생기는것이기때문에좀더정확 하다고할수있습니다. 하지만불행이도, 이방법에는약간의문제가있습니다. 만일 PositivePower(x, -n) 에서언더플로우가발생한다면, 언더플로우예외처리기가 작동하던지, 아니면언더플로우상태플래그가설정이될것입니다. 그런데사실이는 틀렸습니다. 왜냐하면 x n 이언더플로우가된다면, x n 이범위안에있을수있기 때문입니다. 4) 하지만 IEEE 표준에서는사용자로하여금모든플래그에접근가능하게 하므로서브루틴에서이문제를쉽게고칠수있습니다. 먼저오버플오우와언더플로우 예외가능비트들을꺼버리고, 오버플로우와언더플로우상태비트들을저장합니다. 그리고 1/PositivePower(x, -n) 을계산한다음에, 만일두상태비트들모두설정 되어있지않다면이들의예외가능비트들을다시켜버리면됩니다. 만일상태비트들 중하나가설정되어있다면플래그를다시복구한뒤에, PositivePower(1/x, -n) 을 이용하여계산을해서올바른예외가발생하도록하면됩니다. 플래그를사용하는다른예제로, arccos x = 2 arctan 1 x 1+x 를계산할때가있습니 다. 만일 arctan( ) 가정확히 π/2 로계산된다면, arccos 1 은언제나무한대연산 덕택에 2 arctan( ) = π 가되겠지요. 하지만여기에약간의문제가있는데, arccos( 1) 이전혀예외적인상황이아님에도불구하고 1 x 1+x 는 0 으로나눠짐플래그를설정해버리기때문입니다. 이문제를해결하는방법은직관적인데, 이연산이전에 0 으로 나눠짐플래그의값을저장해놓고연산후에원래플래그의값으로복구하면됩니다. 4) 왜냐하면 x < 1, n < 0 일때만일 x n 이언더플로우문턱인 2 e min 보다겨우몇비트차이가난다면, x n = 2 e min < 2 e max 가될수있기때문에오버플로우가나지않을수있습니다. 참고로모든 IEEE 표준에서 e min < e max 입니다 31

38 제 3 장 시스템적측면 컴퓨터시스템을설계함에있어서부동소수점에관련한지식은필수적이라할수있습니다. 대부분의컴퓨터아키텍쳐에는부동소수점연산이포함되어있고, 컴파일러들은반드시이러한부동소수점명령을생성할수있어야하며, 운영체제는부동소수점예외가발생하였을때무엇을해야할지결정할수있어야만합니다. 하지만이러한정보들은주로컴퓨터디자이너들이아닌소프트웨어제작자들에게전달되지, 컴퓨터시스템디자이너들이부동소수점의수치적측면에대해정보를얻을기회가많지않습니다. 아래예제는디자인측면에서타당해보이는 BASIC 코드가어떻게예측하지못한동작을할수있는지살펴보도록합시다. q = 3. 0 / 7. 0 i f q = 3. 0 / 7. 0 then p r i n t Equal : e l s e p r i n t Not Equal 놀랍게도 IBM PC 에서 Borland 사의 Turbo Basic 을이용해서컴파일하면놀랍게도, Not Equal 을출력하게됩니다. 이예제는다음명령어집합 (instruction sets) 에서살펴보도록합시다. 그런데, 어떤사람들은위와같은코드에서발생하는문제가바로부동소수점을어떠한오차범위 E 를가지고비교하는것이아니라직접적으로비교함으로써발생하는오류라고생각합니다. 하지만이러한처방은몇가지궁금증을불러일으키는데, 과연그렇다면 E 의값으로얼마를써야할까요? 만일 x < 0 이고 y > 0 인데 E 범위안에있다면두수의부호자체가다르지만같다고해야할까요. 게다가 ab a b < E 는동치가아닌데, 왜냐하면과 abandbc 라고해서 ac 를보장하지는않기때문입니다. 하드웨어에서지원하는명령어들을모아놓은것을명령어집합이라고부릅니다. 3.1 명령어집합 (instruction sets) 어떤알고리즘에서계산을할때계산중간에특정계산들은높은정밀도로계산해야하는경우가종종있습니다. 예를들어서, 이차방정식의근의공식 ( b± b 2 4ac)/2a 을계산하는것을생각해봅시다. 이전에도이야기하였지만, b 2 4ac 일경우, 반올림오차로인하여제곱근연산시자리수의절반이상을신뢰할수없게만들수있습니다. 32

39 하지만 b 2 4ac 부분만단정밀도로아닌배정밀도로계산을하게된다면, 배정밀도의절반자리수가부정확하더라도, 단정밀도로변환시에전체자리수를신뢰할수있게되는것입니다. 만일두개의단정밀도값을입력받아서곱셈을한후배정밀도로결과를돌려주는명령어가있다면 a, b, c 가단정밀도일때 b 2 4ac 를배정밀도로계산하는일은어렵지않을것입니다. 정확히반올림되는두개의 p 자리수의곱을계산하기위해서, 어차피곱셈기 (multiplier) 는대략 2p 자리의결과값을생성하기떄문에하드웨어로하여금두개의단정밀도값을인자로받아서결과값을배정밀도로계산하는것은결과값을단정밀도로계산하는것보다약간의비용이더들뿐더러, 배정밀도곱셈기를이용하는것보다는훨씬효과적입니다. 그럼에도불구하고, 대부분의명령어집합은피연산자와같은타입으로결과값을반환하는명령만을지원하는경우가많습니다. 만일두개의단정밀도피연산자를가지고배정밀도결과값을반환하는명령이오직근의공식뿐이라면, 이새로운명령을명령어집합에포함시키는일이별로효과적이지않을수도있습니다. 하지만, 이명령어는그외수많은용도로사용가능합니다. 예를들어다음과같은연립일차방정식들을푸는시스템을생각해봅시다. a 11 x 1 + a 12 x a 1n x n = b 1 a 21 x 1 + a 22 x a 2n x n = b 2 a n1 x 1 + a n2 x a nn x n = b n 이는행렬을이용해서 Ax = b 꼴로다시쓸수있는데, 여기서 a 11 a 12 a 1n a 21 a 22 a 2n A = a n1 a n2 a nn 로쓸수있습니다. 이제, 위방정식의한해 x (1) 를가우시안소거법 (Gaussian elimination) 과같은방법을사용해서구했다고생각해봅시다. 그럼, 그해의정확도를높일수있는간단한방법이하나있는데요, 바로반복향상법 (iterative improvement) 라고부르는방법입니다. 먼저 ζ = Ax (1) b (3.1) 를계산합니다. 그리고, 아래와같은시스템의해를구합니다. Ay = ζ (3.2) 만일 x (1) 이정확한해였다면, ζ 는영벡터였을것이고, 따라서 y 도영벡터였을 것입니다. 하지만, 일반적으로 ζ 와 y 를계산하면서반올림오차가생기기때문에 33

40 Ay ζ Ax (1) b = A(x (1) x) 가됩니다. ( 여기서 x 는알려지지않은참값 ) 그렇다면 y x (1) x 가되기때문에, 해의좀더향상된근사치로 x (2) = x (1) y (3.3) 를생각할수있습니다. 여기서 x (1) 을 x (2) 로바꾸고, x (2) 를 x (3) 으로바꾸는식으로해서세단계 3.1, 3.2, 3.3 를계속반복적으로수행할수있습니다. 이때, 구하게된 x (i+1) 은언제나 x (i) 보다정확하게되는데, 자세한내용은 [Golub and Van Loan 1989] 를참조하시기바랍니다. 반복향상법을사용할때, ζ 벡터의원소들은모두근접한부동소수점수들의차이로구성되어있기때문에, 나쁜지워짐이발생할수있습니다. 따라서 ζ = Ax (1) b 가배정밀도로계산되지않는한큰효과를볼수없겠지요. 이는역시, 두개의단정밀도수 (A, x (1) ) 들의곱셈을배정밀도결과값으로반환해야하는경우입니다. 정리하자면, 두개의부동소수점수를곱한다음그결과를피연산자의두배정밀도로반환하는명령을부동소수점명령어집합에추가하는것은매우효율적일것입니다. 이러한명령어를컴파일러에어떻게구현하는냐는다음섹션에서다루고있습니다 언어와컴파일러 (Languages and Compilers) 컴파일러와부동소수점과의관계는 Farnum [1988] 의논문에서잘찾아볼수있고, 이섹션에서다루는대부분의내용은그논문에서따온것입니다. 모호함 (Ambiguity) 이상적으로컴퓨터언어는언어의모든요소들을정확하게정의하여서프로그램의모든문장들을엄밀하게증명할수있어야만합니다. 이러한사실은언어의정수부분에서참인가운데, 부동소수점의경우엄밀하게정의되어있지않은경우가많습니다. 이는많은언어설계자들이부동소수점부분을정의할때, 부동소수점연산은언제나반올림오차를수반하기때문에엄밀하게증명하는것은불가능한것이라판단한것으로생각됩니다. 하지만우리는이전섹션들에서이러한논리는틀렸다는것을보여왔습니다. 모호함 (Ambiguity) 에서는언어정의들에서공통적으로모호하게다루어지는정의부분에대해논의할것이고, 이들을어떻게처리해야하는지도제안할것입니다. 놀랍게도, 몇몇언어에서는 x 가부동소수점변수일때 ( 예를들어값이 3.0/10.0 이라고합시다 ), 프로그램의모든부분에서 10.0 * x 가같은값을가지도록명확히보장하지않는다는점에있습니다. 예를들어서 Brown 의공리를기반으로한 Ada 라는언어의경우, 부동소수점연산은오직 Brown 의공리만을만족시킬필요가있는것처럼작동하는데, 따라서많은가능한값들중에임의로그들중하나의값을가지면된다고설계되어있습니다. 이는모든부동소수점연산이정확하게정의되어있는 IEEE 34

41 모델과정반대입니다. IEEE 모델에서는우리는 3.0/10.0)*10.0 이 3 과동일하다는사실을 ( 정리 1.5.3) 증명할수있습니다. 하지만 Brown 의모델은그렇지않습니다. 대부분의언어의정의에서나타나는모호함으로오버플로우와언더플로우그리고다른예외에서찾을수있습니다. IEEE 표준에서는명확하게이들각각의예외가발생하였을때의처리를규정해놓았기때문에, 언어측면에서이표준을모델로사용한다면모호함을해결할수있습니다. 또다른모호한부분으로, 괄호를해석하는데관련이있습니다. 반올림오차로인하여, 부동소수점수의경우결합법칙이성립하지않을수있습니다. 예를들어서 (x+y)+z 는 x+(y+z) 와전혀다를수있습니다. ( 예를들어서 x = 10 30, y = 10 30, z = 1) 따라서괄호를정확히지키는것의중요성은정말아무리강조해도지나칠수없습니다. 심지어이전에정리 1.4.1, 1.4.2, 에서보았듯이, 이들은모두괄호에기반을하고있습니다. 예를들어정리 에서식 x h = mx (mx x) 는만일괄호가없었더라면 x h = x 가되어서전체알고리즘을망가뜨렸습니다. 언어측면에서괄호를명확히지키지않는다면, 부동소수점계산은무용지물이될것입니다. 많은컴퓨터언어에서부분식계산은부정확하게정의되어있습니다. 예를들어서 ds 는배정밀도이고 x, y 는단정밀도라고해봅시다. 그렇다면 ds + x*y 는단정밀도로계산되어야할가요, 아니면배정밀도로계산되어야할까요. 다른예제로, x + m/n 에서 m, n 이모두정수일때, 나눗셈을정수나눗셈으로해야할까요, 부동소수점나눗셈으로해야할까요. 이러한문제에대해해결하는방으로두가지방법이있는데, 이두방법모두완벽히만족스럽지는않습니다. 첫번째방법은, 어떠한식에서모든변수들이같은타입이되도록강요시키는것입니다. 이는매우간단한해결책이지만몇가지문제점이있습니다. 예를들어서 Pascal 과같은언어에서는부분범위타입 (subrange type) 이란것이있어서, 부분범위변수들과정수변수들을혼합해서사용하는것을허용하고있습니다. 따라서, 단정밀도와배정밀도변수들을혼합하는것을제한한다는것은이상한일이지요. 다른문제로, 상수들을다룰때, 예를들어서 0.1*x 를계산할때대부분의언어에서는 0.1 을단정밀도상수로해석합니다. 만일프로그래머가모든부동소수점변수의정의를단정밀도에서배정밀도로바꾼다고하면, 0.1 은계속단정밀도상수로취급되기때문에이는컴파일오류를발생시키게됩니다. 따라서프로그래머들은이를고치기위해모든부동소수점상수들을찾아나서야하겠지요. 두번째방법으로는혼합연산을허용하는대신에, 부분식들을계산할때반드시규칙들을정의하는것입니다. 예를들어원래 C 언어정의에따르면모든부동소수점수들은배정밀도로계산되게요구하였습니다 [Kernighan and Ritchie 1978]. 이를통해서이섹션맨처음에나온예제같은문제점들을해결할수있었습니다. 3.0/7.0 과같은식들은모두배정밀도로변환되서계산되는데 q 는단정밀도이기때문에그몫은다시단정밀도로반올림되서저장됩니다. 그런데, 3/7 는이진수로무한히반복되는무한소수이기때문에, 배정밀도로계산된값과, 단정밀도로저장된값에는차이가생기게됩니다. 따라서 q == 3/7 과같은문장은성립하지않게됩니다. 결과적으로, 모든식을가능한높은정밀도로계산하는것은좋은해결책이아닐수있습니다. 위문제에대한해답을제시해줄수있는예로내적 (inner product) 계산이있습니 35

42 식을트리로분해해서생각하는것입니다. 자세한것은 http: //en.wikipedia. org/wiki/binary_ expression_tree 을참고하세요 다. 만일천개가넘는항들의내적계산을수행한다면, 합에서발생하는반올림오차는상당히클수있습니다. 이러한반올림오차를줄이는방법으로이합들을배정밀도로수행하는데있습니다. ( 이는최적화기에서이야기할것입니다 ) 만일 d 가배정밀도변수이고, x[] 와 y[] 가단정밀도배열들이라면, 내적계산루프는 d = d+x[i]*y[i] 처럼생겼겠지요. 만일곱셈이단정밀도로계산된다면, 배정밀도로덧셈을해서얻는이점이줄어들것입니다. 왜냐하면계산된단정밀도의곱은배정밀도변수에더해지기직전에단정밀도로잘려나가기때문입니다. 위두예제를아우르는하나의규칙은바로, 그식에서가장높은정밀도의변수로계산을한다는것입니다. 그렇게된다면 q == 3.0/7.0 은그식에서가장높은정밀도가단정밀도이므로 1), 단정밀도로계산되어식이참이되겠고, 마찬가지로 d = d+x[i]*y[i] 는배정밀도로계산되어배정밀도덧셈의모든효과를누릴수있게됩니다. 하지만, 이규칙은모든경우를깔끔하게처리하기에는너무간단합니다. 예를들어서 dx, dy 가배정밀도변수라면, 식 y = x + single(dx-dy) 는배정밀도변수를가지고있지만, 합을배정밀도로하는것은의미가없게됩니다. 왜냐하면두피연산자모두단정밀도이기때문이죠. 조금더상세한규칙은다음과같습니다. 일단, 각각의연산에잠정적인정밀도를부여합니다. 이잠정적인정밀도는피연산자들의최대정밀도로구성됩니다. 이잠정적인정밀도는식트리 (expression tree) 의잎부터뿌리까지모든부분에해당되게됩니다. 그다음으로, 뿌리에서잎으로쭉읽어들이게되는데, 각각의연산에잠정적정밀도와부모로부터기대되는정밀도를중최대값을할당하게됩니다. q == 3.0/7.0 의경우, 모든잎이단정밀도이기때문에, 모든작업은단정밀도에서작동하게됩니다. d = d+x[i]*y[i] 의경우, 곱셈연산에서잠정적정밀도는단정밀도이지만그다음을읽어들일때부모연산에서피연산자를배정밀도로기대하기때문에배정밀도가되버리게됩니다. 그리고마지막으로 y=x+single(dx-dy) 의경우, 덧셈이단정밀도로수행되게됩니다. Farnum [1988] 은이알고리즘을구현하는데큰어려움이없음을보였습니다. 이규칙의단점은로는, 부분식을계산할때정밀도가이식이포함되어있는전체식에의존된다는것입니다. 이러한사실은몇가지짜증나는문제점을야기할수있는데, 예를들어서여러분이프로그램을디버깅하는데중간식의값을알고싶다고합시다. 하짐나여러분은단순히디버거에그부분식을쳐서그값을계산해달라고요청할수없습니다. 왜냐하면부분식의값은그식이포함되어있는전체식에의존하기때문에그값만따로뽑아올수없게되는것입니다. 마지막으로부분식에대해한가지코멘트해보자면, 십진수상수를이진수로변환하는과정에서부분식계산규칙이십진수상수를해석하는데중요한역할을한다는사실입니다. 이는 0.1 과같은이진수로정확히표현할수없는 ( 무한소수 ) 들을처리하는데중요한사실입니다. 언어의또다른가능한모호함은지수함수를내장연산으로계산하면서발생합니다. 1) 여기서 3.0 은흔히단정밀도상수임을가정합니다. 참고로배정밀도상수의경우 3.0D0 로표현됩니다. 36