Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

附加的實現說明(Additional Implementation Notes)

實現說明的主要部分是每個函數和分佈的文檔。這裡的說明是更普遍的性質並反映所使用的更普遍的實現哲學(implementation philosophy)。

實現哲學(Implemention philosophy)

"首先讓其正確,然後讓其快速(First be right, then be fast)."

在速度和精度之間永遠都會有潛在的妥協。可能找到更快的算法,尤其是對於特定範圍的參數,但是對於大多數的使用數學函數和數學分佈的應用程序來說,我們判定(judge):速度極少與精度一樣重要。

所以我們的優先級是精確度。

為了允許特殊函數的計算精度,對生成非常精確的測試數據表做出了巨大的努力。

(同樣也要求使用更多的CPU資源 - JM的筆記本電腦有流出塑料熔漿的危險,所以,PAB的雙核台式電腦在計算一些測試數據表的日子裡的50%的時間處於繁忙的狀態).

對於一個特定的RealType類型,比如說float 或 double,對於一些函數可能可以找到更小因此也更快的逼近方法,但精確度更低一些 (或許因為沒有優化的迭代,例如,當計算反函數的時候)。

如果這些算法「對於滿足要求」是足夠精確的,那麼用戶可以使用他自己的特定實現方法。

例如,從計算代價很昂貴的時代開始就存在逼近算法:

H Goldberg and H Levine, Approximate formulas for percentage points and normalisation of t and chi squared, Ann. Math. Stat., 17(4), 216 - 225 (Dec 1946).

A H Carter, Approximations to percentage points of the z-distribution, Biometrika 34(2), 352 - 358 (Dec 1947).

對於一些速度至關重要的應用程序來說,這些算法仍然可以提供足夠的精度。

精度以及測試數據的表示

為了使足夠多的浮點類型滿足足夠的精度,如果可能的話,常量值被給定為50個十進制數字(雖然許多資料只提供精度為接近64個bit的double精度)。除非它們可以被精確表示,例如,整數,或類似於0.125的分數,否則通過在數值的後面附加一個L使得數值被當作long double類型。這避免在轉換為double(缺省的類型)時發生精度損失。在 static_cast<RealType>之後使用來為抽樣測試提供合適的RealType類型。

返回常量值的函數,類似於函數 kurtosis ,寫作:

static_cast<RealType>(-3) / 5;

來提供編譯器可以為RealType計算的最精確的值。 (分母是一個整數,所以可進行精確的類型提升)。

測試1/3,不能 使用基底(radix)為浮點數的兩個值來表示,而(應當)使用,例如:

static_cast<RealType>(1) / 3;

如果函數對於輸入非常敏感,將一個不精確的值(例如 0.1 )作為輸入可能會使結果偏移很遠:例如 ~1e-7 是「錯誤的」,因為0.1沒有精確的二進製表示。這就是為什麼精確的二進制值-1/2 , 1/4 , 1/8等等-用在測試中,以及偶爾使用的分數a/b,其中,b是2的冪(為了確保結果是一個精確表示的二進制值)。

容許度測試(Tolerance of Tests)

容許度需要設置為最大值:

否則,當long double比測試數據有更多的數字個數時,那麼沒有一個基於epsilon的容許度可以通過(Otherwise when long double has more digits than the test data, then no amount of tweaking an epsilon based tolerance will work)。

一個常見的問題是,當對於類似於 Microsoft VS.NET 的實現,容許度是合適的:double 和 long double 是同樣的大小,在其它的long double比double 更精確的系統上,容許度測試會失敗。首先檢查是否有L後綴,然後檢查容許度是否太高。

處理不合適的參數

數學特殊函數中的錯誤, J. Marraffino & M. Paterno 提議:當參數使得計算結果是一個數學未定義值的時候,產生一個定義域錯誤是強制性的。

一個函數被稱作是在 a = (a1, a2, . . .) 處定義的,如果 x = (x1, x2, . . .)的極限在所有的方向的逼近值都是相同的。定義的值可能是任何值 +infinity,或 -infinity。

如果函數函數定義延伸到 + infinity 然後 '折回' 到 - infinity,那麼函數是未定義的。

這個庫中的函數在數學未定義的值上計算也會產生一個錯誤。

