Integer square root, round 2

Context

The Newton's method we used to compute integer square root is generally a very efficient algorithm to compute integer square roots.
However, using BrainFuck, the algorithm itself takes time to be processed, a lot of instructions are needed in order to iterate on different approximations.

Another intuitive but generally inefficient way to compute integer square root is to iterate on all integers and compute their squares until we go over our target number.
A bit less naive approach is to compute the difference between 2 squares, and subtract this difference to N while it's positive.
Example: 27

• 1: remove 1, result is 26
• 2: remove 3 (2*2 - 1*1), result is 23
• 3: remove 5 (3*3 - 2*2), result is 18
• 4: remove 7 (4*4 - 3*3), result is 11
• 5: remove 9 (5*5 - 4*4), result is 2
• 6: remove 11 (6*6 - 5*5), result is negative, so integer square root is 5

We'll see that, due to the nature of BrainFuck, this algorithm is actually really more efficient than Newton's method.

Initial state

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

Process

• Store 1 (square difference) on third cell and 1 (actual result) on fourth
• While A is not null
• remove 1
• remove 1 to square difference
• If square difference is null
• Copy actual result's double on third cell
• Increase both actual result and cell square difference by 1
• Decrease result by 1 (because the current result is over our target)

Code - try it

start from 1
>>>+>+<<<<

while N is not null
[
  decrease N
  -
  set else bit
  >>+
  decrease square difference
  >-
  if current square difference is null
  [<-]<[-
    reinit square difference
    copy current integer double
    >>[-<++>>+<]>[-<+>]<
    increase both integer and square difference by 1
    +<+
    <
  <]
  <
]
decrease result by 1 and print
>>>[-]>-
[>>>>++++++++++<<<<[->+>>+>-[<-]<[->>+<<<<[->>>+<<<]>]<<]>+[-<+>]>>>[-]>[-<<<<+>>>>]<<<<]<[>++++++[<++++++++>-]<-.[-]<]

Final state

• Memory: 0 0 0 D R with R=isqrt(A) and D the remaining part of the square difference 
• Cursor: first cell
• Input: unchanged
• Output: unchanged

This algorithm takes about 15 times less instructions to be executed...

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