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

Я пытаюсь преобразовать значение двойной точности с плавающей точкой 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)
  }

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

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

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

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

  4. Если он вообще работает, есть ли шанс, что алгоритм продолжит работать, если "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.

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

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 (и интересен только случай, когда>0, потому что 1e98_2>0).

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

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

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

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

  • (x/1e98) > 1000000000183.50000215

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

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

Мои ответы:

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

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

  3. Принимая первый неудачный пример (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 двойников, так что посчитайте это несерьезным ответом...

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

Источник
Другие вопросы по тегам