第十五章　用組語實作對數函數

這一章將介紹用組合語言，配合 SSE 指令集，撰寫對數函數。這兒的對數是指自然對數，也就是以 e 為底的對數，一般記為 ln。這兒的 e 稱為自然常數或是尤拉數（Euler's number，e = limn→∞( 1+1n )ⁿ ），約等於 2.71828182845904523536……，是無限不循環小數。計算出自然對數之後，很容易就能經由換底公式計算出常用對數來，log x = ln xlog_e 10，其中的分母 log_e 10 為一常數，約等於 0.4342944819。

原理

利用基礎對數律與泰勒展開式計算對數

想要算出任意一個正浮點數 x 的自然對數，可以依照下面方式思考。任意正數，x 可以寫成

x = m × 2^p;兩邊取自然對數
ln x = ln ( m×2^p ) = ln m + p×ln 2;此為第⑴式，且 1≤ m < 2，p 為整數

上面的式子中，ln 2 是常數，約等於 0.69314718055994530942。p 是 2 的冪次方，是整數，在 IEEE 754 雙精確度浮點數的表示方式中，位於第 52～62 位元（請參閱第十一章）中。但從這十一個位元得到的是指數偏移值，還須減去 1023 才是真正的指數。到此，上式的第二項已經算出來了。

至於計算第一項，先要取得有效數 m，它存放於雙精確度浮點數的第 0～51 位元，但這部分僅僅是小數部分，還須處理正規化的問題，因此程式如下：

        movq    r9,xmm0
        and     r9,0FFFFFFFFFFFFFh   ;只留有效數
        or      r9,3FF0000000000000h ;正規化，使RDX=m=1.xxxx

這樣 R9 暫存器內就是 m，接下來就可以用泰勒級數計算 ln m 了。但這兒用一個較偷懶的想法：考慮一個無窮等比數列，首項，a₁=1；公比，r= -y；且 | y | < 1，所以這個無窮等比數列是 1、-y、y²、-y³、y⁴、-y⁵……，如果要計算它的總和，代入等比級數公式： a₁-a₁rⁿ 1 - r，因為 a₁=1，r=-y，又因為 n→∞ 故 rⁿ=0，所以此等比數列總和是11+y，

1 - y + y² - y³ + y⁴ - y⁵ +……= 11+y;反過來寫  
11+y = 1 - y + y² - y³ + y⁴ - y⁵ + ……;兩邊取積分  
∫11+y dy = ∫( 1 - y + y² - y³ + y⁴ - y⁵ +……) dy;根據微積分，11+y 的積分等於 ln(1+y)
ln(1+y) = y - y²2 + y³3 - y⁴4 + y⁵5 - y⁶6 + ……+ C;欲求 C，可使 y=0 代入，使等號右邊的 y、y²2……皆為 0，得 C=0
ln(1+y) = y - y²2 + y³3 - y⁴4 + y⁵5 - y⁶6 + ……;第⑵式

令 y = m-1，則 m= 1+y，代入⑵得

ln m = y - y²2 + y³3 - y⁴4 + y⁵5 - y⁶6 + ……

將上式代入第⑴式，就可以算出 ln x 了。與前一章一樣，可以用秦九韶算法把上式轉變成

ln m = y( 1- y2 + y²3 - y³4 + y⁴5 - y⁵6 + ……)
     = y( 1-y( 12 - y3 + y²4 - y³5 + y⁴6 - ……) )
      ⁝
     ≈ y( 1-y( 12 -y( 13 -y( 14 -y( 15 -y( 16 - y7 ) ) ) ) ) ) ;假如只近似到 y⁷

藉由秦九韶算法可以使用迴圈，讓程式編寫變得精簡。

最佳化

用上面的對數律與泰勒展開式，來計算自然對數，是可行的。但是因為正整數從 1 開始，而其倒數，1正整數，收斂得很慢，尤其是當 m 離 2 越近，y 就離 1 越近，收斂得越慢。改良的方法是，當 m > √2 時，把 m 除以 2，然後使 p 增加一，這樣的話能收斂得快一些。

回到問題的源頭，我們要計算的是 ln x，所以把 x 轉換成類似科學記號的表示方式，但並非以 10 為底數，而是以 2 為底數：

x = m × 2^p

