lut.py

from math import sqrt

def generate_lookup_table() -> str:
    roots = [list[int]() for _ in range(32)]
    for i in range(256):
        roots[i >> 3].append(sqrt(i << 8))
    table = [sum(roots[i]) / len(roots[i]) for i in range(32)]

    # Convert to a single 32-byte value
    lookup_value = 0
    for i, value in enumerate(table):
        lookup_value |= round(value) << (8 * (31 - i))

    return f"0x{lookup_value:064x}"

print(generate_lookup_table())

Sqrt.opcodes

// Efficient SQRT method
// Constant gas cost of 293 discounting RETURN and CALLDATALOAD logic, or 311 including everything.
//
// Special thanks to @caironeth for have found a better log2(x) here:
//  - https://github.com/Lohann/openzeppelin-contracts/pull/1
//
// For testing use this tool: https://www.evm.codes/playground?fork=cancun
// Authors:
// - Lohann Paterno Coutinho Ferreira <developer@lohann.dev>
// - Cairo <https://github.com/cairoeth>
// - Tiago de Paula <tiagodepalves@gmail.com>

// LOAD X
PUSH0
CALLDATALOAD

// ------ 2 ** (log2(x) / 2) ------

// If value has upper 128 bits set, log2 result is at least 128
DUP1                                       // stack: x x
PUSH16 0xffffffffffffffffffffffffffffffff  // stack: 2**128-1 x x
LT                                         // stack: x>=2**128 x
PUSH1 7                                    // stack: 7 x>=2**128 x
SHL                                        // stack: ±log2(x) x

// If upper 64 bits of 128-bit half set, add 64 to result
// stack: log2(x) 1 x
DUP2                     // stack: x ±log2(x) x
DUP2                     // stack: ±log2(x) x ±log2(x) x
SHR                      // stack: x/±log2(x) ±log2(x) x
PUSH8 0xffffffffffffffff // stack: 2**64-1 x/±log2(x) ±log2(x) x
LT                       // stack: (x/2**±log2(x) >= 2**64) ±log2(x) x
PUSH1 6                  // stack: 6 (x/2**±log2(x) >= 2**64) ±log2(x) x
SHL                      // stack: ±log2(x) ±log2(x) x
OR                       // stack: ±log2(x) x

// If upper 32 bits of 64-bit half set, add 32 to result
DUP2
DUP2
SHR
PUSH4 0xffffffff
LT
PUSH1 5
SHL
OR

// If upper 16 bits of 32-bit half set, add 16 to result
DUP2
DUP2
SHR
PUSH2 0xffff
LT
PUSH1 4
SHL
OR

// If upper 8 bits of 16-bit half set, add 8 to result
DUP2
DUP2
SHR
PUSH1 0xff
LT
PUSH1 3
SHL
OR // stack: ±log2(x) x

// Table lookup
PUSH32 0x1b3647545f69737b838b9299a0a6acb2b7bdc2c8cdd2d6dbe0e4e9edf1f6fafe
// stack: table ±log2(x) x
DUP3  // stack: x table ±log2(x) x
DUP3  // stack: ±log2(x) x table ±log2(x) x
SHR   // stack: 8*i table ±log2(x) x
PUSH1 3
SHR   // stack: i table ±log2(x) x
BYTE  // stack: table[i] ±log2(x) x

// Initial guess: ±sqrt(x) := table[i] * 2**(±log2(x)/2) / 16
SWAP1 // stack: ±log2(x) table[i] x
PUSH1 1
SHR   // stack: ±log2(x)/2 table[i] x
SHL   // stack: 16*±sqrt(x) x
PUSH1 4
SHR   // stack: ±sqrt(x) x
// ---------------------------------

// Newton's method
// r = (r + x/r) / 2
DUP1
DUP3
DIV
ADD
PUSH1 1
SHR

DUP1
DUP3
DIV
ADD
PUSH1 1
SHR

DUP1
DUP3
DIV
ADD
PUSH1 1
SHR

DUP1
DUP3
DIV
ADD
PUSH1 1
SHR

DUP1
DUP3
DIV
ADD
PUSH1 1
SHR

DUP1
DUP3
DIV
ADD
PUSH1 1
SHR

// y = x/r
DUP1
SWAP2
DIV

// min(r, x/r)
DUP2
GT
SWAP1
SUB

// Return result
PUSH0
MSTORE
PUSH1 32
PUSH0
RETURN

Sqrt.sol

// SPDX-License-Identifier: GPL-2.0-only
// Gas efficient SQRT method, computes floor(sqrt(x)).
// Constant gas cost of 339.
//
// Authors: Lohann Ferreira, Tiago de Paula
pragma solidity >=0.7.0 <0.9.0;

function sqrt(uint256 x) pure returns (uint256 r) {
    assembly ("memory-safe") {
        // r ≈ log2(x)
        r := shl(7, lt(0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF, x))
        r := or(r, shl(6, lt(0xFFFFFFFFFFFFFFFF, shr(r, x))))
        r := or(r, shl(5, lt(0xFFFFFFFF, shr(r, x))))
        r := or(r, shl(4, lt(0xFFFF, shr(r, x))))
        r := or(r, shl(3, lt(0xFF, shr(r, x))))

        let lut := 0x1B3647545F69737B838B9299A0A6ACB2B7BDC2C8CDD2D6DBE0E4E9EDF1F6FAFE
        let i := shr(3, shr(r, x))
        r := shr(4, shl(shr(1, r), byte(i, lut)))

        // Newton's method
        r := shr(1, add(r, div(x, r)))
        r := shr(1, add(r, div(x, r)))
        r := shr(1, add(r, div(x, r)))
        r := shr(1, add(r, div(x, r)))
        r := shr(1, add(r, div(x, r)))
        r := shr(1, add(r, div(x, r)))

        // r = min(r, x/r)
        r := sub(r, gt(r, div(x, r)))
    }
}