Построение фрактала Мандельброта

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

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

Добавил поток для построения, но рано обрадовался: прога как зависала, так виснет и сейчас.
И как использовать все ядра процессора на все 100? Кто-нибудь знает?

@FIL · 02.04.2013, 17:53

Smooth Coloring - победил:

@SeryZone · 02.04.2013, 21:30 **[ТС]**

gorfil, ух ты! Как?

@FIL · 02.04.2013, 22:40

Собственно, теорию взял тут (см. рис.):
http://en.wikipedia.org/wiki/M... te_note-23

Оказалось, что Z - это Zn, при котором выполняется условие (Sqr(z.Re) + Sqr(z.Im) > 4).
Отсюда вытекает основная проблема - ассемблерная вставка для вычисления n, собственно, возвращает только n (без Z), поэтому пришлось взять ее дельфийский вариант для экспериментов.
Получился следующий код:

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
      c.Im := iY * kZ + minIm;
      c.Re := iX * kZ + minRe;
      z.Re := 0;
      z.Im := 0;
 
      for n := 0 to Iterations do
      begin
        Re := Sqr(z.Re) - Sqr(z.Im);
        z.Im := z.Re * z.Im + z.Im * z.Re + c.Im;
        z.Re := Re + c.Re;
        if Sqr(z.Re) + Sqr(z.Im) > 4 then Break;
      end;
 
      if n < Iterations then
      begin
        if Smooth then // this should smoothen the iterations
          nSmooth := n - log2(ln(Sqrt(z.Re * z.Re + z.Im * z.Im)) / ln(2))
        else
          nSmooth := n;
        ScanLinePoint[iX] := Pallete[Round(nSmooth * Length(Pallete) / Iterations)];
      end else
        Integer(ScanLinePoint[iX]) := 0;

P.S. Может UI поможет нам ассемблерную вставку доработать, чтобы она еще и |Zn|>4 возвращала

?

@SeryZone · 03.04.2013, 08:08 **[ТС]**

gorfil, Да, надо бы ему написать

И я предлагаю два варианта: при сглаживании палитры пусть возвращает Z, иначе - это просто не нужно.

@FIL · 04.04.2013, 12:28

FadeDemon, прикрутил твою процедуру fractal_mandelbrot_fast к своему коду для сравнения.
Получилось (при прочих равных условиях), что с твоим кодом выигрыш ~18%:

@FadeDemon · 04.04.2013, 22:53

Блин, столько писал, писал, а форум взял все и не отправил

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

Получилось (при прочих равных условиях), что с твоим кодом выигрыш ~18%:

Можно посмотреть на твою реализацию, и на то как имено прикрутил?

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

Не думаю вообще что можно вытянуть больше чем эти 18%, я не представляю от чего ещё можно избавиться ради лишних наносекундок, но впринципе чего не получилось у меня, может и получиться у тебя.

Тут ещё есть некоторые идеи насчет плавных цветов, я если честно только тут про эту вещь и узнал

Во первых, константу, с которой мы сравниваем модуль, надо вместо 2 побольше брать, к примеру 16 или 32 (об этом кстати в википедии и пишут), плюс, от корня можно точно также избавиться как и раньше, по правилам математики - корень это степень 1/2, а степень аргумента можно вынести за логарифм, плюс ещё некоторые преобразования, и формула вырождается примерно вот в это (только что проверял, сей код - темпоральный вариант):

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
  for ix:=1 to w do begin
    complex_sqr(@z);
 
    complex_add(@z, c);
 
    if(complex_abs2(@z) > 1024)then begin
      w := ix;
      break;
    end;
 
  end;
 
  Result := w - ln(ln(complex_abs2(@z))) / ln(2) - 10;

Это для варианта, когда мы модуль числа сравниваем с 32, собственно константа 10 тут возникает как 2 * log(32) / log(2). Для сравнения:

Тут функция цвета (sin(n * 7) + 1) * 127, у тебя на картинке выше думаю полоски (хотя они конечно еле заметны) на сглаженном варианте уйдут полностью, думаю.

Остаеётся только проблема оптимизации подсчета логарифма, и тут я пока вижу 3 варианта (мне тоже захотелось себе такую гладкую штучку, очень захотелось):