有效數 m 的範圍是 1 ≤ m < 2，而我們使用的泰勒展開式是 y 的函數，而 y=m-1，故 y 的範圍在 0 ≤ y < 1。當 y 越接近 1，泰勒展開是收斂得越慢。這裏說的最佳化是假如 m > √2 時，把 m 除以 2，然後使 p 增加一。

這時會遇到兩種情形：

⑴、如果 1 ≤ m ≤ √2：因為 m 已經夠小了，所以程式什麼也不做，m 的範圍是 1 ≤ m ≤ √2。

⑵、如果√2 < m < 2：原來的數，x，大小仍不變，但把有效數縮小變為一半，縮小的部分挪移到指數去了：

x = m × 2^p = m2 × 2^p+1

m 原來的範圍是 √2 < m < 2，把 m 除以 2 之後，新的 m 範圍在 √22 < 新的 m < 1，整理後 √0.5 < 新的 m < 1。

因為 √0.5 約等於 0.7071，所以綜合上述兩種情形，m 的範圍落在 √0.5 < m ≤ √2。而我們計算 ln m 的泰勒展開式是 y 的函數，而 y= m-1，因此 y 的範圍是 √0.5 -1 < y ≤ √2 -1，整理後得到 y 的範圍約在 -0.2929 < y ≤ 0.4142。這個範圍更接近 0，比未最佳化的範圍，0 ≤ y < 1，好多了。

C 語言標準函式：strtod

strtod 是 C/C++ 標準函式，它能把 ANSI 數字字串轉換為雙精確度浮點數，用法是：

double strtod(
   const char *strSource,
   char **endptr
);

*strSource 是要轉換的 ANSI 字串位址，該字串內含有一連串的阿拉伯數字並以 0 結尾，strtod 會將這字串轉換成雙精確度浮點數，並存於 XMM0 暫存器返回。**endptr 是某個變數位址，這個變數必須是四字組長，strtod 會把結尾的 0 之位址或是第一個無法轉換為數值字元的位址，存入於此變數中。應用程式藉此檢查 endptr 變數內的位址所指向的地方是否為 0，如果是 0，表示正常結束；如果不是 0，表示轉換的 ANSI 字串含有未能轉換成數值的字元。

要轉換的 ANSI 字串可以以數個空白或跳格字元起頭，strtod 會忽略。此字串在阿拉伯數字字元之前，可以添加一個「+」或「-」字元，表示這轉換後的雙精確度浮點數是正數還是負數，但正負號與數字字元之間不能有空白。要轉換的 ANSI 字串也可以包含「E」或「e」，代表時的冪次方。例如底下的例子：

szValue1    DB      9,"  96485.3321",0  ;XMM0=96485.3321，strtod 可以跳過前面的空白或跳格字元，9 是跳格字元
szValue2    DB      "96,485.3321",0     ;XMM0=96，「,」是無效的字元
szValue3    DB      "96485.33 21",0     ;XMM0=96485.33，不在開頭的空白或跳格字元，視為無效的字元
szValue4    DB      "-9.6485E4",0       ;XMM0=-96485，「E」後面的數字字元代表十的冪次方
szValue5    DB      "9.6485E-4",0       ;XMM0=0.00096485，「E」後面的數字字元代表十的冪次方
szValue6    DB      "2.175E 6",0        ;XMM0=2.175，「E」後面不可接空白，如果這樣會導致連「E」也是無效的

如果 strtod 成功地將字串轉換成雙精確度浮點數，那麼回傳值就是 XMM0，其低半部就是轉換後的雙精確度浮點數。如果只有字串的前面部分能轉換成數值，而後有未能轉換成數值的字元，那麼 XMM0 會等於已轉換的數值，而 endptr 是第一個無法轉換為數值之字元的位址。如果字串的第一個字元就是無法轉換成數值，或是字串最前面是數個空白或跳格字元而緊接其後的字元無法轉換成數值，那麼 XMM0 等於 0.0，而且 endptr 內的位址就是字串的位址。

GetCommandLine API

GetCommandLine 可以獲得命令提示列上的字串。它沒有參數，所以呼叫方式可以直接使用 x64 指令，CALL，因為 MASM64 SDK 的 invoke 巨集，都是為了將參數存入暫存器或堆疊的，既然 GetCommandLine 沒有參數，就可以直接使用 CALL 指令。底下是 GetCommandLine 的原型：

LPSTR GetCommandLine();