這個庫中的函數在數學函數計算出的值不是一個實數值的時候也會產生一個錯誤。

這些實現是遵循這些提議的並與ISO/IEC 9899:1999 Programming languages - C and with the Draft Technical Report on C++ Library Extensions, 2005-06-24, section 5.2.1, paragraph 5. See also domain_error.相兼容的。

參考策略指南 瞭解更多錯誤處理策略,應當允許用戶遵從這些建議以及其它的錯誤處理行為的信息。

參考錯誤處理 瞭解這些機制的一個全面的解釋,以及錯誤處理實例error_handling_example.cpp

[Caution] 注意

如果你允許拋出異常但是沒有使用 try & catch ,那麼程序將會以一個未捕獲的異常終止並可能會退出。因此為了從報錯消息中得到幫助,推薦對於所有的應用程序允許拋出所有的異常並使用 try&catch 。然而,為了簡化代碼,大多數的例子程序都沒有這樣做。

處理數學上沒有定義的函數

數學上沒有定義的函數,類似於柯西均值,在缺省情況下不能通過編譯。一個 策略 可以用來控制這些行為。

如果這個函數允許未定義的行為,在缺省情況下,調用這些函數將會拋出定義域錯誤的異常。但是錯誤處理策略可以設置為不拋出異常,取而代之的是返回 一個NaN值。例如,在任何的Boost頭文件#include之前:

#define BOOST_MATH_DOMAIN_ERROR_POLICY ignore_error

那麼如果調用未實現的函數,函數 mean(cauchy<>()) 將返回 std::numeric_limits<T>::quiet_NaN()。

[Warning] 警告

如果std::numeric_limits<T>::has_quiet_NaN 為 false (例如T是一個用戶自定義類型),那麼當定義域錯誤出現時總會拋出一個異常。因此也強烈推薦捕獲異常。

分佈的均值(Median of distributions)

對於許多的分佈類型,我們不能夠找到一個分析方程(analytic formula),這就使得我們不能為這些分佈實現均值函數(median functions),一系列值的中間值。

然而,對於分佈 dist 的有用的均值逼近也許是可用的:

quantile(dist, 0.5).

Mean, Median, 以及 Skew, Paul T von Hippel

Descriptive Statistics,

以及

Mathematica Basic Statistics. 給出了更多信息,尤其是離散分佈。

處理浮點無限值(Handling of Floating-Point Infinity)

有一些函數和分佈可以將正無窮和負無窮作為參數,但是在將無限值作為特殊測試情況的實驗之後 ,我們推斷禁止正無窮和負無窮作為參數會更有用一些,取而代之的是返回定義域錯誤

將無限值作為特殊情況增添了額外的複雜性,因為不像內建的類型-但不是所有-平台-並不是所有的用戶自定義類型都提供了對std::numeric_limits<RealType>::infinity() 的特化並返回是無限值的任何表示。

原理在於,非有限的值可能因為溢出錯誤而出現在代碼中,並且將這個錯誤情況提示出來而不是仍舊繼續下去可能更有用一些。代碼也變得更複雜,更易於出錯,需要更多的測試,可讀性也更差。

然而,在一些情況下,例如正態分佈,我們很明顯地感覺到,我們應當允許參數為無限值,假定在那個實現上RealType實現了無限值的表示方法。

需要對無限值(或其它特殊值)進行特殊處理可以在調用函數和分佈之前攔截這些值並返回他們設置的值,或他們選擇的行為。這可能要比處理異常更簡單一些。

向下溢出,向上溢出,denorm 可以通過使用錯誤處理策略來處理。

我們也嘗試了捕獲邊界情況:產生除數為0或向下溢出並產生一個異常。在pole處的函數的行為可以由錯誤處理策略來控制。

Scale, Shape and Location

我們考慮把 位置(location) 和尺度(scale)添加到函數列表中,例如:

template <class RealType>
inline RealType scale(const triangular_distribution<RealType>& dist)
{
  RealType lower = dist.lower();
  RealType mode = dist.mode();
  RealType upper = dist.upper();
  RealType result;  // of checks.
  if(false == detail::check_triangular(BOOST_CURRENT_FUNCTION, lower, mode, upper, &result))
  {
    return result;
  }
  return (upper - lower);
}

