捨入策略

為了盡可能的通用, 這個庫使用一個類來完成函數所需要的向上捨入和向下捨入功能. 這個類是策略類policies的第一個參數, 它同樣也是區間類interval的策略定義中名為rounding的類型.

默認情況下,這個捨入策略是interval_lib::rounded_math<T>. interval_lib::rounded_math已經針對標準的浮點類型 (float , double 和 long double)進行了特化. 所以，如果你的區間模板類的類型不是這三種之一，一種好的解決方案是提供一個該類的特化類. 但是如果 默認的針對 float, double, 或 long double 的特化不是你想要的, 或者你不想提供一個interval_lib::rounded_math<T>特化類(比如說你想在自己的名字空間中進行工作)，你同樣可以定義你自己的捨入策略並將其直接傳遞給 interval_lib::policies.

要求

下面是這個庫所提供的內容，參數範圍在對應的函數後面(如同你所看到的，這些函數不用擔心非法的參數，但是它們必須處理參數是無限的情況).

/*捨入要求*/
struct rounding {
  // 缺省構造函數, 析構函數
  rounding();
  ~rounding();
  // 數學操作
  T add_down(T, T); // [-∞;+∞][-∞;+∞]
  T add_up  (T, T); // [-∞;+∞][-∞;+∞]
  T sub_down(T, T); // [-∞;+∞][-∞;+∞]
  T sub_up  (T, T); // [-∞;+∞][-∞;+∞]
  T mul_down(T, T); // [-∞;+∞][-∞;+∞]
  T mul_up  (T, T); // [-∞;+∞][-∞;+∞]
  T div_down(T, T); // [-∞;+∞]([-∞;+∞]-{0})
  T div_up  (T, T); // [-∞;+∞]([-∞;+∞]-{0})
  T sqrt_down(T);   // ]0;+∞]
  T sqrt_up  (T);   // ]0;+∞]
  T exp_down(T);    // [-∞;+∞]
  T exp_up  (T);    // [-∞;+∞]
  T log_down(T);    // ]0;+∞]
  T log_up  (T);    // ]0;+∞]
  T cos_down(T);    // [0;2π]
  T cos_up  (T);    // [0;2π]
  T tan_down(T);    // ]-π/2;π/2[
  T tan_up  (T);    // ]-π/2;π/2[
  T asin_down(T);   // [-1;1]
  T asin_up  (T);   // [-1;1]
  T acos_down(T);   // [-1;1]
  T acos_up  (T);   // [-1;1]
  T atan_down(T);   // [-∞;+∞]
  T atan_up  (T);   // [-∞;+∞]
  T sinh_down(T);   // [-∞;+∞]
  T sinh_up  (T);   // [-∞;+∞]
  T cosh_down(T);   // [-∞;+∞]
  T cosh_up  (T);   // [-∞;+∞]
  T tanh_down(T);   // [-∞;+∞]
  T tanh_up  (T);   // [-∞;+∞]
  T asinh_down(T);  // [-∞;+∞]
  T asinh_up  (T);  // [-∞;+∞]
  T acosh_down(T);  // [1;+∞]
  T acosh_up  (T);  // [1;+∞]
  T atanh_down(T);  // [-1;1]
  T atanh_up  (T);  // [-1;1] 
  T median(T, T);   // [-∞;+∞][-∞;+∞]
  T int_down(T);    // [-∞;+∞]
  T int_up  (T);    // [-∞;+∞]
  // 轉換函數
  T conv_down(U);
  T conv_up  (U);
  // 不受保護的捨入類
  typedef ... unprotected_rounding;
};

捨入類的構造函數和析構函數有非常重要的語義要求: 它們負責設置和重置在T類型上運算的捨入方式. 例如，如果T是一個標準的浮點類型，並且浮點運算是依照 IEEE 754標準,構造函數可以保存當前的捨入狀態, each _up (_down) 函數將會向上捨入(向下捨入),析構函數將會 重置已保存的捨入狀態.實際上這是默認的捨入策略。

所有的數學函數的意義都是很清楚的直到函數atanh_up:每個函數返回在類型T中可以表示的相應數學函數運算結果真值的下邊界的值(對於 _down) 或者上邊界值(對於 _up).函數median 計算兩個參數的平均值，並將結果捨入到最鄰近的可以表示的數值. 函數int_down 和 int_up 計算小於或者大於其參數的最小或最大整數. 最後, conv_down 和 conv_up 負責將其它類型的數值轉換到基類型數值: 第一個函數將數值向下捨入，第二個函數將數值向上捨入.