GetCommandLine 的回傳值是 RAX，代表命令字串的位址。所謂命令字串是指使用者在命令提示字元底下所輸入的文字，包含應用程式名稱及之後的參數（亦可稱為選項），最後還有結尾的零。例如下圖中紅色框內，使用者執行「ln.exe」時，在「ln」之後緊接著應用程式能接受的參數「300.889」，那麼命令字串就是「ln 300.899」最後再加上結尾的 0。

如果應用程式在「E:\HomePage\SOURCE\My App\LN.EXE」，而當前目錄是「C:\Users\me\」，那麼要執行 LN.EXE，須以下面方式輸入指令。這是因為「My App」含有空白字元，必須用一對 " " 括住。

C:\Users\me>"E:\HomePage\SOURCE\My App\LN" 300

這時候，命令字串變成「"E:\HomePage\SOURCE\My App\LN" 300」,0。如果使用者輸入時，僅用一對 " " 把含空白字元的目錄名稱括住，也是可以的，例如輸入

C:\Users\me>E:\HomePage\SOURCE\"My App"\LN 300

這樣的話，命令字串就變成「E:\HomePage\SOURCE\"My App"\LN 300」,0。但是不可以像下面那樣：

C:\Users\me>"E:\HomePage\SOURCE\My App\LN 300"

這樣的話，命令提示字元會以為要執行「LN 300.EXE」可執行檔，而非「LN.EXE」，這當然會發生錯誤。

SSE 指令：CVTSI2SD 與 CVTSI2SS

CVTSI2SD 指令全名為 convert signed integer to scalar double-precision floating-point，意思就是把一個 32 位元或 64 位元的有號整數轉換為雙精確度浮點數，並存入 XMM 暫存器的低半部。它的語法是

CVTSI2SD XMM暫存器,來源運算元
CVTSI2SS XMM暫存器,來源運算元

CVTSI2SD 的目的運算元只能是 XMM 暫存器，來源運算元可以是通用暫存器或記憶體變數，但不能是 XMM 暫存器。來源運算元可以是 32 位元或 64 位元，不過都會轉換成雙精確度浮點數，存於目的運算元的低半部，而高半部不會清零。

要注意的是 64 位元有號數的精確度是 63 位元，而雙精確度浮點數的有效數是 53 位元，所以如果是很大或很小的整數（絕對值超過 2⁵³），就會依據 MXCSR 暫存器的 RC 位元，進行捨入，因此這種情形會損失精確度。至於 32 位元整數轉換成雙精確度浮點數，就沒有這個問題。

CVTSI2SS 指令與 CVTSI2SD 類似，把位於來源運算元內的 32 位元或 64 位元有號整數，轉換成單精確度浮點數，並存於目的運算元的 XMM 暫存器第 0～31 位元內。XMM 暫存器的 32～127 位元不受影響。單精確度浮點數的有效位數僅 24 位元，因此如果是很大或很小的整數，也會因為捨入而面臨精確度損失的問題。

CVTSI2SS 與 CVTSI2SD 的來源運算元如果是記憶體變數，不需要對齊節位址；但是如果對齊，會有較好的效率。這兩個指令都不會產生無效操作、溢位、除零等例外；但如果發生捨入，會設定 MXCSR 暫存器的 PE 旗標。

SSE 指令：UCOMISD 與 UCOMISS

UCOMISD 是用來比較兩個雙精度浮點數的 SSE2 指令，全名是 unordered compare scalar double-precision floating-point values and set integer EFLAGS。它的語法是

UCOMISD XMM暫存器,來源運算元
UCOMISS XMM暫存器,來源運算元

UCOMISD 的目的運算元只能是 XMM 暫存器，來源運算元可以是 XMM 暫存器或 64 位元長的記憶體變數。UCOMISD 將目的運算元的低半部與來源運算元比較（如果來源運算元是 XMM 暫存器，那就比較低半部），比較結果並不寫回目的運算元，而是設定旗標暫存器，然後搭配 JA、JB、JAE、JBE、JE、JNE、JP 等條件跳躍指令。

UCOMISS 與 UCOMISD 類似，但是 UCOMISS 的目的運算元只能是 XMM 暫存器，來源運算元可以是 XMM 暫存器或 32 位元長的記憶體變數。UCONISS 只比較目的運算元的 0～31 位元單精確度浮點數與來源運算元比較（如果來源運算元是 XMM 暫存器，那就比較第 0～31 位元），其他性質都與 UCOMISS 相同。

旗標暫存器（EFLAGS）常用於整數（integer）比較，所以 UCOMISD、UCOMISS 中有個 I，就是因為這原因。比較結果如下：