但是發現這些概率對於大多數的分佈來說是未定義的(或者它們的定義太具有爭論性)。因為它們不是成員函數,所以,它們可以依據需要來添加。

特殊函數和分佈的實現說明(Notes on Implementation of Specific Functions & Distributions)
使用的有理逼近(Rational Approximations Used)

這個庫中的一些特殊函數是使用有理逼近來實現的。它們或者來自於文獻資料,或者使用由John Maddock 發明的Remez 代碼

使用有理逼近而不是多項式逼近來確保精度:多項式逼近通常可以很好地達到一個特定的精度,但是不管增加多少項,都不能達到更高的精度。

我們自己的逼近方法是出於增加精度(為了支持128bit 的 long double),或者是因為文獻中沒有可用的逼近方法或者基於 non-BSL compatible license而發明。我們的Remez 代碼與文獻資料十分相符,結果是相當「toy」的情況。所有的逼近針對收斂性進行了檢查並確保它們不是ill-conditioned(這些係數可以給出一個理論上好的解答,但導致有理函數在某些固定精度無法計算)。

使用不同的Remez實現來重新計算會產生不同的係數:在通常情況下,這個問題是 ill conditioned ,並且,對於其它的許多逼近,我們的Remez實現找出了一個寬的並且不是很清楚的極值(當然,對於類似於逼近函數exp的簡單例子,極小值是定義良好的,不管使用誰的實現方法,係數都是一樣的)。通常這不會影響逼近的合理性:有很好的文獻資料支持這種思想。注意:"in error"在這裡有一個特殊的含義,參考"Approximate construction of rational approximations and the effect of error autocorrection.", Grigori Litvinov, eprint arXiv:math/0101042。因此係數仍然需要精確地計算 ,即使與「真正的」極小值相比是錯的。

數學常量的表示(Representation of Mathematical Constants)

在 constants.hpp中的宏BOOST_DEFINE_MATH_CONSTANT用於為數學函數和分佈提供高精度的數值常量,因為,對內建的 float, double 和 long double 類型,以及對於用戶自定義的類型,類似於 NTL::quad_float 和 NTL::RR.提供一致的值是很重要的。

為了允許在數學工具包以及測試程序中使用精確度為100個十進制數字的NTL::RR類型,明顯需要將常量定義到這個精度。

然而,一些編譯器並不接受這樣長的十進制數字串。因此常量被分為三部分,第一部分包含至少是long double的精度,如果是已知的或是不需要的,則第二部分為0,如果需要的話,第三部分允許提供一個指數(如果沒有就使用0)-另兩個參數可能只包含十進制數字(以及正負號和小數點),並不包含一個類似於1.234E99的指數。僅當T是一個用戶定義的類型的時候才使用第二個十進制字符串,當這個常量被轉化為一個字符串或使用lexical_cast 轉換為類型T。 (這是必須的,因為你不能使用一個數值常量,即使是long double也沒有足夠的數字用於存儲這麼高的精度)。

例如, pi 定義為:

BOOST_DEFINE_MATH_CONSTANT(pi,
  3.141592653589793238462643383279502884197169399375105820974944,
  5923078164062862089986280348253421170679821480865132823066470938446095505,
  0)                                              

並這樣使用:

using namespace boost::math::constants;

double diameter = 1.;
double radius = diameter * pi<double>();

or boost::math::constants::pi<NTL::RR>()

注意,明確地指定類型是必須的:

所以,你不能這樣寫

double p = boost::math::constants::pi<>();  // 不能為模板推導參數「T」

你也不能這樣寫

double p = boost::math::constants::pi; // Context does not allow for disambiguation of overloaded function     
double p = boost::math::constants::pi(); // Context does not allow for disambiguation of overloaded function     
線程安全性(Thread safety)

通過設置errno來報告錯誤已經是安全的(否則沒有一個標準庫數學函數是線程安全的?)。如果你通過拋出異常來報告錯誤,errno就不再使用。

除那之外,這些代碼是準備對於內建數據類型 是線程安全的:因此float, double 和 long double 都是線程安全的。

