Square Root (제곱근)

설명

주어진 값의 제곱근을 구하는 코드입니다.

어떠한 값에 1, 3, 5, 7, 9, ... 와 같이 홀수를 순서대로 뺄셈하다 보면 0보다 작은 값이 되기전까지의 뺄셈을 수행한 횟수가 제곱근이 됩니다.

여기서는 단순히 이러한 논리를 좀더 컴퓨터가 계산하기 용이한 속도최적화에 맞춰서 응용한 구현이라고 소개하면 맞을 듯 싶습니다. 그래서 double형 타입을 통해서 간단히 구현하지는 않고 정수화하여 계산하는 방식을 취합니다.

사용방법


코드

  • 일반적인 SquareRoot (실수연산)
    double hwport_fsqrt(double s_value)
    {
        int s_exponent;
        double s_temp_value;
    
        if(isnan(s_value) != 0) {
            errno = EDOM;
            return(s_value);
        }
    
        if(s_value <= 0.0) {
            errno = EDOM;
            return(0.0);
        }
    
        s_temp_value = frexp(s_value, (int *)(&s_exponent));
        if((s_exponent & 1) != 0) {
            --s_exponent;
            s_temp_value *= 2.0;
        }
    
        s_temp_value = ldexp(s_temp_value + 1.0, (s_exponent >> 1) - 1);
        for(s_exponent = 4 /* HOW TO LOOOPING? */; s_exponent >= 0; s_exponent--) {
            s_temp_value = (s_temp_value + (s_value / s_temp_value)) / 2.0;
        }
    
        return(s_temp_value);
    }
    
  • 정수연산화 SquareRoot
    typedef unsigned long long hwport_ulonglong_t;
    typedef unsigned long long hwport_uintmax_t;
    
    hwport_uintmax_t hwport_isqrt(hwport_uintmax_t s_value)
    {
        hwport_uintmax_t s_result;
        hwport_uintmax_t s_one;
    
        if(s_value <= ((hwport_uintmax_t)0u)) { return((hwport_uintmax_t)0u); }
    
        /* s_one: the second-to-top bit is set */
        s_one = ((hwport_uintmax_t)1u) << ((((hwport_uintmax_t)sizeof(hwport_uintmax_t)) << 3) - ((hwport_uintmax_t)2u));
    
        /* s_one: starts at the highest power of four <= than the s_value */
        while(s_one > s_value) { s_one >>= 2; }
    
        for(s_result = (hwport_uintmax_t)0u;s_one > ((hwport_uintmax_t)0u);) {
            if(s_value >= (s_result + s_one)) {
                s_value -= s_result + s_one;
                s_result += s_one << 1;
            }
    
            s_result >>= 1;
            s_one >>= 2;
        }
    
        /* do arithmetic rounding to nearest integer */
        if(s_value > s_result) { ++s_result; }
    
        return(s_result);
    }
    
    int main(void)
    {
        hwport_uintmax_t s_value;
        hwport_uintmax_t s_value_input;
        hwport_uintmax_t s_value_output;
    
        /* 소숫점 이하 3자리를 구하려면 1000을 기점으로 입력값에는 이의 제곱인 1000000을 곱해주고 결과값은 1000으로 나누면 정수부분, 1000의 나머지는 소숫점 이하부분이 되겠죠. */
    
        s_value = (hwport_uintmax_t)987654321u; /* => 31426.968 */
        s_value_input = s_value * ((hwport_uintmax_t)1000000u);
        s_value_output = hwport_isqrt(s_value_input);
    
        (void)hwport_printf("%llu => %llu => %llu.%03llu\n",
            (hwport_ulonglong_t)s_value,
            (hwport_ulonglong_t)s_value_output,
            (hwport_ulonglong_t)(s_value_output / ((hwport_uintmax_t)1000u)),
            (hwport_ulonglong_t)(s_value_output % ((hwport_uintmax_t)1000u))
        );
    
        return(EXIT_SUCCESS);
    }
    



/*

[ FrontPage | PrintView | RawView | RSS ]

Copyright ⓒ MINZKN.COM
All Rights Reserved.

MINZKN

----

*/