unprotected_rounding 允許去掉所有的捨入保護. 至於為什麼有人會這麼做, 參看下面保護章節
.

已提供的類概覽

提供了許多類. 這此類依據層次來組織.底部是rounding_control類. 接下來是 rounded_arith_exact, rounded_arith_std 和 rounded_arith_opp類. 然後是rounded_transc_dummy, rounded_transc_exact, rounded_transc_std 和 rounded_transc_opp類. 最後是save_state 和 save_state_nothing類.每一個類都提供了該類所在層次的下一個層次中的類所要求的函數. 例如 , 一個rounded_transc_... 類需一個rounded_arith_...類的成員.

已提供的捨入策略類有_std和_opp兩種版本,第一個版本每次都會調整捨入方式, 而第二種版本試圖保持向正無窮捨入. _opp 的主要目的是通過使用"取反技巧"(opposite trick)來加速計算(參見性能提示 ). 這一版本要求在進入類的任何運算函數之前捨入方式應當是向上捨入，它保證在函數退出時，捨入方式仍然是向上捨入。

請注意：混合使用_opp和_std是一個非常壞的主意，因為它們沒有兼容的性質.

第三個版本稱作_exact，在計算時不會改變捨入方式. 它是一個精確的版本，因為它針對於產生精確結果的基類型.

最後一個版本是_dummy . 它不進行任何的計算，但仍舊產生可兼容的結果

請注意：可以對一個不精確的基類型使用「精確」版本的捨入策略，例如 . float 或 double. 在這種情況下, 包含性質不再得到保證, 當包含性質沒有受到嚴格要求的時候，這種方式可以用來加速計算.例如，在計算機圖形學中，由於浮點捨入而產生小的計算誤差是可以接受的.

下面是每個類所定義的內容，稍後，將會有更全面的描述，這些成員不會被重複，請轉到這裡來參看這些定義，繼承也被用來避免重複。

template <class T>
struct rounding_control
{
  typedef ... rounding_mode;
  void set_rounding_mode(rounding_mode);
  void get_rounding_mode(rounding_mode&);
  void downward ();
  void upward   ();
  void to_nearest();
  T to_int(T);
  T force_rounding(T);
};

template <class T, class Rounding>
struct rounded_arith_... : Rounding
{
  void init();
  T add_down(T, T);
  T add_up  (T, T);
  T sub_down(T, T);
  T sub_up  (T, T);
  T mul_down(T, T);
  T mul_up  (T, T);
  T div_down(T, T);
  T div_up  (T, T);
  T sqrt_down(T);
  T sqrt_up  (T);
  T median(T, T);
  T int_down(T);
  T int_up  (T);
};

template <class T, class Rounding>
struct rounded_transc_... : Rounding
{
  T exp_down(T);
  T exp_up  (T);
  T log_down(T);
  T log_up  (T);
  T cos_down(T);
  T cos_up  (T);
  T tan_down(T);
  T tan_up  (T);
  T asin_down(T);
  T asin_up  (T);
  T acos_down(T);
  T acos_up  (T);
  T atan_down(T);
  T atan_up  (T);
  T sinh_down(T);
  T sinh_up  (T);
  T cosh_down(T);
  T cosh_up  (T);
  T tanh_down(T);
  T tanh_up  (T);
  T asinh_down(T);
  T asinh_up  (T);
  T acosh_down(T);
  T acosh_up  (T);
  T atanh_down(T);
  T atanh_up  (T);
};

template <class Rounding>
struct save_state_... : Rounding
{
  save_state_...();
  ~save_state_...();
  typedef ... unprotected_rounding;
};

概要

namespace boost {
namespace numeric {
namespace interval_lib {
/*基本的捨入控制*/
template <class T>  struct rounding_control;

/* 算術函數捨入 */
template <class T, class Rounding = rounding_control<T> > struct rounded_arith_exact;
template <class T, class Rounding = rounding_control<T> > struct rounded_arith_std;
template <class T, class Rounding = rounding_control<T> > struct rounded_arith_opp;

/* 超越函數捨入 */
template <class T, class Rounding> struct rounded_transc_dummy;
template <class T, class Rounding = rounded_arith_exact<T> > struct rounded_transc_exact;
template <class T, class Rounding = rounded_arith_std<T> > struct rounded_transc_std;
template <class T, class Rounding = rounded_arith_opp<T> > struct rounded_transc_opp;

/* 捨入狀態保存類 */
template <class Rounding> struct save_state;
template <class Rounding> struct save_state_nothing;

/* 針對於類型T的缺省策略 */
template <class T>  struct rounded_math;
template <>  struct rounded_math<float>;
template <>  struct rounded_math<double>;

/* 一些用於從未保護方式到保護方式轉換的元編程代碼 */
template <class I> struct unprotect;

} // namespace interval_lib
} // namespace numeric
} // namespace boost

