직각삼각형의 빗변의 길이 구하기 / underflow 와 overflow 방지

 

C 언어, Java 언어, Python 언어 등으로 프로그래밍할 때 부동소수점수(floating point number) 를 많이 사용한다. 보동 double 타입이라고 부르는 수치형 데이터인데, 이 데이터의 처리 가능 범위는
대략


      (   -1.79769313486231E308    ~    1.79769313486232E308  )

 

정밀도	비트수	범위	유효숫자 개수(10진수)
단정밀도	32 bits	±1.18×10−38 to ±3.4×1038	약 7개
배정밀도	64 bits	±2.23×10–308 to ±1.80×10308	약 15개

이다.

  그렇다면 1.0E-190 에 해당하는 부동소수점수는 0 이 아니며  이것의 제곱 1.0E-380 도 0이 아니다. 그런데 1.0E-380 을 컴퓨터 계산에서는 어떻게 처리될까? 아마도 0 이라고 하거나, 아주 엉뚱한 값이라고 할 수도 있다. 이와 같이 0 이 아닌 수의 제곱이나 세제곱이 컴퓨터 계산으로는 엉뚱한 값이  될 수 있음을 염두에 두어야 한다. 0은 아니지만 아주 작은 두 수의 곱의 경우에도 이런 현상(언더플로우)을 초래할 수 있다.

  이제 직각삼각형의 빗변의 길이를 구하는 문제를 생가해 보자.
직각을 낀 두 변의 길이가 1 인 직각이등변삼각형의 빗변의 길이는 얼마일까?
당연히 그 답은 sqrt(2) = 1.4142.... 이다,
그러면 직각을 낀 두 변의 길이가 1.0E-190 인 직각이등변삼각형의 빗변의 길이는 얼마일까?
이 때도 당연히 그 답은 sqrt(2)*1.0E-190 = 1.4142....E-190 이다,
컴퓨터 계산 상으로도 그렇게 나올까?
a = 1.0E-190 로 놓고,

        빗변의 길이 = sqrt(a*a + a*a)
                         = sqrt(1.0E-380 + 1.0E-380)
                         = sqrt(0 + 0)   (또는 엉뚱한 값)
                         = 0   (또는 엉뚱한 값)

단순히 수학 시간에 배운대로 피타고리스 정리를 쓰면 위의 계산 처럼 기대하지 않은 엉뚱한 값이 될 수 있다. 그 빗변의 길이 1.4142....E-190 는 배정밀도 부동소수점수로 처리 가능하므로 이런 언더플로우(underflow) 문제를 해결하여야 한다,  두변의 길이가 매우 큰 경우에는 오버플로우가 발생한다, 그래서 직각삼각형의 빗변의 길이를 언도플로우나 오버플로우를 방지하는 다음 알고리즘이 필요하다.

pythagAlgorithm          ; 언더플로우와 오버플로우를 방지
    Input a, b     ; 직각을 낀 두 변의 길이
    a1 <- abs(a)
    b1 <- abs(b)
    if a1 >b1:
        Output a1*sqrt(1.0 + (b1/a1)*(b1/a1))
    else:

if b1 == 0.0:
       Output 0.0
 else:
        Output b1*sqrt(1.0 + (a1/b1)*(a1/b1))  

* 불완전한 C 언어 구현
  (Stabe, c = ...  부분의 출력이 정확하지 못하다. gcc 와 Visual C++ 의 결과가 다르다. )

#include <stdio.h>
#include <math.h>
  
double pythag(double a, double b) {
    double absa = abs(a);
    double absb = abs(b);
    if (absa > absb) {
         return absa*sqrt(1.0 + (absb/absa)*(absb/absa));
    }
    else {
        if (absb == 0.0) {
            return 0.0;
         }
         else {
             return absb*sqrt(1.0 + (absa/absb)*(absa/absb));
         }
    }
}

int main() {
    double a = 1.0E-190;
    double b = 1.0E-190;
    double c = sqrt(a*a + b*b);
    printf("a = %g\n", a);
    printf("b = %g\n", b);
    printf("Unstable, c = sqrt(a*a + b*b) = %g\n", c);
    c = pythag(a, b);
    printf("Stable, c = %g\n", c);
    printf("b*sqrt(2) = %g\n", b*sqrt(2));

    return 0;
}

/*
gcc 로 컴파일하여 실행한 결과:
-----------------------------
a = 1e-190
b = 1e-190
Unstable,  c = sqrt(a, b) = 0
Stable,  c = 1.96204e+09
b*sqrt(2) = 1.41421e-190

Visual C++ 로 컴파일하여 실행한 결과:
-----------------------------------
a = 1e-190
b = 1e-190
Unstable, c = sqrt(a*a + b*b) = 0
Stable, c = 0
b*sqrt(2) = 1.41421e-190
*/

* 수정된 C 언어 소스
  (위의 소스는 절대값을 구하는 abs 함수에 문제가 있었다.
    그래서 절대값 계산하는 ABS 매크로를 만들어 사용하였다 ).

#include <stdio.h>
#include <math.h>
 