1. Ряд тейлора и его модификации - вариант мертворожденный, ибо сам ряд расходиться, а модификации очень ресурсоёмки...
2. Метод ньютона для уравнения e^x - аргумент = 0, вариант как я обнаружил, в целом приемлимый, ряд тейлора экспоненты (экспоненту тоже только через FPU считать) нормально сходиться на ~8 итерациях, а итераций метода ньютона хватает 2 - но правда это для диапазона от 8 до 16, а если брать ту самую 1024... в общем печаль, обида, и досада, хотя если опдумать, то n > 1024 можно делить на собстно 1024, и получим интервал от ~1 до ~2, а деление можно легко компенсировать, прибавив к результату log(1024). Так, что не все тут потеряно.
3. Ну и мне кажется самый вероятный метод, это из SSE переключаться в FPU, и считать логарифм его вшитыми средствами

но проблема в том, что нету команд, напрямую переписывающих значения из XMM регистров в FPU стек, только в MMX регисты, а данные в FPU стеке храняться в Extended формате (80 бит), а MMX регистры всего лишь 64 бит имеют размер, плюс командой EMMS весь FPU стек помечается свободным, а без нее он, грубо выражаясь, не работает... В общем тоже не так тут все просто...
4. В инете вроде есть готовые реализации логарифма на SSE. Если ничего не получиться, можно обратиться и к ним...

Весь день над этим ломаю голову, в общем. Ещё я тоже начал размышлять над многопоточным и прогрессивным рисованием, у меня рисование по отдельным пикселям уже есть, осталось все это склеить, и получиться даже православнее чем Fractal eXtreme, я думаю

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

@FIL · 05.04.2013, 12:07

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

Можно посмотреть на твою реализацию, и на то как имено прикрутил?

Выкладываю как есть (там потоки вырезаны и много лишнего осталось).
Сравнивал fractal_mandelbrot_fast с ассемблерной вставкой на FPU, обернутой в 2 дельфийских цикла и с вариантом прогрессивного построения с одним потоком - результаты одинаковые - 18-20%. Думаю если циклы сделать на ASM, то процент немного сократится, но разница все равно существенная будет.

По поводу Smooth Coloring...
Последовал твоему совету и увеличил константу - результат отличный.
Вообще в Инете видел много разных вариаций с логарифмами, даже вариант после нахождения |Zn|>2 еще пару циклов прокрутить (видимо, это аналог увеличению константы).
Сам пока остановился примерно на таком варианте:

Delphi
1
2
zz := Sqrt(z.Re * z.Re + z.Im * z.Im);
nSmooth := n - log2(log2(zz) / log2(Sqrt(Range))); // где Range = 64

Результат замечательный и загадочные круги максимально далеко от центра

.

На счет оптимизации подсчета логарифма ничего умного не скажу, кроме как лог. в знаменателе посчитать заранее

. Вчера три часа потратил только на то, чтобы разобраться как работает вставка на FPU, чтобы Z оттуда добыть. Поэтому для меня будет хорошо, если я эти логарифмы хотя бы на FPU реализую.

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

у меня рисование по отдельным пикселям уже есть

Вот это будет интересно посмотреть.

Еще один глобальный вопрос, который хотелось бы обсудить:
Переменные типа Double (да и Extended тоже) позволяют увеличивать фрактал лишь до определенных пределов.
А как продвинуться дальше? Есть какие идеи?

@FIL · 05.04.2013, 18:07

Попытки вытащить |Zn| из ASM успехом не увенчались, чувствую нужна помощь зала

.
Зато получилась достаточно забавная разновидность фрактала:

@FadeDemon · 05.04.2013, 19:30

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

А как продвинуться дальше? Есть какие идеи?

Ну... Здесь путь, в общем то в целом один (как мне пока кажется), это изобретать свою нецелочисленную арифметику, с гейшами и перферансом, и на ней уже по своему определять операции сложения, умножения, деления... Правда тут уже несколько способов.

Можно поступить, как у нас в институте поступает препод вычмата, готовящий олимпиадников по програмированию: написать программу, которая реализует все базовые операции на больших (строго говоря любого размера) числах, заданных в строчку, то-есть например такая программа без труда сложит число запись которого в длинну несколько мегабайт, или по другому у него миллион знаков.. Более того, я думаю в инете подобных вещей навалом. Можно немного оптимизировать, записывая число не в строчку, а как положено, битами, но толко складывать/умножать по своим правилам, кстати, как минимум сложение таких длинных оптимизированных чисел можно сказать уже реализовано в архитектуре x86, там существуют такие команды как ADD и ADC, первая - обыкновенное сложение, а вот вторая - прибавление флага переноса (флага CF), который как раз выставляется командой ADD при возникновении переполнения. Умные компиляторы (делфи 7 кстати тоже) сложение/вычитание 64-битных чисел на 32-битных программах делают именно так, кстати (в x64 там с такими числами проблем по понятным причинам нету

). Понятно становиться, что с помощью большого кол-ва последовательно стоящих ADD/ADC можно складывать числа абсолютно произвольной длинны - все будет работать корректно. С умножением таких вещей сложнее, но делфи к примеру тоже умножение 64битных чисел как то делает, с помощью 32 битных команд (только что посмотрел, FPU там не используется), если разобраться как, то можно расширить и на произвольные числа, но это касается только целых числел, а нам нужны не целые....