條件	ZF	CF	PF	條件跳躍指令
目的運算元>來源運算元	0	0	0	JA/JNBE
目的運算元=來源運算元	1	0	0	JE/JZ
目的運算元<來源運算元	0	1	0	JB/JC
任一運算元為 NaN	1	1	1	JP

此外要注意的是，浮點數不分有號數或無號數，所以英特爾規定 UCOMISD 之後的條件跳躍指令都統一使用 JA、JB，不可使用 JL、JG。例如底下的例子，因為－2.1>－5.4，所以會發生跳躍至 greater: 處，故最後 AX=1。程式碼如下：

x        REAL8   -2.1
y        REAL8   -5.4
          ⁝
         movq    xmm6,x
         ucomisd xmm6,y
         ja      greater
         mov     ax,0
         jmp     short ok
greater: mov     ax,1
ok:

UCOMISD、UCOMSS 對於特殊數值的處理，也很聰明。例如它們會把忽略負零的符號，因此負零等於正零。對無限大的處理是：正無限大 > 任何有限數 > 負無限大。對 NaN 的處理是只要比較的運算元中，有任何一個運算元是 NaN，那麼就會設定同位旗標（parity flag，註一）。

NaN 是 not a number 的意思，中文可翻譯成非數值，根據 IEEE 754 只要是指數偏移值所有位元均為 1，有效數字不可全為 0（只要任一位元或數個位元是一即可），這樣的資料稱為 NaN。NaN 有幾個特殊的特性：
①、NaN 與任何數（包括自身）比較都不相等。
②、涉及 NaN 的大部分運算結果仍為 NaN。例如 NaN+5、NaN×0 等結果都是 NaN。

NaN 分為兩種：

⑴、QNaN：Q 是 quiet 的縮寫，QNaN 可翻譯成靜默型非數值，有效數字的最高位元為 1 的 NaN 就是 QNaN。當程式執行了數學上無法定義或超出了浮點數表示範圍的運算時，例如 0÷0、∞－∞、∞×0、負數平方根等等，就會產生 QNaN。當產生 QNaN 時，不會觸發任何硬體異常或例外（exception），而是默默的（quiet）把 NaN 的狀態傳播下去。

⑵、SNaN：S 是 signaling 的縮寫，SNaN 可翻譯成信號型非數值，有效數字的最高位元為 0 的 NaN 就是 SNaN。一般的數學運算幾乎不會產生 SNaN，而是預先把 SNaN 存於記憶體變數裏，由應用程式去讀取這個含有 SNaN 的變數，主要目的是為了除錯。當 CPU 遇到 SNaN 參與運算時，會自動將 MXCSR 暫存器的第零位元，IE 設為 1，表示發生無效運算；再檢查 MXCSR 暫存器的 IM 位元是 0 還是 1：

ⓐ、若 IM 位元為 0，就引發硬體中斷，程式崩潰或進入例外處理程式。

ⓑ、若 IM 位元為 1，不引發中斷，自動將 SNaN 變為 QNaN，繼續執行原來的應用程式。

實作自然對數函數

瞭解了原理之後，就可以用組合語言寫出計算自然對數的原始程式了。原始程式如下：

;LN.ASM用於計算自然對數函數。

INCLUDE E:\masm64\include64\masm64rt.inc

EXTRN strtod:PROC,printf:PROC

;****************************************************************************************

.DATA

number DQ ?

endptr DQ ?

szHint DB "使用方式：於命令提示字元輸入",0Dh,0Ah," ln 正數",0Dh,0Ah,\

"本應用程式就能幫你算出該正數的自然對數（即ln值）。",0

szErr DB "無法計算零或負數的對數值，或輸入非數值。",0

szFmt DB "ln %g=%.15g",0

;****************************************************************************************

.CODE

;----------------------------------------------------------------------------------------

ln PROC x:QWORD

;計算x的自然對數，即計算ln x

;輸入：XMM0－XMM0內存有要計算自然對數的數值(雙精確度浮點數)，且必須是正數

;輸出：若RAX=0，表示輸入的x為0、負數、-∞及NaN，無法計算其自然對數

;　　：若RAX=1，XMM0=ln x(若x=+∞，ln x=+∞，RAX=1)

;原理：ln x=ln(m×2^p)=ln m + p*ln 2

;　　　ln m=y-(y^2)/2+(y^3)/3-(y^4)/4+...

