分享实现 L^β(0<=β<1)

dwuw · 2025年02月22日 · 最后由 shooter 回复于 2025年02月24日 · 56 次阅读

1.原理

目标实现 L^β的计算，其中 0 <=β <1 这里将 L 和β都换成 Q64.64 的形式来计算

$L^β = 2^{β * log_2(L)}$
(用 2 为底会方便很多)

2.计算$log_2(L)$

我们要先计算 log2(L) 因为任何一个数字都可以换成 $L = 2^k * (y)$ 所以有: $$ log_2(x) = log_2(2^k * (y)) = log_2(2^k) + log_2(y) = k*log_2(2) + log_2(y) =k+log_2(y) $$

其中 k 是 msb (Most Significant Bit , 如果最高 bit 数量 - leading_zero 得到 , 如:u128 就是 127-leading_zero)

let leading_zeros = x.leading_zeros();
let scale_msb = 127 - leading_zeros; // 这才是正确的 MSB 位置
let msb = scale_msb - 64;

eg:
假设输入一个数 x = 300
 300 = 256 + 44
     = 2^8 + 44
     = 2^8 * (1 + 44/256)
     = 2^8 * (1 + 0.171875)
因此:
log_2(300) = log_2(2^8 * (1 + 0.171875))
        = 8*log_2(2) + log_2(1 + 0.171875)
        =8+log(1.171875)

同理 如果 x =100
100 = 64 + 36
    = 2^6 ( 1 + 36/64)
    = 2^6( 1 + 0.5625)
因此 
log_2(100) = 6log_22 + log_2(1+0.5625) = 1+ log_2(1.5625)

同理 如果 x =100.5
可以分为两部分100
100 = 64 + 36.5
    = 2^6 ( 1 + (36+0.5)/64)
    = 2^6( 1 + 0.57)
因此 
log_2(100) = 6log_2 2 + log_2(1+0.57)
            = 6 + log2_(1.57)

当然如果 y=1 那么$log_2 y$ =0，也就是说 L 正好就是 2 的幂

if y == 1 << 64 {
    return Ok((scale_msb as u128) << 64);
}

所以接下来我们主要算的就是另一部分 $log_2(y)$ y 该如何得到呢

$x = 2^k * y => y = x / 2^k => y = x * 2^-k$

用10这个例子来看
十进制:10 = 二进制:1010  此时的msb = 3 即k=3 
而 *2^-k 其实就是 >>k
所以 y = 二进制:1.010 = 十进制: 1.125

所以用代码来说就是

let mut res = (msb as u128) << 64;
let y = x >> msb;

其中 $1<= y <2$ , 这部分的原理是 : 如果 $y ∈ [1,2)$，那么 $log₂(y) ∈ [0,1)$

