# 2025q1 一對一討論
contributed by < `thwuedwin` >
[記錄:一對一討論](https://docs.google.com/document/d/1BTUCBCq4-JX5VDAts3B1pLymftVmmnRBlTvKf1pB2Us/edit?tab=t.0)
## Q1: Assume x is a single precision floating point (32bits or more) and always >= 1. Calculate (int) sqrtf(x) .
### 數值範圍分析
輸入:所有正浮點數 `f`
輸出:`f` 的平方根,整數 `x`,其中 $x^2\leq f< (x+1)^2$
當 $0\leq f<1$ 時,`x` 必為 $0$。
考慮 `f` 的最大值為 $1.11...1_{(2)}\times2^{127}$,$\text{f_max}<2^{64}$,可以用 64 位元整數儲存,但在各種實作內都會有超出範圍的情況。
這邊列出 3 種演算法
- digit-by-digit
- 牛頓法
- 二分搜尋
在 digit-by-digit 做法中,會使用變數 `y` 來記錄答案,在初始情況下,`y` 的位元數會和 `f` 相同,若使用大數來儲存,就可以非常簡單的實作。不使用大數的情況下,因為 `f` 的精度有限,在有效位數外會有很多 0,在數值不大時,可以只儲存有效位數,另外處理指數的部分。但在 `f` 很大,大約 $2^{122}$ 以上時,在儲存有效位數後,沒有辦法留下足夠的空間來繼續運算。
牛頓法需要進行除法,但浮點數除以整數的誤差非常大,只要精度不足就會導致錯誤的結果。我認為還是有必要將其轉換成整數在進行運算,這樣無可避免仍需要使用大數。
二分搜尋法也需要進行平方的比較,這部分也會超出 64 位元的範圍。
我的實作中用兩個 64 位元整數來表示 128 位元整數,搭配自定義的加減法、位移運算、比較大小等函式來使用。
```c
typedef struct big_num {
uint64_t high, low;
} bn_t;
```
### 128-bits digit-by-digit 實作
這部分是我目前的實作方法,可以處理所有正浮點數輸入。對於負數、正負無窮大、NaN會直接報錯。
32 位元的浮點數最大值是大約是 $2^{127}$,可以用 64 位元無號整數來儲存。本文中若沒特別說明,則浮點數是指 IEEE 754 定義的 32 位元浮點數。
- 不使用除法
- 除了輸入的參數外,不使用浮點數
- 有些平台沒有支援 GCC [128 bits 整數](https://gcc.gnu.org/onlinedocs/gcc/_005f_005fint128.html),所以我沒有使用。
用 digit-by-digit 作法,將主要計算平方根的變數改為自定義的大數可解決溢位或精度不足的問題。
以三元平方和為例
$$
(a+b+c)^2=a^2+(2a+b)b+(2(a+b)+c)c
$$
$m$ 元平方和可以拆解成 $1$ 到 $m-1$ 元平方和,再加上一項
$$
Y_m=(2P_{m+1}+a_m)a_m
$$
其中 $P_m=a_n+a_{n-1}+...+a_{m+1}$。
實作上我們會先選擇一個足夠大且是 2 的冪,作為初始值,以三元平方和為例就是 $a^2$,接著逐個位元檢查「加入這個位元是否會超過」,不會超過代表應該把該位元設成 1。
```c
uint64_t my_sqrtf_bn(float f)
{
union {
float f;
uint32_t u;
} u;
u.f = f;
uint32_t sign = u.u & 0x80000000;
uint32_t exp = (u.u >> 23) & 0xFF;
uint32_t man = (u.u & 0x007FFFFF) | 0x00800000 ;
assert(!sign);
assert(exp != 255);
if (exp < 127)
return 0;
exp -= 127;
bn_t x;
bn_init(&x);
if (exp <= 23) {
man >>= 23 - exp;
x.low = man;
} else {
x.low = man;
bn_lshift(&x, exp - 23);
}
int shift = exp & ~1;
bn_t m, y;
bn_init(&m);
m.low = 1ULL;
bn_init(&y);
bn_lshift(&m, shift);
while (m.high || m.low) {
bn_t b;
bn_assign(&b, &y);
bn_add(&b, &m);
bn_rshift(&y, 1);
if (bn_geq(&x, &b)) {
bn_sub(&x, &b);
bn_add(&y, &m);
}
bn_rshift(&m, 2);
}
assert(!y.high);
return y.low;
}
```
下圖是 `f` 對 `my_sqrtf_bn` 作圖的結果

### 測試
#### 效能
我目前的實作 `my_sqrtf_bn` 有使用大數,可以預期效能會比較差。我將目前的實作和沒有使用大數,正確但沒辦法處理所有輸入的 `my_sqrtf` 進行比較。我使用 `clock_gettime` 來記錄執行時間,結果如下

可以看到使用大數所需的時間大約是不使用大數的兩倍,或許可以在較小的數值不使用大數。
將很大的離群值(outliers)移除後可以看出執行時間是隨著數值變大而增加,符合預期

讓我~~比較~~不解的是,為何在值很小的時候反而經常出現執行時間特別長的情況,這些離群值要如何解釋。但我目前沒有想到要如何驗證。
#### 正確性
因為已經實作大數,測試就變得非常簡單,僅需檢查 $x^2\leq f < (x+1)^2$。經過測試,我的實作是正確的。
### 實作 & 測試過程紀錄
這部分是我實作過程思路和發現的紀錄
#### 第一版想法
實作上參考了過去 `sqrti` 的實作向下取整數,所以在數值小於 1 時,開平方根後會是 0,因此沒有特別考慮 subnormal number。
IEEE 754 會將浮點數以科學記號表示,可以表示為 $v=(-1)^S\times m \times 2^e$,實作中先將浮點數的三個部分分別取出。因為小數部分的精度有限,以 32 位元為例,包含 implicit leading significand 在內有 24 位元,儲存在 `man` 中,現在 man 的資料代表的意思是 $b_0.b_1b_2...b_23$,我想將它小數點移到最右邊。所以若指數小於 23,第 17 行的部分,我會把指數乘到小數部分,再將乘完後小數點以後的部分省略。若指數大於 23,則只把 $2^{23}$ 乘到小數部分,剩下的部分可以直接除以 2 來計算,只要注意指數的奇偶即可。
```c
uint64_t my_sqrtf(float f)
{
union {
float f;
uint32_t u;
} u;
u.f = f;
uint32_t sign = u.u & 0x80000000;
uint32_t exp = (u.u >> 23) & 0xFF;
uint32_t man = u.u & 0x007FFFFF;
man |= (1 << 23); // add the implicit leading bit
assert(!sign);
if (exp < 127)
return 0;
exp -= 127; // unbiased exponent
if (exp < 23) {
man >>= (23 - exp);
exp = 0;
} else {
exp -= 23;
if (exp & 1) {
exp &= ~1;
man <<= 1;
}
}
uint64_t result = sqrti(man);
result <<= (exp >> 1);
}
```
簡單測試後發現這個做法在 `f` 很大是錯的。若對 $v=m \times 2^e$ 開平方根,會變成 $\sqrt v=\sqrt m\times 2^{e/2}$,假設 $e$ 是偶數,若不是的話塞到 $m$ 裡面即可。接著我會先計算 $\sqrt m$ 並將其捨去到整數,再乘上 $2^{e/2}$,當 $e$ 大時,先捨棄就會造成很大的誤差。所以不能事先捨棄,必須計算到所有位元的貢獻。
#### 第二版想法
必須計算到所有位元的貢獻不等於需要把所有位元都儲存進變數內。科學記號有個很好的性質 $1.234\times10^{100}$,可以被改寫為 $1234\times10^{97}$。我也可以用同樣的概念來改寫原先 IEEE 754 的格式,分別儲存有效位數和指數。有效位數後方接著的 0 都不放進變數裡,而是另外紀錄後面應該接著幾個 0,並使用整數來操作。這樣做的好處是有效位數是一個整數,就可以直接採用 `sqrti` 中 digit-by-digit 的方式計算平方根。所以有效位數後方接著的 0 都不放進變數裡,而是另外紀錄後面應該接著幾個 0。這樣的做法可以節省空間,但在位元運算時需要注意對齊,才能得到正確的結果。
```c=
uint64_t my_sqrtf(float f)
{
union {
float f;
uint32_t u;
} u;
u.f = f;
uint32_t sign = u.u & 0x80000000;
uint32_t exp = (u.u >> 23) & 0xFF;
uint64_t man = u.u & 0x007FFFFF;
man |= (1 << 23); // add the implicit leading bit
assert(!sign);
assert(exp != 255);
if (exp < 127)
return 0;
exp -= 127; // unbiased exponent
// exp is the total bits
uint64_t m;
int total_bits = exp + 1;
int hidden_bits;
int total_shift = (total_bits - 1) & ~1;
if (total_bits <= 24) {
man >>= (24 - total_bits);
m = 1ULL << total_shift;
hidden_bits = 0;
} else {
hidden_bits = total_bits - 24;
if (hidden_bits & 1) {
man <<= 1;
hidden_bits -= 1;
}
int shift = (total_bits - hidden_bits - 1) & ~1;
m = 1ULL << shift;
total_shift -= shift;
}
uint64_t y = 0;
while (m || hidden_bits) {
uint64_t b = y + m;
y >>= 1;
if (man >= b) {
man -= b;
y += m;
}
m >>= 2;
if (m == 1 && hidden_bits) {
int lz = (clz64(man) < clz64(y)? clz64(man): clz64(y));
int shift;
if (lz >= hidden_bits) {
shift = hidden_bits;
m <<= shift;
man <<= shift;
y <<= shift;
} else {
shift = lz & ~1;
m <<= shift;
man <<= shift;
y <<= shift;
}
hidden_bits -= shift;
}
}
return y;
}
```
24-37 行在計算有效位數,因為小數點後的值並不影響開平方根的結果,所以我將小數點後的部分捨去。 39-62 行是 digit-by-digit 計算的主迴圈。47 行開始是當 `m` 變成 1 時,計算無法繼續,若還有隱藏位元的話,就將其放入有效位數中,繼續計算。在加入有效位數時,會先使用 `clz` 來計算有多少空間可以釋放,避免 overflow 造成錯誤。
但這個做法還是有一個致命的問題。在輸入的 `f` 過大時,具體來說是 $2^{122}$ 左右開始,leading zero 的空間不足造成無法新增有效位數。但這個做法在較小的輸入都是正確的。
#### 測試過程紀錄
我一開始是窮舉,列出很多組值,從其中抽幾個點用計算機計算,小至個位數,大至 $10^{36}$ 都抽幾個數字檢查正確。下圖為 $f$ 對 $\sqrt f$ 作圖。

但用眼睛看不夠嚴謹,需要進一步測試。
測試會遇到很大的困難點在於浮點數的誤差。若是整數的平方根,一般會用 `x*x<=n && n<(x+1)*(x+1)` 進行測試。
我用窮舉的方式,發現在帶入的數值 `f=16785408` 會產生錯誤。這個值輸出的平方根為 $4096$,實際值是 $4096.99987...$,其精度超過浮點數的表示範圍,實際上的數值會被轉換為 $4097$,這和 IEEE 754 的捨入模式有關。經測試,將上界的判斷式改為 `n<=(x+1)*(x+1)` 可成功通過測試。
IEEE 754 的捨入模式為 Round Ties To Even,因為 $4096.99987...$ 較接近 $4097$,而不是小於 $4097$ 的最大浮點數,所以實際儲存值為 $4097$
> roundTiesToEven, the floating-point number nearest to the infinitely precise result shall be
delivered; if the two nearest floating-point numbers bracketing an unrepresentable infinitely
precise result are equally near, the one with an even least significant digit shall be delivered;
接著遇到的錯誤產生在 `f=18446744073709551616`,這個值是 $2^{64}$,平方根為 $2^{32}$,很明顯是相乘之後溢位了。
在這之後我嘗試改成將 `x` 轉型成浮點數,但測試錯誤產生的點提早了,我認為是浮點數誤差造成的,這時我沒有深究。~~我決定將小於 $2^{32}$ 的數據用前面的方法測試,超過的部分用科學記號的方式來測試。再將整數轉換為科學記號時,~~
我決定用整數表示科學記號,以此觀察浮點數的誤差來源。我認為 IEEE 754 的規範中,有效位數只有 24 位元,所以我也只需要保留 24 位元。實際測試之後發現,在數值很大像是 $2^{32}$ 時,僅僅數個位元的誤差就非常可觀,當我捨棄第 25 位元時,實際捨棄的值為 $2^7$,平方後會是 $2^{14}$,在考慮進位後會是一個不容忽視的值。我會用以下這個算是來判斷是否小於 $(x+1)^2$,但誤差會導致 $f-x^2$ 數量級遠大於 $x$。
$$
(x+1)^2\>f\leq x^2
\Rightarrow f \leq x^2\ \text{and}\ f - x^2>2x+1
$$
舉個更實際的例子,$2^{80}$ 的平方根為 $2^{40}$。若只保留有效位數 24 位元,平方後有 48 位元有效位數,我們計算出來的結果在 $2^{40}$ 到 $2^{41}-1$ 都會被判斷為合法的解答,這顯然是錯誤的。
### 數學推導
#### digit-by-digit 開平方根
在[作業二](https://hackmd.io/2pHlZxn7SnuPh3gIhH-Z3g?view#程式碼運作原理分析17)有整理過 digit-by-digit 的做法的推論。
#### 128-bits 運算證明
##### 位移運算
根據 C 語言規格書 6.5.7 第 3 段,位移量 `shift` 應介於 0 到 $W - 1$,其中 $W$ 為該變數的位元數。
> If the value of the right operand is negative or is greater than or equal to the width of the promoted left operand, the behavior is undefined.
若 `shift` 大於或等於 128 時,我將所有位元都設成 $0$,這是未定義行為所以我可以自行決定,且在我的實作中有確保不會超過。
在這次討論中,最右邊的位元(MSB)為第 $127$ 位元,最左邊的位元(LSB)為第 $0$ 位元。
在其他情況下讓 `n` 左移 `shift` 位元,將 `n` 的 $0$ 到 $127 - \text{shift}$ 位元,搬移至 $\text{shift}$ 到 $127$ 位元,並將 $0$ 到 $\text{shift}$ 位元補 $0$。
若 `shift` 在 64 到 127 之間,保留的位元會小於 64 位元,於是全部來自 `low`,保留位元為 $\text{shift}-64$,直接存入 `high`,並將 `low` 設為 $0$。
若 `shift` 在 $0$ 到 $63$ 之間,可直接將 `high` 左移,超出的範圍本身就會被捨棄。但直接將 `low` 左移則會導致錯誤的捨棄,所以要將 `low` 的 $\text{shift}$ 到 $63$ 位元先存入 `high`,在進行左移。
右移運算和左移基本上相同。
```c
void bn_lshift(bn_t *n, int shift)
{
if (shift == 0) return;
if (shift >= 128) {
n->high = 0;
n->low = 0;
} else if (shift >= 64) {
n->high = n->low << (shift - 64);
n->low = 0;
} else {
n->high <<= shift;
n->high |= n->low >> (64 - shift);
n->low <<= shift;
}
}
```
##### 加法運算
根據 C 語言規格書 6.2.5 小節第 9 段,無號整數的不會發生溢位(overflow)。當運算結果超出可以表示的最大值時,會自動對可以表示的最大值加一的數取模。以 64 位元無號整數為例,超出可表示範圍會對 $2^{63}+1=2^{64}$ 取模。
> A computation involving unsigned operands can never overflow, because a result that cannot be represented by the resulting unsigned integer type is reduced modulo the number that is one greater than the largest value that can be represented by the resulting type.
加法需要將 `low` 和 `high` 分別相加,但各自可能進位,根據 C 語言的計算規則取模,相當於直接捨棄超出 64 位元的部分。我使用 `carry` 這個變數來記錄 `low` 相加進位導致被捨棄的部分。若 `high` 的部分發生進位代表 128 位元仍然無法表示,按照 C 語言的規定直接捨棄。實作中有確保不會超過 128 位元的表示範圍。
將兩個無號整數 `a` 和 `b` 相加,最大值發生在兩者皆為 $2^{64}-1$,相加結果為 $2^{65} - 2$,需要 65 位元才能表達,可以總結發生溢位時,最多就是在第 65 位元加 1。所以在 `low` 相加時`carry` 只可能是 0 或 1,並將其加入到 `high`。
Claim: 在 `uint64_t c = a + b` 的情況下,`c < a` 若且唯若發生溢位
若在沒有溢位的情況下,`c` 必然大於或等於 `a` 和 `b`。設 64 位元整數可表示的最大值為 `ULLONG_MAX`,值為 $2^{64}-1$。在發生溢位的情況下,因為 `a` 和 `b` 各別都小於 `ULLONG_MAX`,所以 `c` 要表示的真實值 $C$ 最大不會超過 $2\times \text{ULLONG_MAX}$,所以儲存在 `c` 內的值為
$$
\begin{align}
c &= a + b - \text{ULLONG_MAX} \\
&=a - (\text{ULLONG_MAX} - b) < a
\end{align}
$$
得證。因此以下 128 位元整數加法正確。
```c
void bn_add(bn_t *lhs, const bn_t *rhs)
{
uint64_t low = lhs->low + rhs->low;
uint64_t carry = low < lhs->low;
lhs->low = low;
lhs->high = lhs->high + rhs->high + carry;
}
```
##### 減法運算
減法和加法概念類似,也是分別相減,但需注意借位。借位 `borrow` 的判斷很簡單,只需要比較減數和被減數即可。
```c
void bn_sub(bn_t *lhs, const bn_t *rhs)
{
uint64_t borrow = lhs->low < rhs->low;
lhs->low -= rhs->low;
lhs->high = lhs->high - rhs->high - borrow;
}
```
##### 大於或等於
兩數大小是從高位數開使比較,於是先判斷 `high` 的大小,若 `high` 不相同才需要近一步比較 `low`。
```C
bool bn_geq(const bn_t *lhs, const bn_t *rhs)
{
if ((lhs->high > rhs->high) ||
(lhs->high == rhs->high && lhs->low >= rhs->low))
return true;
return false;
}
```
##### 平方運算
此函式是用於測試,不包含在 `my_sqrtf_bn` 內。這個函式會把 64 位元無號整數 `x` 平方並存入 `out` 內。我先將 `x` 拆成高位 `h` 和低位 `l`,此時
$$
\begin{align}
x^2 &= (h\times2^{32}+l)^2 \\
&=(h^2\times2^{64} + 2hl\times 2^{32} + l^2)
\end{align}
$$
因為平方項 $h^2$ 和 $l^2$ 是 32 位元相乘,不可能會溢位,無需特別處理。但交叉項 $2hl$ 有可能溢位,需小心處理。
所以要將 $h^2$ 放入第 64 位元到第 127 位元,$l^2$ 放入第 0 位元到第 63 位元,而 $2hl$ 放入第 32 位元到第 95 或第 96 位元。
我先計算交叉項 `hl`,接著需要將其乘以 $2$,但要先檢查是否會溢位,並將其存在 `cross_carry`,可以注意到溢位最多 1 位元。接著僅需對齊位元,逐個相加即可。
```c
void bn_square(bn_t *out, uint64_t x)
{
uint64_t xh = x >> 32;
uint64_t xl = x & 0xFFFFFFFF;
uint64_t xh_xh = xh * xh;
uint64_t xl_xl = xl * xl;
uint64_t cross = xh * xl;
uint64_t carry, cross_carry;
cross_carry = (cross >> 63);
cross <<= 1;
out->low = xl_xl + (cross << 32);
carry = out->low < xl_xl;
out->high = xh_xh + (cross >> 32) + carry + (cross_carry << 32);
}
```
### 問題紀錄
> Linux 核心使用 digit-by-digit 的理由?Linux 核心有支援浮點數的根號嗎?
我想到可能的原因有
- 牛頓法需要除法
- 牛頓法需要選取初始值,可能會有表現不穩定的情況?根據課程教材[從 √2 的存在談開平方根的快速運算](https://hackmd.io/@sysprog/sqrt)中和牛頓法和 IEEE 754 + 二分搜尋法的比較結果,某些值的所需的收斂時間較長,或許可以用作 side channel attack,目前僅止於猜測,尚未實驗觀察。
#### 牛頓法執行時間跳點實驗
我以[教材](https://hackmd.io/@sysprog/sqrt#牛頓法求平方根)內的程式碼進行測試。
我首先測試了這個現象的可重現性,發現是可重現的,因此更有實驗觀察的價值。

我畫出了牛頓法迭代次數和時間的圖,可以看到迭代次數多的部分對應到執行時間較長。

而迭代次數和初始值的選擇有關,我目前沒有想到要如何預測迭代次數。牛頓法收斂非常快,且理論上只要初始值選取得當,收斂速率不受數值大小影響。實際測試整數範圍內不會有太大的執行時間變動,僅在固定範圍內浮動。因為預測迭代次數沒辦法很精準,迭代次數很少就會造成預測的誤差很大。所以我認為很難用牛頓法的執行時間進行 side channel attack 。
## Q2: 查閱 IEEE 754 規格書關於 exponent 的描述
根據 [IEEE 754](https://ieeexplore.ieee.org/document/8766229) 規格書,浮點數(floating-point number)應可以表示以下三種數值:
> - An Signed zero and non-zero floating-point numbers of the form $(−1)^s\times b^e\times m$, where
> - $s$ is 0 or 1
> - $e$ is any integer $emin\leq e\leq emax$
> - $m$ is a number represented by a digit string of the form $d_0.d_1d_2...d_{p-1}$ where $d_i$ is an integer digit $0\leq d_i<b$ (therefore $0\leq m<b$)
> - Two infinities, $+\infty$ and $-\infty$.
> - Two NaN, qNaN (quiet) and sNaN (signaling).
在 3.4 小節中給出的格式如下:
> 
> - 1-bit sign $S$
> - $w$-bit biased exponent $E=e+\text{bias}$
> - $(t = p − 1)$-bit trailing significand field digit string $T = d_1d_2…d_{p −1}$; the leading bit of the significand,
$d_0$, is implicitly encoded in the biased exponent $E$.
其中根據指數(biased exponent)的值 $E$ 分成三種
- $1$ 到 $2^w-2$ 之間的所有整數用來表示 normal number,此時對應浮點數的值為 $v=(-1)^S\times 2^{E-\text{bias}}\times(1+2^{1-p}\times T)$。normal number 的 implicit leading significand 為 1
- 保留 $0$ 來表示 $\pm0$ 和 subnormal number
- $E = 0$ 且 $T = 0$ 時,對應的值為 $v=(-1)^S\times(+0)$,(區分正負 0)
- $E = 0$ 且 $T \not= 0$時,對應的值為 $v=(-1)^S\times2^{\text{emin}}\times(0+2^{1-p}\times T)$。subnormal number 的 implicit leading significand 為 0
- 保留 $2^w-1$ 來表示 $\pm\infty$ 和 NaNs
- $E = 2^{w-1}$ 且 $T = 0$,表示正負無窮大 $v=(-1)^S\times(+\infty)$
- $E = 2^{w-1}$ 且 $T \not= 0$,表示 NaN,由 $d_1$ 來區分 qNaN 和 sNaN
根據 Table 3.5,可以總結 32 位元浮點數正數的表示範圍:

- normal number
- 最大值: $1.111...1_{(2)}\times 2^{127}$,大約為 $3.4\times10^{38}$
- 最小值: $1.0_{(2)}\times 2^{-126}$ 大約為 $1.17\times10^{-38}$
### 問題紀錄
> 為何 exponent 部分要設計 bias?
根據規格書 2.1 小節,設計 bias 的主要理由是讓指數部分以非負整數表示。
> biased exponent: The sum of the exponent and a constant (bias) chosen to make the biased exponent’s range non-negative.
> significand 和 mantissa 的差別
我過去看到的很多文章都是使用 mantissa 這個字,但我發現在 IEEE 754 的規格書內一次都沒有使用這個字,而是使用 significand。
[significand](https://en.wikipedia.org/wiki/Significand):
> the first (left) part of a number in scientific notation or related concepts in floating-point representation, consisting of its significant digits.
[mantissa](https://en.wikipedia.org/wiki/Common_logarithm#Mantissa_and_characteristic):
> An important property of base-10 logarithms, which makes them so useful in calculations, is that the logarithm of numbers greater than 1 that differ by a factor of a power of 10 all have the same fractional part. The fractional part is known as the mantissa.
significand 是指科學記號的小數部分,mantissa 原意是專指對數的尾數(小數部分)。因此我認為使用 significand 較為精確。但在 glibc 的 float.h 內定義的巨集 `# define FLT32_MANT_DIG FLT_MANT_DIG` 也採用了 mantissa 的說法,或許用 mantissa 也可以。
## 參考資料
- [IEEE 754](https://ieeexplore.ieee.org/document/8766229)
- [從 √2 的存在談開平方根的快速運算](https://hackmd.io/@sysprog/sqrt#:~:text=存在談開-,平方根,-的快速運算)