Решение одной проблемы численного интегрирования при накоплении малых приращений.

В задаче «Движение трех тел» столкнулся с проблемой потери данных во время сложения переменных с плавающей точкой, у которых порядок значительно отличается. Оказывается решить эту проблему можно, разделив переменную на старшую (HIGH) и младшую (LOW) части. Операции сложения и вычитания выполняются так:

1. Сложение LOW и запись во временную переменную WRKL.

2. Сложение HIGH и запись во временную переменную WRKH.

3. Вычисление разности порядков

DifOrd = order(WRKH) – order(WRKL)

4. Цикл, пока DifOrd > 0

5. WRKH = WRKH + 2^DifOrd

6. WRKL = WRKL - 2^DifOrd

7. DifOrd = order(WRKH) – order(WRKL)

8. Конец цикла

Для понимания подпрограммы на языке Ассемблер примем следующие обозначения:

Src1 – первое слагаемое;

Src2 – второе слагаемое;

Wrk – временная переменная;

Ord – временная переменная;

.Low – младшая часть числа;

.High – старшая часть числа;

Теперь несложно выполнить алгоритм сложения:

Fld Src1.Low ; загрузить младшую часть слагаемого 1

Fadd Src2.Low ; прибавить младшую часть слагаемого 1

Fstp Wrk.Low ; сохранить временную переменную

Fld Src1.High ; загрузить старшую часть слагаемого 1

Fadd Src2. High ; прибавить старшую часть слагаемого 1

Fstp Wrk. High ; сохранить временную переменную

Выполнено раздельное сложение младших и старших частей, надо перенести старшие разряды Low, которые подходят для части High.

Cicl:

Movzx eax, word ptr Wrk.Low[6] ; порядок младшей части

and eax, 07FFFh ; убрать знак числа

Movzx edx, word ptr Wrk.High[6] ; порядок старшей части

and edx, 07FFFh ; убрать знак числа

Sub edx, eax ; разность порядков

Jng EndCicl ; выход, если edx <=0

Mov Ord, edx

Fild Ord ; загрузить порядок двойки

Fld1

Fscale ; возвести 2 в степень

Fld Wrk.High

Fadd st, st(1) ; добавить младший разряд к High

Fstp Wrk.High ; сохранить High

Fld Wrk.Low

Fsub st, st(1) ; Вычесть старший разряд из Low

Fstp Wrk.Low ; сохранить Low

Finit ; сбросить загруженные значения

Jmp Cicl ; возврат на начало цикла

EndCicl;

Переписать Wrk по адресу результата.

Приведенный алгоритм правильно работает, если Wrk – положительное число. Если предполагается, что Wrk может быть отрицательным, необходимо сделать ветвление программы. Для отрицательных чисел строка

Fadd st, st(1) ; добавить младший разряд к High

Заменяется на

Fsub st, st(1) ; вычесть младший разряд из High

Все остальное не меняется, так как мантиссы чисел с плавающей точкой всегда положительны.