movq x,xmm0

jmp start

ALIGN 16

SignMsk DQ 8000000000000000h,0 ;變號

coeff DQ 0.05555555555555555556,0.05882352941176470588,\ ;1/18, 1/17

0.0625, 0.06666666666666666667,\ ;1/16, 1/15

0.07142857142857142857,0.07692307692307692308,\ ;1/14, 1/13

0.08333333333333333333,0.09090909090909090909,\ ;1/12, 1/11

0.1, 0.11111111111111111111,0.125, 0.14285714285714285714,\ ;1/10, 1/9, 1/8, 1/7

0.16666666666666666667,0.2, 0.25, 0.33333333333333333333,\ ;1/6, 1/5, 1/4, 1/3

0.5, 1.0 ;1/2, 1

one DQ 1.0 ;1.0=3FF0000000000000h

two DQ 2.0

ln2 DQ 0.69314718055994530942

sqrt2 DQ 1.41421356237309504880

ExpMsk DQ 7FF0000000000000h ;ExpMsk恰好也是正無限大

ManMsk DQ 000FFFFFFFFFFFFFh

HidBit DQ 0010000000000000h ;第52位元，正規化數的隱含位元

start: movq r8,xmm0

movq r9,xmm0

;檢查x是否為0(包含-0)或負數(包含-∞、負NaN)

or r8,r8

jg posi ;若為正數、正無限大、正NaN，跳至posi；

inval: xor rax,rax ;若為負數或0，則使RAX變為0，退出ln副程式

jmp exit2

;到此已過濾掉負值與零，先檢查x是否為正NaN，再檢查是否為+∞

posi: ucomisd xmm0,xmm0

jp inval ;若PF=1，代表x=+NaN

cmp r8,ExpMsk ;檢查x是否為+∞

je exit1

;到此，x為大於零的正規化數或大於零的非正規化數

and r8,ExpMsk ;R8=指數偏移數

jnz normal

;處理非正規化數

mov r8,-1022

nxtbit: test r9,HidBit ;檢查有效數的第52位是否為1

jnz done ;如果開頭變1了，代表正規化完成

shl r9,1 ;有效數左移1位

dec r8 ;指數減 1

jmp nxtbit

;處理正規化數：x=m×2^p，計算指數p，並存於R8

normal: shr r8,52

and r9,ManMsk ;R9=有效數

sub r8,1023 ;R8=p

;計算有效數，m，並存於XMM2

done: or r9,one ;正規化，使R9=m=1.xxxx

movq xmm3,r9

movq xmm2,r9 ;XMM2=m，m的範圍在[1.0, 2.0)

cmpsd xmm3,sqrt2,6;若m大於sqrt2，則m=m/2且p=p+1這樣可以讓

movq rax,xmm3 ;y更接近0，大幅加快泰勒級數收斂速度

or rax,rax

jz lt_sq2

inc r8 ;p=p+1

divsd xmm2,two ;m=m/2

lt_sq2: cvtsi2sd xmm1,r8 ;XMM1=p

subsd xmm2,one ;XMM2=y

mulsd xmm1,ln2 ;XMM1=p*ln 2

;計算ln m=ln(1+y)=y*(1-y*(0.5-y*(0.333..-y*(0.25-y*(0.2-y*(0.166..-y*(0.142....)))))))

movsd xmm3,xmm2

xorpd xmm3,SignMsk ;XMM3=-y

lea r8,coeff

movsd xmm0,QWORD PTR [r8]

mov rcx,LENGTHOF coeff-1

next: add r8,8

mulsd xmm0,xmm3

addsd xmm0,QWORD PTR [r8]

loop next

mulsd xmm0,xmm2 ;XMM0=ln m

addsd xmm0,xmm1 ;XMM0=ln x=ln m + p*ln 2

exit1: mov rax,1

exit2: ret

ln ENDP

;----------------------------------------------------------------------------------------

main PROC

call GetCommandLine

;找到第一個參數所在位址，方法是每個參數前都有空白字元

nxtchr: cmp BYTE PTR [rax],22h

je quot

cmp BYTE PTR [rax]," "

je cvt2d

q_exit: inc rax

jmp nxtchr

quot: inc rax ;處理兩個"之間的字元

cmp BYTE PTR [rax],22h

jne quot

cmp BYTE PTR [rax],0

jz hint

jmp q_exit

cvt2d: invoke strtod,rax,ADDR endptr ;將參數轉換成雙精確度浮點數