我们要找到这个 $[0,1)$ 范围内的值的二进制表示： $log_2(y) = 0.b₁b₂b₃...$ (如果用二进制表示，那么这里$bn$ = 0 or 1

又因为 : $log₂(y²) = 2log₂(y)$

如果 $y² ≥ 2$ -> $log_2(y^2) >1$ -> $2log₂(y) ≥ 1$ -> $log₂(y) ≥ 0.5$ 所以 b₁ = 1
如果 $y² < 2$ -> $log_2(y^2) <1$ -> $2log₂(y) < 1$ -> $log₂(y) < 0.5$ 所以 b₁ = 0

如果 b1=1

当 b₁ = 1： $log₂(y) = 0.1b₂b₃...$ $= 0.5 + 0.b₂b₃...$

接下来，令 $y_{new}=y^2/2$ , 则有 $log₂y_{new} = log₂(y²/2) = 2log₂y - 1$

$将 log₂(y) = 0.5 + 0.b₂b₃...$ 代入等式之后就可以得到

$log₂y_{new}$ $= log₂(y²/2) = 2log₂y - 1$ $=2(0.5+0.b_2b_3...b_n) - 1$ $= 1 + 2(0.b_2b_3...b_n) -1$ $= 2* 0.b₂b₃...$

接着和上面一样

如果 $y_{new}² ≥ 2$ -> $log_2(y_{new}^2) >1$ -> $2log₂(y_{new}) ≥ 1$ -> $log₂(y_{new}) ≥ 0.5$ -> $2* 0.b₂b₃... \geq 0.5$ -> $0.b_2b_3.. \geq 0.25$ -> $b_2 =1$
如果 $y_{new}² < 2$ -> $log_2(y_{new}^2) <1$ -> $2log₂(y_{new}) < 1$ -> $log₂(y_{new}) < 0.5$ -> $2* 0.b₂b₃... <0.5$ -> $0.b_2b_3.. < 0.25$ -> $b_2 =0$ ... 接下来同理

用一个例子来看，第一次迭代

第一轮：
y = 1.5
y² = 2.25 > 2
log₂(1.5) > 0.5
b₁ = 1
res += 0.5
y_new = 2.25/2 = 1.125

第二轮
y = 1.125
y² = 1.265625 < 2
log₂(1.125) < 0.5
b₂ = 0
y_new = 1.265625  (因为 < 2，不需要除2)

第三轮
y = 1.265625
y² = 1.601806640625 < 2
log₂(1.265625) < 0.5
b₃ = 0
y_new = 1.601806640625

第四轮
y = 1.601806640625
y² = 2.5657... > 2
log₂(1.601806640625) > 0.5
b₄ = 1
res += 0.0625 (2⁻⁴)
y_new = 2.5657.../2 ≈ 1.28285...

总的来说，每一轮的操作

1.计算 y²
2.判断是否 ≥ 2
    如果 ≥ 2：当前位为1，加上对应权重，y = y²/2
    如果 < 2：当前位为0，y = y²
3.继续下一轮，权重位退一位

let double_unit = TWO_X64; // 2 * 2^64
let mut delta = HALF_X64; // 0.5 * 2^64
let mut y = y; // 从前面传入的 y

for _ in 0..64 {
    // y = (y * y) / ONE_X64
    y = fixed_mul(y, y).unwrap();

    // 检查 y 是否大于等于 2 * 2^64
    if y >= double_unit {
        // 添加 2^{-m} 因子到结果中
        res = res + delta;

        // 将 y 除以 2
        y >>= 1;
    }

    // 每次迭代将 delta 除以 2
    delta >>= 1;
}

为什么新的 y 对应 b2

第一轮后
y_new = 1.125
log₂(1.125) = log₂(2.25/2)
             = log₂(2.25) - 1
             = 2log₂(1.5) - 1
             = 2(0.5 + 0.b₂b₃...) - 1
             = 2*0.b₂b₃...

判断第二位
y_new² = 1.265625 < 2
意味着 log₂(y_new) < 0.5

而 log₂(y_new) = 2*0.b₂b₃...
所以 2*0.b₂b₃... < 0.5
因此 0.b₂b₃... < 0.25
所以 b₂ = 0

2.计算 $β * log_2L$

这很简单直接通过定点数相乘中间用 u256 过度一下即可

// 定点数乘法 (Q64 * Q64 => Q64)
pub fn fixed_mul_x64(a: u128, b: u128) -> Option<u128> {
    let a = U256::from(a);
    let b = U256::from(b);
    let result = (a * b) >> 64;
    result.try_into().ok()
}

3.计算$2^x$

这就是最后一步因为 $L^β = 2^{β * log_2(L)}$

而 $β * log_2L$ 我们已经算出来了接下来就是算 $2^x$

计算 $2^x$，其中 x 可以表示为： $x = n + f$，n 是整数部分，f 是小数部分所以就有 : $2^x = 2^n * 2^f$ (2^f 可以用二进制分解法计算)

对于 f ,他是一个小数 ,所以就有 $f = 0.b₁b₂b₃...$ $= b₁/2 + b₂/4 + b₃/8 + ...$ $= b₁*2⁻¹ + b₂*2⁻² + b₃*2⁻³ + ...$

所以 $2^f = 2^(b₁*2⁻¹ + b₂*2⁻² + b₃*2⁻³ + ...)$ 又因为指数的分配律： $2^{a+b} = 2^a 2^b$ 所以 $2^f = 2^{b₁*2⁻¹} * 2^{b₂*2⁻²} * 2^{b₃*2⁻³} * ...$

所以可以得到 $2^x = 2 ^ n * 2^{b₁2⁻¹} * 2^{b₂*2⁻²} * 2^{b₃*2⁻³ ...}$
其中可以通过判断 f 的二进制权重来判断 $b_1,b_2,b_3 ...$ 是否为 0 当$b_1 , b_2 , b_3 ...$ 为 0 的时候其实就是 $ 1$ 而我们可以提前算出 $2^{2^{-1}}<<64 , 2^{2^{-2}} <<64 , 2^{2^{-3}} <<64 ...$的值 , 同时因为我们是 Q64.64，因此 f 最多最多之后 64 为小数，所以可以提前将这 64 个值算出，减少 cu

所以我们先判断 $2^ f$ 的部分

这里的核心就是通过 bit 判断 x 的 64-0 位是否有，如果有 1，就直接再 result 上乘上我们与计算的值即可

如这里 x & 0xFF00000000000000

十六进制0xFF00000000000000
二进制  0000 1111 1111 0000 0000
        ---- 63 --  56  -------0

这里的 0xFF00000000000000 就是 63-56 位置都是 1 的二进制如果 x & 0xFF00000000000000 > 0 那么就说明 x 在 63-56 位上至少有一个 1，就可以进入判断

而在 if 的内部就是对 63-56 每一位的单独判断如果有 1 那么就说明这里的$b_i$ 为 1，直接将 res 乘以我们预计算的结果即可 : 如 0x16A09E667F3BCC909 就是 2^(2^-1)

// 从1 <<64 开始
   let mut result: u128 = ONE_X64;

 // 先检查第63位到第56位（共8位）中是否有任何一位是1
   //x & 0xFF00000000000000 > 0 意味着第63-56位 至少有一个1
   //0x16A09E667F3BCC909 是 2^(2^-1) <<64
   //0x1306FE0A31B7152DF 是 2^(2^-2) <<64
   //0x1172B83C7D517ADCE 是 2^(2^-3) <<64
   //0x10B5586CF9890F62A 是 2^(2^-4) <<64
   //0x1059B0D31585743AE 是 2^(2^-5) <<64
   //0x102C9A3E778060EE7 是 2^(2^-6) <<64
   //0x10163DA9FB33356D8 是 2^(2^-7) <<64
   //0x100B1AFA5ABCBED61 是 2^(2^-8) <<64
   if x & 0xFF00000000000000 > 0 {
       if x & 0x8000000000000000 > 0 {
           result = fixed_mul_x64(result, 0x16A09E667F3BCC909).unwrap();
       }
       if x & 0x4000000000000000 > 0 {
           result = fixed_mul_x64(result, 0x1306FE0A31B7152DF).unwrap();
       }
       if x & 0x2000000000000000 > 0 {
           result = fixed_mul_x64(result, 0x1172B83C7D517ADCE).unwrap();
       }
       if x & 0x1000000000000000 > 0 {
           result = fixed_mul_x64(result, 0x10B5586CF9890F62A).unwrap();
       }
       if x & 0x800000000000000 > 0 {
           result = fixed_mul_x64(result, 0x1059B0D31585743AE).unwrap();
       }
       if x & 0x400000000000000 > 0 {
           result = fixed_mul_x64(result, 0x102C9A3E778060EE7).unwrap();
       }
       if x & 0x200000000000000 > 0 {
           result = fixed_mul_x64(result, 0x10163DA9FB33356D8).unwrap();
       }
       if x & 0x100000000000000 > 0 {
           result = fixed_mul_x64(result, 0x100B1AFA5ABCBED61).unwrap();
       }
   }

   // 检查56-48位
   if x & 0xFF000000000000 > 0 {
       if x & 0x80000000000000 > 0 {
           result = fixed_mul_x64(result, 0x10058C86DA1C09EA2).unwrap();
       }
       if x & 0x40000000000000 > 0 {
           result = fixed_mul_x64(result, 0x1002C605E2E8CEC50).unwrap();
       }
       if x & 0x20000000000000 > 0 {
           result = fixed_mul_x64(result, 0x100162F3904051FA1).unwrap();
       }
       if x & 0x10000000000000 > 0 {
           result = fixed_mul_x64(result, 0x1000B175EFFDC76BA).unwrap();
       }
       if x & 0x8000000000000 > 0 {
           result = fixed_mul_x64(result, 0x100058BA01FB9F96D).unwrap();
       }
       if x & 0x4000000000000 > 0 {
           result = fixed_mul_x64(result, 0x10002C5CC37DA9492).unwrap();
       }
       if x & 0x2000000000000 > 0 {
           result = fixed_mul_x64(result, 0x1000162E525EE0547).unwrap();
       }
       if x & 0x1000000000000 > 0 {
           result = fixed_mul_x64(result, 0x10000B17255775C04).unwrap();
       }
   }

   // 检查48-40位
   if x & 0xFF0000000000 > 0 {
       if x & 0x800000000000 > 0 {
           result = fixed_mul_x64(result, 0x1000058B91B5BC9AE).unwrap();
       }
       if x & 0x400000000000 > 0 {
           result = fixed_mul_x64(result, 0x100002C5C89D5EC6D).unwrap();
       }
       if x & 0x200000000000 > 0 {
           result = fixed_mul_x64(result, 0x10000162E43F4F831).unwrap();
       }
       if x & 0x100000000000 > 0 {
           result = fixed_mul_x64(result, 0x100000B1721BCFC9A).unwrap();
       }
       if x & 0x80000000000 > 0 {
           result = fixed_mul_x64(result, 0x10000058B90CF1E6E).unwrap();
       }
       if x & 0x40000000000 > 0 {
           result = fixed_mul_x64(result, 0x1000002C5C863B73F).unwrap();
       }
       if x & 0x20000000000 > 0 {
           result = fixed_mul_x64(result, 0x100000162E430E5A2).unwrap();
       }
       if x & 0x10000000000 > 0 {
           result = fixed_mul_x64(result, 0x1000000B172183551).unwrap();
       }
   }

   // 检查40-32位
   if x & 0xFF00000000 > 0 {
       if x & 0x8000000000 > 0 {
           result = fixed_mul_x64(result, 0x100000058B90C0B49).unwrap();
       }
       if x & 0x4000000000 > 0 {
           result = fixed_mul_x64(result, 0x10000002C5C8601CC).unwrap();
       }
       if x & 0x2000000000 > 0 {
           result = fixed_mul_x64(result, 0x1000000162E42FFF0).unwrap();
       }
       if x & 0x1000000000 > 0 {
           result = fixed_mul_x64(result, 0x10000000B17217FBB).unwrap();
       }
       if x & 0x800000000 > 0 {
           result = fixed_mul_x64(result, 0x1000000058B90BFCE).unwrap();
       }
       if x & 0x400000000 > 0 {
           result = fixed_mul_x64(result, 0x100000002C5C85FE3).unwrap();
       }
       if x & 0x200000000 > 0 {
           result = fixed_mul_x64(result, 0x10000000162E42FF1).unwrap();
       }
       if x & 0x100000000 > 0 {
           result = fixed_mul_x64(result, 0x100000000B17217F8).unwrap();
       }
   }

   // 检查32-24位
   if x & 0xFF000000 > 0 {
       if x & 0x80000000 > 0 {
           result = fixed_mul_x64(result, 0x10000000058B90BFC).unwrap();
       }
       if x & 0x40000000 > 0 {
           result = fixed_mul_x64(result, 0x1000000002C5C85FE).unwrap();
       }
       if x & 0x20000000 > 0 {
           result = fixed_mul_x64(result, 0x100000000162E42FF).unwrap();
       }
       if x & 0x10000000 > 0 {
           result = fixed_mul_x64(result, 0x1000000000B17217F).unwrap();
       }
       if x & 0x8000000 > 0 {
           result = fixed_mul_x64(result, 0x100000000058B90C0).unwrap();
       }
       if x & 0x4000000 > 0 {
           result = fixed_mul_x64(result, 0x10000000002C5C860).unwrap();
       }
       if x & 0x2000000 > 0 {
           result = fixed_mul_x64(result, 0x1000000000162E430).unwrap();
       }
       if x & 0x1000000 > 0 {
           result = fixed_mul_x64(result, 0x10000000000B17218).unwrap();
       }
   }

   // 检查24-16位
   if x & 0xFF0000 > 0 {
       if x & 0x800000 > 0 {
           result = fixed_mul_x64(result, 0x1000000000058B90C).unwrap();
       }
       if x & 0x400000 > 0 {
           result = fixed_mul_x64(result, 0x100000000002C5C86).unwrap();
       }
       if x & 0x200000 > 0 {
           result = fixed_mul_x64(result, 0x10000000000162E43).unwrap();
       }
       if x & 0x100000 > 0 {
           result = fixed_mul_x64(result, 0x100000000000B1721).unwrap();
       }
       if x & 0x80000 > 0 {
           result = fixed_mul_x64(result, 0x10000000000058B91).unwrap();
       }
       if x & 0x40000 > 0 {
           result = fixed_mul_x64(result, 0x1000000000002C5C8).unwrap();
       }
       if x & 0x20000 > 0 {
           result = fixed_mul_x64(result, 0x100000000000162E4).unwrap();
       }
       if x & 0x10000 > 0 {
           result = fixed_mul_x64(result, 0x1000000000000B172).unwrap();
       }
   }

   // 检查16-8位
   if x & 0xFF00 > 0 {
       if x & 0x8000 > 0 {
           result = fixed_mul_x64(result, 0x100000000000058B9).unwrap();
       }
       if x & 0x4000 > 0 {
           result = fixed_mul_x64(result, 0x10000000000002C5D).unwrap();
       }
       if x & 0x2000 > 0 {
           result = fixed_mul_x64(result, 0x1000000000000162E).unwrap();
       }
       if x & 0x1000 > 0 {
           result = fixed_mul_x64(result, 0x10000000000000B17).unwrap();
       }
       if x & 0x800 > 0 {
           result = fixed_mul_x64(result, 0x1000000000000058C).unwrap();
       }
       if x & 0x400 > 0 {
           result = fixed_mul_x64(result, 0x100000000000002C6).unwrap();
       }
       if x & 0x200 > 0 {
           result = fixed_mul_x64(result, 0x10000000000000163).unwrap();
       }
       if x & 0x100 > 0 {
           result = fixed_mul_x64(result, 0x100000000000000B1).unwrap();
       }
   }

   // 检查8-0位
   if x & 0xFF > 0 {
       if x & 0x80 > 0 {
           result = fixed_mul_x64(result, 0x10000000000000059).unwrap();
       }
       if x & 0x40 > 0 {
           result = fixed_mul_x64(result, 0x1000000000000002C).unwrap();
       }
       if x & 0x20 > 0 {
           result = fixed_mul_x64(result, 0x10000000000000016).unwrap();
       }
       if x & 0x10 > 0 {
           result = fixed_mul_x64(result, 0x1000000000000000B).unwrap();
       }
       if x & 0x8 > 0 {
           result = fixed_mul_x64(result, 0x10000000000000006).unwrap();
       }
       if x & 0x4 > 0 {
           result = fixed_mul_x64(result, 0x10000000000000003).unwrap();
       }
       if x & 0x2 > 0 {
           result = fixed_mul_x64(result, 0x10000000000000001).unwrap();
       }
       if x & 0x1 > 0 {
           result = fixed_mul_x64(result, 0x10000000000000001).unwrap();
       }
   }

当小数判断结束之后就可以判断整数，判断整数 n 非常简单直接>>64 之后就可以得到 n 然后将 res 乘以 2^n 即可

let integer_part = x >> 64;
result <<= integer_part;

至此我们就完成了 pow 的计算结果也是 Q64.64 的形式

参考：

https://github.com/raydium-io/raydium-clmm/blob/master/programs/amm/src/libraries/tick_math.rs

https://github.com/PaulRBerg/prb-math

1 个赞

shooter #1 2025年02月24日

需要登录后方可回复, 如果你还没有账号请注册新账号

分享 实现 L^β(0<=β<1)

分享 实现 L^β(0<=β<1)

1.原理

2.计算$log_2(L)$

2.计算 $β * log_2L$

3.计算$2^x$

分享实现 L^β(0<=β<1)

分享实现 L^β(0<=β<1)