#define ABS(x) (((x) < 0) ? (-(x)) : (x))
 
double pythag(double a, double b) {
    double absa = ABS(a);
    double absb = ABS(b);
    if (absa > absb) {
        return absa*sqrt(1.0 + (absb/absa)*(absb/absa));
    }
    else {
        if (absb == 0.0) {
            return 0.0;
        }
        else {
            return absb*sqrt(1.0 + (absa/absb)*(absa/absb));
        }
    }
}

int main() {
    double a = 1.0E-190;
    double b = 1.0E-190;
    double c = sqrt(a*a + b*b);
    printf("a = %g\n", a);
    printf("b = %g\n", b);
    printf("Unstable, c = sqrt(a*a + b*b) = %g\n", c);
    c = pythag(a, b);
    printf("Stable, c = %g\n", c);
    printf("b*sqrt(2) = %g\n", b*sqrt(2));

    return 0;
}

/*
gcc 와 Visual C++ 로 컴파일하여 실행한 결과:
--------------------------------------------
a = 1e-190
b = 1e-190
Unstable, c = sqrt(a*a + b*b) = 0
Stable, c = 1.41421e-190
b*sqrt(2) = 1.41421e-190
*/

* C++ 언어로 직각삼각형의 빗변의 길이 구하기
   ( C++ 언어에서는 위의 C 언어에서의 절대값 문제가 발생하지 않는다. )

#include <iostream>
#include <cmath>
 
using namespace std;
 
double pythag(double a, double b) {
    double absa = abs(a);
    double absb = abs(b);
    if (absa > absb) {
        return absa*sqrt(1.0 + (absb/absa)*(absb/absa));
    }
    else {
        if (absb == 0.0) {
            return 0.0;
        }
        else {
            return absb*sqrt(1.0 + (absa/absb)*(absa/absb));
        }
    }
}

int main() {
    double a = 1.0E-190;
    double b = 1.0E-190;
    double c = sqrt(a*a + b*b);
    cout << "a = " << a << endl;
    cout << "b = " << b << endl;
    cout << "Unstable, c = sqrt(a*a + b*b) = " << c << endl;
    c = pythag(a, b);
    cout << "Stable, c = pythag(a, b) = " << c << endl;
    cout << "b*sqrt(2.0) = " << b*sqrt(2.0) << endl;
   
    return 0;
}

/*
g++와 Visual C++ 로 컴파일하여 실행한 결과:
--------------------------------------------
a = 1e-190
b = 1e-190
Unstable, c = sqrt(a*a + b*b) = 0
Stable, c = 1.41421e-190
b*sqrt(2) = 1.41421e-190
*/

