Написание программы с 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.

Новые блоги и статьи Все статьи Все блоги /
Как надо было делать правильно - вообще не слушать ИИ, его дело - текст. Hrethgir 15.07.2026 А вообще надо было делать не так, чтобы не исправлять ошибку https:/ / www. cyberforum. ru/ blogs/ 223907/ 10996. html , правильно было : удалил кнопку, а потом добавил её. Всё это стресс конечно. Почему?. . .	Знаете, от чего появляется стресс? kumehtar 15.07.2026 Стресс появляется оттого, что вы, зная как поступить правильно, делаете наоборот. Ларри Уингет	Нас бьют, а мы крепчаем? Какая чушь! kumehtar 15.07.2026 Вот это вот бытующее мнение, что человек меняется только под действием проблем, всегда считал вредоносным. Поверь, от проблем у людей появляются не изменения, а нервное расстройство и ПТСР. А. . .	сукцессия 29. Переход от одних деревьев на другие делать более или менее вероятностным? anaschu 12.07.2026 Насколько смена типов микоризы — исключительное событие в двухвековой сукцессии? Оценка вероятности в пространстве параметров В текущей версии модели успешно реализован ключевой механизм. . .
сукцессия 27. Думаю, как переделывать уже написанную статью с планами на сукцессию. anaschu 12.07.2026 Анализ соответствия модели требованиям Реализованные компоненты: Механизм закисления почвы через протонную помпу Конкуренция между типами микориз pH как триггер сукцессии C/ P соотношение. . .	Сукцессия 26. Мат модель создана. anaschu 12.07.2026 Модель смены растительных сукцессий посредством управления грибами работает внутри небольшой ячейки почвы, восстанавливающейся после пожара, где ненадолго бывшее царство хвойных снова захватили. . .	Решил проблему с ошибкой пагинации сообщений с сервера на алгоритме обхода дерева "Эстафета хвоста". Hrethgir 12.07.2026 Проблема была в том, что удалялась именно новая кнопка, а не старая. Ни один ИИ не обнаружил это, а сам я смог только когда с работой стало попроще и когда заставил работать будущее автономное. . .	сукцессия 25. Хронология ошибок anaschu 12.07.2026 # От 50-тонного гриба до устойчивого леса: хроника ошибок при построении модели вековой сукцессии микоризы ## О чём эта статья В процессе построения ОДУ-модели (система дифференциальных. . .

@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 3582 / 1712 / 236 Регистрация: 26.02.2009 Сообщений: 8,638 Записей в блоге: 6
	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