已提供的類的描述

現在我們以它們出現在捨入策略定義中的順序來描述每一個類 (從最外到最內的順序是概要中順序的反序).

保護控制

保護基於這樣一個事實：區間操作將被捨入方式控制所包圍. 對一個類去保護意味著去掉所有的捨入控制，每一種捨入策略都提供了 unprotected_rounding類型. 當嵌套在另一個捨入類中的時候，要求的 unprotected_rounding類型使得另一個捨入類可以工作. 例如，下面最開始的三行將產生同樣的結果 (因為第一個操作是捨入構造函數，而最後一個是析構函數，確保設置捨入方式); 並且最後一行允許有未定義的行為(因為沒有調用構造函數或析構函數).

T c; { rounding rnd; c = rnd.add_down(a, b); }
T c; { rounding rnd1; { rounding rnd2; c = rnd2.add_down(a, b); } }
T c; { rounding rnd1; { rounding::unprotected_rounding rnd2; c = rnd2.add_down(a, b); } }
T d; { rounding::unprotected_rounding rnd; d = rnd.add_down(a, b); }

自然地， rounding::unprotected_rounding 可以僅僅是 rounding 本身.但是如果它是一個構造函數和析構函數為空的簡化版本，它就可以提高性能. 為了避免未定義的行為,在這個庫中, 僅當一個rounding類型的對象已經存在時，rounding::unprotected_rounding 類型的對象才會被創建，參見 性能提示瞭解其它的信息.

支持庫定義了帶有一個模板參數為區間類型I的元編程類模板unprotect，它返回一個已經去保護的區間類型unprotect <I>::type.關於類型的一些信息: interval<T, interval_lib::policies<Rounding, _> >::traits_type::rounding與Rounding是同樣的類型, unprotect<interval<T, interval_lib::policies<Rounding, _> > >::type 與interval<T, interval_lib::policies<Rounding::unprotected, _> >是同樣的類型.

狀態保存

首先是save_state類. 這個類負責保存當前的捨入方式並在構造函數中調用函數init,在析構函數中重置已保存的捨入方式. 這個類也定義了unprotected_rounding 類型 .

如果捨入方式不需要任何的狀態保存或者初始化, save_state_nothing 可以用來代替save_state.

超越函數

rounded_transc_exact, rounded_transc_std 和 rounded_transc_opp 期望名字空間std提供函數exp log cos tan acos asin atan cosh sinh tanh acosh asinh atanh. 對於_std和_opp版本, 所有的這些函數應當對調用函數downward或upward所改變的當前捨入方式作出反應 .

請注意: 不幸的是, 後者(_opp)極少是這種情況. 這就是為什麼提供了一個不依賴於名字空間std中函數的類rounded_transc_dummy . 當然，這是沒有什麼神奇的地方. rounded_transc_dummy 類的成員函數並不進行任何的計算. 它們僅返回合法的數值. 例如, cos_down 永遠返回-1，在這種情況下，我們為默認實現確保包含性質, 即使這些函數完全沒有可以返回用戶的返回值.為了獲得有用的值，應當明確地另一個策略, 它極有可能導致違反包含性質.在這種情況下,我們確保向用戶清楚地指明這些違例，用戶也就知道了他違反了什麼， 這個類或許應當被用作默認的超越捨入類, 但是我們決定由於缺少聲明而編譯失敗將會比那些合法的但是無用的函數更有意義一些。

基本的算術函數

rounded_arith_std類和rounded_arith_opp類期望運算符 + - * /和函數std::sqrt顧及當前的捨入方式.

rounded_arith_exact類要求定義函數std::floor 和 std::ceil，因為它不能依賴於to_int

捨入控制

前面的每個類所定義的函數不需要作任何的解釋.例如, add_down 的作用是計算兩個數的和並將結果向下捨入. 對於 rounding_control, 情況稍有點複雜.

基礎函數是force_rounding，如果它的參數不是依據於當前的捨入方式進行捨入，那麼它返回依據於當前捨入方式進行正確捨入的函數參數.這個函數必須處理延遲捨入，實際上， 依賴於計算所完成的方式, 運算的中間結果可能在內部被保存為更精確的形式，並且可能導致錯誤的捨入. 這裡 是一個當捨入沒有強制實施時所發生的情況的例子。

