Цілочисельний квадратний корінь

Цілочисельний квадратний корінь (isqrt) із невід’ємного цілого числа n — це невід’ємне ціле число m, яке є найбільшим цілим числом, меншим або рівним квадратному кореню з n ,

${\displaystyle {\mbox{isqrt}}(n)=\lfloor {\sqrt {n}}\rfloor .}$

Наприклад, ${\displaystyle {\mbox{isqrt}}(27)=\lfloor {\sqrt {27}}\rfloor =\lfloor 5.19615242270663...\rfloor =5.}$

Вступне зауваження

Нехай ${\displaystyle y}$ та ${\displaystyle k}$ є цілими невід’ємними числами.

Алгоритми, які обчислюють (^[en]) ${\displaystyle {\sqrt {y}}}$ працюють вічно на кожному вході ${\displaystyle y}$ який не є ідеальним квадратом.

Алгоритми, що обчислюють ${\displaystyle \lfloor {\sqrt {y}}\rfloor }$ працюють не вічно. Проте вони можуть обчислювати ${\displaystyle {\sqrt {y}}}$ з будь-якою бажаною точністю ${\displaystyle k}$ .

Можна обрати будь-яке ${\displaystyle k}$ і обчислити ${\displaystyle \lfloor {\sqrt {y\times 100^{k}}}\rfloor }$ .

Наприклад (нехай ${\displaystyle y=2}$ ):

${\displaystyle {\begin{aligned}&k=0:\lfloor {\sqrt {2\times 100^{0}}}\rfloor =\lfloor {\sqrt {2}}\rfloor =1\\&k=1:\lfloor {\sqrt {2\times 100^{1}}}\rfloor =\lfloor {\sqrt {200}}\rfloor =14\\&k=2:\lfloor {\sqrt {2\times 100^{2}}}\rfloor =\lfloor {\sqrt {20000}}\rfloor =141\\&k=3:\lfloor {\sqrt {2\times 100^{3}}}\rfloor =\lfloor {\sqrt {2000000}}\rfloor =1414\\&...\\&k=8:\lfloor {\sqrt {2\times 100^{8}}}\rfloor =\lfloor {\sqrt {20000000000000000}}\rfloor =141421356\\&...\\\end{aligned}}}$

Порівняйте результати з ${\displaystyle {\sqrt {2}}=1.41421356237309504880168872420969807856967187537694...}$

Виявляється, що множення входу алгоритму на ${\displaystyle 100^{k}}$ дає точність ${\displaystyle k}$ десяткових цифр.

Щоб обчислити (повне) десяткове представлення ${\displaystyle {\sqrt {y}}}$, можна виконати ${\displaystyle {\mbox{isqrt}}(y)}$ нескінченну кількість разів, збільшуючи ${\displaystyle y}$ в 100 разів при кожному проході.

Припустимо, що в наступній програмі ( ${\displaystyle {\mbox{sqrtForever}}}$ ) процедура ${\displaystyle {\mbox{isqrt}}(y)}$ вже визначена, і — заради аргументації — що всі змінні можуть містити цілі числа необмеженої величини.

Тоді ${\displaystyle {\mbox{sqrtForever}}(y)}$ надрукує повне десяткове представлення ${\displaystyle {\sqrt {y}}}$.

// Друкує sqrt(y) без зупинки void sqrtForever( unsigned int y ) { unsigned int result = isqrt( y ); printf( "%d.", result );// print result, followed by a decimal point while( true )// повторюється вічно ... { y = y * 100;// теоретичний приклад: переповнення ігнорується result = isqrt( y ); printf( "%d", result % 10 );// надрукувати останню цифру результата } } 

Висновок полягає в тому, що алгоритми, які обчислюють isqrt(), обчислювально еквівалентні алгоритмам, які обчислюють sqrt().

Основні алгоритми

Алгоритм з використанням лінійного пошуку

Наступні програми мовою C є простими реалізаціями.  

Лінійний пошук

За допомогою додавання

У наведеній вище програмі (лінійний пошук, за зростанням) можна замінити множення на додавання, використовуючи еквівалентність

${\displaystyle (L+1)^{2}=L^{2}+2L+1=L^{2}+1+\sum _{i=1}^{L}2}$ .

// Цілочислельний квадратний корінь // (лінійний пошук, за зростанням) використовуючи додавання unsigned int isqrt( unsigned int y ) { unsigned int L = 0; unsigned int a = 1; unsigned int d = 3; while( a <= y ) { a = a + d;// (a + 1)^2 d = d + 2; L = L + 1; } return L; } 

За допомогою віднімання

Jarvis, (2006) опублікував метод отримання цілочисельного квадратного кореня з використанням виключно віднімання.