mov r8,endptr

cmp BYTE PTR [r8],0 ;若[R8]不為0，表示使用者輸入了非數字字元

jz calc

hint: lea rcx,szHint

jmp pnt_it

calc: movq number,xmm0

call ln

or rax,rax ;若RAX=0，表示輸入零或負值

jnz ok

lea rcx,szErr

jmp pnt_it

ok: movq r8,xmm0

mov rdx,number

lea rcx,szFmt

pnt_it: call printf

invoke ExitProcess,0

main ENDP

;****************************************************************************************

END

將它儲存成 LN.ASM 之後，進入命令提示字元，依據前面提過的方法組譯，並在命令提示字元下執行。底下說明 LN.ASM。

ln 副程式

ln 副程式只有一個參數，而且此參數為雙精確度浮點數，依照 x64 呼叫慣例，主程式須將此參數存於 XMM0 傳遞給 ln。進入 ln 之後，首先把存於 XMM0 的參數存入堆疊，即影子空間（第 23 行）。

接下來是跳躍指令，跳至 start: 標記處。在 JMP 與 start: 標記處，存放的是 ln 需要用到的常數資料。也須有人會問，為何不把這些資料放在 .DATA 或 .CONST 等資料區段，卻要這樣安排？這是因為，往後小木偶打算將 ln 副程式獨立出來，製作成 DLL 檔，供其他應用程式呼叫。為了將其封裝起來做成「黑箱」，意思是只要給 ln 參數它就能計算出自然對數，與其他外在環境無關，所以才把這些常數放在 .CODE 區段內。

因為對數的定義域不可是負數或零，所以接下來要檢查參數是否為負數或是零，如果是就應該使 RAX 變為零並退出 ln 副程式。再者，IEEE 754 定義了幾種特殊的資料，例如非數值（NaN）、正負無限大、負零與正零，也需特別處理，分成底下幾種情況：

⑴、第 43 行的 OR 指令檢查存放在 R8 的參數是否為負數或零（前兩行已將使用者輸入的參數搬移至 R8）

        or      r8,r8
        jg      posi    ;若為正數，跳至pass；

這兩行指令用得相當精妙，僅兩行就能檢查參數是否為正值。還記的 JG 指令要發生跳躍的條件，必須是 ZF=0 且 SF=OF。如果參數是 0，執行「OR R8,R8」後，ZF=1，因此不會跳躍；如果參數是負數（包含負無限大、負零、負數、-NaN），那麼符號位元是一，執行 OR 後 SF 為一，但是執行 OR 指令後不論結果如何，必定清除 OF，故 SF 必不等於 OF，不會跳躍。所以綜合來看，只有正的浮點數（包含正無限大、正 NaN、一般正數），才會跳躍至 posi: 處。

綜合上述說明，當參數是零（其實零和正零的 IEEE 754 編碼相同）、負無限大、負零、負數、-NaN 都會過濾掉。這幾種情形，都應把 RAX 設為 0，表示輸入的參數是無效的。

⑵、通過⑴之後的參數，還有正的非數值（+NaN）或正無限大必須過濾掉，但這兩種情形不同。先說 +NaN，可以用 UCOMISD 除掉 +NaN（第 48 行）。如果參數是 +NaN，應把 RAX 設為 0，表示輸入的參數是無效的。

⑶、第 50 行，檢查參數是否為正無限大，這道指令沒有特別之處，就只是把參數與代表正無限大的十六進位數值（此數值是 7FF0000000000000h）比較，如果相同，就代表參數為正無限大。因為正無限大取對數值，也是正無限大，並非無效的參數，但卻也不必去計算，應直接將 RAX 設為一，且使 XMM0 保持原來的正無限大返回即可。

⑷、如果能執行到第 53 行，表示輸入的參數為一般正數，而且不是 +NaN，也不是正無限大。一般正數可能是正規化數，也可能是非正規化數，計算它們的對數值過程不太一樣，因此第 54 行判斷是正規化數還是非正規化數，如果是正規化數，就跳躍至第 66 行，依據上面提過的對數律、最佳化、以秦九韶算法計算泰勒展開式，便可求出自然對數。

如果是非正規化數，過程稍稍不同，這兒就說說如何計算非正規化數的自然對數。

非正規化數