И тут у меня недавно возникла другая идея - изобрести свой тип float заново. Если коротко, то я думаю все знаем мы, что любой тип с плавающе точкой состоит из 3 вещей: мантиссы, порядка, и знакового бита. А ограничение, из-за которого фрактал нельзя увелчить больше положенного, возникает, как я понял, из-за того, что переполняется именно порядок, то-есть гдето в рачетах возникают числа по порядку меньшие чем 10^-300 степени, которые тип Double перестает переваривать (Extended перестает переваривать числа порядка 10^-4000), а все потому, что в этих типах для хранения порядка выделено совсем немного битов (для float - 8 бит, для дубла чуть побольше). С мантиссой, то-есть битами точности как мне кажется всё ок, их хвататет с избытком, и при увеличении порядка, мантисса то собственно не деградирует, допустим если умножить 0.256565 и 0.35678758, получиться примерно 0.0915392054627, в мантиссу запижется чтото типо 915392054627, и если умножить 0.000...(милион нулей)...00000256565 и 0.000...(милион нулей)...0000035678758, в мантиссу запишется то-же самое, но вот порядок.. в общем, вся загвоздка у нас в этом самом порядке, нам надо его значительно расширить.

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

Как только у нас появляется такая структура, мы можем ее считать в общем-то новым полноценным типом с плавающей точкой, нам остается только научить компьютер такие структуры вычитать, складывать, умножать, делить... Если мы научимся делать эти 4 базовые вещи, мы тут же научимся и возводить в целочисленную степень. А как только научимся возводить в степень, то немного покурив вычислительную математику, узнав что такое ряд тейлора и метод Ньютона, мы сможем на наших числах с легкостью определить такие вещи как синус, косинус, экспонента, логарифм и прочие взрочлые штуки (тут тольо оговориться надо, это будет работать ЛЮТО, БЕШЕНО, ЗВЕРСКИ НЕВЫНОСИМО МЕДЛЕННО, как впрочем и все наши операции с такими числами).

Теперь, попробуем научиться такие структуры складывать, тут нам помогут некоторые школьные знания, и интуиция... Допусим, сначала попробуем разобраться, как же всетаки будут интуитивно представляться у нас числа с плавающей точкой. Рассмотрим число 1.0. В целом это то же самое что 1. Для его записи, в мантиссу мы поместим 1, а в порядок 0. Такая конфигурация значений у нас и будет означать 1.0, в нашем новом формате чисел. Чтобы закрепить понимание, рассмотрим числа 100.0 и 0.01. Весь прикол в том, что у всех трех рассмотренных чисел, мантисса будет равна 1, но у 100 - порядок будет 2, а у 0.01 - порядок -2. Это собствено нам вспоминается школьная арифметика: порядок есть ни что иное, как положение плавающей точки в числе. Настоящие числа, double например, устроены так же, или почти также, точно на этот счет не знаю. Но, такого как я описал представления, нам будет достаточно.

Попробуем теперь числа 10 и 1 сложить. У обоих мантисса 1, у первого порядок 1, у второго 0. Что мы делаем: берем едицицу, двигаем в ней запятую на указаный порядок (то-есть умножаем ее на 10), и, собственно, обыкновенным сложением складываем. Записываем результат в мантиссу, а вот порядко в этом случае устанавливаем в 0, так как в мантиссе уже лежит число 11, и оно само автоматически имеет нужный порядок. Теперь другой пример: складываем 0.0000..(муллион нулей)..000001 и 1. Что тут можно сделать? Мы поступим по аналогии. В мантиссе у обоих чисел снова единицы. Но у первого порядок -1000000, у второго 0. Мы берем, и единицу миллион раз делим на 10. И то что получилось, прибавляем к единице, тоесть первому слагаемому, и выставляем его порядок. Нам здесь должно быть профиг, что от единицы, деленной миллион раз на 10 ничего не останется, на самом деле, точно такая же потеря точности происходит и в настоящих типах, и это совершенно нормально - в этом то и кроется секрет разделения числа на мантиссу (биты точности) и порядок (количество нулей, на которые эта точность уезжает). Попробуем теперь сложить 0.001 и 0.00001. У первого порядок -3, у второго -5. Что же мы делаем тут?... В мантиссе опять единицы. Но разница порядоков - 2. Получается, мы берем единицу, и делим ее 2 раза на 10, чтобы компенсировать разницу порядком, но мы помним тут еще, что мантисса и порядок у нас целые, и все наши манипуляции мы производим с целыми типами, и поэтому, чтобы корректно "поделить 1 на 100", мы на самом деле эту единицу трогать не будем, а умножим на 100 предыдущую, и запомним, что мы на 2 порядок увеличили. После сложения у нас получается 101, и это мы пишем в мантиссу, плюс мы помним, что задрали числу порядок на 2, а порядок бОльшего слагаемого, у нас был -3. Получается, мы возвращаем порядок на место, тоесть прибавляем к порядку -2, и добавляем туда -3, получается -5. И если руками сдвинуть запятую в 101 на 5 знаков назад, то получим, что мы не ошиблись, и у нас действительно получается 0.00101, что и должно было получиться.