Розрахунок ${\displaystyle {\mbox{isqrt}}(y)}$, де ${\displaystyle y}$ є невід’ємним цілим числом, є окремим випадком, який не потребує повного алгоритму. Досить наступної програми.

// Цілочислельний квадратний корінь (використовуючи віднімання) unsigned int isqrt( unsigned int y ) { unsigned int a = 5 * y; unsigned int b = 0; while( a >= b ) { a = a - b; b = b + 10; } return b / 10; } 

Числовий приклад

Наприклад, якщо обчислити ${\displaystyle {\mbox{isqrt}}(200)}$ використовуючи метод «віднімання», отримуємо ${\displaystyle [a,b]}$ послідовність

${\displaystyle {\begin{aligned}&[1000,0]\rightarrow [990,10]\rightarrow [980,20]\rightarrow [960,30]\rightarrow [930,40]\rightarrow [890,50]\rightarrow [840,60]\rightarrow [780,70]\rightarrow [710,80]\\&\rightarrow [630,90]\rightarrow [540,100]\rightarrow [440,110]\rightarrow [330,120]\rightarrow [210,130]\rightarrow [80,140].\end{aligned}}}$

${\displaystyle 140/10=14=\lfloor {\sqrt {200}}\rfloor }$

Це обчислення займає 14 кроків ітерації, як і лінійний пошук (за зростанням, починаючи з ${\displaystyle 0}$ ).

Алгоритм із використанням бінарного пошуку

Лінійний пошук послідовно перевіряє кожне значення, доки воно не досягне найменшого ${\displaystyle x}$ де ${\displaystyle x^{2}>y}$ .

Прискорення досягається за допомогою двійкового пошуку. Наступна програма мово C є реалізацією.

// Цілочисильний квадратний корінь (використовуючи бінарний пошук) unsigned int isqrt( unsigned int y ) { unsigned int L = 0; unsigned int M; unsigned int R = y + 1;  while( L != R - 1 )  {  M = (L + R) / 2; if( M * M <= y ) L = M; else R = M; }  return L; } 

Числовий приклад

Наприклад, якщо обчислити ${\displaystyle {\mbox{isqrt}}(2000000)}$ використовуючи двійковий пошук, можна отримати ${\displaystyle [L,R]}$ послідовність

${\displaystyle {\begin{aligned}&[0,2000001]\rightarrow [0,1000000]\rightarrow [0,500000]\rightarrow [0,250000]\rightarrow [0,125000]\rightarrow [0,62500]\rightarrow [0,31250]\rightarrow [0,15625]\\&\rightarrow [0,7812]\rightarrow [0,3906]\rightarrow [0,1953]\rightarrow [976,1953]\rightarrow [976,1464]\rightarrow [1220,1464]\rightarrow [1342,1464]\rightarrow [1403,1464]\\&\rightarrow [1403,1433]\rightarrow [1403,1418]\rightarrow [1410,1418]\rightarrow [1414,1418]\rightarrow [1414,1416]\rightarrow [1414,1415]\end{aligned}}}$

Це обчислення займає 21 крок ітерації, тоді як лінійний пошук (за зростанням, починаючи з ${\displaystyle 0}$ ) потрібно 1414 кроків.

Алгоритм з використанням методу Ньютона

Один спосіб розрахунку ${\displaystyle {\sqrt {n}}}$ і ${\displaystyle {\mbox{isqrt}}(n)}$ полягає у використанні ^[en], який є окремим випадком методу Ньютона, щоб знайти розв’язок рівняння ${\displaystyle x^{2}-n=0}$, за даною ітераційною формулою

${\displaystyle {x}_{k+1}={\frac {1}{2}}\left(x_{k}+{\frac {n}{x_{k}}}\right),\quad k\geq 0,\quad x_{0}>0.}$

Послідовність ${\displaystyle \{x_{k}\}}$ сходиться квадратично до ${\displaystyle {\sqrt {n}}}$ коли${\displaystyle k\to \infty }$ .

Критерій зупинки

Можна довести  це ${\displaystyle c=1}$ є найбільшим можливим числом, для якого критерій зупинки

${\displaystyle |x_{k+1}-x_{k}|<c}$

забезпечує ${\displaystyle \lfloor x_{k+1}\rfloor =\lfloor {\sqrt {n}}\rfloor }$ в наведеному вище алгоритмі.

У реалізаціях, які використовують числові формати, які не можуть точно представити всі раціональні числа (наприклад, з плаваючою комою), для захисту від помилок округлення слід використовувати константу зупинки, меншу за одиницю.

Область обчислень