第十一章提過，不論是單精確度還是雙精確度浮點數，其格式均分為三部分：最高的一個位元是符號位元（sign）、其次是指數偏移值（exponent bias）、最後是有效數字（significand）。對雙精確度浮點數來說，其格式為浮點數的數值大小為

 (-1)^sign×1.significand×2^exponent

上面的 exponent 等於指數偏移值減 3FFh 或 7Fh（對雙精確度而言是 3FFh，即十進位 1023；對單精確度而言是 7Fh，即十進位 127）。此外 1.significand 中的 1 稱為隱含整數，它可使精確度多一位。採用這樣的方法表示的數值，稱為正規化數（normal number），所以正規化數的定義是指數偏移值既不是全部都是零、也不是全部都是一，並且具有一位隱含整數，其值是 1。

對雙精確度浮點數而言，可以推論出最大的正規化數是 7FEF FFFF FFFF FFFFh，換算成實際數值是 1.111……1×2^7FE-3FF ≈ 1.798×10³⁰⁸（1.111……1 是二進位，小數點後共 52 個 1）；而最小值是 0010 0000 0000 0000h，換算成實際數值是 1.0×2^1-1023 ≈ 2.225×10^-308。

假如把這個最小的正規化數除以 2，結果指數偏移值會變為零，無法用正規化數表示，如果捨棄變成零，就變得很奇怪，像這樣的情形叫突然下溢（abrupt underflow）。非正規化數就是為了解決這樣的問題，應運而生的。非正規化數是指有效數字不全是零，而且指數偏移值為零，並且隱含一位整數，其值是 0。非正規化數的指數是固定值，對雙精確度浮點數而言，指數是－1022；對單精確度浮點數而言，指數是－126。除了指數固定之外，非正規化數的隱含整數是零，不是一。

根據非正規化數的規則，對雙精確度浮點數而言，最大的非正規化數是 000F FFFF FFFF FFFFh，換算成實際數值是 0.111……1×2^-1022 ≈ 2.225×10^-308（0.111……1 是二進位，小數點後共 52 個 1）。它比最小的正規化數稍微小一點，也就是說，它接續著最小的正規化數往零的方向繼續變小。當計算結果小於最小正規化數時，不是瞬間變為 0，而是逐漸「流」向 0，非正規化數來保留盡可能多的精確度，這樣的情形稱為漸進下溢（gradual underflow）。

順帶提一提，對雙精確度浮點數而言，最小的非正規化數是 0000 0000 0000 0001h，換算成實際數值是 0.000……01×2^-1022 = 2^-52×2^-1022 ≈ 4.941×10^-324。其中 0.000……01 是二進位，小數點後共 51 個 0，只有第零位元是 1。

瞭解了非正規化數之後，就能計算了。參考右圖，假設有個非正規化雙精確度浮點數，其編碼為 0001 C000 0000 0000h，圖中有效數字是以黃色和紫色表示，黃色代表 0，紫色代表 1。圖中有效數以二進位表示是 0.000111000……0（1 後面有 46 個 0），為了要符合泰勒展開式計算 ln m 中，m 須在 1 與 2 之間，所以必須把圖中的有效數往左移 4 個位元。還記得每往左移一個位元，相當於乘以二，但又得使 m 維持不變；如果使 2 的冪次方減一，就能達成目的。原始碼第 57～63 行便是處理這事。

經由第 57～63 行得到非正規化數的指數與 m 後，就可以併入計算正規化數的程式中，執行最佳化、再以秦九韶算法計算泰勒展開式，便可求出自然對數。

得到參數

LN.EXE 執行方式與前面的例子都不相同。當使用者在命令提示字元中執行 LN.EXE 時，必須要在應用程式之後輸入一個正數，像下面那樣，所輸入的正數一般稱為參數獲選項：

E:\HomePage\SOURCE\Win64\CONSOLE>ln 57.3

這樣的使用方式，就好像數學家在紙上計算 57.3 的自然對數一樣，寫「ln 57.3」。「ln 57.3」稱為命令字串，而此處的 57.3 並非以浮點數格式存在記憶體內，而是以數字字串方式存於記憶體中。那麼，要怎麼樣才能把數字字串的「57.3」，變成浮點數的「57.3」呢？大約有底下三個步驟：

呼叫 GetCommandLine 得到命令字串位址。
找到數字字串位址（第 103～115 行）。
呼叫 strtod 將數字字串轉換成浮點數。

