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);
}

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);
}