對於非內建的數據類型-例如, NTL::RR -初始化這個庫中的各種常量是潛在的非線程安全的 這是最令人不愉快的,但修改這一點是一個不小的挑戰。一些編譯器可能提供以線程安全的方式來初始化常量的方法(Commeau, 或者其它的?),如果是這種情況,那麼這個問題就解決了。這是下一個C++標準的討論熱點,所以,希望所有的編譯器都被要求在這裡做出同樣正確的事情。

測試數據(Sources of Test Data)

我們找到了大量的測試數據。如果這些數據與我們的測試的數據結果相符,那麼我們假定它們是「已知良好的(known good)」,僅在出現嚴重的不相符時我們才咨詢其它的參考資料。精度,實際上並且聲明,是變化非常廣泛的。僅有Wolfram Mathematica functions 提供高於C++double精度的數據(64-bit 浮點數),並作為到目前為止最值得信賴的數據源。

數據源的一個有用的索引:Web-oriented Teaching Resources in Probability and Statistics

Statlet: 是一個Javascript 應用程序,用於計算和為統計分佈繪圖,並提供最完整的統計分佈範圍:

Bernoulli, Binomial, discrete uniform, geometric, hypergeometric, negative binomial, Poisson, beta, Cauchy-Lorentz, chi-sequared, Erlang, exponential, extreme value, Fisher, gamma, Laplace, logistic, lognormal, normal, Parteo, Student's t, triangular, uniform, and Weibull.

它計算 pdf, cdf, survivor, log survivor, hazard, tail areas, & critical values for 5 tail values.

它也是為韋伯分佈找到的唯一的獨立資源;然而,不幸的是,在那些特殊函數所驗證的實現區域,它會產生的精度很差。

生成和操作方程(Creating and Managing the Equations)

這些數學方程的原始資源現在是MathML: 參看 libs/math/doc/sf_and_dist/equations/中的.mml文件。

這些方程大都使用類似於Mathcast的GUI編譯器進行編譯,請注意, Open Office 當前提供的方程編譯器會損壞這些文件,因此現在不應該使用這個編譯器。

使用SVGMath 來轉換為SVG格式,命令行為:

$for file in *.mml; do 
>/cygdrive/c/Python25/python.exe 'C:\download\open\SVGMath-0.3.1\math2svg.py' \
>>$file > $(basename $file .mml).svg
>done

注意,SVGMath 要求這些mml文件沒有被包裝在XHTML XML中-這些缺省由Mathcast 添加 - 一種繞道方式是拷貝已存在的mml文件,然後使用 Mathcast編輯:存在的文件應當被保留。在SVGMath所使用的XML解析器中有一個bug,這個工具的作者也注意到這一點了。

如果需要的話,這些XHTML包裝可以使用下面的方法來去掉:

cat filename | tr -d "\r\n" | sed -e 's/.*\(<math[^>]*>.*</math>\).*/\1/' > newfile

為SVGMath設置字體現在是一件非常棘手的事情,在WindowsXP系統上,JM的字體設置與與SVGMath 提供了樣本配置文件是相同 的,但是:

<!-- Double-struck -->
    <mathvariant name="double-struck" family="Mathematica7, Lucida Sans Unicode"/>

改為:

<!-- Double-struck -->
    <mathvariant name="double-struck" family="Lucida Sans Unicode"/>

注意,不像SVGMath提供的樣本配置文件,這裡沒有使用Mathematica 7 字體,因為這種字體缺少足夠的Unicode信息使得它可以與SVGMath或XEP一起使用。

同樣也請請注意,在這裡的SVG 文件也大部分是Windows特定的,因為它們使用多種Windows字體。

可以使用Batik 來從SVG文件生成PNG文件,命令行為:

java -jar 'C:\download\open\batik-1.7\batik-rasterizer.jar' -dpi 120 *.svg

或者使用Inkscape ,命令行為:

for file in *.svg; do 
  /cygdrive/c/progra~1/Inkscape/inkscape -d 120 -e $(cygpath -a -w $(basename $file .svg).png) $(cygpath -a -w $file); 
done

當前, Inkscape 似乎生成更漂亮的png圖片。

在當前目錄\math_toolkit\libs\math\doc\sf_and_dist下使用一個命令將PDF 文件生成到 \pdf\math.pdf,典型的做法是:

bjam -a pdf