函數get_rounding_mode返回當前的捨入方式 , 函數set_rounding_mode將捨入方式設置為由函數get_rounding_mode返回的捨入方式. 函數downward, upward 和to_nearest將捨入方式設置為三種捨入方向中的一種. 這種捨入方式對於所有的使用類型T的函數應當是全局的，例如 ,在調用函數downward之後, 函數force_rounding(x+y)返回的和是向 -∞捨入的.

函數to_int依據當前的捨入方式來計算最鄰近的整數值.

非特化版本的rounding_control 不做任何事情. 對應的函數實現是空的，並且函數to_int和函數force_rounding 是鑒定函數 . 函數pi_ constant返回合適的整數(例如, pi_up返回T(4)).

類模板rounding_control針對於 float, double 和 long double進行了特化，為的是更好的使用計算機浮點運算單元.

模板類rounded_math

默認的策略(也稱作 rounded_math<T>)簡單地定義如下:

template <class T> struct rounded_math<T> : save_state_nothing<rounded_arith_exact<T> > {};

並且針對於float, double 和 long double進行了特化，使用rounded_arith_opp如下:

template <> struct rounded_math<float>       : save_state<rounded_arith_opp<float> >       {};
template <> struct rounded_math<double>      : save_state<rounded_arith_opp<double> >      {};
template <> struct rounded_math<long double> : save_state<rounded_arith_opp<long double> > {};

性能事項

這一章主要處理使用計算機的浮點運算單元時類庫的性能，讓我們將[a,b] 與 [c,d] 的和作為一例子，結果是 [down(a+c), up(b+d)], 其中 down 和up 表明需要使用捨入方方式 .

捨入方式轉換

如果FPU(計算機浮點運算單元)可以針對每一個操作使用不同的捨入方式，這裡就不會有什麼問題了，例如，對於Alpha處理器是這種情況: 每一個浮點指令可以指定一個不同的捨入方式，然而，IEEE-754 標準並沒有要求這種行為，所以，大多數的FPU僅僅提供一些指令用於設置接下來的所有的浮點運算的捨入方式，通常情況下，這些指令需要刷新FPU的指令流水線

在這種情況下，計算 [a,b] 與 [c,d]之和的時間遠比計算a+b 和 c+d 的時間要多。因為這兩個加法不能並行執行，因此，客觀的方法是減少捨入方式調整的數量。

如果這個庫不是用來提供精確的計算，而僅僅是針對成對算術( pair arithmetic)，那麼解決方案就相當簡單:不使用捨入，在這各情況下，計算[a,b]與[c,d]的和將會與計算a+b和c+d一樣快 . 這樣，每一件事都很完美。

然而，如果要求精確的計算，這樣一種解決方案是完全不予考慮的，難道我們無計可施(penniless)了嗎，不，這裡仍然有一個技巧，實際上，如果一元取負運算是一種精確的運算，那麼down(a+c) = -up(-a-c)， 現在就可以使用同樣的捨入方式來計算整個和，通常，捨入方式改變的代價遠比符號改變的代價要大。

加速連續運算

區間加法並不是唯一的運算，大多數的區間運算可以只設置一次FPU的捨入方式而得以完成，因此，浮點捨入策略的操作假定捨入方向被正確的設置，在一個程序中這種假設通常是不成立的 (用戶和標準庫期望捨入方式是向最近的值捨入), 所以，這些操作應當被包含在一個用於設置浮點運算環境的外殼中，這種保護由捨入策略的構造函數和析構函數來完成，

現在讓我們考慮兩個連續的區間加法的情況[a,b] + [c,d] + [e,f]， 生成的代碼應當是這樣的:

init_rounding_mode();    // 在第一次加法中構造捨入對像
t1 = -(-a - c);
t2 = b + d;
restore_rounding_mode(); // 捨入對像析構
init_rounding_mode();    // 在第二次加法中構造捨入對像
x = -(-t1 - e);
y = t2 + f;
restore_rounding_mode(); // 捨入對像析構
// 結果為區間 [x,y]

在兩步操作之間，捨入方式被重置了，然後再一次被初始化，理想情況下，編譯器應當可以優化掉這些無用的代碼，但不幸的是，它們無法做到這一點，並且這會減慢代碼的執行速度 . 為了避免這些瓶頸， 用戶可以告訴區間運算他們不再需要捨入保護，然後就會由用戶來保護區間運算，編譯器將能夠生成下面的代碼 :

init_rounding_mode();    // 用戶完成
x = -(-a - c - e);
y = b + d + f;
restore_rounding_mode(); // 用戶完成

