Какая математическая теория лежит в основе алгоритма квадратного корня Apple с фиксированной точкой

Я работаю над проектом, в котором мне нужно вычислить квадратный корень из 16-битного числа. Я использую микроконтроллер NXP, который не поддерживает операции с плавающей запятой. Поискав в Интернете, я нашел следующий подход: en&sa=X&ved=0ahUKEwiH9b7X6ZXXAhWrhlQKHebPDpoQ6AEIMDAB#v=onepage&q=apple%20fixed%20point%20square%20root&f=false" rel="nofollow noreferrer">Квадратный корень с фиксированной точкой Кена Турковски из Apple. Я также нашел тот же алгоритм с некоторыми изменениями в WEB, в том числе здесь: arm-thumb2">Ищу эффективный алгоритм целочисленного квадратного корня для ARM Thumb2 Я копирую код сюда, чтобы сделать его доступным в случае неработающей ссылки:

    void put_space(unsigned int number)
{
    int i;
    char spaces = 12;
    char buffer[15];
    char number_of_chars = sprintf(buffer,"%u",number);
    for(i=0;i<spaces-number_of_chars;i++)
    {
        printf(" ");    
    }   
}
void print_binary(unsigned int number)
{
    int i;
    unsigned char bit;
    for(i=0;i<32;i++)
    {
        if(i%8 == 0 && i!=0)
        {
            printf(" ");
        }
        bit = (number & (2147483648 >> i)) == (2147483648 >> i);
        printf("%u",bit);
    }
}
void print_number(unsigned int number)
{
    print_binary(number);
    printf(" ");
    printf("(%u)",number);
    put_space(number);
}


typedef signed int TFract; /* 2 integer bits, 30 fractional bits */
    TFract FFracSqrt(TFract x)
    {
        register unsigned int root, remHi, remLo, testDiv, count;
        root = 0; /* Clear root */ printf("root: "); print_number(root);        
        remHi = 0; /* Clear high part of partial remainder */  printf("remHi: "); print_number(remHi);
        remLo = x; /* Get argument into low part of partial remainder */ printf("remLo: "); print_number(remLo);
        count = 30; /* Load loop counter */

        printf("\n\n");
    do 
    {
        remHi = (remHi<<2) | (remLo>>30); printf("remHi: "); print_number(remHi);
        remLo <<= 2; /* get 2 bits of arg */ printf("remLo: "); print_number(remLo);
        root <<= 1; /* Get ready for the next bit in the root */ printf("root: "); print_number(root);
        testDiv = (root << 1) + 1; /* Test radical */ printf("testDiv: "); print_number(testDiv);       
        if (remHi >= testDiv) 
        {
            remHi -= testDiv; printf("remHi: "); print_number(remHi);
            root++; printf("root: "); print_number(root);
        }
        printf("\n");
    } while (count-- != 0);

    printf("\n\nResult: %u\n\n",root);  
    return(root);
}

Я могу понять, что происходит в алгоритме, включая сдвиги и операцию ИЛИ. Однако я не могу понять математику этого кода. FFracSqrt получает целое число «x» в качестве параметра и возвращает пол (sqrt (x) * 32768). Как это происходит? Пожалуйста, помогите мне и другим людям понять это, потому что я пока ничего не нашел о теории. Искренне,

Я добавляю здесь вывод, когда вход 2: Перед тем как сделать, в то время как корень петли: 00000000 00000000 00000000 00000000 (0)
remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000000 00000010 (2)

В то время как сделай петли remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000000 00001000 (8)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000000 00100000 (32)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000000 10000000 (128)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00000010 00000000 (512)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00001000 00000000 (2048)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 00100000 00000000 (8192)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000000 10000000 00000000 (32768)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00000010 00000000 00000000 (+131072)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00001000 00000000 00000000 (524288)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 00100000 00000000 00000000 (2097152)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000000 10000000 00000000 00000000 (8388608)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00000010 00000000 00000000 00000000 (33554432)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0) remLo: 00001000 00000000 00000000 00000000 (134217728)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 00100000 00000000 00000000 00000000 (536870912)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000000 (0)
remLo: 10000000 00000000 00000000 00000000 (2147483648)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000010 (2)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000000 00000000 (0)
testDiv: 00000000 00000000 00000000 00000001 (1) Если условие ИСТИНА <бр > remHi: 00000000 00000000 00000000 00000001 (1)
корень: 00000000 00000000 00000000 00000001 (1)

remHi: 00000000 00000000 00000000 00000100 (4)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000000 00000010 (2)
testDiv: 00000000 00000000 00000000 00000101 (5)

remHi: 00000000 00000000 00000000 00010000 (16)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000000 00000100 (4)
testDiv: 00000000 00000000 00000000 00001001 (9) Если условие ИСТИНА remHi: 00000000 00000000 00000000 00000111 (7)
корень: 00000000 00000000 00000000 00000101 (5)

