[devel] Минутка математики (Was: Re: [#238855] [test-only] FAILED (try 4) openblas.git=0.3.9-alt1 R-base.git=4.0.0-alt1)

[Alexey Sheplyakov]

Здравствуйте!

On Tue, Apr 28, 2020 at 10:53:03PM +0300, Kirill Maslinsky wrote:
 
> $(R_EXE) на ppc64le зависает не в специфическом месте при сборке, а
> вообще всегда при запуске. Похоже, это происходит вот в этом цикле:
> 
> http://git.altlinux.org/people/kirill/packages/R-base.git?p=R-base.git;a=blob;f=src/main/machar.c;h=8db54f1350fe85b0204c7affb4ef79180fc3e5ad;hb=edfd434c0f879f16a9a8f124648df4b4eb23a7e3#l156
> 
> src/main/machar.c
> 
> 155         for(;;) {
> 156                 temp = one - a;
> 157                 if (temp - one != zero)
> 158                         break;
> 159                 a = a * beta;
> 160                 *negep = *negep - 1;
> 161         }
> 
> Что бы это могло быть?

Исходя из базы (beta) и длины мантиссы (it) пытаются вычислить минимальное
число eps, такое, что 1.0 - eps != 1.0 (точнее, его логарифм по основанию
beta). Когда этот код писали, еще не было стандарта IEEE 754, <float.h>,
макросов {FLT,DBL,LDBL}_EPS. Приходилось вот так изворачиваться.

Но на ppc64 (и некоторых SPARCv9) вместо нормального (IEEE 754) long doulbe
почему-то используется double double, т.е. представление числа в виде пары
(суммы) двух double. При таком способе большИе числа можно представить
точнее, чем double. Например, число 2^53 1/2 можно точно (без округления)
представить в виде cуммы "большого" числа 2^53 и "маленького" 1/2. А в double
такое число не влезет, 53 бит не хватит, придется округлять. Причем диапазон
"маленького" числа гораздо шире, чем у целого той же разрядности (53 бита),
ведь есть 11-битный показатель, и им можно масштабировать! Но не выйдет
представить с большей (чем double) точностью положительное число, меньшее
2^{-1021}, по той же причине -- показатель все равно 11 бит. Таким образом,
у double double длина мантиссы не постоянна и зависит от показателя.
А функция MACH_NAME исходит из обратного, т.е. из (в общем-то разумного)
предположения, что длина мантисы фиксирована.


А теперь вернемся к 1 + eps != 1.0. С помощью "честного" числа с плавающей
точкой с длиной мантисы 106 бит (что соответствует паре double) можно
представить b1.0000{100 нулей}1, или 1 + 1/2^105. А с помощью double double
можно представить даже b1.0{1019 нулей}1 в виде пары (1, 1/2^(-1021)).
Потому попытка искать с 

149         *negep = *it + 3;

или -109, и в дальнейшем *увеличивать* показатель

155         for(;;) {
156                 temp = one - a;
157                 if (temp - one != zero)
158                         break;
159                 a = a * beta;
160                 *negep = *negep - 1;
161         }

заведомо обречена. Поэтому этот цикл никогда не закончится (ну, может
когда int со знаком переполнится).

К счастью, "подарок" от IBM (double double) легко отключть с помощью флагов
компиляции

-mabi=ieeelongdouble -mlong-double-128

Такие дела. Простите за много букв.

Всего доброго,
    Алексей