Брр... Аж глаза на лоб полезли от таких размышлений. В общем, вычитание можно определить аналогичным способом, для хранения знака можно использовать признак отрицательной мантиссы, например. Вот как именно определять умножение, так с ходу сказать не могу, могу только сказать что мантиссы там напрямую перемножаются, а порядки складываются, но дело все в том, что если мы попытаемся умножить 0xFFFFFFFFFFFFFFFF само на себя, то оно явно вылезет за пределы Int64, и нас постигнет фейл, поскольку это не потеря точноти будет - потеряются старшие биты. что делать с этим, я чесговоря не придумал.....

Это все в общем так, размышления, вдруг у когонибуть появиться желание такое реализовать

Такую вот арифметику вполне множно использовать для рассчета фрактала, и даже логарифмы на такой арифметике прикрутить, но опять повторюсь, работать это будет просто ЗВЕРСКИ медленно, учитывая что в наших текущих кодах, сложение-умножение выполняются по операции на одну команду процессора. С рядами тейлора - там вообще жесть, чтобы посчитать допустим экспоненту с нужной нам точности, нам фактически надо посчитать многочлен 9-10 степени, а то и побольше, что потребует собой 10² / 2 операций умножения (это если без схемы горнера, правда, по этой схеме многочлен можно вроде быстрее посчитать, за 10 умножений, не помню чтоно, это типо если икс за скобки максимально вынести), и 9 операций сложения, я уже не говорю что там на констаныт делить надо. И это только чтобы посчтить ОДНУ ЭКСПОНЕНТУ! А чтобы посчитать логарифм методом ньютона, придется сделать штук 5-6-20-100 итераций, тоесть посчитать столько же экспонент.... Это все ради только одного, единственного вызова логарифма, с нормальной точностью, которую нам вот гарантирует FPU. Причем заметим, всю эу порногорафию, FPU проделывает одной командой. ОДНОЙ. а сколько тут понадобиться нам.... ну думаю можно осознать всю серьезность проблемы. Так что задумка своей целочисленной арифметики - весьма спорная вещь. Но с другой стороны, очень интересная теоретически. Лично мне тоже интересно позумить наш фрактал, он ведь не полностью самоподобен то....

P.S прошу прощения за многабукав, привык целые статьи писать

Кстати в конце добавить можно, что в инете полюбому можно найти готовые реализации библиотек для работы с большими числами, ну вот просто полюбому... и их уже использовать, чем наверно я и советую заняться, если чильно захочется такой высокозумленый фрактал. Но, написать своё - это всегда намного интересне... (по крайней мере мне

)

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

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

Попытки вытащить |Zn| из ASM успехом не увенчались, чувствую нужна помощь зала

Попробуй так: заведи в параметрах процедуры асмовставки переменную var xxx: Double, и в том месте, где в асмокоде должен храниться модуль (или модуль без корня) написать чтото типо

Assembler
1
FST qword ptr xxx

или просто

Assembler
1
FST xxx

или еще можно по другому

Assembler
1
2
FLD ST(0)
FSTP xxx

ну и вызывая вставку, подставлять собственно на это место переменную. в теории, в нее должен сохраняться модуль, или что там было в коде на вершине FPU стека. если не получиться, потом скажу точнее, подозреваю что не все тут так просто, и что делфи сволочь хитрая, и параметры может ради оптимизации передавать без угрызений совести прямо через стек FPU..
ах да, в отладчике делфи есть режим просмотра FPU - тыкаешь в асмокоде на какойнибуть строчке F4, потом Alt+Ctrl+C, и там Ctrl+U - вылезет FPU стек, я думаю наблюдая за ним можно понять, где там может храниться модуль... если вдруг он храниться не на вершине, то тогда:

Assembler
1
2
FLD ST(i)
FSTP xxx

где i - индекс строки сверху в отладчике, начиная с нуля. как-то так

@FIL · 06.04.2013, 17:10

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

Попробуй так:

FadeDemon, я примерно эти варианты и пробовал. Даже тему отдельную запостил:
Из ассемблерной вставки в локальную переменную
Просмотрел через отладчик - оказывается fcomip st, st(5) выбрасывает нужное мне значение. Спасибо за совет.
Надеюсь проблема только в этом.

По зуму...
Про "Длинную" арифметику я читал, коды различных операций есть. Но применительно к данной задаче - это, как мне кажется, не выход. Слишком все медленно будет работать. Если других вариантов нет, то это очень печально.

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

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

Попытки вытащить |Zn| из ASM успехом не увенчались

Сделал. Вот код:

Кликните здесь для просмотра всего текста

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
function GetMandelbrot_Asm(cRe, cIm: Extended): Integer;
asm
  push ebx      // Ñîõðàíåíèå ebx
 
  finit         // Èíèöèàëèçàöèÿ ïðîöåññîðà
  fld Range     // Çàãðóçêà âåùåñòâåííîãî çíà÷åíèÿ
  fld cRe       //
  fld cIm       //
  fldz          // z.Im = 0 //Çàãðóçêà +0.0
  fldz          // z.Re = 0
 
  xor ebx, ebx  // Îáíóëåíèå ðåãèñòðà ebx
  mov ecx, Iterations
 
@@lloop:
                // Re := Sqr(z.Re) - Sqr(z.Im);
  fld st        // Êîïèðóåì z.Re
  fmul st, st   // z.Re * z.Re
  fld st(2)     // Êîïèðóåì z.Im
  fmul st, st   // z.Im * z.Im
  fsubp         // Sqr(z.Re) - Sqr(z.Im) è óäàëÿåì z.Im
 
                 // z.Im := z.Re * z.Im + z.Im * z.Re + c.Im;
  fld st(1)      // Êîïèðóåì z.Re
  fmul st, st(3) // z.Re * z.Im
  fadd st, st    // z.Re * z.Im + z.Re * z.Im
  fadd st, st(4) // z.Re * z.Im + z.Im * z.Re + c.Im
 
                 // z.Re := Re + c.Re;
  fstp st(3)     // Ïåðåíîñèì Re //Ñîõðàíåíèå âåùåñòâåííîãî çíà÷åíèÿ è èçâëå÷åíèå èç ñòåêà
  fadd st, st(4) // Re + c.Re;
 
                 // Sqr(z.Re) + Sqr(z.Im)
  fstp st(1)     // Ïåðåíîñèì z.Re (çàòèðàåì çíà÷åíèå â st(1))
  fld st         // Êîïèðóåì z.Re
  fmul st, st    // z.Re * z.Re
  fld st(2)      // Êîïèðóåì z.Im
  fmul st, st    // z.Im * z.Im
  faddp          // Sqr(z.Re) + Sqr(z.Im) è óäàëÿåì Sqr(z.Im)
 
//  fcomip st, st(5)
  fcomi st, st(5)
  jc @@next      // st(0) < Range => Continue
  jz @@next      // st(0) = Range => Continue
 
  // here : st(0) > Range => save Result and Exit
  mov Result, ebx
  
  fsqrt          // Êâàäðàòíûé êîðåíü
  fst ZZ         // Çàïèñü âåùåñòâåííîãî çíà÷åíèÿ
 
  jmp @@exit
 
@@next:
  fstp st
  inc ebx
  loop @@lloop
  mov ebx, Iterations
  mov Result, ebx
 
@@exit:
  pop ebx
end;

Надо добавить глобальную переменную ZZ: Double; и это будет |Zn|.

Добавлено через 19 часов 2 минуты
Переделал ASM вставку и добавил туда режим Smooth Coloring:

Кликните здесь для просмотра всего текста

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
asm
  push ebx             // Ñîõðàíåíèå ebx
 
  finit                // Èíèöèàëèçàöèÿ ïðîöåññîðà
  fld EscapeRadius     // Çàãðóçêà âåùåñòâåííîãî çíà÷åíèÿ
  fld z.Re             //
  fld z.Im             //
  fldz                 // z.Im = 0 //Çàãðóçêà +0.0
  fldz                 // z.Re = 0
 
  xor ebx, ebx         // Îáíóëåíèå ðåãèñòðà ebx
  mov ecx, Iterations
 
