Воробьёвы

(-:

№29.5. Экстрим программированние для дZенствующих. Бесконечные просторы чисел REAL

Последние события, происходящие на WASM.RU поражают трагичностью и одновременно комизмом... Демон коннекта потерял покой и сон, от чего не спится Aquilе. На COMPO пришли потрясающие предложения и примеры. Так держать!!! Serrgio пробует роль Тутанхамона. А Edmond купил четыре плитки шеколада и Дzenствует по полной программе, компилируя шеколад с коньяком!!! Все желающие могут к ниму присоеденится по DialUp.

Опыты с CSS закончились плачевно: Windows обиделась и не желает показывать синий экран :(.

Вообщем так: "Танцуют ВСЁ!!!"

После просвещения FPU мы продолжим математическое Дzenствование по бесконечным просторам чисел REAL...

Chingachguk :: HI-TECH

Сплайн. Аппроксимация данных.

Вводная часть.

В данном документе речь пойдет об интерпретации экспериментальных данных.

Под экспериментальными данными понимается некий набор измерений, являющийся результатом опыта.

Так, например вы стоите на выходе из метро и считаете число людей, проходящих мимо вас за каждую минуту. Допустим, в результате у вас получились следующие измерения: 45,23,55,87,53,48,61... Обобщив измерения в виде гистограммы, вы получаете некоторое распределение:

    |Количество измерений данного числа людей в минуту
    |30        ----
    |29       |    |
    |28     --      --
    |27    |          |
    |26  --           |
    |25 |              --
    |24 |                |
    |23 |                 --
    | --                    |    Число людей за минуту
    ----------------------------------->
      44 45 46 47 48 49 50 51 52

Для экспериментатора может быть интересно понять, не может ли данное распределение быть описано некоторой аналитической формулой. Иначе говоря, определить природу процесса: с какой вероятностью можно ожидать выхода определенного числа людей в произвольном измерении, среднее число людей, выходящих из метро, наиболее вероятное и т.д.

У вас может возникнуть вопрос, зачем такие сложности, если все эти вопросы можно снять прямым измерением. Т.е. пусть перед нами стоит задача определить, сколько людей проходит за 8 часов через выход за рабочий день. Да, можно стоять все эти восемь часов и старательно записывать каждого отдельного человека. Именно так и поступает наше гос-во, когда хочет переписать всех мемберов нашей большой тусовки. При этом на исходе восьмого часа у вас уже слезятся глаза, у вас трижды проверяли документы и дважды подходили какие-то амбалы...

Однако анализируя результаты получасового опыта, в течении которого вы замеряли число человек за минуту, вы имеете другую возможность решить поставленную задачу. Допустим, вы предполагаете, что данное распределение (вероятность выхода определенного числа людей в минуту) описывается определенной формулой, проверяете свою гипотезу, убеждаетесь, что она более-менее достоверна. Формула дает Вам среднее значение распределения. Вам стало известно среднее число людей в минуту - достаточно умножить его на число минут в восьми часах для получения результата. Более того, вы также можете оценить ошибку ваших расчетов, поскольку подразумевается кратковременное (неполное) измерение, очевидно содержащее некоторую неопределенность результата в сравнении с абсолютно точным измерением, т.е. формула также может дать точность оценки, например:

    Число человек в минуту = 47 ± 5;

где 5 - точность оценки за минуту, или примерно 10%. С такой же точностью можно будет определить число людей за восемь часов:

    Число человек за 8 часов = (47 ± 5)?(8?60) = 22560  ± 2400;

Пример аппроксимации данных опыта.

Перейдем от абстрактного опыта к более жизненному. Предположим, вы анализируете свойства генератора 8-ми байтных ключей. Насколько случайными получаются ключи? Т.е. насколько близка псевдослучайная последовательность к истинно случайной.

Критерием качества псевдослучайной последовательности может служить число нулевых/единичных бит, пар бит 00-01-10-11, троек бит... Все они должны быть равновероятны, т.е. вероятность выпадения в каждом ключе определенного числа нулевых бит должно описывается так называемым биномиальным распределением:

              k
     f(k) = C  / 2n;
              n

где n - полное число бит в ключе (в нашем случае 8?8 = 64), k - определенное число бит, Ckn - число вариантов размещения k нулевых бит по n свободным местам:

       k
     C   = n! / k! ? (n-k)!;
       n

2n - число всех битовых комбинаций в числе из n бит, f(k) же вероятность выпадения k бит в числе из n бит.

Таким образом, появляется возможность проверки генератора ключей на уровне числа битов, пар битов, троек... Для этого необходимо получить достаточное число ключей (провести опыт), получить данные опыта - число нулевых бит в каждом ключе и попытаться аппроксимировать, т.е. описать полученные данные некоторой зависимостью (формулой). Если распределение числа нулевых бит недостоверно описывается биномиальным распределением, что будет наблюдается, например, в случае следующего алгоритма генерации ключа:

     LOOP INDEX=1 TO 8
       KEY[INDEX] = PASSWORD[INDEX] XOR MAGIC_NUMBER
     END LOOP

то можно будет сделать определенные выводы о свойствах генератора.

В некоторых случаях аппроксимацию данных называют "сплайном". Данное наименование можно встретить в таких популярных программах, как Excel, в котором можно попытаться выполнить сглаживание (splain) зависимости, например полиномом степени от 2 до 6, правда весьма примитивно. Также можно упомянуть визуальное отображение звука в некоторых CD/mp3 плеерах в виде огибающей кривой. В качестве более-менее профессиональных средств анализа численных данных рекомендую Origin.

Еще один пример аппроксимации (реальный опыт).

Приведу также пример одного реального опыта, выполненного учеными во льдах Антарктики.

Как известно, полярная шапка нашей планеты представляет собой огромную толщу льда (до нескольких километров), лежащую на твердой (каменистой) поверхности Земли. Она образовывалась тысячи лет: вода замерзала и увеличивала эту толщу.

Ученые провели следующий опыт: они бурили небольшую скважину в толще льда и исследовали особенности намерзания (вспоминается файл "Лёд" из "Секретных материалов"). Оказалось, что намерзание происходило этапами: например, на некоторой глубине обнаружилось намерзание толщиной 10 метров, потом 50, потом 30 метров - подобно годичным кольцам деревьев.

Как же ученые интерпретировали полученные ими результаты толщин слоев намерзшего льда? Они сделали следующее достаточно простое предположение: допустим, толщина намерзания льда в зависимости от времени описывается суммой гармоник (синусоид), т.е.:

    Намерзло льда (t) = S(ai×sin(wi×t + bi));

где S означает сумму синусоид, ai - амплитуда i-ой синусоиды, wi - ее период, bi - ее начальная фаза.

Проведя соответствующие расчеты (называемые разложением в ряд Фурье, когда некоторую зависимость f(x) представляют как сумму синусоид с различными частотами и амплитудами), ученые определили, что процесс намерзания хорошо описывается суммой всего трех гармоник с тремя различными периодами ! Это означает, что на нашей планете происходили (и происходят) некие периодические процессы, влияющие в том числе и на намерзание льда на полярных шапках. Так, один из периодов составил около 10 000 лет - следовательно, с такой периодичностью на Земле наступает некоторый процесс похолодания.

Приближение экспериментальных данных аналитической зависимостью.

Критерий Хи-квадрат.

Как же осуществить проверку достоверности выбранной аппроксимирующей функции? Допустим, у вас есть следующая зависимость:

    X  Y
    1  2
    2  4
    3  6
    4  7
    5  10
    6  13
     ...

Разумно полагая, что это что-то вроде y=ax, причем скорее всего a=1, только точки x=4 и x=6 не совсем вписываются в эту гипотезу, вы хотели бы проверить свою теорию и получить что-то вроде примерного значения "a" и некоторую оценку достоверности своего предположения.

Каким образом можно осуществить такую проверку? Сформулируем задачу в наиболее общем виде: обозначим аргументы (входные данные) за Xi, значения зависимости (выходные данные) за Yi, погрешность значений за Eyi. Аппроксимирующую функцию (допускаемый аналитический вид зависимости) за YiТеор. Тогда назовем критерием достоверности аппроксимации экспериментальной зависимости Yi = function (Xi) теоретической зависимостью YiТеор = function (Xi) функцию:

    Cr = function(Xi,(Yi,Eyi),YiТеор);

значение которой говорит о том, насколько точно хорошо теоретическое распределение YiТеор = function (Xi) описывает реальное распределение Yi = function (Xi).

В нашем простом случае у входных данных (X=1,2,3,4...) нет никаких погрешностей, равно как нет их и у выходных данных (Y=2,4,6,7,10,13). Теоретическое распределение имеет простой вид YТеор = aX с одним неизвестным параметром (a). Предполагаем a=2 и получаем ряд значений теоретического распределения:

    YТеор(X) = 2X = 2,4,6,8,10,12

Полученные значения YТеор необходимо подставить в критерий-функцию Cr и по ее значению оценить достоверность выбранной аппроксимации.

Какой же вид должна иметь критерий-функция Cr? Допустим, чем меньше ее значение, тем достовернее аппроксимация. Очевидно, она должна возрастать тем больше, чем больше получается отклонение экспериментальных значений Y от аппроксимирующих их YТеор, учитывая все точки распределения.

Таких критерий-функций существует несколько. Приведу пример функции Xq ("Хи-квадрат"). Вот ее вид:

    Xq = S( (Yi-YТеорi)2 / Eyi2 )i=1..N

где S означает суммирование квадратов разницы между экспериментальными и теоретическими значениями, деленными на квадраты погрешностей измерения.

Оценим поведение этой функции. Чем больше отклонение теоретического приближения от экспериментального, тем больше значение Xq, т.е. чем больше значение Xq, тем менее достоверным является приближение.

А вот как будет выглядеть запись Xq для нашего случая с линейной (y=ax) зависимостью:

    X  Y  YТеор=ax YТеор-Y (Y-YТеор)2
    1  2     2       0        0
    2  4     4       0        0
    3  6     6       0        0
    4  7     8       1        1
    5  10    10      0        0
    6  13    12      -1       1
     ...

Итого в этом случае Xq=1+1=2.

Тонкий момент. Обратите внимание на то, что каждый квадрат отклонения делится на квадрат погрешности измерения Y - Eyi. Смысл в том, что разные точки распределения могут иметь различные ошибки, другими словами - давать различное количество информации о характере искомого распределения: чем больше погрешность точки, тем меньше ее вклад в оценку достоверности приближения. В нашем случае Eyi отсутствуют у всех точек, что эквивалентно тому, что все точки вносят одинаковый вклад в искомый вид теоретического распределения.

Оценку достоверности аппроксимации при помощи Xq можно определить по так называемым "таблицам Стьюдента", в которых приводится максимально допустимое значение Xq для аппроксимаций распределений с различным числом входных данных и различным числом параметров аппроксимирующего распределения. Однако должно быть ясно, что чем больше параметров имеет приближение и чем больше входных данных, тем большим будет минимально допустимое значение Xq, определяющее достоверность аппроксимации.

В нашем случае говорят, что зависимость Y(x) = 2,4,6,7,10,13.. аппроксимируется функцией y=ax,a=2 с точностью Xq = 2 для 6 точек, число параметров аппроксимирующей функции равно 1 (это a=2).

Применение критерия Хи-квадрат для аппроксимации.

Разберем более детально пример с аппроксимацией распределения линейной функцией y=ax и покажем, как можно найти искомый параметр распределения "a" не предполагая его значение, а вычислив его с помощью функции Xq.

Способ первый. Аналитическое решение.

Замечаем, что в случае приближения y=ax Xq - это функция от одного аргумента:

     Xq = function(a);

поскольку множества Y и X фиксированы на момент расчета, а YТеор вычисляется как функция от X. Тогда для того, чтобы найти минимум функции Xq достаточно решить уравнение:

 dXq
 ---  = S(2×(Yi-aXi)×(-Xi) ) = 0;
 da

найдя неизвестное "a".

В более общем случае теоретическое распределение может быть записано в виде:

    YТеор = function( a0 , a1 , a2... ),

где a0..am - параметры искомого теоретического приближения. Очевидно, уравнение нахождения неизвестных коэффициентов aj можно записать в виде:

     dXq
     --- = 0,
     da1
     dXq
     --- = 0,
     da2
     ...
     dXq
     --- = 0,
     dam

т.е. получается система из m уравнений, образованных записью частной производной по каждому из a1..am параметру. Решая эти уравнения совместно, можно получить неизвестные параметры a1..am приближающего распределения.

Способ второй. Численные методы.

Несмотря на то, что Способ 1 может дать однозначное и наиболее точное значение параметров приближающего распределения aj, существует и другой способ определить неизвестные aj. Он заключается в поиске на некотором поле значений aj таких, при которых значение Xq минимально. Поясню мысль рисунком:

     a1 minimum -- a1 minimum+step_a1 -- a1 minimum+step_a1×2 -- ... a1 max
     a2 minimum -- a2 minimum+step_a2 -- a2 minimum+step_a2×2 -- ... a2 max
     ...
     am minimum -- am minimum+step_am -- am minimum+step_am×2 -- ... am max

Т.е. вначале задаются некие граничные (максимальные и минимальные) значения каждого из неизвестных параметров aj и шаг разбиения интервала (aj max - aj min) - step_aj. Каждое значение aj пробегает значения от минимума до максимума. Для каждой комбинации параметров aj (очевидно, всего возможно ((число шагов разбиения))число параметров aj вариантов) рассчитывается значение Xq и ищется такая комбинация, которая отвечает минимальному Xq

Позволю себе назвать такой способ поиском "на решетке"(lattice) значений неизвестных параметров aj. Значение step_aj назовем "шагом решетки".

После того, как на данной "решетке" значений параметров будут найдены параметры, отвечающие минимальному Xq, необходимо будет выполнить модификацию решетки с целью локализации значения каждого искомого параметра в области истинного значения параметра:

- Если значение параметра, отвечающего минимальному Xq будет обнаружено внутри решетки (т.е. он будет между минимальным и максимальным значением для данного параметра), то новые границы будут определены следующим образом:

Новое значение минимума/максимума = Значение параметра ± Шаг решетки

- Если значение параметра, отвечающего минимальному Xq будет обнаружено на краю решетки (т.е. он будет равен либо минимальному либо максимальному значению для данного параметра), то новые границы будут определены следующим образом:

Новое значение минимума/максимума = Значение параметра ± ширина решетки.

Рассмотрим пример того, как это бы происходило в случае нашей простой линейной зависимости Y(x) = 2,4,6,7,10,13.. Мы предполагаем, что у аппроксимирующей зависимости один параметр: a. Найдем его. Сначала зададим максимальное и минимальное начальное значение a:

     amax = 100, amin = -100;

число разбиений решетки пусть будет равно 5, тогда в первом приближении (на первой решетке) будет проверены следующие значения a:

    -100 -50 0 +50 +100,

при этом шаг решетки равен:

    (amax - amin) / (5-1) = ( 100 - (-100) ) / 4 = 50,

и мы вычислим пять значений функции Xq для всех значений a:

     a     -100       -50        0          +50        +100
     Xq(a) 9.47E+0005 2.46E+0005 3.74E+0002 2.09E+0005 8.73E+0005

Минимум Xq в данном случае получается если a=0. Поскольку точка лежит внутри решетки (a <> -100 и a <> 100), то новые максимальные и минимальные значения "a" необходимо выбрать следующим образом:

     amax = 0 + 50, amin = 0 - 50;

и перейти к следующему расчету Xq на новом интервале.

Чем больше будет выполнено таких приближений, тем точнее будут найдены искомые параметры, но до известного предела, заложенного в самом критерии Xq.

Программная реализация именно такого способа будет рассмотрена ниже. В чем его достоинства/недостатки по сравнению с Способом 1? Достоинством является универсальность относительно вида приближающей функции YТеор: нет необходимости решать сложную систему уравнений в частных производных по aj, недостатком является некоторая неточность результатов определения параметров aj, что будет показано ниже.

Реализация аппроксимации данных полиномом численными методами.

Реализуем программно решение следующей задачи: пусть имеется набор значений X,Y, предположительно образованных полиномом вида:

    Y(x) = anxn + an-1x(n-1) + ... a0

Попытаемся построить алгоритм, способный аппроксимировать заданную зависимость Y(X) полиномом вида Y(x) = anxn + .. +a0 : найти неизвестные коэффициенты an,an-1,..,a0 и дать оценку достоверности приближения с помощью критерия Xq.

Реализация алгоритма будет проведена на win32-ассемблере с использованием сопроцессора, все особенности работы с которым будут по возможности прокомментированы.

Заголовок программы на MASM: стандартное подключение библиотек и т.д.

     .386
     .model flat,stdcall
    option casemap:none
    include \masm32\include\windows.inc
    include \masm32\include\kernel32.inc
    includelib \masm32\lib\kernel32.lib
    include \masm32\include\user32.inc
    includelib \masm32\lib\user32.lib

Описываем данные и константы программы.

    ; Флаги вывода промежуточных результатов (Debug):
; выводить значение степени полинома из файла
ShowPolinomPower_NOT_DEFINED equ TRUE
; выводить значения X,Y из файла
ShowXYAtFileReading_NOT_DEFINED equ TRUE
FloatMask equ 14 ; Маска для вывода чисел с плавающей точкой:
                        ; 14 цифр после десятичной точки
MaxPower equ 10         ; Максимальная степень входного полинома:
                        ; y(x)=a0+a1×x+a2×x2+...aMaxPower×xMaxPower

; Размер решетки для вычислений:
LatticeSize equ 4 ; (число разбиений интервала поиска)
ImpLattices equ 500     ; Максимальное число итераций
MaxPoints equ 40        ; Максимальное число входных точек
 .data
; Максимальное (по модулю) начальное значение коэффициентов
MaxStartCoffValue dd 1.0E+10
BreakValueOfXq dd 0.1   ; Величина, контролирующая момент, когда
                        ; приближение считается достоверным
coff struc              ; Параметры полинома: его коэффициенты
  Cofficient dd (MaxPower+1) dup (?)
coff ends
Point struc             ; Описание одной точки: пара (X,Y)
  X dd?
  Y dd?
Point ends
 .data?
Power dd?              ; Предполагаемая степень из файла
EntryData db MaxPoints×sizeof(Point) dup(?) ; Входные точки (X,Y)
NumberDataPoints dd?   ; Число точек в файле
CurrLattice dd?        ; Текущий номер решетки
TempReal dd?           ; Временная переменная
                        ; Массивы коэффициентов полинома:
amin coff <>            ; Набор минимальных коэффициентов
amax coff <>            ; Набор максимальных коэффициентов
awork coff <>           ; Набор промежуточных коэффициентов
afoundmin coff <>       ; Набор коэффициентов, отвечающих минимальному Xq
Xq dd?                 ; Значение Xq
Xqmin dd?              ; значение минимального Xq на очередном шаге вычислений
; Набор индексов для прохождения узлов решетки
LatticeIndexes dd (MaxPower+1) dup(?)
; Индексы, отвечающие меньшему Xq
NewIndexes dd (MaxPower+1) dup(?)

Константы ShowPolinomPower и ShowXYAtFileReading определяются для вывода отладочных сообщений при чтении параметров из файла. Они начнут работать, если удалить "_NOT_DEFINED".

FloatMask управляет выводом чисел с плавающей точкой: в данном случае эта величина определяет вывод 14 знаков после десятичной точки.

MaxPower определяет максимально допустимую степень аппроксимирующего полинома, поддерживаемого программой - 10. MaxPoints определяет максимально возможный размер входных данных (точек x,y)

LatticeSize является количеством разбиений интервала [максимальное значение - минимальное значение] каждого коэффициента полинома. ImpLattices определяет число итераций в численном поиске минимального Xq (число анализируемых решеток). MaxStartCoffValue определяет начальные значения всех коэффициентов полинома: вначале полагаем макисмально/минимально возможными значения коэффициентов равными ± MaxStartCoffValue. BreakValueOfXq определяет точность определения коэффициентов полинома: чем больше этот порог, тем хуже будет определены коэффициенты полинома: на самом деле эта величина должна быть согласована с таблицей Стьюдента.

.data
ErrTitle  db "Error? ;)",0 ; Заголовок сообщения отладки/ошибки
 .code
start:
 .data
  AppName db "Splain test program",0
  Welcome db "Example of Approximate Data by Polinom. Coded by Chingachguk, 2002.",10,10
          db "Entry Data(must be in file 'splain.dat'):",10
          db "<polinom power>",10
          db "<X1>,<Y1>",10
          db "...",10
          db "<Xn>,<Yn>",0
 .code
  invoke MessageBox,NULL,addr Welcome,addr AppName,MB_OK+MB_ICONINFORMATION
 .data
IniFileName db "splain.dat",0 ; Имя файла данных
 .data?
IniDescr dd?           ; Дескриптор файла
 .code
; ×××××××××××××××× Read Entry File ×××××××××××××××××××××××××××××

; Open entry file (not Create !)
  invoke CreateFile,offset IniFileName,GENERIC_WRITE or GENERIC_READ,\
    FILE_SHARE_WRITE or FILE_SHARE_READ,0,OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,0
  mov  IniDescr,eax     ; Save file descriptor
  cmp  eax,0ffffffffh
  jz   FileNotOpen
 .data?
CurrString dd?
FileStringSize equ 80
FileString db FileStringSize dup(?)
 .code
  ; Ведем счет строкам входного файла
  mov  dword ptr CurrString,1
  ; Читаем предполагаемую степень полинома (max: MaxPower)
  push FileStringSize   ; размер приемника
  push offset FileString ; Адрес приемника строки
  push IniDescr         ; Дескриптор файла
  call ReadFileString
  ; Переводим строку со степенью в число (real)
  push offset Power
  push offset FileString
  call StringToVal
  jc   InvalidNumericFormat
  ; Переводим степень из real в int
  fld  dword ptr Power
  ; Задаем режим округления (биты 11 - отбросить дробную часть)
  call Set_FPUTruncType
  frndint
  fistp dword ptr Power
  ; Проверяем степень на валидность
  cmp  dword ptr Power,0
  jl   InvalidPowerValue
  cmp  dword ptr Power,MaxPower
  ja   InvalidPowerValue
; ××× Отладочное сообщение: значение предполагаемой степени полинома ×××
ifdef ShowPolinomPower
 .data
PolinomPowerName db "Power: "
 .code
  mov  edi,offset FileString
  mov  esi,offset PolinomPowerName
  mov  ecx,sizeof PolinomPowerName
  cld
  rep  movsb
  mov  eax,Power
  mov  ebx,10000000
  call DecChar
  mov  byte ptr [edi+8],0
  invoke MessageBox,0,offset FileString,offset ErrTitle,MB_OK
endif
  ; Читаем пары точек (X,Y) из файла
  mov  dword ptr NumberDataPoints,0
  ; Стартовое значение базы массива EntryData (ebp)
  xor  ebp,ebp
@@ReadPoints:
  ; Ведем счет строкам входного файла
  inc  dword ptr CurrString
  ; Читаем X
  push FileStringSize   ; размер приемника
  push offset FileString ; Адрес приемника строки
  push IniDescr         ; Дескриптор файла
  call ReadFileString
  test eax,eax
  jz   @@ReadPointsDone
  ; Переводим строку в число (real)
  lea  eax,EntryData.Point[ebp].X
  push eax
  push offset FileString ; Адрес строки
  call StringToVal
  jc   InvalidNumericFormat
  ; Читаем Y
  push FileStringSize   ; размер приемника
  push offset FileString
  push IniDescr
  call ReadFileString
  test eax,eax
  jz   @@ReadPointsDone
  ; Переводим строку в число (real)
  lea  eax,EntryData.Point[ebp].Y
  push eax
  push offset FileString
  call StringToVal
  jc   InvalidNumericFormat
; ××× Отладочное сообщение: значения X,Y ×××
ifdef ShowXYAtFileReading
  cld
  mov  edi,offset FileString
  mov  eax,' :X '
  stosd
  push FloatMask  ; Маска вывода
  lea  eax,EntryData.Point[ebp].X
  push eax   ; Число
  push edi   ; Изображение числа
  call ValToString
  push FloatMask
  lea  eax,EntryData.Point[ebp].Y
  push eax ; Число
  mov  esi,offset FileString
  call Str_Len
  lea  edi,[esi+ecx]
  mov  eax,' :Y '
  stosd
  push edi   ; Изображение числа
  call ValToString
  invoke MessageBox,0,offset FileString,offset ErrTitle,MB_OK
endif
  ; Максимальное число входных точек не должно быть больше MaxPoints
  inc  dword ptr NumberDataPoints
  cmp  dword ptr NumberDataPoints,MaxPoints
  ja   InvalidPointsValue
  ; Индексируем следующие точки в массиве
  add  ebp,sizeof Point
  jmp  @@ReadPoints
@@ReadPointsDone:
; ×××××××××× approximate Data ×××××××××××××
; Прочитав все точки из файла, задаем начальные значения коэффициентов
; set initial coefficient's values (равные ±MaxStartCoffValue)
  mov  ecx,Power
  inc  ecx
  xor  esi,esi
@@SetStartCoffs:
  fld  MaxStartCoffValue
  fchs                  ; Сделать отрицательным
  fstp amin.Cofficient[esi]
  fld  MaxStartCoffValue
  fstp amax.Cofficient[esi]
  ; Следующий коэффициент
  add  esi,sizeof coff / (MaxPower+1)
  loop @@SetStartCoffs
; ××××××× Основной цикл: поиск коэффициентов, отвечающих минимальному Xq ×××××××

; В этом цикле мы вычисляем Xq на всех узлах решетки,
; находим минимальный Xq и запоминаем отвечающие ему
; значения коэффициентов решетки,
; затем сжимаем/сдвигаем решетку минимального Xq в зависимости от того,
; лежат ли коэффициенты на краю (внутри) решетки.
; Если значение Xq показывает, что приближение достоверно,
; мы заканчиваем основной цикл поиска
  ; Номер решетки приближения
  mov  dword ptr CurrLattice,1
NextLattice:
  ; find Xq minimum at the current Lattice
  ; Получить стартовое значение Xq
  ; (Получить значение полинома с минимальными коэффициентами (массив amin))
  fldz
  fstp Xqmin
  xor  esi,esi          ; Индекс - esi - номер точки из файла
  mov  ecx,NumberDataPoints ; Число точек в файле
@@GetStartXq:
  push dword ptr EntryData.Point[esi].X ; значение аргумента (X)
  push Power            ; размер (степень) полинома
  push offset amin      ; адрес массива коэффициентов полинома
  call polinom
  ; Результат: значение полинома в ST(0)
  ; вычислить разницу -(Yi-polinom(Xi))
  fsub dword ptr EntryData.Point[esi].Y
  ; Возвести в квадрат разницу -(Yi-polinom(Xi)), она в ST(0)
  mov  eax,2                            ; степень - 2
  call get_pow          ; вернется в ST(0) квадрат
  fadd dword ptr Xqmin  ; Прибавить к накопителю Xqmin
  fstp dword ptr Xqmin  ; Запомнить накопитель
  add  esi,sizeof Point
  loop @@GetStartXq
  ; Обнуляем стартовые коэффициенты для прохода по значениям коэффициентов
  mov  ecx,Power
  inc  ecx
  mov  edi,offset LatticeIndexes
  cld
  xor  eax,eax
  rep  stosd
  ; выполняем поиск Xq на текущих значениях коэффициентов
FindXqMinimum:
  ; Определяем коэффициенты текущего приближенного полинома (awork)
  xor  esi,esi          ; Индексы - esi, edi
  xor  edi,edi
  mov  ecx,Power        ; По всем степеням
  inc  ecx
@@GetCurrentCoffs:
  ; вычисляем разницу между максимальным и минимальным значением коэффициента
  fld  amax.Cofficient[esi]
  fsub amin.Cofficient[esi]
  ; Делим разницу на число (узлов решетки-1), получая шаг решетки
  mov  TempReal,LatticeSize-1
  fidiv dword ptr TempReal
  ; Умножаем на номер узла текущего коэффициента
  fimul dword ptr LatticeIndexes[edi]
  ; Получаем значение коэффициента в узле, добавляем минимальное значение
  fadd dword ptr amin.Cofficient[esi]
  ; вычисляем текущие коэффициенты полинома (awork)
  fstp dword ptr awork.Cofficient[esi]
  ; Следующий коэффициент
  add  esi,sizeof coff / (MaxPower+1)
  add  edi,4            ; Следующий номер узла
  loop @@GetCurrentCoffs
  ; вычисляем Xq для промежуточных значений коэффициентов (awork)
  fldz
  fstp Xq
  xor  esi,esi          ; Индекс - esi - номер точки из файла
  mov  ecx,NumberDataPoints ; Число точек в файле
@@GetCurrentXq:
  push dword ptr EntryData.Point[esi].X
  push Power            ; размер (степень) полинома
  push offset awork     ; адрес массива коэффициентов полинома
  call polinom
  ; Результат: значение полинома в ST(0)
  ; вычислить разницу -(Yi-polinom(Xi))
  fsub dword ptr EntryData.Point[esi].Y
  ; Возвести в квадрат разницу -(Yi-polinom(Xi)), она в ST(0)
  mov  eax,2            ; степень - 2
  call get_pow          ; вернется в ST(0) квадрат
  fadd dword ptr Xq     ; Прибавить к накопителю Xqmin
  fstp dword ptr Xq
  add  esi,sizeof Point
  loop @@GetCurrentXq
  ; Сравниваем Xq и Xqmin
  fld  Xq
  fcomp dword ptr Xqmin ; Сравнить Xq с Xqmin и вытолкнуть из стека ST(0)
  fstsw ax              ; Перенести флаги в регистр ax
  sahf                  ; Дублировать их в регистре флагов
  jnc  @@CurrXqNoLow
  ; Если Xq < Xqmin, то присвоить новое значение Xqmin
  ; Т.е. мы нашли коэффициенты с меньшим Xq, запомним их
  fld  dword ptr Xq
  fstp dword ptr Xqmin
  ; Также запомнить текущие коэффициенты в afoundmin
  mov  ecx,Power
  inc  ecx
  mov  esi,offset awork
  mov  edi,offset afoundmin
  cld
  push ecx
  rep  movsd
  pop  ecx
  ; Запомнить текущие номера узлов решетки
  mov  esi,offset LatticeIndexes
  mov  edi,offset NewIndexes
  rep  movsd
@@CurrXqNoLow:
  ; Инкрементировать номера узлов решетки,
  ; перебирая все узлы решетки
  xor  esi,esi
@@NextApprox:
  inc  dword ptr LatticeIndexes[esi*4]
  cmp  dword ptr LatticeIndexes[esi*4],LatticeSize
  jb   FindXqMinimum    ; Следующая решетка
  mov  dword ptr LatticeIndexes[esi*4],0
  inc  esi
  cmp  esi,dword ptr Power
  jbe  @@NextApprox
  ; Мы прошли по всем узлам решетки
  ; Анализируем, где был найден минимум Xq
  mov  ecx,Power
  inc  ecx
  xor  esi,esi
@@AnaliseLattice:
  mov  eax,dword ptr NewIndexes[esi]
  test eax,eax
  jz   @@XqMinLieOnRange
  cmp  eax,LatticeSize-1
  jnz  @@XqMinNoLieOnRange
@@XqMinLieOnRange:
  ; point lie on the rage of lattice: no change size of lattice !
  ; Строим решетку симметрично относительно найденного узла
  ; Новый минимум/максимум=найденный узел±ширина решетки
  ; Получаем в ST(0) новое значение минимального коэффициента
  fld  dword ptr afoundmin.Cofficient[esi]
  fsub dword ptr amax.Cofficient[esi]
  fadd dword ptr amin.Cofficient[esi]
  ; Получаем в ST(0) новое значение максимального коэффициента, мин.->ST(1)
  fld  dword ptr afoundmin.Cofficient[esi]
  fadd dword ptr amax.Cofficient[esi]
  fsub dword ptr amin.Cofficient[esi]
  jmp  @@NextLatticePoint
@@XqMinNoLieOnRange:
  ; point lie inside of lattice: minimize size of lattice
  ; Строим решетку симметрично относительно найденного узла
  ; Новый минимум/максимум=найденный узел±шаг решетки
  ; Получаем в ST(0) новое значение минимального коэффициента
  fld  dword ptr amin.Cofficient[esi]
  fsub dword ptr amax.Cofficient[esi]
  ; Делим разницу на число (узлов решетки-1), получая шаг решетки
  mov  TempReal,LatticeSize-1
  fidiv dword ptr TempReal
  fadd dword ptr afoundmin.Cofficient[esi]
  ; Получаем в ST(1) новое значение максимального коэффициента, ST(0)->ST(1)
  fld  dword ptr amax.Cofficient[esi]
  fsub dword ptr amin.Cofficient[esi]
  ; Делим разницу на число (узлов решетки-1), получая шаг решетки
  fidiv dword ptr TempReal
  fadd dword ptr afoundmin.Cofficient[esi]
@@NextLatticePoint:
  ; ST(0) - новый максимум, ST(1) - новый минимум
  fstp dword ptr amax[esi] ; <-ST(0)
  fstp dword ptr amin[esi] ; <-ST(0)
  add  esi,4
  loop @@AnaliseLattice
  ; Переходим к следующей решетке, выполняя не более ImpLattices приближений
  inc  dword ptr CurrLattice
  ; Проверяем величину текущего минимума Xq
  ; Возможно, мы уже нашли коэффициенты
  fld  dword ptr BreakValueOfXq
  fmul dword ptr Power
  fadd dword ptr BreakValueOfXq
  fimul dword ptr NumberDataPoints
  ; ST(0)=NumberDataPoints × BreakValueOfXq × (Power+1)
  ; Если величина Xqmin =< NumberDataPoints×BreakValueOfXq×(Power+1),
  ; то приближение считается достоверным
  ; Сравнить NumberDataPoints×BreakValueOfXq×(Power+1) с Xqmin и вытолкнуть из стека ST(0)
  fcomp dword ptr Xqmin
  fstsw ax              ; Перенести флаги в регистр ax
  sahf                  ; Дублировать их в регистре флагов
  jnc  @@ApproximateValid
  cmp  dword ptr CurrLattice,ImpLattices
  jbe  NextLattice
@@ApproximateValid:
; ×××××××××× Цикл поиска коэффициентов закончен ×××××××××××××××××

; выводим полученные коэффициенты, значение Xq для анализа достоверности
; нашего приближения данных из файла полиномом
 .data
ResultTitle db "Calculations resume",0
ResultHeader0 db "Assume that you data was approximated by polinom with power????.",10,10
ResultHeaderD db "a[?]=?.????????????????????  "
ResultHeaderE db "Xq =?.???????????????????? for????? points.",10
 .data?
ResultHeader db 1000 dup(?)
 .code
  cld
  mov  edi,offset ResultHeader
  ; выводим степень полинома
  push edi
  mov  eax,Power
  mov  ebx,1000
  mov  edi,offset ResultHeader0+60
  call DecChar
  pop  edi
  mov  esi,offset ResultHeader0
  mov  ecx,sizeof ResultHeader0
  rep  movsb
  ; выводим полученные коэффициенты
  mov  ecx,Power        ; По всем степеням
  inc  ecx
  xor  esi,esi
@@OutCoffs:
  pushad
  mov  eax,Power
  inc  eax
  sub  eax,ecx
  mov  ebx,1
  mov  edi,offset ResultHeaderD+2
  call DecChar
  push FloatMask        ; Маска вывода
  lea  eax,afoundmin.Cofficient[esi]
  push eax
  push offset ResultHeaderD+6 ; Изображение числа
  call ValToString
  popad
  push esi
  push ecx
  mov  esi,offset ResultHeaderD
  mov  ecx,sizeof ResultHeaderD
@@CopyCoff:
  lodsb
  test al,al
  jnz  @@UsalSym
  mov  al,10
  mov  ecx,1
@@UsalSym:
  stosb
  loop @@CopyCoff
  pop  ecx
  pop  esi
  ; Следующий коэффициент
  add  esi,sizeof coff / (MaxPower+1)
  loop @@OutCoffs
  mov  al,10
  stosb
  ; выводим Xqmin
  push edi
  push FloatMask ; Маска вывода
  push offset Xqmin ; Число
  push offset ResultHeaderE+5 ; Изображение числа
  call ValToString
  mov  byte ptr ResultHeaderE+27,' '
  mov  eax,NumberDataPoints
  mov  ebx,10000
  mov  edi,offset ResultHeaderE+32
  call DecChar
  pop  edi
  mov  esi,offset ResultHeaderE
  mov  ecx,sizeof ResultHeaderE
  rep  movsb
  xor  al,al
  stosb
  invoke MessageBox,0,offset ResultHeader,offset ResultTitle,MB_OK+MB_ICONINFORMATION

Exit:
  invoke CloseHandle,IniDescr
  invoke ExitProcess,NULL
;×××××××××××× Подпрограмма вычисления значения полинома ××××××××××××××××××××××××
; Get value of LocCoff[i]×xi, i=0..power
; [ebp+8] - адрес массива коэффициентов полинома
; [ebp+0Ch] - размер (степень) полинома
; [ebp+10h] - значение аргумента (X)
; Результат: значение полинома в ST(0)
polinom proc
  push ebp
  mov  ebp,esp
  ; result: [ebp-4], локальная степень полинома: [ebp-8]
  sub  esp,2*4
  push esi
  ; Определяем начальный результат
  fldz
  fstp dword ptr [ebp-4]
  ; Загружаем адрес коэффициентов полинома
  mov  esi,[ebp+8]
  ; Степень полинома
  fldz
  fstp dword ptr [ebp-8]
@@polinom:
  ; Получить XТекущая степень (0..[ebp+0Ch])

  fld  dword ptr [ebp+10h] ; Загрузить X
  mov  eax,dword ptr [ebp-8] ; степень
  call get_pow          ; вернется в ST(0)
  ; Умножить на текущий коэффициент
  fmul dword ptr [esi]
  ; Добавить сумматор
  fadd dword ptr [ebp-4]
  ; Запомнит новый результат
  fstp dword ptr [ebp-4]
  add  esi,sizeof coff / (MaxPower+1) ; Следующий коэффициент
  inc  dword ptr [ebp-8]
  mov  eax,dword ptr [ebp-8]
  cmp  eax,[ebp+0Ch]
  jbe  @@polinom
  ; Возвращаем в ST(0) результат
  fld  dword ptr [ebp-4]
  pop  esi
  leave
  ret  3*4
polinom endp
;×××××××××××××××××××××××××××××××× Сообщения об ошибках ×××××××××××××××××××××××××
FileNotOpen:
 .data
IniOpOK   db "Entry Data File (splain.dat) not open (not exist?) !",0
 .code
  invoke MessageBox,0,offset IniOpOK,offset ErrTitle,MB_OK or MB_ICONERROR
  jmp  Exit
InvalidNumericFormat:
 .data
InvalidNumericFormatMsg db "Invalid numeric format in string "
ErrStrNum               db "00000000",0
 .code
  mov  eax,CurrString
  mov  ebx,10000000
  mov  edi,offset ErrStrNum
  call DecChar
  invoke MessageBox,0,offset InvalidNumericFormatMsg,offset ErrTitle,MB_OK or MB_ICONERROR
  jmp  Exit
InvalidPowerValue:
 .data
InvalidPowerValueMsg db "Invalid Power Value",0
 .code
  invoke MessageBox,0,offset InvalidPowerValueMsg,offset ErrTitle,MB_OK or MB_ICONERROR
  jmp  Exit
InvalidPointsValue:
 .data
InvalidPointsValueMsg db "Too much data points in file !",0
 .code
  invoke MessageBox,0,offset InvalidPointsValueMsg,offset ErrTitle,MB_OK or MB_ICONERROR
  jmp  Exit

Кратко о вспомогательных программах. Для реализации алгоритма нам потребуется несколько подпрограмм:

  • семейство Set_FPU..Type. Подпрограммы для задания режима округления сопроцессора при операциях конвертирования чисел с плавающей точкой в числа типа inteher.
  • get_pow. Подпрограмма возведения в целую степень числа с плавающей точкой. Реализована примитивно умножением аргумента самого на себя степень раз.
  • StringToVal: Подпрограмма деформатирования строки в число типа Real. Позволяет читать из файла числа с плавающей запятой в виде строки и преобразовывать их в формат real (4 байта). Реализована следующим образом: накопитель инициализируется числом 0, затем читается цифра из строки, добавляется к приемнику, приемник предварительно доумножается на 10. Десятичная точка учитывается подсчетом числа цифр после нее и финальным делением на 10 в степени "число цифр после десятичной точки".
  • ValToString: Подпрограмма форматирования числа типа Real в строку. Позволяет выполнять форматный вывод чисел с плавающей точкой с заданным числом знаков после десятичной точки. Алгоритм вывода довольно прост и, возможно, допускает некоторую погрешность: определяется максимальная степень десяти в числе при помощи команды fyl2x сопроцессора: фактически вычисляется max_power = int( log10(x) ). Далее последовательным делением на 10 в степени max_power, max_power-1, ... определяются коэффициенты разложения числа по степеням десяти.
  • ReadFileString: парсер входного файла. Позволяет разделять строки с числовыми данными из входного файла, отделенные стандартными разделителями (10 или 13+10 или пробелами) и читать их в заданный буфер.
; ××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××
; Вспомогательные подпрограммы
; ××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××
 .data
  DecNumber dd 10.0
; ×××××××××××× Set_FPU...Type: Задать режим округления сопроцессора ××××××××××××
; Set_FPURoundType: округление до ближайшего целого (00)
; Set_FPUTruncType: отбрасывание дробной части (11)
 .data?
control dw?
 .code
Set_FPURoundType proc
  fstcw word ptr control
  and word ptr control,0F3FFh
  fldcw word ptr control
  ret
Set_FPURoundType endp
Set_FPUTruncType proc
  fstcw word ptr control
  or  word ptr control,0C00h
  fldcw word ptr control
  ret
Set_FPUTruncType endp
; ×××××××××××× get_pow: Подпрограмма возведения в целую степень числа ××××××××××

; ST(0) - аргумент, eax - степень
; result: ST(0)
; Для возведения в произвольную стпень числа следует воспользоваться
; более сложным алгоритмом с использованием fyl2x:
; см., например, П.И. Рудаков, К.Г. Финогенов
; "Программируем на языке ассемблера IBM PC"
 .code
get_pow proc
  test eax,eax
  jnz  @@NoZeroPower
  fstp ST(0)           ; вытолкнуть аргумент из стека
  fld1                  ; Загрузить 1.0
  ret                   ; Вернуть 1.0
@@NoZeroPower:
  push ecx
  mov  ecx,eax
  ; Сделать степень положительной
@@NegLoop:
  neg  ecx
  js   @@NegLoop
  fld1                  ; Загрузить 1.0, ST(0) - накопитель
@@get_pow:
  fmul ST(0),ST(1)      ; Накопитель=Накопитель x arg
  loop @@get_pow
  fxch ST(1)            ; Поменять накопитель и аргумент
  fstp ST(0)            ; вытолкнуть аргумент из стека
  ; Если степень отрицательна, необходимо Накопитель=1/Накопитель
  cmp  eax,0
  jge  @@NoNegPow
  fld1                  ; Загрузить 1
  fxch ST(1)            ; Обменять ST(1) и ST(0)
  fdivp ST(1),ST(0)     ; Разделить 1 на Накопитель
@@NoNegPow:
  pop  ecx
  ret
get_pow endp
; ×××× StringToVal: Подпрограмма деформатирования строки в число типа Real ×××××
; [ebp+8] - Адрес строки(Zero-terminated)
; [ebp+0Сh] - Адрес приемника (real или dd)
; Десятичным разделителем считается символ "."
; Пробелы перед цифрами пропускаются
; Ошибка - флаг CF установлен
 .code
StringToVal proc
  push ebp
  mov  ebp,esp
  sub  esp,1*4          ; Число из строки [ebp-4]
  mov  ebx,[ebp+0Ch]    ; Адрес приемника храним в ebx
  fldz                  ; Обнуляем приемник
  fstp dword ptr [ebx]
  cld
  ; Запоминаем знак числа
  xor  edi,edi
  mov  esi,dword ptr [ebp+8]
  lodsb
  dec  esi
  cmp  al,'-'
  jnz  @@NoSignBefore
  ; Устанавливаем 31-ый бит знака числа
  inc  edi
  ror  edi,1
  inc  esi
@@NoSignBefore:
  xor  edx,edx          ; Флаг десятичной точки и число знаков после нее
@@ParseString:
  xor  eax,eax
  lodsb
  test al,al
  jz   @@EndParse
  cmp  al,' '           ; Пропускаем пробелы слева
  jz   @@ParseString
  cmp  al,'.'
  jnz  @@NoDecPoint
  test dh,dh
  jnz  @@ErrorSign
  mov  dx,256
  jmp  @@ParseString
@@NoDecPoint:
  sub  al,'0'
  js   @@ErrorSign
  cmp  al,9
  ja   @@ErrorSign
  mov  dword ptr [ebp-4],eax
  fld  dword ptr [ebx]  ; Загружаем текущее значение результата
  fmul DecNumber        ; Умножаем результат на 10
  fiadd dword ptr [ebp-4] ; Добавляем к результату текущее число
  fstp dword ptr [ebx]  ; Запоминаем получившееся значение
  ; Учитываем число знаков после десятичной точки
  inc  dl
  jmp  @@ParseString
@@EndParse:
  ; Учесть десятичную точку
  test dh,dh
  jz   @@ThereAintPoint
  mov  dh,0
  ; Разделить результат на 10edx
  fld  dword ptr DecNumber ; Загрузить аргумент (10)
  mov  eax,edx          ; Целая степень
  call get_pow          ; 10edx вернется в ST(0)
  fld  dword ptr [ebx]  ; Загрузить результат без точки
  fxch ST(1)            ; Обменять ST(1) и ST(0)
  fdivp ST(1),ST(0)
  fstp dword ptr [ebx]  ; Res'=Res/10edx

@@ThereAintPoint:
  ; Учесть знак числа
  or   dword ptr [ebx],edi
  clc
  leave
  ret  2*4
@@ErrorSign:
  stc
  leave
  ret  2*4
StringToVal endp
; ××××× ValToString: Подпрограмма форматирования числа типа Real в строку ××××××
; [ebp+8] - Адрес строки(будет Zero-terminated) в виде 12345.89
; [ebp+0Сh] - Адрес числа типа real (или dd)
; [ebp+10h] - Маска вывода числа (byte)
ValToString proc
 .code
  push ebp
  mov  ebp,esp
  ; Локальные переменные:
  ; [ebp-4] - модуль числа
  ; [ebp-08] - максимальная степень 10-ти в числе
  ; [ebp-0Ch] - временная
  sub  esp,4*3
  cld
  mov  edi,dword ptr [ebp+8] ; Адрес выходной строки
  mov  ebx,[ebp+0Ch]    ; Загружаем адрес числа в ebx
  mov  edx,[ebx]        ; Читаем его в edx
  ; Определяем знак, удаляя его одновременно (31 бит)
  shl  edx,1
  jnc  @@NoSign
  mov  al,'-'
  stosb
@@NoSign:
  shr  edx,1
  mov  [ebp-4],edx      ; Запоминаем модуль числа в локальной переменной
  ; Проверка на ноль
  test edx,edx
  jnz  @@NonZero
  mov  al,'0'
  stosb
  jmp  @@DoneValToStr
@@NonZero:
  ; Определяем старшую степень десяти - Log10(x) = Log2(x)/Log2(10)
  ; выполняется в два приема, поскольку нет операции Log10 в FPU,
  ; а есть y×Log2(x)
  ; Получаем 1/Log2(10)
  fld1                  ; Загрузить 1.0
  fldl2t                ; Загрузить log2(10)
  fdivp ST(1),ST(0)     ; Разделить 1/log2(10), результат в ST(0)
  fld  dword ptr [ebp-4] ; Загрузить наше число, 1/log2(10) протолкнуть в ST(1)
  fyl2x ; ST(0)=ST(1)×log2(ST(0))
  ; Задаем режим округления (биты 11 - отбросить дробную часть)
  call Set_FPUTruncType
  frndint
  fistp dword ptr [ebp-8] ; выгрузить степень как целое со знаком
  ; Задаем режим округления (биты 00 - округлить до ближайшего целого)
  call Set_FPURoundType
  ; Сохранить степень числа
  push dword ptr [ebp-8]
  ; вывести число, деля его на 10...
  ; вывести целую часть числа, точку, дробную часть
  mov  ecx,dword ptr [ebp+10h] ; Сколько цифр после десятичной точки выводить
  inc  ecx              ; Цифра перед десятичной точкой
  xor  edx,edx          ; Флаг вывода точки
@@OutNumbers:
  fld  dword ptr DecNumber ; Загрузить 10 - основание
  mov  eax,dword ptr [ebp-8] ; Степень
  call get_pow          ; ST(0) = 10Степень
  fld  dword ptr [ebp-4] ; Загрузить число
  ; ST(0)=Число, ST(1) = 10Степень

  fdiv ST(0),ST(1)      ; ST(0)=ST(0)/ST(1)
  ; Задаем режим округления (биты 11 - отбросить дробную часть)
  call Set_FPUTruncType
  frndint
  fistp dword ptr [ebp-0Ch] ; выгружаем как целое без округления
  ; Задаем режим округления (биты 00 - округлить до ближайшего целого)
  call Set_FPURoundType
  ; выводим цифру
  mov  eax,dword ptr [ebp-0Ch]
  mov  ebx,1
  call DecChar          ; В формате - @n1
  inc  edi
  ; выводим точку, если это необходимо
  test edx,edx
  jnz  @@PointAlreadyOut
  mov  al,'.'
  stosb
  inc  edx
@@PointAlreadyOut:
  ; вычесть из числа 10Степень(сейчас в ST(0)) x Текущая цифра
  fimul dword ptr [ebp-0Ch]
  fld  dword ptr [ebp-4] ; Загрузить число
  fsub ST(0),ST(1)
  fstp dword ptr [ebp-4] ; выгрузить число с вычетом
  fstp ST(0)             ; Очистить стек
  ; Уменьшить степень на 1
  dec  dword ptr [ebp-8]
  loop @@OutNumbers
  ; Восстановить степень числа
  pop  dword ptr [ebp-8]
  ; вывести степень 10-ти
  mov  ax,'+E'
  mov  edx,[ebp-8]
  cmp  edx,0
  jge  @@NoNegPower
  neg  edx
  mov  ah,'-'
@@NoNegPower:
  stosw
  mov  eax,edx
  mov  ebx,1000
  call DecChar           ; В формате @n4
  add  edi,4
@@DoneValToStr:
  mov  byte ptr [edi],0
  leave
  ret  3*4
ValToString endp
; ××××××××××××× ReadFileString: Подпрограмма чтения строки из файла ××××××××××××
; [ebp+8] - Дескриптор файла
; [ebp+0Ch] - Адрес приемника строки
; [ebp+10h] - размер приемника
; -> eax реальный размер строки в приемнике и ноль в конце приемника
; Разделителями считаются символы 0A или 0D,0A, пробел.
; Пробелы слева читаются до любого другого символа,
; затем считываются символы числа, пробел служит знаком конца строки
ReadFileString proc
 .code
  push ebp
  mov  ebp,esp
  ; Переменная для байта из файла размером в 1 байт [ebp-8]
  ; Число байт, прочтенных ReadFile [ebp-4]
  ; Число байт в строке [ebp-0Ch]
  sub  esp,3*4
  cld
  mov  edi,[ebp+0Ch]
  mov  dword ptr [ebp-0Ch],0
  ; Флаг чтения любого символа, кроме пробела
  xor  esi,esi
@@NextByte:
  push NULL
  lea  eax,[ebp-4]
  push eax
  push 1
  lea  eax,[ebp-8]
  push eax
  push dword ptr [ebp+8]
  call ReadFile
  cmp  dword ptr [ebp-4],0
  jz   @@EndRead
  mov  al,byte ptr [ebp-8]
  cmp  al,0Ah
  jz   @@EndRead
  cmp  al,0Dh
  jz   @@NextByte        ; Пропускаем байт 0Ah
  cmp  al,' '            ; Пробел?
  jnz  @@NoFiller
  test esi,esi
  jnz  @@EndRead         ; Пробел после строки - break
  jmp  @@AfterFiller     ; Пробелы перед строкой заносим в буффер
@@NoFiller:
  inc  esi               ; Установим флаг, что был символ <> пробел
@@AfterFiller:
  mov  ebx,dword ptr [ebp-0Ch]
  inc  ebx
  cmp  ebx,dword ptr [ebp+10h]
  jae  @@EndRead
  inc  dword ptr [ebp-0Ch]
  stosb
  jmp  @@NextByte
@@EndRead:
  mov  byte ptr [edi],0
  mov  eax,dword ptr [ebp-0Ch]
  leave
  ret  3*4
ReadFileString endp
; ××××××××××××× DecChar: Подпрограмма форматирования строки из числа ×××××××××××
; eax=number to digit, edi=offset result string in format 00000000(@n[ebx])
; ebx=начальный делитель
DecChar proc
  pushad
  pushfd
  cld
@GetDec:
  xor  edx,edx
  div  ebx
  add  al,'0'
  stosb
  push edx
  mov  eax,ebx
  xor  edx,edx
  mov  ebx,10
  div  ebx
  mov  ebx,eax
  pop  eax
  test ebx,ebx
  jnz  @GetDec
  popfd
  popad
  ret
DecChar endp
; ××××××××××××××××××××××× Str_Len: Get lenght of string ×××××××××××××××××××××××××
; esi=addr of string; result: ecx=length
Str_Len proc
  xor  ecx,ecx
  push eax
  push esi
@@GetStrLen:
  lodsb
  test al,al
  jz   @@ExitStrLen
  inc  ecx
  jmp  @@GetStrLen
@@ExitStrLen:
  pop  esi
  pop  eax
  ret
Str_Len endp
end start

Оценим работоспособность алгоритма на его способности аппроксимировать Y(X) полиномом вида Y(x) = a2x2 + a1x+a0 (параболой): найти неизвестные коэффициенты a2, a1, a0.

Для примера рассмотрим простую параболу:

    Y(x) = 3x2 - 5x -8;

на интервале значений x={1..10}:

      X     Y
     1.0 -10.00
     2.0  -6.00
     3.0   4.00
     4.0  20.00
     5.0  42.00
     6.0  70.00
     7.0 104.00
     8.0 144.00
     9.0 190.00
    10.0 242.00

Поместим эту зависимость Y(x) в файл входных данных программы (splain.dat), указав предполагаемую степень полинома первой строкой:

    Содержимое входного файла splain.dat:
    2
     1.0 -10.00
     2.0  -6.00
     ...
    10.0 242.00

В результате программа должна выдать следующий результат (содержимое MessageBox'а):

    Assume that you data was approximated by polinom with power 0002

    a[0]= -7.282...
    a[1]= -5.288...
    a[2]= 3.0247...

    Xq = 0.413... for 00010 points.

Как видно, программа достаточно точно смогла определить коэффициенты параболы. Заметим, однако, что несмотря на то, что входные данные являлись абсолютно точно параболой с опредленными коэффициентами, они не были определены на 100%. Почему так происходит? Дело в том, что критерий Xq выполняет подбор коэффициентов по величине квадратов отклонения точек экспериментальных данных от теоретических. Очевидно, что для небольших значений Y приближение будет менее точным, нежели чем для больших: см. например, то, что значение a[2] определено гораздо точнее, чем a[1] или тем более a[0].

Оценка предложенного метода численной аппроксимации с помощью Xq.

В результате проведенного теста оказывается, что предложенный способ достаточно точно позволяет определять параметры аппроксимирующей зависимости, по крайней мере позволяет выполнять предварительные оценки свойств экспериментальных зависимостей. Вместе с этим он претендует на некоторую универсальность алгоритма относительно выбора вида и параметров аппроксимирующей зависимости: она может быть практически сколь угодно сложной и иметь широкий спектр параметров; алгоритм не требует сложных математических вычислений.

(C) Chingachguk /HI-TECH

#################################################################################

masquer

Фильтр для Photoshop

  1. Предисловие
  2. Пишем фильтр
    1. Ресурсы
    2. Структура FilterRecord и написание фильтра
  3. Заключение

Предисловие

Довольно часто слышу мнение, что ассемблер мертвый язык, что он никому не нужен, за исключением особых ситуаций и т.д. и т.п. Хочу надеяться, что эта статья будет неплохим гвоздем в гроб этих мненийеще одним гвоздем в крышку гроба этих мнений, т.к. в этом рассмотренном ниже случае его применение более чем оправдано. Как говорил Остап Бендер: «Слухи о моей смерти явно преувеличены». Тоже самое может сказать и за ассемблер, который периодически закапывают не уставая все кому не лень со времени его появления. И все-таки он жив. Жив и активно развивается. Убедитесь сами!

Пишем фильтр

Являясь большим поклонником как замечательного продукта Adobe Photoshop, так и не меньшим (если не большим) поклонником ассемблера, мне подумалось - почему бы не усилить эффект и соединить два в одном - так, собственно, и родилась идея написания плагина-дополнения для Photoshop, и непременно на ассемблере. Порывшись в интернете и скачав Photoshop 6.0 SDK, я обнаружил практически полное отсутствие информации об этом, а платить Adobe за участие в их форумах не хотелось. Но кое-какую информацию найти все же удалось - совершенно случайно на сайте WhizKid-а я обнаружил пример фильтра для Photoshop, но написан он был достаточно громоздко с использованием синтаксиса nasm. Пришлось заняться рассмотрением примеров и углубиться в чтение документации, содержащимися в SDK, благо - все оказалась достаточно информативным.

Что удалось выяснить. Photoshop изначально создавался как продукт для Apple Macintosh, соответственно первые плагины и интерфейс взаимодействия программы-хоста (Photoshop) и плагина создавался, учитывая специфику яблочных компьютеров. Отголоски этого при написании плагина в среде Windows состоят в специфическом формате ресурсов PiPL. Но об этом формате немного позже. Еще одним наследием Apple является то, что данные и указатели для Photoshop формируются в формате big endian (в отличие от little endian у процессоров Intel), но я пока не замечал этого влияния, так что об этом пока достаточно просто упомянуть (с этим столкнутся те, кто будет писать плагины импорта и экспорта).

2.1 Ресурсы

Формат ресурсов для фильтров Photoshop отличается от формата ресурсов программ для Windows. Остановимся на его основных элементах подробнее.

Kind { Filter } - тип плагина, возможные варианты - "Filter", "Parser", "ImageFormat", "Extension", "Acquire", "Export". Соответственно у нас тип плагина - фильтр.
Name { "masquer" } - имя в меню Plugins, здесь только фантазия разработчика.
Category { "masquer test" } - имя в субменю, здесь будет masquer->masquer test.
Version { (1 << 16) | 0 } - версия плагина, здесь 1.0.
CodeWin32X86 { "PluginMain" } - точка входа в наш плагин.
SupportedModes - режимы изображений, которые мы будем поддерживать, у нас будет так:
noBitmap, doesSupportGrayScale, noIndexedColor, noRGBColor, noCMYKColor, noHSLColor, noHSBColor, noMultichannel, noDuotone, noLABColor
EnableInfo - строка, при выполнении условий которой, наш плагин будет находиться в положении Enable.
Например,
"in (PSHOP_ImageMode, GrayScaleMode)"
Функция in возвращает истинное значение только тогда, когда первый параметр соответствует хотя бы одному из следующих. В данном случае мы имеем истинное значение, если изображение находится в режиме градаций серого.
FilterCaseInfo - массив из 7 элементов, каждый из которых состоит из четырех байт. Первые два значения определяют действия хоста для пре- и постпроцессинга. Т.е. определив в нашем случае inStraightData, outStraightData мы говорим хосту (Adobe Photoshop) что мы хотим получить наши данные немаскированными, и отдадим их с тем же условием, что хост не будет их демаскировать (dematte). Оставшиеся элементы составляют комбинацию флагов:
  • doNotWriteOutsideSelection - должен ли фильтр писать за границы выделения
  • filtersLayerMasks - можно ли обрабатывать маски слоев
  • worksWithBlankData - может ли фильтр работать с абсолютно прозрачным выделением/изображением
  • copySourceToDestination - будут ли копироваться исходные данные для фильтрования в другое место и уже там обрабатываться

Полный список функций и параметров можно найти в файле Plug-in Resource Guide.pdf.

Теперь, когда у нас готов файл с ресурсами (расширение *.r), для пишущих под Windows необходимо привести этот формат к стандартному. Для этого с Photoshop 6.0 SDK идет утилита Cnvtpipl.exe. В качестве параметра нужно просто указать наш файл с ресурсом и все - у нас уже есть файл с "нормальными" ресурсами (*.rc).

2.2 Структура FilterRecord и написание фильтра

Собственно теперь мы можем приступать к непосредственному написанию кода нашего плагина. Но перед этим рассмотрим одну очень важную для нас структуру. А именно речь пойдет о структуре FilterRecord. Всю структуру я рассматривать не буду, она слишком большая, да и нет в этом смысла, я начну, а те кому это нужно будет, уже сами смогут разобрать остальные поля. Итак:

serialNumber dd 0 В предыдущих версиях здесь был серийный номер Photoshop-а и плагин мог применять его для copy-protection. Не знаю как в 6, а в 7 версии здесь 0
abortProc dd 0 Адрес процедуры TestAbort. Не интересно
progressProc dd 0 Адрес процедуры UpdateProgress. Для создания ProgressBar. Обойдемся без него.
parameters dd 0 Параметры, задаваемые пользователем. Для нас неинтересно. В начале содержит 0
imageSize pPOINT<> Размеры изображения. Структура из двух word-ов. Отсюда вытекает, что изображение не может быть больше, чем 65535х65535 пикселей
planes dw 0 количество каналов, например для CMYK - 4, для RGB - 3, мы упростим себе задачу и выберем GrayScale - 1.
filterRect pRECT<> Общий фильтруемый прямоугольник. Даже если мы фильтруем неровное выделение, все равно нам будет передан граничащий прямоугольник, внутри которого располагается выделение. Структура из 4 word-ов: top, left, bottom и right соответственно.
background RGBColor<> Цвет подложки (background)
foreground RGBColor<> Цвет переднего плана (foreground)
  dw 0 Для выравнивания внутри структуры
maxSpace dd 0 Максимально возможный объем данных...
bufferSpace dd 0 ...и буфера
inRect pRECT<> Прямоугольник, который будет передан фильтру в текущей итерации. Для уменьшения расхода памяти рекомендуется условно разбить изображение на прямоугольники типа 128х128, или 256х256 и "кусками" обрабатывать изображение. Мы этим пренебрежем сейчас и закажем все изображение.
inLoPlane dw 0 Первый запрашиваемый канал...
inHiPlane dw 0 ...и последний
outRect pRECT<> Прямоугольник для результирующего изображения. У нас inRect и outRect будут совпадать
outLoPlane dw 0 так же, как и для inRect
outHiPlane dw 0  
inData dd 0 адрес, где находится исходное изображение для обработки
inRowBytes dd 0 смещение между строками изображения. Дело в том, что если ширина изображения не кратна 32, а адрес каждой новой строки должен быть кратен 32. Соответственно в таких случаях остаток до границы строки заполняется 0 и не учитывается при выводе результата.
outData dd 0 адрес, где мы можем разместить результат
outRowBytes dd 0 то же смещение

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

Основной экспортируемой нашим плагином функцией является PluginMain. Экспортируемой, потому что наш плагин не что иное, как библиотека dll, имеющая расширение 8bf. В качестве входящих параметров мы имеем selector, адрес filterRecord, data, result. Начнем с конца, как с наименее важных параметров:

  1. result - результат выполнения, некоторые ошибки могут нами обрабатываться, некоторые - хостом. Будем считать, что мы все делаем без ошибок и памяти у нас достаточно, поэтому принимать во внимание этот параметр не будем.
  2. data - здесь можно хранить хендл к нашим структурам данных. Обойдемся.
  3. filterRecord - адрес нашей главной структуры.
  4. selector - здесь остановимся поподробнее. Именно через этот селектор хост (Photoshop) сообщает нам, что нам нужно делать, передавая нам одно из пяти значений. Если селектор равен 0, то о нашем фильтре хотят получить больше информации, а именно: выбрали наш плагин в Help->About Plug-In. Adobe, правда, хочет, чтобы мы выводили не просто месадж бокс, а создавали диалоговое окно без кнопок, но с надписями, ну да переживут как-нибудь.
    align 4
    DoAbout proc
      invoke MessageBox, 0, ADDR szText, ADDR szText, MB_OK
      ret
    DoAbout endp

Параметром номер 1 является filterSelectorParameters, где по идее мы должны обработать передаваемые нам параметры и отразить диалог, мы упростим все до одного ret.

    align 4
    DoParameters proc
      ret
    DoParameters endp

2 - filterSelectorPrepare - здесь мы должны посчитать и выделить память, а можно и нули занести, тогда Photoshop сам это за нас сделает. Так и решим.

    align 4
    DoPrepare proc
      xor eax, eax
      mov fr.bufferSpace, eax
      mov fr.maxSpace, eax
      ret
    DoPrepare endp

3 - filterSelectorStart - здесь мы уже обрабатываем изображение

    align 4
    DoStart proc uses edi esi
      mov edx, dword ptr fr
      assume edx:PTR FilterRecord ; edx указывает на FilterRecord
      movzx eax, [edx].imageSize.h ; получаем размеры изображения
      movzx ebx, [edx].imageSize.v

      mov word ptr [edx].inRect.top, 0 ; укажем хосту что именно мы хотим обрабатывать
      mov word ptr [edx].inRect.left, 0
      mov word ptr [edx].inRect.bottom, bx
      mov word ptr [edx].inRect.right, ax

      mov word ptr [edx].outRect.top, 0 ; здесь inRect = outRect
      mov word ptr [edx].outRect.left, 0
      mov word ptr [edx].outRect.bottom, bx
      mov word ptr [edx].outRect.right, ax

      mov eax, [edx].advanceState ; пускай хост обновит параметры,
      call eax  ; хотя нам в данном случае это и не нужно. Это необходимо при
          ; обработке изображения "кусками"

      mov edx, dword ptr fr
      assume edx:PTR FilterRecord
      mov esi, dword ptr [edx].outData ; так как у нас inRect = outRect,
               ; то esi указывает на данные для обработки
      mov edi, esi ; и edi тоже

      movzx eax, [edx].imageSize.h ; рассчитаем количество итераций цикла обработки
      mov ebx, eax
      shr eax, 4
      test ebx, 0Fh
      jz @@_1
      inc eax
    @@_1:
      shl eax, 4
      movzx ebx, [edx].imageSize.v
      mul ebx
      mov ecx, eax
      shr ecx, 4 ; наконец-то рассчитали
    ; дальше идет основной цикл вычислений.
    ; Чтобы особо не выдумывать я сдвигаю каждое значение яркости на 1 бит влево используя MMX.
    ; Но об этом в другой раз.
    @@next:
      movq mm0, [esi]
      movq mm1, [esi+8]

      movq mm2, mm0
      movq mm3, mm1

      psllw mm0, 9
      psrlw mm2, 8
      psllw mm2, 9
      psrlw mm2, 8

      psllw mm1, 9
      psrlw mm3, 8
      psllw mm3, 9
      psrlw mm3, 8

      por mm0, mm2
      por mm1, mm3

      movq [esi], mm0
      movq [esi+8], mm1
      add esi, 16
      dec ecx
      jz @@fin
      jmp @@next
    ;end of main
    @@fin:
      emms ; вдруг кому-то еще с сопроцессором работать нужно
      mov edx, dword ptr fr
      assume edx:PTR FilterRecord
      mov dword ptr [edx].inRect.top, 0 ; скажем хосту - все готово
      mov dword ptr [edx].inRect.bottom, 0
      mov dword ptr [edx].outRect.top, 0
      mov dword ptr [edx].outRect.bottom, 0
      mov dword ptr [edx].maskRect.top, 0
      mov dword ptr [edx].maskRect.bottom, 0
      ret
    DoStart endp

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

    align 4
    DoContinue proc
      ret
    DoContinue endp

5 - filterSelectorFinish - финиш - он и есть финиш - очистка ресурсов, а так как никаких ресурсов мы не занимали, то и чистить нам нечего.

    align 4
    DoFinish proc
      ret
    DoFinish endp

Теперь посмотрим, как все это вызвать, чтобы работало. Собственно, процедура PluginMain

     .data
      szText db "masquer's Photoshop Plugin",0
      procs dd DoAbout, DoParameters, DoPrepare, DoStart, DoContinue, DoFinish 
          ;LUT (LookUp Table)
      align 4
      fr FilterRecord <>
      outRect pRECT <>
      imageSize pPOINT <>
     .code
     ...
    align 4
    PluginMain proc uses ebp edi ecx, selector:DWORD,
         filterRecord:DWORD, data:DWORD, result:DWORD
      mov ecx, result ; нулевой результат нам не нужен
      jcxz @@exit
      mov eax, filterRecord
      lea eax, [eax]
      mov dword ptr [fr], eax ; настроим указатели
      mov eax, selector
      cmp eax, 5
      ja @@exit
      call [procs+eax*4] ; вызовет процедуру в соответствии с селектором
    @@exit:
      mov esp, ebp ; подчистим стек
      pop ebp
      retn
    PluginMain endp
     ...

Ну, и, как и у любой нормальной dll, пропишем DllEntry

    align 4
    DllEntry proc hInstDLL:DWORD, reason:DWORD, unused:DWORD
      mov eax, 1
      ret
    DllEntry Endp

Для отладки я использовал незаменимый во всех случаях жизни SoftIce, а точнее я цеплялся к MessageBoxA, смотрел хелп, ну а дальше уже дело техники. Напоследок приведу только содержимое .def файла

    LIBRARY test.8bf
    EXPORTS
    PluginMain

Да, чуть не забыл. Для того чтобы все это заработало, достаточно скопировать полученный test.8bf в папку к плагинам [Здесь путь к каталогу Photoshop]\Plug-Ins\test.8bf.

Заключение

Как говорится - вуаля! Заготовка полностью написана, теперь только знание математики и ее применение к обработке изображений сдерживает наш творческий порыв. Но к следующему разу я что-нибудь обязательно придумаю, про MMX, например, попробуем рассказать, что это такое и где его можно применить.

© masquer, 2002, masquer@pochta.ws

##########################################################################

Подготовленно к выпуску 20.11.2002
HI-TECH group

Выпуск компилировали и билдили:

  • Составители: Serrgio, Edmond
  • Главный комментатор: FatMoon
  • HTML вёрстка: Masquer
  • Редактор рассылки: страшный маньяк CiberManiac
  • Техническая поддержка Aquila
  • WASM.RU