用戶需要生成一個捨入對像 ，只要這個捨入對像依然存在，就可以繼續使用區間運算未受保護的版本，通過使用特定捨入策略的區間類型來選定，如果初始的區間類型是I，那麼，I::traits_type::rounding就是捨入對象的類型，並且interval_lib::unprotect<I>::type是未受保護的區間類型.

因為FPU的捨入方式在捨入對象的生存期內被改變，任何不涉及這個區間庫的算術浮點運算將會導致不可預期的結果，相應的，當沒有捨入對像存在而使用未受保護的區間運算將不保證產生真正的結果區間。

例子

下面是一個使用霍納法則(Horner's scheme)來計算多項式值的例子. 在整個計算過程中捨入方式的改變都被禁止

// I 是一個區間類， 多項式是一個簡單的數組

template<class I>
I horner(const I& x, const I p[], int n) {

  // 保存和初始化捨入方式
  typename I::traits_type::rounding rnd;

  // 定義區間類型的未受保護版本
  typedef typename boost::numeric::interval_lib::unprotect<I>::type R;

  const R& a = x;
  R y = p[n - 1];
  for(int i = n - 2; i >= 0; i--) {
    y = y * a + (const R&)(p[i]);
  }
  return y;

  // 由 rnd 的析構函數來重置捨入方式 
}

請注意：為了保護所有的區間計算，一個捨入對像被特地構造出來。在任何操作之前，每一個I類型的區間在一個R類型的區間中被轉化， 如果這種轉化沒有完成，結果仍然是正確的，但整個性能優化上的獲益都沒有了，盡可能的情況下，轉換為const R&而不是R會更好一些。實際上，這個函數可能已經在一個未保護的代碼塊中被調用，那麼類型R和I將會是相同的區間,因而也就沒有轉換的必要了。

不那麼有趣的備註(Uninteresting remark)

在最開始的時候提到過，Alpha處理器可以針對每個操作使用不同的捨入方式， 然而，由於指令的形式，向正無窮的捨入是不可用的，只有向負無窮的捨入是可用的，所以，最基本的技巧是改變符號，但是沒有必要在一個操作的兩個方面保存和恢復捨入方式。

擴展寄存器

除了捨入方式切換的開銷之外，還有另外一個問題,一些FPU使用擴展的寄存器(例如 , 浮點運算(float computations)將會使用雙精度浮點寄存器(double registers), 或者雙精度浮點運算(double computations)使用長雙精度浮點寄存器(long double registers)來完成). 這必然會產生很多的問題。

第一個問題是由於尾數的擴充精度，捨入同樣在這個擴充精度上完成，必然地， 在擴展寄存器中我們仍然可以得到:down(a+b) = -up(-a-b)，但回到標準精度，現在我們得到: down(a+b) < -up(-a-b)而不是相等. 一種解決方案是不要使用這種方法，但仍然還有其它的問題，例如，數字之間的比較.

自然的，對於指數的擴展精度仍然有一個問題，為了舉例說明這個問題，假定 m 為+inf之前最大的數，如果我們計算 2*[m,m], 結果將會是 [m,inf]. 但是由於擴展寄存器，FPU將會首先存儲 [2m,2m]，然後在運算的最後將[2m,2m]轉化為 [inf,inf](當捨入方式向正無窮捨入時). 所以，結果不再是精確的.

這裡僅有一個解決方案迫使FPU在每個計算之後將擴展精度值轉換為標準精度值，一些FPU提供了一個指令來完成這種轉化(例如PowerPC處理器). 但是對於沒有提供這種指令的處理器 (x86 處理器), 唯一的方法就是將這些值先寫入內存中，然後再將它們從內存中取回來，很明顯，這種操作的開銷很大。

例子

下面列舉了一些情況：

• 如果你需要進行float或double類型的精確運算，使用默認的 rounded_math<T>;
• 對於快速而廣泛的且不需要任何捨入和精度的區間運算使用 save_state_nothing<rounded_transc_exact<T> >;
• 對於不支持超越函數的精確類型(像int或是對無限(infinite )和NaN提供處理方法的有理數) 解決方案為 save_state_nothing<rounded_transc_dummy<T, rounded_arith_exact<T> > > 或都直接使用 save_state_nothing<rounded_arith_exact<T> >;
• 如果比前面的情況更複雜，最好的方法或許是直接定義你自己的策略類.

Revised 2006-12-24

Copyright © 2002 Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, Polytechnic University
Copyright © 2004-2005 Guillaume Melquiond, ENS Lyon

Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)