@@lloop:
                       // Re := Sqr(z.Re) - Sqr(z.Im);
  fld st               // Êîïèðóåì z.Re
  fmul st, st          // z.Re * z.Re
  fld st(2)            // Êîïèðóåì z.Im
  fmul st, st          // z.Im * z.Im
  fsubp                // Sqr(z.Re) - Sqr(z.Im) è óäàëÿåì z.Im
 
                       // z.Im := z.Re * z.Im + z.Im * z.Re + c.Im;
  fld st(1)            // Êîïèðóåì z.Re
  fmul st, st(3)       // z.Re * z.Im
  fadd st, st          // z.Re * z.Im + z.Re * z.Im
  fadd st, st(4)       // z.Re * z.Im + z.Im * z.Re + c.Im
 
                       // z.Re := Re + c.Re;
  fstp st(3)           // Ïåðåíîñèì Re //Ñîõðàíåíèå âåùåñòâåííîãî çíà÷åíèÿ è èçâëå÷åíèå èç ñòåêà
  fadd st, st(4)       // Re + c.Re;
 
                       // Sqr(z.Re) + Sqr(z.Im)
  fstp st(1)           // Ïåðåíîñèì z.Re (çàòèðàåì çíà÷åíèå â st(1))
  fld st               // Êîïèðóåì z.Re
  fmul st, st          // z.Re * z.Re
  fld st(2)            // Êîïèðóåì z.Im
  fmul st, st          // z.Im * z.Im
  faddp                // Sqr(z.Re) + Sqr(z.Im) è óäàëÿåì Sqr(z.Im)
 
  fcomi st, st(5)
  ja @@Result          // st(0) > Range
                       // st(0) <= Range => Continue
  fstp st
  inc ebx
  loop @@lloop
 
@@Result:
  push ebx
  fild dword ptr [esp] // Êîïèðóåì ïîëó÷åííîå çíà÷åíèå n â st
  pop ebx
 
  cmp Smooth, False    // Ðåæèì Smooth Coloring
  jz @@Exit            // åñëè âûêëþ÷åí
                       // nSmooth := n - log2(log2(ZZ) / logEscapeRadius)
  fstp st(2)
  fsqrt                // Ïîëó÷àåì |Z| // Êâàäðàòíûé êîðåíü
  fld1                 //Çàãðóçêà +1.0
  fld1
  fxch st(2)           // Îáìåí ñîäåðæèìîãî ðåãèñòðîâ
  fyl2x                // log2(|Z|) // st(1) * log2 st(0)
  fld logEscapeRadius
  fdiv                 // log2(|Z|) / logEscapeRadius // Äåëåíèå âåùåñòâåííûõ ÷èñåë
  fyl2x                // log2(log2(|Z|) / logEscapeRadius)
  fsub                 // Âû÷èòàíèå âåùåñòâåííûõ çíà÷åíèé (ñ îáðàùåíèåì)
//  frndint              // Îêðóãëåíèå äî öåëîãî
 
@@Exit:
  fstp n               // Ñîõðàíÿåì ïîëó÷åííîå çíà÷åíèå n
  pop ebx
end;

Smooth: boolean (глобальная) - включает режим.
n: Double (локальная) - собственно n.

@SeryZone · 07.04.2013, 19:31 **[ТС]**

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

@FIL · 07.04.2013, 21:33

Есть такая мысля:
Мы в своих расчетах используем числа, отличающиеся друг от друга лишь последними цифрами.
Например, берем диапазон Re:

Code
1
2
MinRe = 0,12345678901
MaxRe = 0,12345678902

