/*******************************************************************************
 *  isqrt64.c - Square root of a 64 bit number.
 *
 *  Copyright (C) 2011 Texas Instruments Incorporated - http://www.ti.com/ 
 *  ALL RIGHTS RESERVED 
 * 
*/

#include <inttypes.h>
#include "emeter-toolkit.h"

#if defined(EMETER_TOOLKIT_SUPPORT_64BIT)
uint64_t isqrt64(uint64_t h)
{
    uint64_t x;
    uint64_t y;
    int i;

    /* The answer is calculated as a 64 bit value, where the last
       32 bits are fractional.
       Calling with negative numbers is not a good idea :-) */
#if defined(__GNUC__)  &&  defined(__MSP430__)
    x = 0x8000000000000000ULL;
    y = 0;
    i = 64;
    __asm__ (
        "1: \n"
        " setc \n"
        " rlc   %A[x] \n"
        " rlc   %B[x] \n"
        " rlc   %C[x] \n"
        " rlc   %D[x] \n"
        " sub   %A[x],%A[y] \n"
        " subc  %B[x],%B[y] \n"
        " subc  %C[x],%C[y] \n"
        " subc  %D[x],%D[y] \n"
        " jhs   2f \n"
        " add   %A[x],%A[y] \n"
        " addc  %B[x],%B[y] \n"
        " addc  %C[x],%C[y] \n"
        " addc  %D[x],%D[y] \n"
        " sub   #2,%A[x] \n"
        "2: \n"
        " inc   %A[x] \n"
        " rla   %A[h] \n"
        " rlc   %B[h] \n"
        " rlc   %C[h] \n"
        " rlc   %D[h] \n"
        " rlc   %A[y] \n"
        " rlc   %B[y] \n"
        " rlc   %C[y] \n"
        " rlc   %D[y] \n"
        " rla   %A[h] \n"
        " rlc   %B[h] \n"
        " rlc   %C[h] \n"
        " rlc   %D[h] \n"
        " rlc   %A[y] \n"
        " rlc   %B[y] \n"
        " rlc   %C[y] \n"
        " rlc   %D[y] \n"
        " dec   %[i] \n"
        " jne   1b \n"
        : [x] "+r"(x), [y] "+r"(y)
        : [h] "r"(h), [i] "r"(i));
#else
    x =
    y = 0;
    for (i = 0;  i < 64;  i++)
    {
        x = (x << 1) | 1;
        if (y < x)
            x -= 2;
        else
            y -= x;
        x++;
        y <<= 1;
        if ((h & 0x8000000000000000ULL))
            y |= 1;
        h <<= 1;
        y <<= 1;
        if ((h & 0x8000000000000000ULL))
            y |= 1;
        h <<= 1;
    }
#endif
    return  x;
}
#endif