remHi: 00000000 00000000 00000000 00011100 (28)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000000 00001010 (10)
testDiv: 00000000 00000000 00000000 00010101 (21) Если условие ИСТИНА <бр > remHi: 00000000 00000000 00000000 00000111 (7)
корень: 00000000 00000000 00000000 00001011 (11)

remHi: 00000000 00000000 00000000 00011100 (28)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000000 00010110 (22)
testDiv: 00000000 00000000 00000000 00101101 (45)

remHi: 00000000 00000000 00000000 01110000 (112)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000000 00101100 (44)
testDiv: 00000000 00000000 00000000 01011001 (89) Если условие ИСТИНА remHi: 00000000 00000000 00000000 00010111 (23)
корень: 00000000 00000000 00000000 00101101 (45)

remHi: 00000000 00000000 00000000 +01011100 (92)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000000 01011010 (90)
testDiv: 00000000 00000000 00000000 10110101 (181)

remHi: 00000000 00000000 00000001 01110000 (368)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000000 10110100 (180)
testDiv: 00000000 00000000 00000001 01101001 (361) Если условие ИСТИНА <бр > remHi: 00000000 00000000 00000000 00000111 (7)
корень: 00000000 00000000 00000000 10110101 (181)

remHi: 00000000 00000000 00000000 00011100 (28)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000001 01101010 (362)
testDiv: 00000000 00000000 00000010 11010101 (725)

remHi: 00000000 00000000 00000000 +01110000 (112)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000010 11010100 (724)
testDiv: 00000000 00000000 00000101 10101001 (1449)

remHi: 00000000 00000000 00000001 11000000 (448)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00000101 10101000 (1448)
testDiv: 00000000 00000000 00001011 01010001 (2897)

remHi: 00000000 00000000 00000111 00000000 (1792)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00001011 01010000 (2896)
testDiv: 00000000 00000000 00010110 10100001 (5793)

remHi: 00000000 00000000 00011100 00000000 (7168)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00010110 10100000 (5792)
testDiv: 00000000 00000000 00101101 01000001 (11585)

remHi: 00000000 00000000 01110000 00000000 (28672)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 00101101 01000000 (11584)
testDiv: 00000000 00000000 01011010 10000001 (23169) Если условие ИСТИНА <бр > remHi: 00000000 00000000 00010101 01111111 (5503)
корень: 00000000 00000000 00101101 01000001 (11585)

remHi: 00000000 00000000 01010101 11111100 (22012)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 01011010 10000010 (23170)
testDiv: 00000000 00000000 10110101 00000101 (46341)

remHi: 00000000 00000001 01010111 11110000 (88048)
remLo: 00000000 00000000 00000000 00000000 (0)
корень: 00000000 00000000 10110101 00000100 (46340)
testDiv: 00000000 00000001 01101010 00001001 (92681)

Результат: 46340


microcontroller fixed-point square-root
person plinioandrade    schedule 29.10.2017    source источник
comment
ты пройди по математике с карандашом и бумагой?   -  person old_timer    schedule 29.10.2017
comment
Я сделал еще лучше. Я изменил алгоритм, чтобы показать числа в двоичной системе. Как я уже сказал, я знаю, что происходит, но я не могу понять математику.   -  person plinioandrade    schedule 29.10.2017
comment
В Википедии есть хорошая коллекция методов извлечения квадратного корня. Выше выглядит очень много, как двоичная цифра, по-цифр подходить подробно здесь .   -  person John McFarlane    schedule 30.10.2017
comment
Как указал Джон Макфарлейн, этот код реализует многовековой алгоритм длинного квадратного корня, который преподавали в школах до появления карманных калькуляторов. Этот вариант работает с основанием 2, таким образом генерируя один бит результата за итерацию.   -  person njuffa    schedule 31.10.2017

Ответы (1)


arrow_upward
1
arrow_downward

Вот моя попытка с комментариями

TFract FFracSqrt(TFract x)
{
   unsigned root = 0; /* Clear root */

   //remHi & remLo forming a 64bit unsigned integer
   //remHiLo containing the remainder from the radicant not yet transfered into root
   unsigned remHi = 0;
   unsigned remLo = x;

   unsigned bitcount = 30; /* Load loop counter */

   do 
   {
       //shift left the remHiLo package by two bit positions
       remHi = (remHi<<2) | (remLo>>30);
       remLo <<= 2;

       //As we multiplied remHiLo with 4, it is fair to multiply the root with sqrt(4)
       root <<= 1;

       //We just shifted in a zero bit in root. Should it be set?

       //We should check whether (root+1) is too large

       //(root+1) * (root+1)
       //root^2 + 2*root + 1

       //We always deduct root^2 from remHi, so we need to check only the delta
       //between root^2 and (root+1)^2

       unsigned testDiv = (root << 1) + 1;

       if (remHi < testDiv)
       {
          //root+1 would overshoot the radicant, so leave everything as it is
       } 
       else
       {
           root++; //set LSB
           remHi -= testDiv; //update remHi to reflect the new bit in root
       }
    } while (bitcount-- != 0);
    return root;
person user5329483    schedule 31.10.2017