Eigen: Что такое вычислительное ограничение для log()? (Windows vs. Linux)

Question

У меня странная проблема с кодом Eigen. Он хорошо работает на Linux (64-разрядный) и Mac OSX (64-разрядный), но он не работает в Windows (32 бит) с 1.#INF значением.

Я отследил вниз проблему к этой функции, где я реализовал в Box-Мюллер-преобразование:

Eigen::MatrixXd box_muller (const Eigen::VectorXd vRand) 
{ 
    unsigned long n = vRand.rows(); 
    unsigned long m = n/2; 

    Eigen::ArrayXd rand1 = vRand.head (m); 
    Eigen::ArrayXd rand2 = vRand.tail (m); 

    /* Implemented according to 
    * http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform 
    */ 

    for (unsigned long i=0; i<rand1.cols(); i++) 
    { 
     if (rand1 (i) < 1e-8) 
     { 
      rand1 (i) = 1e-8; 
     } 
    } 

    rand1 = (-2 * rand1.log()).sqrt(); // something must be wrong here 

    std::cout << rand1.mean() << std::endl; // prints 1.#INF 

    rand2 = rand2*2*M_PI; 

    Eigen::MatrixXd result (2*m, 1); 
    Eigen::MatrixXd res1 = (rand1 *rand2.cos()).matrix(); 
    Eigen::MatrixXd res2 = (rand1 *rand2.sin()).matrix(); 
    result << res1, res2; 

    return result; 
} 

Этот фрагмент кода не удается, но только на Windows, и только для больших входных векторов vRand. Он не работает, потому что в rand1 есть непредвиденные значения, которые генерируют значения 1.#INF, -1.#INF, 1.#IND или 1.#QNAN в последующих вычислениях. Когда я устанавливаю количество элементов в vRand на небольшое число, скажем, 10000, он отлично работает. Но когда число велико (например, 100000), он терпит неудачу.

Я пробовал все, что мог придумать, и теперь у меня нет идей. Что я могу попытаться устранить эту проблему?

Подробнее:

• Мы строим это с MS Visual Studio 12.0. Он работает на сервере buildserver, поэтому я не могу легко найти точную версию.
• Мы используем C++ 0x, насколько я знаю.

Answer 1

Я предлагаю заменить этот цикл:

for (unsigned long i=0; i<rand1.cols(); i++) 
{ 
    if (rand1 (i) < 1e-8) 
    { 
     rand1 (i) = 1e-8; 
    } 
} 

с:

rand1 = rand1.max(1e-8); // limit min value to 1e-8 
rand1 = rand1.min(1.0); // limit max value to 1.0 

Это позволило бы гарантировать, что значения rand1.log() будет <= 0 и последующее sqrt() не должны затем потерпеть неудачу. Если это исправляет проблему, вы можете захотеть вернуться назад, чтобы определить, откуда поступают входные значения вне диапазона.