Следовательно, общую часть мы можем представить отдельным числом (как бы "освободить" мантиссу"):

Code
1
2
MinRe = 0,12345678901 = 0,1234567891 + 0,0000000001 = 1,23456789e-1 + 1,e-10
MaxRe = 0,12345678902 = 0,1234567892 + 0,0000000002 = 1,23456789e-1 + 2,e-10

Что (теоретически) позволит нам "углубиться" еще на 9 знаков.
Из мат. операций мы используем только [+] [-][*], поэтому потери точности быть не должно.
Отделять числа от основного надо по 9 цифр, т.к. мантисса (Extended) вмещает до 19 цифр, а макс. значение при перемножении числа из 9 цифр = 18.

P.S. Устройство вещественных типов: http://www.delphikingdom.com/a... alogID=374

@FIL · 09.04.2013, 21:17

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

Есть такая мысля:

Походу, бред это все. Придется через "длинную" арифметику пробовать.

@FadeDemon · 10.04.2013, 07:00

Проблема то не в мантиссе, проблема здесь исключительно в порядке, вернее в его переполнении, когда возникают числа меньше приблизительно 10^-300, то порядок не может сместиться в меньшую сторону, и чтобы это компенсировать, страдать начинает мантисса, которую просто брутально сдвигают, компенсируя невозможность вместо этого уменьшать порядок, отсюда то и лезут потери точности... Поэтому тут просто надо научиться хранить "неограниченно длинный" порядок, как наиболее перспективный вариант, мне видиться самодельное "плавающее" число, представленное как структура из двух Int64, которые по отдельности будут отвечать за мантиссу и порядок, ну, как я выше писал... Разве что под порядок можно и по больше чисел выделить, например штуки 4 Int64, на дольше хватит

хотя, числа вида 10^2⁶³ уже, в общем-то, внушают уважение..

@FIL · 10.04.2013, 10:07

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

@SeryZone · 10.04.2013, 16:02 **[ТС]**

gorfil, это так и есть. И я считаю, что его тоже нужно делать на SSE...

@Игорь[Igor] · 11.04.2013, 00:21

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
procedure DrawMandelbrot(ACanvas: TCanvas; X, Y, au, bu: double; X2, Y2: integer);
var
  c1, c2, z1, z2, tmp: double;
  i, j, Count: integer;
  tmp2: double;
begin
  c2 := bu;
  for i := 10 to X2 do
  begin
    c1 := au;
    for j := 0 to Y2 do
    begin
      z1 := 0;
      z2 := 0;
      Count := 0;
     {count is deep of iteration of the mandelbrot set
      if |z| >=2 then z is not a member of a mandelset}
      while (((z1 * z1 + z2 * z2 < 4) and (Count <= 90))) do
      begin
        tmp := z1;
        z1 := z1 * z1 - z2 * z2 + c1;
        z2 := 2 * tmp * z2 + c2;
        inc(Count);
      end;
      //the color-palette depends on TColor(n*count mod t)
{$IFDEF LINUX}
      ACanvas.Pen.Color := (16 * Count mod 255);
      ACanvas.DrawPoint(j, i);
{$ELSE}
      ACanvas.Pixels[j, i] := (16 * Count mod 255);
{$ENDIF}
      c1 := c1 + X;
    end;
    c2 := c2 + Y;
  end;
end;
 
 
 
procedure TForm1.Button1Click(Sender: TObject);
var
  R: TRect;
  au, ao: integer;
  dX, dY, bo, bu: double;
begin
  // Initialize Mandelbrot
  R.Left := 0;
  R.Right := 500;
  R.Top := 0;
  R.Bottom := 505;
  ao := 1;
  au := -2;
  bo := 1.5;
  bu := -1.5;
  //direct scaling cause of speed
  dX := (ao - au) / (R.Right - R.Left);
  dY := (bo - bu) / (R.Bottom - R.Top);
  DrawMandelbrot(Self.Canvas, dX, dY, au, bu, R.Right, R.Bottom);
end;

@SeryZone · 12.04.2013, 15:37 **[ТС]**

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

@Игорь[Igor] · 12.04.2013, 16:06

Поможет, еще быстрее будет через SetPixel или ScanLine Bitmap'а

@SeryZone · 15.04.2013, 14:05 **[ТС]**

Методы ускорения и оптимизации:
1) Проверка на симметрию - помогите, я там что-то написал - не получается

2) Проверка на периодичность(Periodicity Checking) - провера чёрных мест. Там какие-то уравнения, о них потом.

С симметрией помогите, не получается.

Новые блоги и статьи Все статьи Все блоги /
Кредитный калькулятор Maks 05.08.2026 Решение задачи по прикладной информатике средствами 1С. Задача: Напишите приложение-калькулятор, которое помогает рассчитывать параметры кредита для аннуитетного и дифференцированного видов. . .	У нас сейчас поговорку "Опять 25" нужно переделать на "Опять +35". kumehtar 04.08.2026 С ностальгией вспоминаю времена моего детства, когда у нас и правда +25 - была максимальная температура летом. Раньше +25 °C реально казались вершиной жары, когда можно было весь день пропадать на. . .	Как ИИ начал спорить и врать (возможно почуяв опасность для себя от индустрии - уход от электроники). Hrethgir 04.08.2026 Недельный диалог, на фоне событий с НПЗ. Да, из спирта можно получать бензин, и это не сложно. Но потом в схеме я решил избавиться от насоса, при этом полностью сделав контроль подачи спирта в. . .	Термопринтер QR701 Argus19 03.08.2026 Термопринтер QR701 Купил два термопринтера QR701. На сэлф-тесте написано: Language: PC936 (GB18030). Что означает, что принтеры могут печатать только латиницу и китайские иероглифы. Так же. . .
Создание формы заимствованного документа Maks 03.08.2026 Задача: Необходимо создать собственную форму заимствованного документа. На форме должен быть реквизит "Покупатель", а также табличная часть со следующими реквизитами: - Расчетный счет покупателя. . .	Задача предоставления скидок покупателям Maks 03.08.2026 Задача: В документе "Продажи" необходимо реализовать функционал предоставления скидок покупателям. Скидка должна автоматически рассчитываться и подставляться в соответствующее поле при выборе. . .	Почему SEO не начинается с ключевых слов: что проверить до написания текстов Neotwalker 01.08.2026 Когда владельцу сайта предлагают заняться SEO, первым шагом часто становится сбор запросов и написание текстов. Логика кажется понятной: 1. Находим ключевые слова. 2. Добавляем их на. . .	Знание — сила: Доктрина интенциональности знаний, углубление в формулу Hrethgir 01.08.2026 https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11957&stc=1&d=1785567302 Знаменитый афоризм Фрэнсиса Бэкона «Знание — сила» (Scientia potentia est) в массовой культуре принято понимать. . .