Хоча ${\displaystyle {\sqrt {n}}}$ для багатьох ${\displaystyle n}$ є ірраціональним, послідовність ${\displaystyle \{x_{k}\}}$ містить лише раціональні значення, коли ${\displaystyle x_{0}}$ є раціональним. Таким чином, за допомогою цього методу немає необхідності виходити з поля раціональних чисел для обчислення ${\displaystyle {\mbox{isqrt}}(n)}$, факт, який має деякі теоретичні переваги.

Використання тільки цілочисельного ділення

Для обчислень ${\displaystyle \lfloor {\sqrt {n}}\rfloor }$ для дуже великих цілих чисел n можна використовувати евклідове ділення для обох операцій ділення. Це має перевагу в тому, що для кожного проміжного значення використовуються лише цілі числа, що робить непотрібним представлення великих чисел з плаваючою комою. Це еквівалентно використанню ітераційної формули

${\displaystyle {x}_{k+1}=\left\lfloor {\frac {1}{2}}\left(x_{k}+\left\lfloor {\frac {n}{x_{k}}}\right\rfloor \right)\right\rfloor ,\quad k\geq 0,\quad x_{0}>0,\quad x_{0}\in \mathbb {Z} .}$

Використовуючи той факт, що

${\displaystyle \left\lfloor {\frac {1}{2}}\left(x_{k}+\left\lfloor {\frac {n}{x_{k}}}\right\rfloor \right)\right\rfloor =\left\lfloor {\frac {1}{2}}\left(x_{k}+{\frac {n}{x_{k}}}\right)\right\rfloor ,}$

можна показати, що буде досягнуто ${\displaystyle \lfloor {\sqrt {n}}\rfloor }$ протягом скінченної кількості ітерацій.

У оригінальній версії дано ${\displaystyle {x}_{k}\geq {\sqrt {n}}}$ для ${\displaystyle k\geq 1}$, і ${\displaystyle {x}_{k}>{x}_{k+1}}$ для ${\displaystyle {x}_{k}>{\sqrt {n}}}$ . Отже, у цілочисельній версії є ${\displaystyle \lfloor {x}_{k}\rfloor \geq \lfloor {\sqrt {n}}\rfloor }$ і ${\displaystyle {x}_{k}\geq \lfloor {x}_{k}\rfloor >{x}_{k+1}\geq \lfloor {x}_{k+1}\rfloor }$ поки остаточне рішення ${\displaystyle {x}_{s}}$ не буде досягнуто. Для остаточного рішення ${\displaystyle {x}_{s}}$ виконується ${\displaystyle \lfloor {\sqrt {n}}\rfloor \leq \lfloor {x}_{s}\rfloor \leq {\sqrt {n}}}$ і ${\displaystyle \lfloor {x}_{s+1}\rfloor \geq \lfloor {x}_{s}\rfloor }$, тому критерій зупинки такий ${\displaystyle \lfloor {x}_{k+1}\rfloor \geq \lfloor {x}_{k}\rfloor }$.

Однак, ${\displaystyle \lfloor {\sqrt {n}}\rfloor }$ не обов’язково є нерухомою точкою наведеної вище ітераційної формули. Дійсно, можна показати, що ${\displaystyle \lfloor {\sqrt {n}}\rfloor }$ є фіксованою точкою тоді і тільки тоді, коли ${\displaystyle n+1}$ не є ідеальним квадратом. Якщо ${\displaystyle n+1}$ є ідеальним квадратом, послідовність закінчується циклом з періодом два між ${\displaystyle \lfloor {\sqrt {n}}\rfloor }$ і ${\displaystyle \lfloor {\sqrt {n}}\rfloor +1}$ замість того, щоб сходитись.

Приклад реалізації на C

// Ціллочисельний квадратний корінь unsigned int int_sqrt ( unsigned int s ) { unsigned int x0 = s / 2;// Початкова оцінка  // Уникнути переповнгення коли s є максимальним значенням, яке можна представити value // Перевірка if ( x0 != 0 ) { unsigned int x1 = ( x0 + s / x0 ) / 2;// Оновлення  while ( x1 < x0 )// Це також перевірка цикла { x0 = x1; x1 = ( x0 + s / x0 ) / 2; }  return x0; } else { return s; } } 

Числовий приклад

Наприклад, якщо обчислити цілочисельний квадратний корінь з ${\displaystyle 2000000}$ використовуючи наведений вище алгоритм, можна отримати послідовність

${\displaystyle {\begin{aligned}&1000000\rightarrow 500001\rightarrow 250002\rightarrow 125004\rightarrow 62509\rightarrow 31270\rightarrow 15666\rightarrow 7896\\&\rightarrow 4074\rightarrow 2282\rightarrow 1579\rightarrow 1422\rightarrow 1414\rightarrow 1414\end{aligned}}}$

