Написание программы с SSE. Перевод

@SeryZone · Регистрация: 09.03.2012

Студворк — интернет-сервис помощи студентам

Недавно строил классическое множество Мандельброта, использовал ранее созданные алгоритмы другими людьми, всё строило хорошо. Но т.к. Я вывел собственную (несложную) формулу для построения аналога, передо мной встала большая проблема: на SSE я не профи, да и с алгоритмом парится - ушло на построение очень много времени, а код я порвал в клочья.
Ближе к теме. Формула классического фрактала - Z = Z²+C
Моя формула - Z = Z² + C + CENTER
На простейшем языке Delphi код выглядел так:

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
function cAbs(z: TComplex): extended; register;
// z :--> |z|
asm
         FLD TComplex.im [EAX]
         FMUL ST, ST
         FLD TComplex.re [EAX]
         FMUL ST, ST
         FADD
         FSQRT
end;
 
function cadd(z1, z2: TComplex): TComplex; register;  //сложение
asm
         FLD TComplex.re [EAX]
         FLD TComplex.re [EDX]
         FADD
         FSTP TComplex.re [ECX]
         FLD TComplex.im [EAX]
         FLD TComplex.im [EDX]
         FADD
         FSTP TComplex.im [ECX]
end;
 
function cSqr(z: TComplex): TComplex; register; //возведение в квадрат комплексного числа
// z :--> z^2
asm
         FLD TComplex.im [EAX]
         FLD ST
         FMUL ST, ST
         FLD TComplex.re [EAX]
         FLD ST
         FMUL ST, ST
         FXCH ST(3)
         FMUL
         FADD ST, ST
         FSTP TComplex.im [EDX]
         FSUB
         FSTP TComplex.re [EDX]
end;
 
function fractal_mandelbrot(c,center:tComplex):Longint;//вычисление
var w:Longint;
z:TComplex;
izf:Double;
begin
  z.re:=0;
  z.im:=0;
  for Result:=0 to iter do begin
    //возводим в квадрат
    z:=cSqr(z);
    //прибавляем оригинал
    z:=cadd(z, c);
    z:=Cadd(z, center);
    //выясняем, больше ли модуль, чем нужно?
    if(cabs(z) > 4)then exit;
  end;
end;
 
procedure Build_fractal()
begin
for iy := 0 to 1919 do
    begin
    for ix := 0 to 12400 do
    begin
      n := fractal_mandelbrot(zData[ix,iy],center);
      if(n >= iter)then
      begin
        NumIter[ix,iy]:=-1;
      end
      else
      begin
        NumIter[ix,iy]:=n;
      end;
end;

Суть задачи - включая циклы переписать весь этот код с функциями с SSE конструкциями (если возможно). Господа, подскажите, где можно почитать о построении алгоритма, и вообще, что можно сделать?

murderer · 13.09.2013, 09:39

extended имеет 112-битную мантиссу

extended - 80 бит, но Delphi округляет размер до 12 байт. 1 бит знака, 15 бит экспоненты и 64 бита мантиссы.

Добавлено через 3 минуты

эта мантисса должна позволять увеличивать до 2^112!

Вобще то так как порядок 15-битный, то диапазон значений [2^-16383;2¹⁶³⁸³], но точно выражаются только 2⁸⁰ чисел, раскиданных по всему диапазону.

Добавлено через 5 минут
Кстати если ещё не видел - загугли "Фрактальный игровой мир", там san base такие вещи вытворяет...

Добавлено через 16 минут

Delphi округляет размер до 12 байт.

Сейчас проверил - sizeof(extended)=10. Но точно помню, что где-то читал о 12 байтах, наверное это относится к Turbo Pascal.

но точно выражаются только 2⁸⁰ чисел

Тут присутствуют ещё и "особые числа" +0, -0, +∞, -∞, NAN.

@SeryZone · 13.09.2013, 15:38 **[ТС]**

murderer, так почему так... А SanBase молодцы, крутые вещи выделывают

Но... Как программировать новый тип? Фрактал екстрим - не самая быстрая программа тоже, но там побитовая точность... Блин... Как же программируют эти типы??? А если через длинную арифметику делать и массивы, максимум, чего мы сможем достичь - это скорости программы Ultra Fractal, который на произвольной точности вычисляет так медленно, что быстрее - об стену убицца, чем дождаться

@SeryZone · 13.09.2013, 23:30 **[ТС]**

