Преобразование с плавающей точкой в ​​десятичную с вычислениями с плавающей запятой

Я пытаюсь преобразовать значение двойной точности с плавающей запятой x в десятичное с 12 (правильно округленными) значащими цифрами. Я предполагаю, что x находится между 10 ^ 110 и 10 ^ 111, так что его десятичное представление будет иметь вид x.xxxxxxxxxxxE110. И, просто для удовольствия, я пытаюсь использовать только арифметику с плавающей запятой.

Я пришел к псевдокоду ниже, где все операции являются операциями с двойной точностью. Обозначение 1e98 - для ближайшего к математике 10 ^ 98, а 1e98_2 является двойным, ближайшим к результату математическое вычитание 10 ^ 98- 1e98. Обозначение fmadd(X * Y + Z) относится к операции с объединенным многократным добавлением с операндами x, Y, Z.

  y = x * 2^-1074;    // exact
  q = y / 1e98;       // q is denormal and the significand of q interpreted
                      // as an integer is our candidate for the 12 decimal
                      // digits of x

  r = fmadd(q * 1e98 - y);  // close to 1e98 * (error made during the division)

  // If 1e98_2 >= 0, we divided by a number that was smaller than we wished
  // The correct answer may be q or q+1.

  if (r and 1e98_2 have opposite signs)
  {
    return the significand of q;
  }

  s = copysign(2^-1074, r);
  r1 = abs(r);
  r2 = abs(1e98_2);

  h = 1e98 * 0.5 * 2^-1074;

  Set rounding mode to downwards

  r3 = fmadd(r2 * q + r1);

  if (r3 < h)
  {
    return the significand of q;
  }
  else
  {
    return significand of (q + s)
  }

Я извиняюсь за путаницу, которая пронизывает вышеупомянутый псевдокод, но мне пока это не совсем понятно, поэтому следующие вопросы:

  • Работает ли первый fmadd по назначению (для вычисления 1e98 * (ошибка, сделанная во время деления))?

  • Знаки. Я не могу убедить себя, что они правы. Но я не могу убедить себя, что они ошибаются.

  • Любая идея, возможно, аргументы о частоте, с которой этот алгоритм может произвести неверный результат?

  • Если это вообще работает, есть ли вероятность, что алгоритм продолжит работу, если "q = y/1e98" изменено на "q = y * 1e-98" (оставив все остальные инструкции одинаковыми )?

Я не тестировал этот алгоритм. У меня нет компьютера с инструкцией fmadd, хотя я надеюсь найти его, чтобы выполнить выше.

Ответы

Ответ 1

Пусть y/d - точная операция, а q=rnd(y/d) - результат, округленный до ближайшего поплавка.
Тогда истинная ошибка, умноженная на d, равна rt=(rnd(y/d)-y/d)*d=q*d-y, а операция, которую мы выполнили с fmadd, равна r=rnd(q*d-y)
Почему q*d-y является точным (fmadd не делает окончательного округления) менее понятно для объяснения, но скажем, что q*d имеет ограниченное количество бит (<nbits(q)+nbits(d)), показатель степени y равен показателю q*d ( +/- 1), а поскольку ошибка |rt|<0.5*ulp(q)*d, это означает, что первые nbits(q) исчезают... Это отвечает на вопрос 1.

So q*1e98 - y = r, где |r|*2^1074 <= 0.5e98 < 5*10^98 (второе неравенство повезло)

q*(10^98) - y = r + (10^98-1e98)*q где |10^98-1e98|*q*2^1074 <= 0.5e95 (предполагая точность не менее 15 цифр, log(2^53)/log(10) > 15)

Итак, вы спрашиваете, |q*(10^98)-y|*2^1074>5*10^97

У вас есть приближение |q*(10^98)-y|, которое r+1e98_2*q

Так как |r| < 5*10^98 и |r+(10^98-1e98)*q|<|r|, если знаки противоположны, я думаю, что это положительно отвечает на вопрос 2. Но я не был бы так уверен, если бы 1e98_2 были < 0.

Если r и 1e98_2 имеют один и тот же знак, он может превышать 5*10^97, поэтому ваша дальнейшая обработка с обсуждением r3 = 1e98_2*q + r по сравнению с h=0.5e98*2^-1074

Для вопроса 3, на первый взгляд, я бы сказал, что две вещи могут привести к сбою алгоритма:

  • 1e98_2 не является точным (10^98-1e98-1e98_2 = -3.6e63 прибл.)

  • и h не ht=0.5*10^98*2^-1074, но немного меньше, как мы видели выше.

Истинная ошибка r3t примерно равна (1e98_2-3e63)*q + r < r3 (и нам интересен только случай, когдa > 0, так как 1e98_2 > 0).

Таким образом, аппроксимация ошибки r3, падающей выше аппроксимированной связи h, когда истинная ошибка r3t ниже истинной связи ht, может привести к неправильному округлению. Возможно ли, и если да, как часто возникает ваш вопрос 3?

Чтобы смягчить выше риск неравенства, вы попытались усечь величину r3, таким образом r3 <= 1e98_2*q + r. Я чувствовал себя немного уставшим, чтобы выполнить истинный анализ границ ошибок...

Итак, я просмотрел ошибку, и первый неудачный пример, который я нашел, был 1.0000000001835e110 (я предполагаю, что правильно округлено до ближайшего двойника, но на самом деле это 1000000000183.49999984153799821120915424942630528225695526491963291846957919215885146546696544423465444842668032e98).

В этом случае r и 1e98_2 имеют один и тот же знак, а

  • (x/1e98) > 1000000000183.50000215

  • q Значение таким образом округляется до 1000000000184

  • r3>h (r3*2^1074 составляет около 5.000001584620017e97), и мы неправильно увеличили q+s, когда он должен был q-s, определенно ошибка.

Мои ответы:

  • да, r=fmadd(q * 1e98 - y) составляет ровно 1e98 * (ошибка во время деления), но мы не заботимся о делении, это просто дает предположение, что считается, что вычитание является точным.

  • да, знак правильный, потому что |r| < 5*10^98 и |r+(10^98-1e98)*q|<|r|, если знаки противоположны. Но я не был бы так уверен, если бы 1e98_2 были < 0.

  • Взяв первый неудачный пример (1.0000000001835e110 - 1.0e110)/1.0e110 ulp -> 1.099632e6, очень наивная гипотеза заключалась бы в том, чтобы сказать, что 1 случай из миллиона, r3 падает на h... Так как раз q + s исправлено в qs, возникновение r3>h, в то время как r3t<ht намного меньше 1/1 000 000 в любом случае... в интересующем диапазоне больше 10 ^ 15 удвоений, так что считайте, что это не серьезный ответ...

  • Да, обсуждение выше касается исключительно предположения q, независимо от способа его создания, и вычитание в 1. будет по-прежнему точным...