Intop-SquareRoot impl

This program is free software.  You can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.

Author

© 2012 Jörg Seebohn

This explains the mathematical background of the high school algorithm.

Definition

The square of a number n is the sum of the first n uneven numbers

n² = 1 + 3 + ... + 2n-1

Proof by Induction

1² = 1
2² = 1 + 3
Now (n²=1+3+...+2n-1) is true
==> (n+1)² = n² + 2n + 1
           = (1+3+...+2n-1) + 2(n+1)-1
==> (n+1)² = (1+3+...+2(n+1)-3 + 2(n+1)-1)

Definition

1.  B is the base of the number (10 for decimal, 2 for binary numbers) 2. a,b,c are numbers in base B and d₀,d₁,d₂... digits in base B 3. sqrt(N) = a * B⁽ⁿ⁻¹⁾ + b * B⁽ⁿ⁻²⁾ + c 4. c = d₍ₙ₋₃₎B⁽ⁿ⁻³⁾ + ...  + d₁B + d₀

Now with these definitions N is

N = (a * B⁽ⁿ⁻¹⁾ + b * B⁽ⁿ⁻²⁾ + c)²
  = a²B⁽²ⁿ⁻²⁾ + 2abB⁽²ⁿ⁻³⁾ + b²B⁽²ⁿ⁻⁴⁾ + 2acB⁽ⁿ⁻¹⁾ + 2bcB⁽ⁿ⁻²⁾ + c²

Partitioning N into pairs of digits and determine a

    ╭─────────┬───────────┬─────────┬─────╮
N:  │d₂ₙ|d₂ₙ₋₁│d₂ₙ₋₂|d₂ₙ₋₃│   ...   │d₁|d₀│
    ╰─────────┴───────────┴─────────┴─────╯
    ╭─────────┬───────────────────────────╮
    │   a²    │     B⁽²ⁿ⁻²⁾               │
    ╰─────────┴───────────────────────────╯

You can determine a by adding the first a uneven numbers until (a+1)² > (d₂ₙ* B + d₂ₙ₋₁)

The reason that you can do so is that the values of b and c does not determine the value of a.  Choosing instead (a+1) and b=0 and c=0 would result in

(a+1)² * B⁽²ⁿ⁻²⁾ > N

You can proof this by setting b and c to their possible maximum values and than add all terms and simplify them so that (a+1)² * B⁽²ⁿ⁻²⁾ is bigger than N.

b = B-1
c = B⁽ⁿ⁻⁴⁾ - 1
N = a²B⁽²ⁿ⁻²⁾ + 2abB⁽²ⁿ⁻³⁾ + b²B⁽²ⁿ⁻⁴⁾ + 2acB⁽ⁿ⁻¹⁾ + 2bcB⁽ⁿ⁻²⁾ + c²
  = a²B⁽²ⁿ⁻²⁾ + 2a(B-1)B⁽²ⁿ⁻³⁾ + (B-1)²B⁽²ⁿ⁻⁴⁾ + 2a(B⁽ⁿ⁻⁴⁾ - 1)B⁽ⁿ⁻¹⁾ + 2(B-1)(B⁽ⁿ⁻⁴⁾-1)B⁽ⁿ⁻²⁾ + (B⁽ⁿ⁻⁴⁾ - 1)²
  = a²B⁽²ⁿ⁻²⁾ + 2aB⁽²ⁿ⁻²⁾-2aB⁽²ⁿ⁻³⁾ +(B²-2B+1)B⁽²ⁿ⁻⁴⁾ + 2aB²ⁿ⁻⁵-2aB⁽ⁿ⁻¹⁾ + 2(B-1)(B⁽ⁿ⁻⁴⁾-1)B⁽ⁿ⁻²⁾ + (B⁽²ⁿ⁻⁸⁾ - 2B⁽ⁿ⁻⁴⁾ + 1)
With B >= 2
==> N < a²B⁽²ⁿ⁻²⁾ + 2aB⁽²ⁿ⁻²⁾-2aB⁽²ⁿ⁻³⁾ +(B²-3)B⁽²ⁿ⁻⁴⁾+ 2aB²ⁿ⁻⁵-2aB⁽ⁿ⁻¹⁾ + 2(B⁽ⁿ⁻⁴⁾-1)B⁽ⁿ⁻²⁾ + (B⁽²ⁿ⁻⁸⁾ - 2B⁽ⁿ⁻⁴⁾ + 1)
      < a²B⁽²ⁿ⁻²⁾ + 2aB⁽²ⁿ⁻²⁾ + B⁽²ⁿ⁻²⁾ -2aB⁽²ⁿ⁻³⁾ -3B⁽²ⁿ⁻⁴⁾ + 2aB²ⁿ⁻⁵-2aB⁽ⁿ⁻¹⁾+ 2B⁽²ⁿ⁻⁶⁾ + (B⁽²ⁿ⁻⁸⁾ - 2B⁽ⁿ⁻⁴⁾ + 1)
With -2aB⁽²ⁿ⁻³⁾ - 2aB⁽ⁿ⁻¹⁾ <= 0 <= 2aB²ⁿ⁻⁵
and  -3B⁽²ⁿ⁻⁴⁾ - 2B⁽ⁿ⁻⁴⁾ < 0 < 2B⁽²ⁿ⁻⁶⁾ + B⁽²ⁿ⁻⁸⁾ + 1
==>   N < a²B⁽²ⁿ⁻²⁾ + 2aB⁽²ⁿ⁻²⁾ + B⁽²ⁿ⁻²⁾
      N < (a+1)² * B⁽²ⁿ⁻²⁾

Determine b

      ╭─────────┬───────────┬─────────┬─────╮
N-a²: │r₂ₙ|r₂ₙ₋₁│d₂ₙ₋₂|d₂ₙ₋₃│   ...   │d₁|d₀│
      ╰─────────┴───────────┴─────────┴─────╯
      ╭─────────────────────┬───────────────╮
      │           2ab |  0  │   B⁽²ⁿ⁻⁴⁾     │
      │         │   + b²    │               │
      ╰─────────────────────┴───────────────╯

You can determine b by adding the first b uneven numbers + b times 2aB until (b+1)² + 2a(b+1)B > (r₂ₙB³+r₂ₙ₋₁B²+d₂ₙ₋₂B+d₂ₙ₋₃)

The reason that you can do so is that the value of c does not determine the value of b.  Choosing instead (b+1) and c=0 would result in

a²B⁽²ⁿ⁻²⁾ + 2a(b+1)B⁽²ⁿ⁻³⁾ + (b+1)²B⁽²ⁿ⁻⁴⁾ > N

Shift

       ╭─────────┬───────────┬─────────┬─────╮
N-...: │r₂ₙ|r₂ₙ₋₁│r₂ₙ₋₂|r₂ₙ₋₃│   ...   │d₁|d₀│
       ╰─────────┴───────────┴─────────┴─────╯

Then combine the digits a and b into a₂=aB+b and compute N₂

N₂ = (a₂ * B⁽ⁿ⁻²⁾ + b₂ * B⁽ⁿ⁻³⁾ + c₂)²

The number a₂ is known so you can compute b₂ as done for the number b.  Then combine the digits a₂ and b₂ into a₃=a₂B+b₂ and compute N₃ and so on.  Repeat this shift step n-2 times.

aₙ is the value of sqrt(N).

If you do the whole computation in base 2 with binary numbers a computed digit is either 0 or 1.  Therefore the test becomes

1 + 4a > (r₂ₙB³+r₂ₙ₋₁B²+d₂ₙ₋₂B+d₂ₙ₋₃)