第一步及第三步很容易做到，但第二步就有些麻煩。如果仔細想想，便可以想得到，參數與應用程式名稱之間以空白字元隔開，所以從第一步得到命令字串位址後，便可經由搜尋空白字元而得到參數所在位址。但不幸的是，理想很豐滿，現實很骨感，凡事總有例外。在這裏的例外，發生在假如 LN.EXE 不在當前目錄裏，而在有空白字元的目錄中。舉個例子說明，假如 LN.EXE 在「E:\MY APP\LN.EXE」，而當前目錄是「D:\MATH」，那麼要執行 LN.EXE 時，可以像下面那樣：

D:\MATH>
D:\MATH>E:\"MY APP"\LN 57.3

目錄名含有空白字元時，在使用該目錄名稱時，必須用一對 " " 將目錄名稱括住，就像上面那樣。這時候就不能單單只要是找到空白字元，就認為後面接著的就是參數，但前面也說了，如果是在一對 " " 內的空白字元，就不是用來隔開參數的而是目錄名稱的一部份。原始程式的第 103～115 行就是藉由這樣的分析，而撰寫出來的。

註釋

註一：同位旗標（Parity Flag）與 x64 指令：JP/JPE、JNP/JPO

同位旗標也稱為奇偶旗標，它位於延伸旗標暫存器的第二個位元，簡稱 PF，早在 8086/8088 CPU 就已使用此位元：運算的結果換成二進位後，最低的八個位元中，若有偶數個 1，則此位元設為 1，反之為 0。

例如底下的例子：

        mov     cx,103h
        or      cx,cx

運算後 CX=0000 0001 0000 00011b，最低的八個位元中只有兩個位元是一，因此 CPU 會把同位旗標設為一。

JP/JPE、JNP/JPO 是 x64 指令，JP 與 JPE 是相同的指令，前者是 Jump short if parity，後者是 Jump short if parity even，兩者都是當同位旗標為一時，發生跳躍。JNP 與 JPO 也是相同的指令，前者是 Jump short if not parity，後者是 Jump short if parity odd，兩者都是當同位旗標為零時，發生跳躍。JP 與 JNP 的語法是：

JP      標記
JNP     標記

標記（label）是一串英文字組成，也可以包含「_」、阿拉伯數字，必須以「:」結尾，它代表著位址。例如底下的例子：

        mov     r8d,103h
        inc     r8d
        jnp     ok
         ⁝
ok:

執行完 INC 後，R8D 的 32 個位元裏有兩個位元是一，但第 0～7 位元只有第二個位元是一（奇數個一），其他七個位元都是零，因此 INC 會把同位旗標清除，故執行 JNP 時，會發生跳躍至 ok: 標記處。

在某些除錯器（例如 DOS 的 DEBUG 或 SYMDEB）中，以縮寫表示同位旗標的狀態：若 PF 被設定變為 1，顯示 PE（parity even 的縮寫，表示偶校驗）；若 PF 被清除變為 0，顯示 PO（parity odd 的縮寫，表示奇校驗）。

同位旗標的設計初衷非常簡單，是為了反映運算結果的「最低八個位元」中，1 的個數是否為偶數。這種情形最常用於早期的資料通訊（如 RS-232 串列埠）中，可用來快速計算資料的校驗位元，以檢查傳輸的正確性。但如今，大多數資料傳輸都由硬體檢驗（如 USB、PCIe、SATA）；而網路通訊也已內建更強大的偵錯機制（如 CRC、Checksum），並且也交由專用硬體完成。因此，一般應用程式開發者幾乎不需要再用同位旗標來檢驗資料正確與否。

儘管同位旗標的原始用途已式微，但是英特爾卻賦予了它一個重要的「第二專長」，那就是專門用來偵測浮點數運算中的 NaN。當 UCOMISD 或 UCOMISS 的運算元中，有任何一個是 NaN 那麼就會設定同位旗標，之後應用程式可以利用 JP 或 JNP 跳躍至某處，去處理 NaN 的情形或不是 NaN 的情形。

第十五章 用組語實作對數函數

原理

利用基礎對數律與泰勒展開式計算對數

最佳化

C 語言標準函式：strtod

GetCommandLine API

SSE 指令：CVTSI2SD 與 CVTSI2SS

SSE 指令：UCOMISD 與 UCOMISS

實作自然對數函數

ln 副程式

非正規化數

得到參數

註釋

註一：同位旗標（Parity Flag）與 x64 指令：JP/JPE、JNP/JPO

第十五章　用組語實作對數函數