murderer, сделал реверс этого "файлика" - библиотеки (в архиве) и... ничего не понял... Имено этот плагинчик позволяет Fractal eXtrem'у вычислять на высокой точности... ~~Чего эт я привязался к Fractal eXtreme...~~
Можно ли что-то выяснить?

@SeryZone · 14.09.2013, 08:25 **[ТС]**

Пользуюсь IDA Disassembler'ом. Там функции по именам sub. Более 100. Кстати, я посмотрел, что FX пользует около 120 кодов для разной точности. Блин.... Где бы найти толковую литературу по АСМ???

Добавлено через 8 часов 41 минуту

Assembler
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
sub_202045D6 proc near
mov     rax, [rsi]
mul     qword ptr [rsi+40h]
mov     r8, rdx
xor     r9, r9
mov     rax, [rsi+8]
mul     qword ptr [rsi+38h]
add     r8, rdx
adc     r9, 0
mov     rax, [rsi+10h]
mul     qword ptr [rsi+30h]
add     r8, rdx
adc     r9, 0
mov     rax, [rsi+18h]
mul     qword ptr [rsi+28h]
add     r8, rdx
adc     r9, 0
add     r8, r8
adc     r9, r9
mov     rax, [rsi+20h]
mul     rax
add     r8, rdx
adc     r9, 0
mov     rax, [rsi+8]
mul     qword ptr [rsi+40h]
mov     r11, rax
mov     r12, rdx
xor     r13, r13
mov     rax, [rsi+10h]
mul     qword ptr [rsi+38h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
mov     rax, [rsi+18h]
mul     qword ptr [rsi+30h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
mov     rax, [rsi+20h]
mul     qword ptr [rsi+28h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
xor     r10, r10
add     r8, r11
adc     r9, r12
adc     r10, r13
add     r8, r11
adc     r9, r12
adc     r10, r13
mov     [rbp+0], r8
mov     rax, [rsi+10h]
mul     qword ptr [rsi+40h]
mov     r11, rax
mov     r12, rdx
xor     r13, r13
mov     rax, [rsi+18h]
mul     qword ptr [rsi+38h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
mov     rax, [rsi+20h]
mul     qword ptr [rsi+30h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
xor     r8, r8
add     r9, r11
adc     r10, r12
adc     r8, r13
add     r9, r11
adc     r10, r12
adc     r8, r13
mov     rax, [rsi+28h]
mul     rax
add     r9, rax
adc     r10, rdx
adc     r8, 0
mov     [rbp+8], r9
mov     rax, [rsi+18h]
mul     qword ptr [rsi+40h]
mov     r11, rax
mov     r12, rdx
xor     r13, r13
mov     rax, [rsi+20h]
mul     qword ptr [rsi+38h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
mov     rax, [rsi+28h]
mul     qword ptr [rsi+30h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
xor     r9, r9
add     r10, r11
adc     r8, r12
adc     r9, r13
add     r10, r11
adc     r8, r12
adc     r9, r13
mov     [rbp+10h], r10
mov     rax, [rsi+20h]
mul     qword ptr [rsi+40h]
mov     r11, rax
mov     r12, rdx
xor     r13, r13
mov     rax, [rsi+28h]
mul     qword ptr [rsi+38h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
xor     r10, r10
add     r8, r11
adc     r9, r12
adc     r10, r13
add     r8, r11
adc     r9, r12
adc     r10, r13
mov     rax, [rsi+30h]
mul     rax
add     r8, rax
adc     r9, rdx
adc     r10, 0
mov     [rbp+18h], r8
mov     rax, [rsi+28h]
mul     qword ptr [rsi+40h]
mov     r11, rax
mov     r12, rdx
xor     r13, r13
mov     rax, [rsi+30h]
mul     qword ptr [rsi+38h]
add     r11, rax
adc     r12, rdx
adc     r13, 0
xor     r8, r8
add     r9, r11
adc     r10, r12
adc     r8, r13
add     r9, r11
adc     r10, r12
adc     r8, r13
mov     [rbp+20h], r9
mov     rax, [rsi+30h]
mul     qword ptr [rsi+40h]
xor     r9, r9
add     r10, rax
adc     r8, rdx
adc     r9, 0
add     r10, rax
adc     r8, rdx
adc     r9, 0
mov     rax, [rsi+38h]
mul     rax
add     r10, rax
adc     r8, rdx
adc     r9, 0
mov     [rbp+28h], r10
mov     rax, [rsi+38h]
mul     qword ptr [rsi+40h]
xor     r10, r10
add     r8, rax
adc     r9, rdx
adc     r10, 0
add     r8, rax
adc     r9, rdx
adc     r10, 0
mov     [rbp+30h], r8
mov     rax, [rsi+40h]
mul     rax
add     r9, rax
adc     r10, rdx
adc     r8, 0
mov     [rbp+38h], r9
mov     [rbp+40h], r10
retn

@snake32 · 16.09.2013, 04:51

Вот выкладываю свой очередной недоделок с исходниками. Мальденброт GLSL. Написан на Delphi XE2.
Мышой(левая кнопка) передвигается комплексная плоскость. Также можно увеличивать/уменьшать масштаб колесом мыши. Только увеличивает он центр экрана, а не место под мышой. Если кто знает подскажите как сделать нормальный зум?
Сначала делал расчёт на CPU без оптимизаций(справа осталась панель, но она не работает), но потом решил задействовать шейдеры так как понял что алгоритм отлично параллелится. На моей HD 7970 GHz Edition всё летает, несмотря на то что во фрагментном шейдере цикл. Как будет на других картах не знаю. Возможно печалька.

@SeryZone · 16.09.2013, 08:26 **[ТС]**

snake32, у меня 5870, но я на делфи разбирал эту тему. Посмотреть бы...

@snake32 · 16.09.2013, 08:55

Сообщение от SeryZone

snake32, у меня 5870, но я на делфи разбирал эту тему. Посмотреть бы...

Дык скачивай и смотри, там как раз на делфи. В чём проблема?

Сейчас проверил на рабочем компе древняя nVidia GTS 250 когда на весь экран разворачиваешь, чувствуется, что туго ей даётся отрисовка.

Добавлено через 9 минут
Или архива не видно вам? Может потому что я с нарушением?

@FIL · 16.09.2013, 08:58

Сообщение от snake32

Если кто знает подскажите как сделать нормальный зум?

Вот реализация на SSE с зумом.

murderer · 16.09.2013, 09:46

snake32
А почему "радуга.bmp" загружается как 2D-текстура? Правильнее будет 1D

@snake32 · 16.09.2013, 09:48

Сообщение от FIL

Вот реализация на SSE с зумом.

У меня конечно не айс рабочий комп Intel Q9450( 4core; 2.66GHz ) но мне кажется что даже GPU GTS 250 опережает в разы CPU. В оправдание этому могу сказать что на GPU всего навсего обычный float а не double.
К сожалению эта карта вообще не способна вычислять Double Precision.

murderer · 16.09.2013, 09:48

Вершинный шейдер не нужен.

@snake32 · 16.09.2013, 10:04

Сообщение от murderer

А почему "радуга.bmp" загружается как 2D-текстура? Правильнее будет 1D

Оно и так не плохо работает

Добавлено через 14 минут

Сообщение от murderer

Вершинный шейдер не нужен.

Как так? Я не умею по другому

@FIL · 16.09.2013, 15:26

Сообщение от snake32

У меня конечно не айс

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

@snake32 · 16.09.2013, 17:40

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

murderer · 17.09.2013, 09:40

А вот и моя версия шейдера. Теоретически формула Z²+C должна скомпилироваться в 2 инструкции MAD. Масштабирование колёсиком, можно таскать изображение левой кнопкой мыши.

@snake32 · 17.09.2013, 12:29

Сообщение от murderer

А вот и моя версия шейдера.

У меня стабильно валится после нескольких изменений камеры (HD 7970)

Добавлено через 2 часа 10 минут
Проверил на рабочем компе - всё круто!( GTS 250 )

murderer · 17.09.2013, 15:21

Думаю к вылетам может приводить этот код (при переполнении экспоненты).

Assembler
1
2
3
4
5
6
7
8
fild      dword[esp-4]
fld       [Scale]
fscale
fstp      [Scale]
fld       [Scale+8]
fscale
fstp      [Scale+8]
fstp      st0

Добавлено через 36 секунд
То есть при слишком большом приближении или отдалении.

@snake32 · 17.09.2013, 15:26

Сообщение от murderer

То есть при слишком большом приближении или отдалении.

Нет, буквально 2-4 шага зума, иногда просто при изменении(лево/право/верх/низ).
Может можно как-то настроить чтобы ошибки не генерировали вылет приложения? Типа как пустой try except end в delphi.

murderer · 18.09.2013, 08:26

snake32
А что если радугу генерить прямо в шейдере?

Вот что у меня получилось.

Code
1
2
3
4
5
void main(void)
{
    float i = gl_FragCoord.x/640.0;
    gl_FragColor = vec4(sin(i*3.14*1.75+3.14*0.5), sin(i*3.14*1.5), sin(i*3.14*1.75-3.14), 0.0)*4.8;
}

Но через синус результат получается слишком размытым - нужна пилообразная функция. Есть идеи?

murderer · 18.09.2013, 08:47

Всё разобрался - там множитель слишком большой был. Вместо 4,8 надо 1,5.

Новые блоги и статьи Все статьи Все блоги /
Саморегулирующийся социальный контракт для сервера cross-section. Hrethgir 14.08.2026 С кодом конечно таких глубоких размышлений пока не было, впрочем я уже привык к алгоритмизации. Суть предмета записи: снова в диалоге с нейросетью (я взял пока себе ник для учётки админа - Rector). . . .	Часы электронные Uhbif79 12.08.2026 Выкладываю программу часов. Программа позволяет: 1. Использовать системное время и дату, 2. Есть возможность вводить время и дату вручную. 3. Реализованы 2 будильника: начало и конец рабочего дня. . . .	Часы с будильником на основе класса QLCDNumber Uhbif79 12.08.2026 Всем добрый день, выкладываю программу часов с будильником на основе класса QLCDNumber. Здесь я пробовал самостоятельно создавал классы, впервые столкнулся с видимостью переменной одного класса из. . .	Установка MinGW GCC 16.2 и CMake 8Observer8 10.08.2026 VK Видео: https:/ / vkvideo. ru/ video-240781534_456239017 YouTube: eY5-5PyI9NM Текстовая версия
Неделя из жизни имитационной модели склада: мои кривые руки растут, откуда надо anaschu 10.08.2026 Неделя из жизни имитационной модели склада: как я почти написал неправильную логику и что с этим делать Работаю сейчас над учебно-рабочим проектом: строю в AnyLogic имитационную модель процессов. . .	Калькулятор для расчета родства russiannick 07.08.2026 1. Задача: Создать калькулятор для расчета родства. Родственных связей существует 8 ступеней, такие как: p - отец P - мать q - муж Q - жена b - брат B - сестра s - сын S - дочь	Мир по моей воле kumehtar 07.08.2026 Когда-то кажется, что всё просто. Ты весь такой светлый. Причиняешь добро. Борешься за справедливость в этом тёмном мире. Потом начинаешь замечать одну неприятную вещь. Почти каждый хороший. . .	Кредитный калькулятор Maks 05.08.2026 Решение задачи по прикладной информатике средствами 1С. Задача: Напишите приложение-калькулятор, которое помогает рассчитывать параметры кредита для аннуитетного и дифференцированного видов. . .

@SeryZone 56 / 28 / 18 Регистрация: 09.03.2012 Сообщений: 726 Записей в блоге: 1
	13.09.2013, 15:38 [ТС]
	murderer, так почему так... А SanBase молодцы, крутые вещи выделывают Но... Как программировать новый тип? Фрактал екстрим - не самая быстрая программа тоже, но там побитовая точность... Блин... Как же программируют эти типы??? А если через длинную арифметику делать и массивы, максимум, чего мы сможем достичь - это скорости программы Ultra Fractal, который на произвольной точности вычисляет так медленно, что быстрее - об стену убицца, чем дождаться 0

@SeryZone 56 / 28 / 18 Регистрация: 09.03.2012 Сообщений: 726 Записей в блоге: 1
	16.09.2013, 08:26 [ТС]
	snake32, у меня 5870, но я на делфи разбирал эту тему. Посмотреть бы... 0

murderer 4195 / 1841 / 223 Регистрация: 06.10.2010 Сообщений: 4,127
	16.09.2013, 09:46
	snake32 А почему "радуга.bmp" загружается как 2D-текстура? Правильнее будет 1D 0

murderer 4195 / 1841 / 223 Регистрация: 06.10.2010 Сообщений: 4,127
	16.09.2013, 09:48
	Вершинный шейдер не нужен. 0

@snake32 3590 / 1720 / 236 Регистрация: 26.02.2009 Сообщений: 8,714 Записей в блоге: 5
	16.09.2013, 17:40
	FIL, спасибо, это я понял. Позже может быть приведу прогу в божеский вид и добавлю double... 0

murderer 4195 / 1841 / 223 Регистрация: 06.10.2010 Сообщений: 4,127
	18.09.2013, 08:47
	Всё разобрался - там множитель слишком большой был. Вместо 4,8 надо 1,5. 0