Integer square root using Newton's method

Context

The integer square root is the operation that returns the integer part of a square root. Namely, isqrt(n) is the biggest integer k so that (k+1)^2 > n.
In most programming languages, there is a fast implementation for this operation, using Newton's method.

Considering a function f, so that f(z) = 0 and let's find z: Newton's method states that, starting with a 'good' value x0, we can define x1 = x0 - f(x0)/f'(x0), and use this formula as a recursive definition, to finally converge to z.

Applied to square root, we have a function f(x) = x^2-n, which is null for x = sqrt(n). So, starting with x0=n, and using the recurrence formula x(k+1) = (x(k) + n/(x(k)) / 2, we can converge to n's square root.

And applied to integers, this gives us:

• Start with X = N
• Define Y = (X+N/X) / 2
• N/X is euclidian quotient
• division by 2 can be done using a shift...
• If X > Y then assign value Y to variable X and loop
• Returns X

Example: sqrt(57)

• X = 57, Y = (57 + 57/57) / 2 = 29
• X = 29, Y = (29 + 57/29) / 2 = 15
• X = 15, Y = (15 + 57/15) / 2 = 9
• X = 9, Y = (9 + 57/9 ) / 2 = 7
• X = 7, Y = (7 + 57/7) / 2 = 7
• Result is 7, and actually 7x7 = 49 < 57 < 8*8 = 64

Example (limit case): sqrt(63)

• X = 63, Y = (63 + 63/63) / 2 = 32
• X = 32, Y = (32 + 63/32) / 2 = 16
• X = 16, Y = (16 + 63/16) / 2 = 9
• X = 9, Y = (9 + 63/9) / 2 = 8
• X = 8, Y = (8 + 63/8) / 2 = 7
• X = 7, Y = (7 + 63/7) / 2 = 8
• Result is 7, and we can see here that X = Y is not a valid stop condition. It has to be X ≤ Y !

Initial state

• Memory: N 0 0 0 0 0 0 0 0 0 0
• Cursor: first cell
• Input: any

Process

• Here is what our memory should looks like when starting a loop
• N [loop indicator] [current value X] 0 0 0 0 0 0 0 0
• Set loop indicator to 1 and current value to N
• Then, while loop indicator is 1
• Reset loop indicator
• Note
• X will be needed
• once to compute Y
• once to be compared to Y
• once to be the final result, if needed
• N will be needed
• once to compute Y
• Hence, copy N once and X twice to get a memory like this
• N 0 X 0 X N 0 0 0 X 0
• Divide N by X. Division also produces R and R'
• R being the remainder
• and R' so that X = R+R'
• N 0 X 0 X 0 R 0 0 R' N/X
• Compute 2Y, by adding X to N/X, so by adding R and R' to N/X
• N 0 X 0 X 0 0 0 0 0 2Y
• Note: Y will be needed
• once to be compared to X
• once to replace X if needed
• Hence, divide 2Y by 2 and store the results at the right places
• N 0 X Y X 0 0 0 Y 0 0
• Compare Y and X
• N 0 X res 0 0 0 0 Y 0 0
• If res = 1
• set loop indicator to 1
• reset X
• copy Y where X stands

Code - try it

initialize memory
[->+>+<<]>[-<+>]+
Memory : N 1 X 0 0 0 0 0 0 0 0
[
  reset loop indicator
  -
  copy X twice
  >[->+>+>>>>>+<<<<<<<]>[-<+>]
  copy N
  <<<[->+>>>>+<<<<<]>[-<+>]
  Memory: N 0 X 0 X N 0 0 0 X 0

  divide N by X
  >>>>[->+>>+>-[<-]<[->>+<<<<[->>>+<<<]>]<<]
  Memory: N 0 X 0 X 0 R 0 0 R' N/X

  add X to quotient
  >[->>>+<<<]>>>[->+<]
  Memory: N 0 X 0 X 0 0 0 0 0 2Y

  divide 2Y by 2 and store result twice)
  >[-[-<<+<<<<<+>>>>>>]<[>]>]
  Memory: N 0 X Y X 0 0 0 Y 0 0

  compare Y less than X
  <<<<<+<<[->[->]>[<<[-]>>->>]<<+<<]>[[-]<+>]>-<<
  Memory: N 0 X res 0 0 0 0 Y 0 0

  if true replace X by Y and loop
  [
    set loop indicator
    -<<+
    reset X
    >[-]
    replace by Y
    >>>>>>[-<<<<<<+>>>>>>]<<<<<
  ]
  go back to loop indicator
  <<
]

Final state

• Memory: N 0 R 0 0 0 0 0 R' 0 0 where R = isqrt(N) and R' = R unless N = R^2 + 1 (then, R' = R+1)
• Cursor: second cell
• Input: unchanged
• Output: unchanged