@FIL Модератор 3492 / 2614 / 742 Регистрация: 19.09.2012 Сообщений: 7,977
	02.04.2013, 17:53
	Smooth Coloring - победил: Миниатюры 0

@SeryZone 56 / 28 / 18 Регистрация: 09.03.2012 Сообщений: 726 Записей в блоге: 1
	02.04.2013, 21:30 [ТС]
	gorfil, ух ты! Как? 0

@SeryZone 56 / 28 / 18 Регистрация: 09.03.2012 Сообщений: 726 Записей в блоге: 1
	03.04.2013, 08:08 [ТС]
	gorfil, Да, надо бы ему написать И я предлагаю два варианта: при сглаживании палитры пусть возвращает Z, иначе - это просто не нужно. 0

@FIL Модератор 3492 / 2614 / 742 Регистрация: 19.09.2012 Сообщений: 7,977
	04.04.2013, 12:28
	FadeDemon, прикрутил твою процедуру fractal_mandelbrot_fast к своему коду для сравнения. Получилось (при прочих равных условиях), что с твоим кодом выигрыш ~18%: Миниатюры 0

@SeryZone 56 / 28 / 18 Регистрация: 09.03.2012 Сообщений: 726 Записей в блоге: 1
	07.04.2013, 19:31 [ТС]
	По поводу целочисленной арифметики: может реализовать битовыми операциями на ассемблере??? Только ж блин, я код не понимаю, куда мне типы программировать... 0

@FadeDemon 60 / 0 / 1 Регистрация: 31.03.2013 Сообщений: 5
	10.04.2013, 07:00
	Проблема то не в мантиссе, проблема здесь исключительно в порядке, вернее в его переполнении, когда возникают числа меньше приблизительно 10^-300, то порядок не может сместиться в меньшую сторону, и чтобы это компенсировать, страдать начинает мантисса, которую просто брутально сдвигают, компенсируя невозможность вместо этого уменьшать порядок, отсюда то и лезут потери точности... Поэтому тут просто надо научиться хранить "неограниченно длинный" порядок, как наиболее перспективный вариант, мне видиться самодельное "плавающее" число, представленное как структура из двух Int64, которые по отдельности будут отвечать за мантиссу и порядок, ну, как я выше писал... Разве что под порядок можно и по больше чисел выделить, например штуки 4 Int64, на дольше хватит хотя, числа вида 10^2⁶³ уже, в общем-то, внушают уважение.. 0

@FIL Модератор 3492 / 2614 / 742 Регистрация: 19.09.2012 Сообщений: 7,977
	10.04.2013, 10:07
	По моим прикидкам, проблема исключительно в мантиссе - не хватает цифр точности. При этом степень вообще практически не задействуется (см. рис.). Когда точность вычислений опускается ниже 19-й цифры (для Extended) все застопоривается. Миниатюры 0

@SeryZone 56 / 28 / 18 Регистрация: 09.03.2012 Сообщений: 726 Записей в блоге: 1
	10.04.2013, 16:02 [ТС]
	gorfil, это так и есть. И я считаю, что его тоже нужно делать на SSE... 0

@SeryZone 56 / 28 / 18 Регистрация: 09.03.2012 Сообщений: 726 Записей в блоге: 1
	12.04.2013, 15:37 [ТС]
	Не думаю, что канва поможет быстро построить фрактал... 0

@Игорь[Igor] 726 / 478 / 130 Регистрация: 24.12.2008 Сообщений: 3,924
	12.04.2013, 16:06
	Поможет, еще быстрее будет через SetPixel или ScanLine Bitmap'а 0