* C# 언어로 직각삼각형의 빗변의 길이 구하기
   ( C# 언어에서는 위의 C 언어에서의 절대값 문제가 발생하지 않는다. )

// Filename: TestPythag.cs
//
//  Purpose:
//        Stable Pythagoras without underflow or overflow.
//
//  Compile: csc TestPythag.cs
//  Execute: TestPythag
//
// Date: 2011. 10. 22 (Sat)

using System;

public class TestPythag {
    public double Pythag(double a, double b) {
  double absa = Math.Abs(a);
  double absb = Math.Abs(b);
  if (absa > absb) {
   return absa*Math.Sqrt(1.0 + (absb/absa)*(absb/absa));
  }
     else {
         if (absb == 0.0) {
          return 0.0;
         }
         else {
          return absb*Math.Sqrt(1.0 + (absa/absb)*(absa/absb));
         }
        }
 }

    public static void Main(string[] args) {
        double a = 1.0E-190;
        double b = 1.0E-190;
        double c = Math.Sqrt(a*a + b*b);
        Console.WriteLine("a = {0}", a);
        Console.WriteLine("b = {0}", b);
        Console.WriteLine("Unstable, c = sqrt(a*a + b*b) = {0}", c);
       
        TestPythag app = new TestPythag();
        c = app.Pythag(a, b);
        Console.WriteLine("Stable, c = Pythag(a, b) = {0}", c);
        Console.WriteLine("b*Math.Sqrt(2) = {0}", b*Math.Sqrt(2));
    }
}

/*
실행 결과:
---------
a = 1E-190
b = 1E-190
Unstable, c = sqrt(a*a + b*b) = 0
Stable, c = Pythag(a, b) = 1.4142135623731E-190
b*Math.Sqrt(2) = 1.4142135623731E-190
*/

* Java 언어로 직각삼각형의 빗변의 길이 안정하게 계산하기.
  ( Java 언어에서도 절대값 문제는 발생하지 않는다. )

// Filename: TestPythag.java
//
//  Purpose:
//        Stable Pythagoras without underflow or overflow.
//
//  Compile: javac TestPythag.java
//  Execute: java TestPythag
//
// Date: 2011. 10. 22 (Sat)

public class TestPythag {
    public double pythag(double a, double b) {
  double absa = Math.abs(a);
  double absb = Math.abs(b);
  if (absa > absb) {
   return absa*Math.sqrt(1.0 + (absb/absa)*(absb/absa));
  }
     else {
         if (absb == 0.0) {
          return 0.0;
         }
         else {
          return absb*Math.sqrt(1.0 + (absa/absb)*(absa/absb));
         }
        }
 }

    public static void main(String[] args) {
        double a = 1.0E-190;
        double b = 1.0E-190;
        double c = Math.sqrt(a*a + b*b);
        System.out.printf("a = %g\n", a);
        System.out.printf("b = %g\n", b);
        System.out.printf("Unstable, c = sqrt(a*a + b*b) = %g\n", c);
       
        TestPythag app = new TestPythag();
        c = app.pythag(a, b);
        System.out.printf("Stable, c = pythag(a, b) = %g\n", c);
        System.out.printf("b*Math.Sqrt(2) = %g\n", b*Math.sqrt(2));
    }
}

/*
Java 6 및 Java 7 에서의 실행 결과:
----------------------------------
a = 1.00000e-190
b = 1.00000e-190
Unstable, c = sqrt(a*a + b*b) = 0.00000
Stable, c = pythag(a, b) = 1.41421e-190
b*Math.Sqrt(2) = 1.41421e-190
*/

* Python 언어로 직각삼각형의 빗변의 길이 안정하게 계산하기

# Filename: testPythag.py
#
#  Purpose:
#        Stable Pythagoras without underflow or overflow.
#
#  Execute: python testPythag.py
#
# Date: 2011. 10. 22 (Sat)

import math

def pythag(a, b):
    absa = abs(a)
    absb = abs(b)
    if absa > absb:
        return absa*math.sqrt(1.0 + (absb/absa)*(absb/absa))
    else:
        if absb == 0.0:
            return 0.0
        else:
            return absb*math.sqrt(1.0 + (absa/absb)*(absa/absb))

a = 1.0E-190
b = 1.0E-190
c = math.sqrt(a*a + b*b)
print "a = %g" % a
print "b = %g" % b
print "Unstable, c = sqrt(a*a + b*b) = %g" % c
       
c = pythag(a, b)
print "Stable, c = pythag(a, b) = %g" % c
print "b*math.sqrt(2) = %g" % (b*math.sqrt(2))

"""
Execution Result:
----------------
a = 1e-190
b = 1e-190
Unstable, c = sqrt(a*a + b*b) = 0
Stable, c = pythag(a, b) = 1.41421e-190
b*math.sqrt(2) = 1.41421e-190
"""

* Ruby 언어로 직각삼감형의 빗변의 길이 안정하게 계산하기

# Filename: testPythag.rb
#
#  Purpose:
#        Stable Pythagoras without underflow or overflow.
#
#  Execute: ruby testPythag.rb
#
# Date: 2011. 10. 22 (Sat)

def pythag(a, b)
    absa = a.abs
    absb = b.abs
    if absa > absb then
        return absa*Math.sqrt(1.0 + (absb/absa)*(absb/absa))
    else
        if absb == 0.0 then
            return 0.0
        else
            return absb*Math.sqrt(1.0 + (absa/absb)*(absa/absb))
        end
    end
end

a = 1.0E-190
b = 1.0E-190
c = Math.sqrt(a*a + b*b)
print "a = %g\n" % a
print "b = %g\n" % b
print "Unstable, c = sqrt(a*a + b*b) = %g\n" % c
       
c = pythag(a, b)
print "Stable, c = pythag(a, b) = %g\n" % c
print "b*Math.sqrt(2) = %g\n" % (b*Math.sqrt(2))

"""
Execution Result:
----------------
a = 1e-190
b = 1e-190
Unstable, c = sqrt(a*a + b*b) = 0
Stable, c = pythag(a, b) = 1.41421e-190
b*Math.sqrt(2) = 1.41421e-190
"""

* Lua 언어로 직각삼각형의 빗변의 길이 안전하게 계산하기

-- Filename: testPythag.lua
--
--  Purpose:
--        Stable Pythagoras without underflow or overflow.
--
--  Execute: lua testPythag.lua
--
-- Date: 2011. 10. 22 (Sat)

function pythag(a, b)
    absa = math.abs(a)
    absb = math.abs(b)
    if absa > absb then
        return absa*math.sqrt(1.0 + (absb/absa)*(absb/absa))
    else
        if absb == 0.0 then
            return 0.0
        else
            return absb*math.sqrt(1.0 + (absa/absb)*(absa/absb))
        end
    end
end

a = 1.0E-190
b = 1.0E-190
c = math.sqrt(a*a + b*b)
print( "a = " .. a )
print( "b = " .. b )
print( "Unstable, c = sqrt(a*a + b*b) = " .. c )
       
c = pythag(a, b)
print( "Stable, c = pythag(a, b) = " .. c )
print( "b*Math.sqrt(2) = " .. (b*math.sqrt(2)) )

--[[
Execution Result:
----------------
a = 1e-190
b = 1e-190
Unstable, c = sqrt(a*a + b*b) = 0
Stable, c = pythag(a, b) = 1.4142135623731e-190
b*Math.sqrt(2) = 1.4142135623731e-190
--]]