注意,XEP 需要配置為使用和嵌入SVG方程使用的任何字體(如果需要的話,編輯由XEP安裝程序提供的樣本xep.xml文件)。

(使用bjam -a 生成的html 文件在math_toolkit\libs\math\doc\sf_and_dist\html\index.html ).

JM's XEP 的配置文件中增加了下面的配置部分:

<font-group xml:base="file:/C:/Windows/Fonts/" label="Windows TrueType" embed="true" subset="true"> 
      <font-family name="Arial">
        <font><font-data ttf="arial.ttf"/></font>
        <font style="oblique"><font-data ttf="ariali.ttf"/></font>
        <font weight="bold"><font-data ttf="arialbd.ttf"/></font>
        <font weight="bold" style="oblique"><font-data ttf="arialbi.ttf"/></font>
      </font-family>

      <font-family name="Times New Roman" ligatures="&#xFB01; &#xFB02;">
        <font><font-data ttf="times.ttf"/></font>
        <font style="italic"><font-data ttf="timesi.ttf"/></font>
        <font weight="bold"><font-data ttf="timesbd.ttf"/></font>
        <font weight="bold" style="italic"><font-data ttf="timesbi.ttf"/></font>
      </font-family>

      <font-family name="Courier New">
        <font><font-data ttf="cour.ttf"/></font>
        <font style="oblique"><font-data ttf="couri.ttf"/></font>
        <font weight="bold"><font-data ttf="courbd.ttf"/></font>
        <font weight="bold" style="oblique"><font-data ttf="courbi.ttf"/></font>
      </font-family>

      <font-family name="Tahoma" embed="true">
        <font><font-data ttf="tahoma.ttf"/></font>
        <font weight="bold"><font-data ttf="tahomabd.ttf"/></font>
      </font-family>

      <font-family name="Verdana" embed="true">
        <font><font-data ttf="verdana.ttf"/></font>
        <font style="oblique"><font-data ttf="verdanai.ttf"/></font>
        <font weight="bold"><font-data ttf="verdanab.ttf"/></font>
        <font weight="bold" style="oblique"><font-data ttf="verdanaz.ttf"/></font>
      </font-family>

      <font-family name="Palatino" embed="true" ligatures="&#xFB00; &#xFB01; &#xFB02; &#xFB03; &#xFB04;">
        <font><font-data ttf="pala.ttf"/></font>
        <font style="italic"><font-data ttf="palai.ttf"/></font>
        <font weight="bold"><font-data ttf="palab.ttf"/></font>
        <font weight="bold" style="italic"><font-data ttf="palabi.ttf"/></font>
      </font-family>
      
    <font-family name="Lucida Sans Unicode">
         <font><font-data ttf="lsansuni.ttf"/></font>
    </font-family>

PAB 不得不修改它的配置文件 ,因為 Lucida Sans Unicode 字體有一個不同的名字。如果你不是使用Windows,那麼修改是非常類似的。

XZ 使用 venerable Latex來生成他的方程, JM 使用mxlatex將其轉換為MathML,這個過程當前是不可靠的並需要一些手工的干預:因此 Latex 文件不是一個自動生成SVG版本的方程的可靠的方法。

使用定義在 math.qbk 中的方程模板來將這些方程嵌入到 quickbook 文件中。生成的 Docbook XML 文件類似於下面這樣:

<inlinemediaobject>
<imageobject role"html">
<imagedata fileref"../equations/myfile.png"></imagedata>
</imageobject>
<imageobject role"print">
<imagedata fileref"../equations/myfile.svg"></imagedata>
</imageobject>
</inlinemediaobject>

MathML 並沒有出現在 Docbook 輸出中,也沒有在生成的HTML文件中:這需要以後的調查研究。

生成圖像(Producing Graphs)

圖像以SVG格式生成,然後用與方程圖像轉換相同的方式將SBG文件轉換為PNG文件。

/libs/math/doc/sf_and_dist/graphs/dist_graphs.cpp 程序和 /libs/math/doc/sf_and_dist/graphssf_graphs.cpp 程序使用[@http:/code.google.com/soc/2007/boost/about.html Google Summer of Code 2007] project of Jacob Voytko (whose work so far is at .\boost-sandbox\SOC\2007\visualization)直接生成SVG文件。


PrevUpHomeNext