Загалом потрібно 13 кроків ітерації. Хоча метод Герона сходиться квадратично близько до розв’язку, на початку досягається менше ніж один біт точності на ітерацію. Це означає, що вибір початкової оцінки є критичним для продуктивності алгоритму.

Якщо доступне швидке обчислення для цілої частини двійкового логарифма або ^[en] (наприклад, std::bit_width у ), краще починати з

${\displaystyle x_{0}=2^{\lfloor (\log _{2}n)/2\rfloor +1}}$ ,

що є найменшим ступенем двійки, більшим за ${\displaystyle {\sqrt {n}}}$. У прикладі цілочисельного квадратного кореня з ${\displaystyle 2000000}$, ${\displaystyle \lfloor \log _{2}n\rfloor =20}$, ${\displaystyle x_{0}=2^{11}=2048}$, а результуюча послідовність дорівнює

${\displaystyle 2048\rightarrow 1512\rightarrow 1417\rightarrow 1414\rightarrow 1414}$ .

У цьому випадку потрібні лише 4 кроки ітерації.

Поцифровий алгоритм

Традиційний ^[en] для обчислення квадратного кореня ${\displaystyle {\sqrt {n}}}$ заснований на роботі від старших розрядів до нижчих, і з кожною новою цифрою вибирається найбільша, яка все ще дасть квадрат ${\displaystyle \leq n}$ . Якщо зупинитися після одиниці, обчисленим результатом буде цілочисельний квадратний корінь.

Використання побітових операцій

При роботі в за основою 2 вибір цифри спрощується до числа між 0 («маленький кандидат») і 1 («великий кандидат»), а маніпуляції з цифрами можна виразити в термінах операцій двійкового зсуву. З * — це множення, << — зсув вліво, а >> — логічний зсув вправо, рекурсивний алгоритм для знаходження цілочисельного квадратного кореня з будь-якого натурального числа виглядає так:

def integer_sqrt(n: int) -> int: assert n >= 0, "sqrt працює тільки для невід'ємних вхідних даних" if n < 2: return n # Рекрсивний виклик: small_cand = integer_sqrt(n >> 2) << 1 large_cand = small_cand + 1 if large_cand * large_cand > n: return small_cand else: return large_cand # еквівалентно: def integer_sqrt_iter(n: int) -> int: assert n >= 0, "sqrt працює тільки для невід'ємних вхідних даних" if n < 2: return n # Find the shift amount. See also [[find first set]], # shift = ceil(log2(n) * 0.5) * 2 = ceil(ffs(n) * 0.5) * 2 shift = 2 while (n >> shift) != 0: shift += 2 # Unroll the bit-setting loop. result = 0 while shift >= 0: result = result << 1 large_cand = ( result + 1 ) # Same as result ^ 1 (xor), because the last bit is always 0. if large_cand * large_cand <= n >> shift: result = large_cand shift -= 2 return result 

Представлення традиційних ручних порозрядних алгоритмів включають різні оптимізації, яких немає в коді вище, зокрема трюк попереднього віднімання квадрата попередніх цифр, що робить непотрібним загальний крок множення. Див. для прикладу.

У мовах програмування

Деякі мови програмування мають явну операцію обчисленню цілочисельного квадратного кореня на додаток до загального випадку або можуть бути розширені бібліотеками для цієї мети.

• (isqrt x) : Common Lisp.
• math.isqrt(x) : Python.

Див. також

• ^[en]

Примітки

1. ^[en]
2. It is no surprise that the repeated multiplication by ${\displaystyle 100}$ is a feature in Jarvis, (2006)
3. The fractional part of square roots of perfect squares is rendered as ${\displaystyle 000...}$.
4. see ''.
5. The method was defined for the positive integers. Te program below computes ${\displaystyle {\mbox{isqrt}}(0)=0}$ as well
6. Woo, 1985.
7. LispWorks Ltd, 1996.
8. Python Software Foundation, 2001.

Посилання

 

• Jarvis, Ashley Frazer (2006). Square roots by subtraction (PDF). Mathematical Spectrum. 37: 119—122.
• LispWorks Ltd (1996). CLHS: Function SQRT, ISQRT. www.lispworks.com.
• Minsky, Marvin (1967). 9. The Computable Real Numbers. Computation: Finite and Infinite Machines. Prentice-Hall. ISBN . OCLC 0131655639.
• Python Software Foundation (2001). Mathematical functions. Python Standard Library documentation (since version 3.8)
• Woo, C (June 1985). .
• A geometric view of the square root algorithm.{{}}: Обслуговування CS1: Сторінки з параметром url-status, але без параметра archive-url ()