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

[Alexey Sheplyakov <asheplyakov@basealt.ru>]

On Wed, Apr 29, 2020 at 04:56:16PM +0400, Alexey Sheplyakov 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)).

Имеется в виду пара (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
>