分形/曼德爾布羅特數值計算
< 分形
mandelbrot-numerics 是一個用於 Claude Heiland-Allen 建立的與曼德爾布羅特集相關的先進數值計算的庫。
這是一個關於它的非官方維基，包含了大多數 maximus/book 的 "lib" 和 "bin"（但不包括渲染）。
檢視
- 字典
- 符號
- 庫的分支
- 來自 Claude Heiland-Allen 的庫
- exrtact 用於操作 EXR 影像的各種工具 與 OpenExr 進行比較。
- kf-extras 用於操作 Kalles Fraktaler 2 輸出的程式。
- mandelbrot-symbolics - 與曼德爾布羅特集相關的符號演算法。
- mandelbrot-graphics - 基於 CPU 的曼德爾布羅特集視覺化。
- mandelbrot-text - 與曼德爾布羅特集相關的解析和漂亮列印。
- bin 用於程式
- lib 用於函式，有兩個版本
- d 用於雙精度
- r 用於任意精度
- 包含
這裡的精度指的是 數值計算中浮點數的精度。
它是一個正整數，表示二進位制數的位數（
- 雙精度：m_d_*()
- 接受指向 double _Complex 的指標的函式使用它們作為輸出。
- 任意精度：m_r_*()
- 從最小精度開始（double = 53 位）。
- 計算新的精度 = 53 + f(引數)，其中引數是過程的一個特徵（例如，對於 m-interior，元件的最小大小）。換句話說，新的精度大於 53。
"你需要 "make install" mandelbrot-symbolics 庫，然後再嘗試使用 mandelbrot-numerics"。
- pkg-config
- gmp (libgmp)
- mpfr (libmpfr)
- mpc (libmpc)
- haskell - ghc
- sndfile 使用：sudo apt-get install libsndfile1-dev[1]
#include <complex.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdlib.h>
#include <gmp.h>
#include <mpfr.h>
#include <mpc.h>
Gcc 標誌
gcc -lmpc -lmpfr -lgmp -lm
從控制檯
sudo aptitude install libmpc-dev
從你想要克隆 git 倉庫的目錄
git clone https://code.mathr.co.uk/mandelbrot-numerics.git
然後從同一個目錄
make -C mandelbrot-numerics/c/lib prefix=${HOME}/opt install
make -C mandelbrot-numerics/c/bin prefix=${HOME}/opt install
然後執行：[2]
export LD_LIBRARY_PATH=${HOME}/opt/lib
檢查
echo $LD_LIBRARY_PATH
結果
/home/a/opt/lib
或者
export PATH=${HOME}/opt/bin:${PATH}
檢查
echo $PATH
要永久設定它，更改 .profile 檔案[3]
sudo gedit ~/.profile
m-describe 53 20 100 0 0 4
LD_LIBRARY_PATH=${HOME}/opt/lib
m-primary-separators
從在 mandelbrot-numerics 目錄中開啟的控制檯
git pull
如果你進行了一些本地更改，你可以撤銷它們
git checkout -f
然後
git pull
現在重新安裝
//mandelbrot-numerics.h
struct m_d_mat2 {
double _Complex a, b, c, d;
};
typedef struct m_d_mat2 m_d_mat2;
- 要包含庫，c 原始碼應該 *只* 包含 #include <mandelbrot-numerics.h>
- 使用 pkg-config 編譯和連結：檢視 mandelbrot-numerics/c/bin/Makefile 中的示例。
- 最快的入門方法是將你的檔案放到 mandelbrot-numerics/c/bin 中並執行 make。
- m-feature-database
- m-misiurewicz
用法
m-misiurewicz precision guess-re guess-im preperiod period maxsteps
描述
- 使用 臨界點的預週期
雙精度示例
m-misiurewicz double -2.1 0 2 1 100 -2.0000000000000000e+00 0.0000000000000000e+00 m-misiurewicz double -1.6 0 3 1 100 -1.5436890126920764e+00 0.0000000000000000e+00 m-misiurewicz double -1.5 0 5 2 100 -1.4303576324513072e+00 0.0000000000000000e+00 m-misiurewicz double -1.4 0 9 4 100 -1.4074051181647020e+00 0.0000000000000000e+00
使用位
./m-misiurewicz 200 -1.4 0 9 4 100 -1.4074051181647020225078282291990509777838059260945479350908287e+00 0.0000000000000000000000000000000000000000000000000000000000000e+00
m-feature-database
usage: ./m-feature-database preperiod period
m-feature-database 1 1
.1(0) 1/2 -2.00 0.00
m-describe
- 官方頁面
- 此過程提供引數平面中點 c 的描述。
無引數使用
m-describe
提供使用說明
usage: m-describe precision maxperiod maxiters re im ncpus
檢查
- 點 c = 0 = 0 +0*I= re+ im*I
- 精度 = 53 位
- 最大週期= 20
- 最大迭代次數 = 100
m-describe 53 20 100 0 0 4
輸出為
the input point was +0e+00 + +0e+00 i the point didn't escape after 100 iterations nearby hyperbolic components to the input point: - a period 1 cardioid with nucleus at +0e+00 + +0e+00 i the component has size 1.00000e+00 and is pointing west the atom domain has size 0.00000e+00 the atom domain coordinates of the input point are -nan + -nan i the atom domain coordinates in polar form are nan to the east the nucleus is 0.00000e+00 to the east of the input point the input point is interior to this component at radius 0.00000e+00 and angle 0.000000000000000000 (in turns) the multiplier is +0.00000e+00 + +0.00000e+00 i a point in the attractor is +0e+00 + +0e+00 i external angles of this component are: .(0) .(1)
另一個例子
m-describe 1000 1000 10000000 -1.749512 0.0 8 the input point was -1.74951199999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999e+00 + +0e+00 i the point didn't escape after 10000000 iterations nearby hyperbolic components to the input point: - a period 1 cardioid with nucleus at +0e+00 + +0e+00 i the component has size 1.00000e+00 and is pointing west the atom domain has size 0.00000e+00 the atom domain coordinates of the input point are -nan + -nan i the atom domain coordinates in polar form are nan to the east the nucleus is 1.74951e+00 to the east of the input point the input point is exterior to this component at radius 1.82808e+00 and angle 0.500000000000000000 (in turns) the multiplier is -1.82808e+00 + +0.00000e+00 i a point in the attractor is -9.14032e-01 + +0e+00 i external angles of this component are: .(0) .(1) - a period 2 circle with nucleus at -1e+00 + +0e+00 i the component has size 5.00000e-01 and is pointing west the atom domain has size 1.00000e+00 the atom domain coordinates of the input point are -0.74951 + -0 i the atom domain coordinates in polar form are 0.74951 to the west the nucleus is 7.49512e-01 to the east of the input point the input point is exterior to this component at radius 2.99805e+00 and angle 0.500000000000000000 (in turns) the multiplier is -2.99805e+00 + +0.00000e+00 i a point in the attractor is -1.49976e+00 + +0e+00 i - a period 3 cardioid with nucleus at -1.754878e+00 + +0e+00 i the component has size 1.90355e-02 and is pointing west the atom domain has size 2.34487e-01 the atom domain coordinates of the input point are -0.022921 + +0 i the atom domain coordinates in polar form are 0.022921 to the west the nucleus is 5.36567e-03 to the west of the input point the input point is exterior to this component at radius 1.24010e+00 and angle 0.000000000000000000 (in turns) the multiplier is +1.24010e+00 + +0.00000e+00 i a point in the attractor is -6.8613231e-02 + +0e+00 i external angles of this component are: .(011) .(100) - a period 237 cardioid with nucleus at -1.7495120000000000000000000000000116053893024686343718029506670414e+00 + +0e+00 i the component has size 5.10303e-60 and is pointing west the atom domain has size 1.34480e-32 the atom domain coordinates of the input point are -0.86149 + -0 i the atom domain coordinates in polar form are 0.86149 to the west the nucleus is 1.16054e-32 to the west of the input point the input point is exterior to this component at radius 5.87695e+33 and angle 0.500000000000000000 (in turns) the multiplier is -5.87695e+33 + +0.00000e+00 i a point in the attractor is +2.5892947214253619092848904520135803680047966975848204640792969525e-02 + +0e+00 i - a period 756 cardioid with nucleus at -1.74951200000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000437443253498601125336414441820009561885311706078736706424485680024027767522121183394207260374919267049e+00 + +0e+00 i the component has size 2.81459e-198 and is pointing west the atom domain has size 8.60400e-102 the atom domain coordinates of the input point are -0.51039 + +0 i the atom domain coordinates in polar form are 0.51039 to the west the nucleus is 4.37443e-102 to the west of the input point the input point is exterior to this component at radius 1.99523e+68 and angle 0.500000000000000000 (in turns) the multiplier is -1.99523e+68 + +0.00000e+00 i a point in the attractor is -1.32154302717027851153258720563184483110950266513429522776342823541752294814605666085090894143487158886716259121852162061940587472967643029955949630872182624723698360235112807868903937707485681031295828661e-02 + +0e+00 i
版本
- m_d.shape.c
- m_r_shape.c
描述[4]
http://ibiblio.org/e-notes/MSet/windows.htm <-- where I got the formula from
檔案
- m-size.c
- m_d_size.c
- 輸出：複數雙精度大小
- m_r_size.c
- 輸出 : mpc_t 大小
// mandelbrot-numerics -- 與 Mandelbrot 集相關的數值演算法 // 版權所有 (C) 2015-2018 Claude Heiland-Allen // 許可證 GPL3+ http://www.gnu.org/licenses/gpl.html
#include <mandelbrot-numerics.h>
extern double _Complex m_d_size(double _Complex nucleus, int period) {
double _Complex l = 1;
double _Complex b = 1;
double _Complex z = 0;
for (int i = 1; i < period; ++i) {
z = z * z + nucleus;
l = 2 * z * l;
b = b + 1 / l;
}
return 1 / (b * l * l);
}
無引數使用
m-size
提供使用說明
usage: ./m-size precision nucleus-re nucleus-im period
引數
結果是 2 個數字
- cabs(size) = size 和 carg(size) = 方向估計
- mpc_abs(r, size, MPFR_RNDN); mpc_arg(t, size, MPFR_RNDN);
It gives the size as a complex number, the magnitude is the size relative to the period 1 continent and the phase / argument is the angle of rotation.
程式碼：$ m-size 53 -8.6915874972342078e-01 2.5565708568620021e-01 48 1.1864737161932281e-08 -5.9691937849660148e-01
我的程式為了方便將其拆分為幅度和相位（以弧度為單位），所以微型布魯瓦是週期 1 大陸大小的 1.1865e-8 倍，所以將初始檢視的檢視大小乘以該數字以使微型布魯瓦顯示為大致相同的大小。
示例用法
m-size double 0 0 1
輸出
1.0000000000000000e+00 0.0000000000000000e+00
換句話說
the component has size 1.00000e+00 and is pointing west
週期為 48 的微型布魯瓦[8]
m-size 53 -8.6915874972342078e-01 2.5565708568620021e-01 48 1.1864737161932281e-08 -5.9691937849660148e-01
計算 複數二次多項式下的週期
"但是，盒週期方法檢測盒內 *核* 的週期 - 如果那裡沒有核，它將失敗。"
Find the period of the central minibrot. You can do this by iterating the 4 c corners of a box around the view, and seeing at which iteration count the 4 z points surround the origin. See http://www.mrob.com/pub/muency/period.html for a more detailed description of the method, I have an implementation which you can read here: http://code.mathr.co.uk/mandelbrot-numerics/blob/HEAD:/c/lib/m_d_box_period.c (there's also an arbitrary precision version in the same directory)
檢視
m-box-period
結果
usage: m-box-period precision center-re center-im radius maxperiod
示例
m-box-period 53 -0.8691524744 0.2556487868 1.25e-5 1000 48
這是一個使用從 Claude 的庫中提取的程式碼的單檔案程式。
/*
http://code.mathr.co.uk/mandelbrot-numerics/blob/HEAD:/c/bin/m-box-period.c
numerical algorithms related to the Mandelbrot set
by Claude Heiland-Allen
https://mathr.co.uk/blog/
GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007
http://code.mathr.co.uk/mandelbrot-numerics/blob/HEAD:/COPYING
mandelbrot-numerics
Numerical algorithms related to the Mandelbrot set: ray tracing, nucleus location, bond points, etc.
---------------------
Dependencies
sudo aptitude install libmpc-dev
---------------------------------
https://gitlab.com/adammajewski/mandelbrot-numerics-box-periood
Existing folder or Git repository
cd existing_folder
git init
git remote add origin git@gitlab.com:adammajewski/mandelbrot-numerics-box-periood.git
git add p.c
git commit -m "sasa"
git push -u origin master
-------------------
gcc p.c -lmpc -lmpfr -lgmp -Wall
./a.out precision center-re center-im radius maxperiod
./a.out 100 0.37496784 0.21687214 1.0 20
*/
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include <math.h>
#include <complex.h>
#include <stdbool.h>
#include <gmp.h>
#include <mpfr.h>
#include <mpc.h>
// static const double twopi = 6.283185307179586;
/* --------------- http://code.mathr.co.uk/mandelbrot-numerics/blob_plain/HEAD:/c/lib/m_d_util.h ---- */
static inline int sgn(double z) {
if (z > 0) { return 1; }
if (z < 0) { return -1; }
return 0;
}
static inline bool odd(int a) {
return a & 1;
}
static inline double cabs2(double _Complex z) {
return creal(z) * creal(z) + cimag(z) * cimag(z);
}
static inline bool cisfinite(double _Complex z) {
return isfinite(creal(z)) && isfinite(cimag(z));
}
static const double pi = 3.141592653589793;
static const double twopi = 6.283185307179586;
// last . takeWhile (\x -> 2 /= 2 + x) . iterate (/2) $ 1 :: Double
static const double epsilon = 4.440892098500626e-16;
// epsilon^2
static const double epsilon2 = 1.9721522630525295e-31;
/* --------------- utils.h -------------------------------------- */
static inline bool arg_precision(const char *arg, bool *native, int *bits) {
if (0 == strcmp("double", arg)) {
*native = true;
*bits = 53;
return true;
} else {
char *check = 0;
errno = 0;
long int li = strtol(arg, &check, 10);
bool valid = ! errno && arg != check && ! *check;
int i = li;
if (valid && i > 1) {
*native = false;
*bits = i;
return true;
}
}
return false;
}
static inline bool arg_double(const char *arg, double *x) {
char *check = 0;
errno = 0;
double d = strtod(arg, &check);
if (! errno && arg != check && ! *check) {
*x = d;
return true;
}
return false;
}
static inline bool arg_int(const char *arg, int *x) {
char *check = 0;
errno = 0;
long int li = strtol(arg, &check, 10);
if (! errno && arg != check && ! *check) {
*x = li;
return true;
}
return false;
}
static inline bool arg_rational(const char *arg, mpq_t x) {
int ok = mpq_set_str(x, arg, 10);
mpq_canonicalize(x);
return ok == 0;
}
static inline bool arg_mpfr(const char *arg, mpfr_t x) {
return 0 == mpfr_set_str(x, arg, 10, MPFR_RNDN);
}
static inline bool arg_mpc(const char *re, const char *im, mpc_t x) {
int ok
= mpfr_set_str(mpc_realref(x), re, 10, MPFR_RNDN)
+ mpfr_set_str(mpc_imagref(x), im, 10, MPFR_RNDN);
return ok == 0;
}
/* -----------c/lib/m_d_box_period.c -----------------------*/
static double cross(double _Complex a, double _Complex b) {
return cimag(a) * creal(b) - creal(a) * cimag(b);
}
static bool crosses_positive_real_axis(double _Complex a, double _Complex b) {
if (sgn(cimag(a)) != sgn(cimag(b))) {
double _Complex d = b - a;
int s = sgn(cimag(d));
int t = sgn(cross(d, a));
return s == t;
}
return false;
}
static bool surrounds_origin(double _Complex a, double _Complex b, double _Complex c, double _Complex d) {
return odd
( crosses_positive_real_axis(a, b)
+ crosses_positive_real_axis(b, c)
+ crosses_positive_real_axis(c, d)
+ crosses_positive_real_axis(d, a)
);
}
typedef struct m_d_box_period {
double _Complex c[4];
double _Complex z[4];
int p;
} m_d_box_period_t ;
m_d_box_period_t *m_d_box_period_new(double _Complex center, double radius) {
m_d_box_period_t *box = (m_d_box_period_t *) malloc(sizeof(*box));
if (! box) {
return 0;
}
box->z[0] = box->c[0] = center + ((-radius) + I * (-radius));
box->z[1] = box->c[1] = center + (( radius) + I * (-radius));
box->z[2] = box->c[2] = center + (( radius) + I * ( radius));
box->z[3] = box->c[3] = center + ((-radius) + I * ( radius));
box->p = 1;
return box;
}
void m_d_box_period_delete(m_d_box_period_t *box) {
if (box) {
free(box);
}
}
extern bool m_d_box_period_step(m_d_box_period_t *box) {
if (! box) {
return false;
}
bool ok = true;
for (int i = 0; i < 4; ++i) {
box->z[i] = box->z[i] * box->z[i] + box->c[i];
ok = ok && cisfinite(box->z[i]);
}
box->p = box->p + 1;
return ok;
}
extern bool m_d_box_period_have_period(const m_d_box_period_t *box) {
if (! box) {
return true;
}
return surrounds_origin(box->z[0], box->z[1], box->z[2], box->z[3]);
}
extern int m_d_box_period_get_period(const m_d_box_period_t *box) {
if (! box) {
return 0;
}
return box->p;
}
extern int m_d_box_period_do(double _Complex center, double radius, int maxperiod) {
m_d_box_period_t *box = m_d_box_period_new(center, radius);
if (! box) {
return 0;
}
int period = 0;
for (int i = 0; i < maxperiod; ++i) {
if (m_d_box_period_have_period(box)) {
period = m_d_box_period_get_period(box);
break;
}
if (! m_d_box_period_step(box)) {
break;
}
}
m_d_box_period_delete(box);
return period;
}
/* -----------http://code.mathr.co.uk/mandelbrot-numerics/blob_plain/HEAD:/c/lib/m_r_box_period.c -----------------------*/
static void cross_r(mpfr_t out, const mpc_t a, const mpc_t b, mpfr_t t) {
mpfr_mul(out, mpc_imagref(a), mpc_realref(b), MPFR_RNDN);
mpfr_mul(t, mpc_realref(a), mpc_imagref(b), MPFR_RNDN);
mpfr_sub(out, out, t, MPFR_RNDN);
}
static bool crosses_positive_real_axis_r(const mpc_t a, const mpc_t b, mpc_t d, mpfr_t t1, mpfr_t t2) {
if (mpfr_sgn(mpc_imagref(a)) != mpfr_sgn(mpc_imagref(b))) {
// d = b - a;
mpc_sub(d, b, a, MPC_RNDNN);
// s = sgn(cimag(d));
int s = mpfr_sgn(mpc_imagref(d));
// t = sgn(cross(d, a));
cross_r(t1, d, a, t2);
int t = mpfr_sgn(t1);
return s == t;
}
return false;
}
static bool surrounds_origin_r(const mpc_t a, const mpc_t b, const mpc_t c, const mpc_t d, mpc_t e, mpfr_t t1, mpfr_t t2) {
return odd
( crosses_positive_real_axis_r(a, b, e, t1, t2)
+ crosses_positive_real_axis_r(b, c, e, t1, t2)
+ crosses_positive_real_axis_r(c, d, e, t1, t2)
+ crosses_positive_real_axis_r(d, a, e, t1, t2)
);
}
typedef struct m_r_box_period {
mpc_t c[4];
mpc_t z[4];
int p;
mpc_t e;
mpfr_t t1;
mpfr_t t2;
} m_r_box_period_t;
extern m_r_box_period_t *m_r_box_period_new(const mpc_t center, const mpfr_t radius) {
m_r_box_period_t *box = (m_r_box_period_t *) malloc(sizeof(*box));
if (! box) {
return 0;
}
// prec
mpfr_prec_t precr, preci, prec;
mpc_get_prec2(&precr, &preci, center);
prec = precr > preci ? precr : preci;
// init
for (int i = 0; i < 4; ++i) {
mpc_init2(box->c[i], prec);
mpc_init2(box->z[i], prec);
}
mpc_init2(box->e, prec);
mpfr_init2(box->t1, prec);
mpfr_init2(box->t2, prec);
// box
mpfr_sub(mpc_realref(box->c[0]), mpc_realref(center), radius, MPFR_RNDN);
mpfr_sub(mpc_imagref(box->c[0]), mpc_imagref(center), radius, MPFR_RNDN);
mpfr_add(mpc_realref(box->c[1]), mpc_realref(center), radius, MPFR_RNDN);
mpfr_sub(mpc_imagref(box->c[1]), mpc_imagref(center), radius, MPFR_RNDN);
mpfr_add(mpc_realref(box->c[2]), mpc_realref(center), radius, MPFR_RNDN);
mpfr_add(mpc_imagref(box->c[2]), mpc_imagref(center), radius, MPFR_RNDN);
mpfr_sub(mpc_realref(box->c[3]), mpc_realref(center), radius, MPFR_RNDN);
mpfr_add(mpc_imagref(box->c[3]), mpc_imagref(center), radius, MPFR_RNDN);
mpc_set(box->z[0], box->c[0], MPC_RNDNN);
mpc_set(box->z[1], box->c[1], MPC_RNDNN);
mpc_set(box->z[2], box->c[2], MPC_RNDNN);
mpc_set(box->z[3], box->c[3], MPC_RNDNN);
box->p = 1;
return box;
}
extern void m_r_box_period_delete(m_r_box_period_t *box) {
if (box) {
for (int i = 0; i < 4; ++i) {
mpc_clear(box->c[i]);
mpc_clear(box->z[i]);
}
mpc_clear(box->e);
mpfr_clear(box->t1);
mpfr_clear(box->t2);
free(box);
}
}
extern bool m_r_box_period_step(m_r_box_period_t *box) {
if (! box) {
return false;
}
bool ok = true;
for (int i = 0; i < 4; ++i) {
// box->z[i] = box->z[i] * box->z[i] + box->c[i];
mpc_sqr(box->z[i], box->z[i], MPC_RNDNN);
mpc_add(box->z[i], box->z[i], box->c[i], MPC_RNDNN);
ok = ok && mpfr_number_p(mpc_realref(box->z[i])) && mpfr_number_p(mpc_imagref(box->z[i]));
}
box->p = box->p + 1;
return ok;
}
extern bool m_r_box_period_have_period(const m_r_box_period_t *box) {
if (! box) {
return true;
}
m_r_box_period_t *ubox = (m_r_box_period_t *) box; // const cast
return surrounds_origin_r(box->z[0], box->z[1], box->z[2], box->z[3], ubox->e, ubox->t1, ubox->t2);
}
extern int m_r_box_period_get_period(const m_r_box_period_t *box) {
if (! box) {
return 0;
}
return box->p;
}
extern int m_r_box_period_do(const mpc_t center, const mpfr_t radius, int maxperiod) {
m_r_box_period_t *box = m_r_box_period_new(center, radius);
if (! box) {
return 0;
}
int period = 0;
for (int i = 0; i < maxperiod; ++i) {
if (m_r_box_period_have_period(box)) {
period = m_r_box_period_get_period(box);
break;
}
if (! m_r_box_period_step(box)) {
break;
}
}
m_r_box_period_delete(box);
return period;
}
/*--------------- */
/*
c:0.37496784+%i*0.21687214;
./a.out 100 0.37496784 0.21687214 1.0 20
./a.out 200 -1.121550281113895 +0.265176187855967 0.005 40
so radius = domain size
*/
static void usage(const char *progname) {
fprintf
( stderr
, "usage: %s precision center-re center-im radius maxperiod\n"
, progname
);
}
/* ========================= mAIN =========================== */
extern int main(int argc, char **argv) {
if (argc != 6) {
usage(argv[0]);
return 1;
}
bool native = true;
int bits = 0;
if (! arg_precision(argv[1], &native, &bits)) { return 1; }
if (native) {
double cre = 0;
double cim = 0;
double radius = 0;
int maxperiod = 0;
if (! arg_double(argv[2], &cre)) { return 1; }
if (! arg_double(argv[3], &cim)) { return 1; }
if (! arg_double(argv[4], &radius)) { return 1; }
if (! arg_int(argv[5], &maxperiod)) { return 1; }
int period = m_d_box_period_do(cre + I * cim, radius, maxperiod);
if (period > 0) {
printf("%d\n", period);
return 0;
}
} else {
mpc_t center;
mpfr_t radius;
int maxperiod = 0;
mpc_init2(center, bits);
mpfr_init2(radius, bits);
if (! arg_mpc(argv[2], argv[3], center)) { return 1; }
if (! arg_mpfr(argv[4], radius)) { return 1; }
if (! arg_int(argv[5], &maxperiod)) { return 1; }
int period = m_r_box_period_do(center, radius, maxperiod);
if (period > 0) {
printf("%d\n", period);
}
mpc_clear(center);
mpfr_clear(radius);
return period <= 0;
}
return 1;
}
Mandelbrot 集中的內部座標
步驟：C 在整個過程中都是固定的 - 它決定了我們要著色的畫素。
第一步是找到一個特殊的 z（也許叫它 w），使得從 w 開始，對 z 進行 p 次 z^2 + c 迭代，可以得到 w。這是 c 的極限迴圈吸引子。這是透過 m_d_attractor 完成的，使用牛頓法求解 f_c^p(z) - z = 0。
然後，要找到 c 的內部座標，取你的 w，並將其迭代 p 次，計算關於 z 的導數
z = w; dz = 1; for (int i = 0; i < p; ++i)
dz = 2 * z * dz z = z * z + c
然後內部座標就是迴圈後的 dz。如果 |dz| < 1，則 c 在週期為 p 的元件內部
您可以根據 dz 對畫素進行著色（例如，角度 -> 色調，半徑 -> 飽和度）
為了突出邊緣，您可以選擇使用 w 找到內部距離估計（m_d_interior_de）。
一個演算法，在給定特定雙曲分量和所需內部角度的情況下，找到 Mandelbrot 集邊界上的點。它涉及二元牛頓法求解 2 個方程組
在哪裡
- F 是函式（複數二次多項式）
- 是迭代多項式
- p 是目標分量的週期
- 其中 |r|≤1。
- θ 是所需的內部角度。
- 第一個方程描述週期點
- 第二個方程描述吸引週期點（乘數 b 的幅度 r 不大於 1）
得到的 c 是邊界上點的座標。它也可以修改以找到內部的點，只需將 b=re2πiθ 設定為 |r|≤1。
Claude 用來數值求解內部座標 b 的演算法
- 當 c 在 Mandelbrot 集的外部時，現在放棄；
- 對於每個週期 p，從 1 開始增加
- 使用一元複數牛頓法求解 Fp(z0,c)=z0 的 z0；
- 透過在 z0 處計算關於 z 的一階導數來求解 b；
- 如果 |b|≤1，則返回 b，否則繼續下一個 p。
complex double z = 0;
complex double c = 0;
m_d_interior(&z, &c, zre + I * zim, cre + I * cim, ir * cexp(I * twopi * it), period, maxsteps);
printf("z = %.16e%+.16e\t c = %.16e%+.16e\n", creal(z), cimag(z), creal(c), cimag(c));
無引數使用
m-interior
給出用法
usage: ./m-interior precision z-guess-re z-guess-im c-guess-re c-guess-im interior-r interior-t period maxsteps
| 引數 | 變數 | 描述 |
|---|---|---|
| 精度 | int bits | 二進位制數字位數 = 有效數字的寬度 = 精度。這裡可以使用
|
| z-guess-re | double zre | 牛頓法初始猜測 z_0 的實部，使用 z = 0 |
| z-guess-im | double zim | |
| c-guess-re | double cre | 牛頓法初始猜測 c_0 的實部，使用給定雙曲分量的核 |
| c-guess-im | double cim | c_0 的虛部 |
| interior-r | double ir | 示例 |
| interior-t | double it | 示例 |
| 週期 | int period | 給定雙曲分量的週期 |
| maxsteps | int maxsteps | 示例 |
z-guess = z-guess-re+z-guess-im*I c-guess = c-guess-re + c-guess=im*I
輸出
- 由週期和核描述的雙曲分量中，具有給定內部座標的點 c
- 週期點 z（核）
示例
./a.out double 0 0 0 0 1 0 1 100 5.0000000000000000e-01 0.0000000000000000e+00 2.5000000000000000e-01 0.0000000000000000e+00
如何閱讀它
- z = 0.5（對於 c = 1/4 的不動點）
- c = 0.25
通用演算法
- 選擇週期
- 從週期計算/選擇中心 c = 初始近似值
- 計算以 c 為中心的元件的大小。從大小計算精度
- 選擇 t = p/q 和 r
- 計算
來自 bin 目錄的程式 m-interior
原始碼
- mandelbrot-numerics/c/bin/ m-interior.c
- mandelbrot-numerics/c/lib/m_d_interior.c
- mandelbrot-numerics/c/lib/m_r_interior.c
函式
- m_d_interior(&z, &c, zre + I * zim, cre + I * cim, ir * cexp(I * twopi * it), period, maxsteps);
- m_r_interior(z, c, z, c, interior, period, maxsteps);
#include <mandelbrot-numerics.h>
#include "m_d_util.h"
extern m_newton m_d_interior_step(double _Complex *z_out, double _Complex *c_out, double _Complex z_guess, double _Complex c_guess, double _Complex interior, int period) {
double _Complex c = c_guess;
double _Complex z = z_guess;
double _Complex dz = 1;
double _Complex dc = 0;
double _Complex dzdz = 0;
double _Complex dcdz = 0;
for (int p = 0; p < period; ++p) {
dcdz = 2 * (z * dcdz + dc * dz);
dzdz = 2 * (z * dzdz + dz * dz);
dc = 2 * z * dc + 1;
dz = 2 * z * dz;
z = z * z + c;
}
double _Complex det = (dz - 1) * dcdz - dc * dzdz;
double _Complex z_new = z_guess - (dcdz * (z - z_guess) - dc * (dz - interior)) / det;
double _Complex c_new = c_guess - ((dz - 1) * (dz - interior) - dzdz * (z - z_guess)) / det;
if (cisfinite(z_new) && cisfinite(c_new)) {
*z_out = z_new;
*c_out = c_new;
if (cabs2(z_new - z_guess) <= epsilon2 && cabs2(c_new - c_guess) <= epsilon2) {
return m_converged;
} else {
return m_stepped;
}
} else {
*z_out = z_guess;
*c_out = c_guess;
return m_failed;
}
}
extern m_newton m_d_interior(double _Complex *z_out, double _Complex *c_out, double _Complex z_guess, double _Complex c_guess, double _Complex interior, int period, int maxsteps) {
m_newton result = m_failed;
double _Complex z = z_guess;
double _Complex c = c_guess;
for (int i = 0; i < maxsteps; ++i) {
if (m_stepped != (result = m_d_interior_step(&z, &c, z, c, interior, period))) {
break;
}
}
*z_out = z;
*c_out = c;
return result;
}
描述[12]
找到一個特殊的 z（也許叫做 w），使得從 w 開始，對 z 進行 p 次迭代，每次迭代為 z^2 + c，最終回到 w。這是 c 的極限迴圈吸引子。這是使用牛頓法來求解 f_c^p(z) - z = 0。
usage: %s precision z-guess-re z-guess-im c-re c-im period maxsteps
找到一個週期性 z 點（動力學平面上的點）在 z-guess 附近
- 對於
- 週期
- c= c-re + c-im* i
- 使用 牛頓法，其中
- z-guess = z-guess-re + z-guess-im* i（根的初始近似值）
- 最大步數 = maxsteps
- 使用有效數字寬度 = 精度 的數字。這裡可以使用
- double 或 53
- 正整數 = MPFR 的位數
./m-attractor double 0 0 -0.965 0.085 2 100
或者
./m-attractor 53 0 0 -0.965 0.085 2 100
結果
-2.7669310528133408e-02 -8.9979332165608591e-02
精度在 m-attractor 中手動設定，讓使用者進行實驗。
引數外部射線:[13]
- 向內追蹤（從無窮大到邊界）
- 使用 牛頓法
When tracing inwards, one peels off the most-significant bit (aka angle doubling) each time the ray crosses a dwell band (integer part of normalized iteration count increases by 1).
m-exray-in usage: m-exray-in precision angle sharpness maxsteps
m-exray-in double 1/3 4 200
m-exray-in double 0 4 200 2.5405429038323935e-01 0.0000000000000000e+00
原始碼
- c/bin/m-exray-in.c = 展示庫函式用法的程式
- c/lib/m_d_exray_in.c = 雙精度（機器雙精度）版本的程式
- c/lib/m_r_exray_in.c = mpfr 版本的程式 =（根據需要動態更改）精度
變數
- native（如果為真，則使用 double = 本機精度）
- precision（可以是 double，參見 arg_precision）
雙精度版本
m_newton m_d_exray_in_step(m_d_exray_in *ray, int maxsteps) {
if (ray->j >= ray->sharpness) {
mpq_mul_2exp(ray->angle, ray->angle, 1);
if (mpq_cmp_ui(ray->angle, 1, 1) >= 0) {
mpq_sub(ray->angle, ray->angle, ray->one);
}
ray->k = ray->k + 1;
ray->j = 0;
}
double r = pow(ray->er, pow(0.5, (ray->j + 0.5) / ray->sharpness));
double a = twopi * mpq_get_d(ray->angle);
double _Complex target = r * (cos(a) + I * sin(a));
double _Complex c = ray->c;
for (int i = 0; i < maxsteps; ++i) {
double _Complex z = 0;
double _Complex dc = 0;
for (int p = 0; p <= ray->k; ++p) {
dc = 2 * z * dc + 1;
z = z * z + c;
}
double _Complex c_new = c - (z - target) / dc;
double d2 = cabs2(c_new - c);
if (cisfinite(c_new)) {
c = c_new;
} else {
break;
}
if (d2 <= epsilon2) {
break;
}
}
ray->j = ray->j + 1;
double d2 = cabs2(c - ray->c);
if (d2 <= epsilon2) {
ray->c = c;
return m_converged;
}
if (cisfinite(c)) {
ray->c = c;
return m_stepped;
} else {
return m_failed;
}
}
mpfr 版本
/*
it is a c console program which
- computes external parameter ray
- for complex quadratic polynomial fc(z) = z^2+c
- uses arbitrary precision ( mpfr) with dynamic precision adjustment
- uses Newton method ( described by T Kawahira) http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf
gcc -o rayin e.c -std=c99 -Wall -Wextra -pedantic -lmpfr -lgmp -lm
./rayin 1/3 25
./rayin 1/3 25 >ray13.txt
./rayin 1/3 25 | ./rescale 53 53 -0.75 0 1.5 0 >rray13.txt
program uses the code by Claude Heiland-Allen
from https://gitorious.org/maximus/book/
" ray_in computes 4 points per dwell band, and ray_out currently computes 16.
Moreover ray_in has dynamic precision adjustment and
ray_out doesn't (so you need a high precision to get through the thin
gaps) - fixing both of these is on my mental todo list. ray_out will
always be slower though, as it needs to compute extra points when
crossing dwell bands to accumulate angle bits (and to work correctly)."
What algorithm do you use for drawing external ray ?
"essentially the Newton's method one in
http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf
precision is automatic, effectively it checks how close dwell bands are
together and uses higher precision when they are narrow and lower
precision when they are wide. in general precision will increase along
the ray as it gets closer to the boundary. "
"periodic rays stop near the landing point at the current zoom level,
manually retrace if you zoom in and find the gap annoying. tracing
periodic rays is a lot slower than tracing preperiodic rays, especially
if you zoom in to the cusp - suggestions welcome for improvements here
(perhaps a larger radius than a single pixel for landing point nearness
threshold?). "
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <mpfr.h> // arbitrary precision
const float pif = 3.14159265358979323846264338327950288419716939937510f;
const double pi = 3.14159265358979323846264338327950288419716939937510;
const long double pil = 3.14159265358979323846264338327950288419716939937510l;
struct mandelbrot_external_ray_in {
mpq_t angle;
mpq_t one;
unsigned int sharpness; // number of steps to take within each dwell band
unsigned int precision; // delta needs this many bits of effective precision
unsigned int accuracy; // epsilon is scaled relative to magnitude of delta
double escape_radius;
mpfr_t epsilon2;
mpfr_t cx;
mpfr_t cy;
unsigned int j;
unsigned int k;
// temporaries
mpfr_t zx, zy, dzx, dzy, ux, uy, vx, vy, ccx, ccy, cccx, cccy;
};
// new : allocates memory, inits and sets variables
// input angle = external angle in turns
extern struct mandelbrot_external_ray_in *mandelbrot_external_ray_in_new(mpq_t angle) {
struct mandelbrot_external_ray_in *r = malloc(sizeof(struct mandelbrot_external_ray_in));
mpq_init(r->angle);
mpq_set(r->angle, angle);
mpq_init(r->one);
mpq_set_ui(r->one, 1, 1);
r->sharpness = 4;
r->precision = 4;
r->accuracy = 4;
r->escape_radius = 65536.0;
mpfr_init2(r->epsilon2, 53);
mpfr_set_ui(r->epsilon2, 1, GMP_RNDN);
// external angle in radions
double a = 2.0 * pi * mpq_get_d(r->angle);
// initial value c =
mpfr_init2(r->cx, 53);
mpfr_init2(r->cy, 53);
mpfr_set_d(r->cx, r->escape_radius * cos(a), GMP_RNDN);
mpfr_set_d(r->cy, r->escape_radius * sin(a), GMP_RNDN);
r->k = 0;
r->j = 0;
// initialize temporaries
mpfr_inits2(53, r->zx, r->zy, r->dzx, r->dzy, r->ux, r->uy, r->vx, r->vy, r->ccx, r->ccy, r->cccx, r->cccy, (mpfr_ptr) 0);
return r;
}
// delete
extern void mandelbrot_external_ray_in_delete(struct mandelbrot_external_ray_in *r) {
mpfr_clear(r->epsilon2);
mpfr_clear(r->cx);
mpfr_clear(r->cy);
mpq_clear(r->angle);
mpq_clear(r->one);
mpfr_clears(r->zx, r->zy, r->dzx, r->dzy, r->ux, r->uy, r->vx, r->vy, r->ccx, r->ccy, r->cccx, r->cccy, (mpfr_ptr) 0);
free(r);
}
// step
extern int mandelbrot_external_ray_in_step(struct mandelbrot_external_ray_in *r) {
if (r->j >= r->sharpness) {
mpq_mul_2exp(r->angle, r->angle, 1);
if (mpq_cmp_ui(r->angle, 1, 1) >= 0) {
mpq_sub(r->angle, r->angle, r->one);
}
r->k++;
r->j = 0;
}
// r0 <- er ** ((1/2) ** ((j + 0.5)/sharpness))
// t0 <- r0 cis(2 * pi * angle)
double r0 = pow(r->escape_radius, pow(0.5, (r->j + 0.5) / r->sharpness));
double a0 = 2.0 * pi * mpq_get_d(r->angle);
double t0x = r0 * cos(a0);
double t0y = r0 * sin(a0);
// c <- r->c
mpfr_set(r->ccx, r->cx, GMP_RNDN);
mpfr_set(r->ccy, r->cy, GMP_RNDN);
for (unsigned int i = 0; i < 64; ++i) { // FIXME arbitrary limit
// z <- 0
// dz <- 0
mpfr_set_ui(r->zx, 0, GMP_RNDN);
mpfr_set_ui(r->zy, 0, GMP_RNDN);
mpfr_set_ui(r->dzx, 0, GMP_RNDN);
mpfr_set_ui(r->dzy, 0, GMP_RNDN);
// iterate
for (unsigned int p = 0; p <= r->k; ++p) {
// dz <- 2 z dz + 1
mpfr_mul(r->ux, r->zx, r->dzx, GMP_RNDN);
mpfr_mul(r->uy, r->zy, r->dzy, GMP_RNDN);
mpfr_mul(r->vx, r->zx, r->dzy, GMP_RNDN);
mpfr_mul(r->vy, r->zy, r->dzx, GMP_RNDN);
mpfr_sub(r->dzx, r->ux, r->uy, GMP_RNDN);
mpfr_add(r->dzy, r->vx, r->vy, GMP_RNDN);
mpfr_mul_2ui(r->dzx, r->dzx, 1, GMP_RNDN);
mpfr_mul_2ui(r->dzy, r->dzy, 1, GMP_RNDN);
mpfr_add_ui(r->dzx, r->dzx, 1, GMP_RNDN);
// z <- z^2 + c
mpfr_sqr(r->ux, r->zx, GMP_RNDN);
mpfr_sqr(r->uy, r->zy, GMP_RNDN);
mpfr_sub(r->vx, r->ux, r->uy, GMP_RNDN);
mpfr_mul(r->vy, r->zx, r->zy, GMP_RNDN);
mpfr_mul_2ui(r->vy, r->vy, 1, GMP_RNDN);
mpfr_add(r->zx, r->vx, r->ccx, GMP_RNDN);
mpfr_add(r->zy, r->vy, r->ccy, GMP_RNDN);
}
// c' <- c - (z - t0) / dz
mpfr_sqr(r->ux, r->dzx, GMP_RNDN);
mpfr_sqr(r->uy, r->dzy, GMP_RNDN);
mpfr_add(r->vy, r->ux, r->uy, GMP_RNDN);
mpfr_sub_d(r->zx, r->zx, t0x, GMP_RNDN);
mpfr_sub_d(r->zy, r->zy, t0y, GMP_RNDN);
mpfr_mul(r->ux, r->zx, r->dzx, GMP_RNDN);
mpfr_mul(r->uy, r->zy, r->dzy, GMP_RNDN);
mpfr_add(r->vx, r->ux, r->uy, GMP_RNDN);
mpfr_div(r->ux, r->vx, r->vy, GMP_RNDN);
mpfr_sub(r->cccx, r->ccx, r->ux, GMP_RNDN);
mpfr_mul(r->ux, r->zy, r->dzx, GMP_RNDN);
mpfr_mul(r->uy, r->zx, r->dzy, GMP_RNDN);
mpfr_sub(r->vx, r->ux, r->uy, GMP_RNDN);
mpfr_div(r->uy, r->vx, r->vy, GMP_RNDN);
mpfr_sub(r->cccy, r->ccy, r->uy, GMP_RNDN);
// delta^2 = |c' - c|^2
mpfr_sub(r->ux, r->cccx, r->ccx, GMP_RNDN);
mpfr_sub(r->uy, r->cccy, r->ccy, GMP_RNDN);
mpfr_sqr(r->vx, r->ux, GMP_RNDN);
mpfr_sqr(r->vy, r->uy, GMP_RNDN);
mpfr_add(r->ux, r->vx, r->vy, GMP_RNDN);
// enough_bits = 0 < 2 * prec(eps^2) + exponent eps^2 - 2 * precision
int enough_bits = 0 < 2 * (mpfr_get_prec(r->epsilon2) - r->precision) + mpfr_get_exp(r->epsilon2);
if (enough_bits) {
// converged = delta^2 < eps^2
int converged = mpfr_less_p(r->ux, r->epsilon2);
if (converged) {
// eps^2 <- |c' - r->c|^2 >> (2 * accuracy)
mpfr_sub(r->ux, r->cccx, r->cx, GMP_RNDN);
mpfr_sub(r->uy, r->cccy, r->cy, GMP_RNDN);
mpfr_sqr(r->vx, r->ux, GMP_RNDN);
mpfr_sqr(r->vy, r->uy, GMP_RNDN);
mpfr_add(r->ux, r->vx, r->vy, GMP_RNDN);
mpfr_div_2ui(r->epsilon2, r->ux, 2 * r->accuracy, GMP_RNDN);
// j <- j + 1
r->j = r->j + 1;
// r->c <- c'
mpfr_set(r->cx, r->cccx, GMP_RNDN);
mpfr_set(r->cy, r->cccy, GMP_RNDN);
return 1;
}
} else {
// bump precision
mpfr_prec_t prec = mpfr_get_prec(r->cx) + 32;
mpfr_prec_round(r->cx, prec, GMP_RNDN);
mpfr_prec_round(r->cy, prec, GMP_RNDN);
mpfr_prec_round(r->epsilon2, prec, GMP_RNDN);
mpfr_set_prec(r->ccx, prec);
mpfr_set_prec(r->ccy, prec);
mpfr_set_prec(r->cccx, prec);
mpfr_set_prec(r->cccy, prec);
mpfr_set_prec(r->zx, prec);
mpfr_set_prec(r->zy, prec);
mpfr_set_prec(r->dzx, prec);
mpfr_set_prec(r->dzy, prec);
mpfr_set_prec(r->ux, prec);
mpfr_set_prec(r->uy, prec);
mpfr_set_prec(r->vx, prec);
mpfr_set_prec(r->vy, prec);
i = 0;
}
mpfr_set(r->ccx, r->cccx, GMP_RNDN);
mpfr_set(r->ccy, r->cccy, GMP_RNDN);
}
return 0;
}
// get
extern void mandelbrot_external_ray_in_get(struct mandelbrot_external_ray_in *r, mpfr_t x, mpfr_t y) {
mpfr_set_prec(x, mpfr_get_prec(r->cx));
mpfr_set(x, r->cx, GMP_RNDN);
mpfr_set_prec(y, mpfr_get_prec(r->cy));
mpfr_set(y, r->cy, GMP_RNDN);
}
void usage(const char *progname) {
fprintf(stderr,
"usage: %s angle depth\n"
"outputs space-separated complex numbers on stdout\n"
"example : \n %s 1/3 25\n"
"or \n %s 1/3 25 > r.txt\n",
progname, progname, progname);
}
// needs input
int main(int argc, char **argv) {
// read and check input
if (argc < 3) { usage(argv[0]); return 1; }
mpq_t angle;
mpq_init(angle);
mpq_set_str(angle, argv[1], 0);
mpq_canonicalize(angle);
int depth = atoi(argv[2]);
struct mandelbrot_external_ray_in *ray = mandelbrot_external_ray_in_new(angle);
mpfr_t cre, cim;
mpfr_init2(cre, 53);
mpfr_init2(cim, 53);
// compute and send output
for (int i = 0; i < depth * 4; ++i) {
mandelbrot_external_ray_in_step(ray);
mandelbrot_external_ray_in_get(ray, cre, cim);
// send new c to output
mpfr_out_str(0, 10, 0, cre, GMP_RNDN); putchar(' ');
mpfr_out_str(0, 10, 0, cim, GMP_RNDN); putchar('\n');
}
// clear
mpfr_clear(cre);
mpfr_clear(cim);
mandelbrot_external_ray_in_delete(ray);
mpq_clear(angle);
return 0;
}
The trick when tracing outwards is to prepend bits when crossing dwell bands, depending if the outer cell was entered from its left or right inner cell. A picture may make it clearer:
引數外部射線
- 向外追蹤（從點 c（從曼德勃羅集的外部）到無窮大）
- 使用 牛頓法
"there are m_r_exray_out_* functions that can be used to get the external angle by tracing an external ray. but it is O(dwell^2), too slow to be practical beyond a few 1000, and it does not yet have adaptive precision so it cam get stuck in tight gaps between spirals..."
檢查哪個射線
- 穿過給定的點 c
- 落在該點 c 上
引數
- precision：輸入字串被轉換為整數，單位為位“用 1000 替換 double 以使用 1000 位而不是本機 53 位”。它可以是字串“double”= 本機精度。
- c-re = 複數 c 的實部
- c-im = 複數 c 的虛部
- int sharpness（用 4 替換 8 以減少銳度，但不要太低，否則可能會導致程式崩潰）
- int maxdwell
- 僅用於格式化輸出（“如果您有所需的（預）週期，請使用它來美化角度格式”）
- preperiod（int）
- period（int）
用法
m-exray-out precision c-re c-im sharpness maxdwell preperiod period
示例
m-exray-out double -0.1 0.651 8 10000 100000 1 | tail -n 1
輸出（僅外部角的二進位制展開式）
.00100100100100100100100100100100100100100100100100100100100100100100100100()
或完整輸出
m-exray-out double -1.7904997268969142e-01 1.0856197533070304e+00 4 20 100 1 -1.7904997268969142e-01 1.0856197533070304e+00 -1.7935443511056054e-01 1.0854961235734493e+00 -1.7969015435633140e-01 1.0853555826931407e+00 -1.8051917357961708e-01 1.0850119514784207e+00 1 -1.8269681168554189e-01 1.0840986080730519e+00 -1.8930978756749961e-01 1.0809967205891908e+00 0 -1.9360134634853554e-01 1.0737664210007076e+00 1 -2.0801776265667046e-01 1.0704315884719302e+00 0 -2.4370145596832851e-01 1.0627226337492655e+00 1 -2.7278498707172455e-01 1.0400463032627758e+00 1 -3.6142246174842857e-01 9.5627327660306660e-01 -4.0875942404913329e-01 9.4318016425218687e-01 -4.5809426255241681e-01 9.9407738931502432e-01 1 -6.0223759082842288e-01 9.3674887784190841e-01 -8.0664626534435835e-01 8.8455840042031386e-01 0 -8.3743553612789834e-01 9.3674304615806725e-01 -8.7287981833135664e-01 9.9952074265609137e-01 -9.1381642842183941e-01 1.0750773818725967e+00 -9.6138941844723436e-01 1.1665050307228815e+00 0 -1.0172147464760770e+00 1.2782611560409487e+00 -1.0836516873794186e+00 1.4168346080502183e+00 -1.1642543320399381e+00 1.5917680433570660e+00 -1.2645279518271537e+00 1.8173272090846146e+00 0 -1.3932089082214374e+00 2.1153654197070240e+00 -1.5644798765083927e+00 2.5204286897782238e+00 -1.8019543103478257e+00 3.0891768094384502e+00 -2.1462162170666716e+00 3.9184694463323830e+00 1 -2.6699272827383060e+00 5.1817748128576184e+00 -3.5099957053533952e+00 7.2066815445348764e+00 -4.9407140183321250e+00 1.0650884836446950e+01 -7.5526171722332034e+00 1.6932133257514206e+01 0 -1.2727799761990893e+01 2.9370322256332429e+01 -2.4030360674257757e+01 5.6528056606204473e+01 -5.1753810887543764e+01 1.2313571132413860e+02 -1.2986129881788275e+02 3.1079020626913712e+02 1 -3.8951035577778663e+02 9.3459826885312373e+02 -1.4412515103402823e+03 3.4614095793299703e+03 -6.8365543576278487e+03 1.6423640211807055e+04 -4.3547867632624519e+04 1.0462387704049867e+05 0 -3.9380510085478675e+05 9.4611788566395245e+05 -5.4017847803638447e+06 1.2977803446720989e+07 -1.2161023621086851e+08 2.9216894171553040e+08 .01010001110101()
輸入
./ray_in ".(0101101011010110101101011010110101100)" 1000
示例
m-exray-out 100 -0.7432918908524301 0.1312405523087976 8 1000 24 4
結果
.010101010101010101010100(1010)
即
.01010101010101010101010(01)
external ray tracing implementation, that switches between: * perturbation methods (including acceleration with bilinear approximation) * arbitrary precision iterations depending on how close neighbouring points on the ray are. Trying to be correct in all cases but faster than just using arbitrary precision for the whole thing.
“使用週期 P 島的大小估計週期 2P 原子的位置，使用牛頓法更精確地找到其核和大小，重複”
計算實軸上的週期加倍級聯的雙曲分量的中心 c（因此虛部為 0）
- 週期 1 的中心：c = 0（省略）
- 週期 2 的中心：c = -1
- 週期 2^2= 4 的中心：c = -1.310702641336833008
- ...
- 週期 2^n 的中心
./m-feigenbaum
-1.000000000000000000e+00
-1.310702641336833008e+00
-1.381547484432061657e+00
-1.396945359704560685e+00
-1.400253081214782869e+00
-1.400961962944840655e+00
-1.401113804939776664e+00
-1.401146325826946537e+00
-1.401153290849922906e+00
-1.401154782546613964e+00
-1.401155102022467736e+00
-1.401155170444417619e+00
-1.401155185098302614e+00
-1.401155188236698823e+00
-1.401155188908866478e+00
-1.401155189052811778e+00
-1.401155189083657104e+00
-1.401155189090263598e+00
-1.401155189091694675e+00
-1.401155189093319819e+00
-1.401155189093314712e+00
-1.401155189093314712e+00
-1.401155189093314712e+00
-1.401155189093314712e+00
-1.401155189093314712e+00
-1.401155189093314712e+00
-1.401155189093314712e+00
-1.401155189093314712e+00
-1.401155189093314712e+00
小修改
printf("period = %d \t\t\t nucleus = %.18e \tsize = %.18e\n", period, creal(nucleus), size);
m-feigenbaum
period = 4 nucleus = -1.310702641336833008e+00 size = 1.179602276312223114e-01
period = 8 nucleus = -1.381547484432061657e+00 size = 2.590331788620404280e-02
period = 16 nucleus = -1.396945359704560685e+00 size = 5.574904503817590395e-03
period = 32 nucleus = -1.400253081214783091e+00 size = 1.195288587163239984e-03
period = 64 nucleus = -1.400961962944842210e+00 size = 2.560528805309172169e-04
period = 128 nucleus = -1.401113804939776664e+00 size = 5.484142131458955168e-05
period = 256 nucleus = -1.401146325826946537e+00 size = 1.174547734653405196e-05
period = 512 nucleus = -1.401153290849925570e+00 size = 2.515527294624666474e-06
period = 1024 nucleus = -1.401154782546613964e+00 size = 5.387491809938387059e-07
period = 2048 nucleus = -1.401155102022462628e+00 size = 1.153835913118543914e-07
period = 4096 nucleus = -1.401155170444413400e+00 size = 2.471163334462234791e-08
period = 8192 nucleus = -1.401155185098302614e+00 size = 5.292465190668938379e-09
period = 16384 nucleus = -1.401155188236713700e+00 size = 1.133481216591179550e-09
period = 32768 nucleus = -1.401155188908850047e+00 size = 2.427895317585669274e-10
period = 65536 nucleus = -1.401155189052770256e+00 size = 5.204848444681121846e-11
period = 131072 nucleus = -1.401155189083645558e+00 size = 1.116488091172383448e-11
period = 262144 nucleus = -1.401155189090176112e+00 size = 2.401481531760724511e-12
period = 524288 nucleus = -1.401155189091663589e+00 size = 6.463020532652650290e-13
period = 1048576 nucleus = -1.401155189092033737e+00 size = 3.635753686080000682e-14
period = 2097152 nucleus = -1.401155189093066022e+00 size = 9.366255511607845319e-19
period = 4194304 nucleus = -1.401155189093066022e+00 size = 2.943937947459073103e-26
period = 8388608 nucleus = -1.401155189093066022e+00 size = 5.438985544182134629e-40
period = 16777216 nucleus = -1.401155189093066022e+00 size = 5.967915741322236566e-68
period = 33554432 nucleus = -1.401155189093066022e+00 size = 8.736600707368986815e-124
period = 67108864 nucleus = -1.401155189093066022e+00 size = 1.222084981000085251e-235
period = 134217728 nucleus = -1.401155189093066022e+00 size = 0.000000000000000000e+00
period = 268435456 nucleus = -1.401155189093066022e+00 size = 0.000000000000000000e+00
雙精度似乎不夠
“使用大小估計，從週期 P 島的觸角尖端（從 -2+0i 重歸一化）開始估計位置，然後使用牛頓法從那裡找到週期 3P 島的核，重複” [14]
用法
m-feigenbaum3 re(location) im(location) periodFactor
m-feigenbaum3 -2 0 3 #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.659992713573920042e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.492879663528381684e+01,0.000000000000000000e+00 T = 43 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.526412200825988208e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.524573293871185342e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.524447443629173904e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.510191709584042741e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 4.825684575389948350e+01,0.000000000000000000e+00 T = 42 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 6.901913875598086001e+00,0.000000000000000000e+00 T = 21 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = inf,-nan T = -2147483648 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = -nan,-nan T = -2147483648 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = -nan,-nan T = -2147483648 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = -nan,-nan T = -2147483648 #endpreset
m-feigenbaum3 -2 0 3 period = 3 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.659992713573920042e+01 +0.000000000000000000e+00 T = 44 size = 1.000000000000000000e+00 period = 9 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.492879663528381684e+01 +0.000000000000000000e+00 T = 43 size = 1.903551591313245098e-02 period = 27 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.526412200825988208e+01 +0.000000000000000000e+00 T = 44 size = 3.464641937910663597e-04 period = 81 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.524573293871185342e+01 +0.000000000000000000e+00 T = 44 size = 6.269877899222824167e-06 period = 243 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.524447443629173904e+01 +0.000000000000000000e+00 T = 44 size = 1.134900148763894982e-07 period = 729 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.510191709584042741e+01 +0.000000000000000000e+00 T = 44 size = 2.054226124560898465e-09 period = 2187 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 4.825684575389948350e+01 +0.000000000000000000e+00 T = 42 size = 3.718249148787694584e-11 period = 6561 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 6.901913875598086001e+00 +0.000000000000000000e+00 T = 21 size = 6.727899706332398247e-13 period = 19683 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = inf -nan T = -2147483648 size = 1.225497697338603077e-14
列印
- 週期 + 1
- 距離
m-nearest-roots
5 -6.0077592739339081e-01
9 -1.2345690574077126e+00
17 -1.5948261888043576e+00
33 -1.7912740981919455e+00
65 -1.8940314652359109e+00
129 -1.9466050018630559e+00
257 -1.9731986196121953e+00
513 -1.9865731876029906e+00
1025 -1.9932800441369114e+00
2049 -1.9966383822549216e+00
4097 -1.9983187808738916e+00
8193 -1.9991592878359983e+00
complex double *c_out;
const complex double c_guess = 0.0;
m_d_nucleus(c_out, c_guess, period, maxsteps);
printf("period = %d \t c = %.16f using double\n", period, creal(*c_out)); //
示例
m-nucleus 53 -0.8691524744 0.2556487868 48 64 -8.6915874972342078e-01 2.5565708568620021e-01
使用 mandelbrot-numerics 程式碼在 c=i 附近找到一個週期為 20000 的小曼德勃羅集
time ( m-size 50000 `m-nucleus 50000 0 1 20000 32 | tee /dev/stderr` 20000 )
結果
Real: -7.9437163502071985489173816041278335565557165593422097986415443718666526638183975269482741283979864796150039328820600171824357494796904100157145276457381327038677304389065791544759044068226298284798485074795465139753011290574940547565632295268309965890736856044225816065002060784840401957436001592831073232116289224835376766024837019729052593614702740064116283928435871912023279827933594666780735810543887516779374159354732526673506932730255777860864410617962251986612293528344539369634355763385689632183660361344255846109415319987535013048216370406305683033823649556367343903980453528078550807179972805086726177286676432146600418321802012705947183289186915513177714879362378374485500397855616627594598759238285491434600163661656097093569019464247435079968347293144424227930353208249468172472673021579306852412165125553964864288143879390896150333653457265029413832025627226640151187411034335007580482643346655743497049637725677688986423224436940553132460564725885956480187539465814934618297444604422558707577928894218074677873014006311238296593575712979281521376321522028222864880290047385477309492930859941138933523459312865570596181676682065664641625195711959850159891958399655264891695837255722367048386224829361833770848853678382189241483407917363077785688704033915945145445141382032098503113386121141821618976435909826878012810311315238838927942886145795725258557067351714858363660176807474230411907332988496361646952820878046993327965167592505974342069923145251125405511298210735382002351444475111195883871537375357019225324966837329585172700415064929059719568669415835536305907605195015067551532717806048959958720723351540772045136711267008091434857778994767951971042965880533276128160757945671341024842631669739934842190421628745823057485321876759312698540929225608881625608367714997633793882137116252066659684867523320525998235559519150039569822804344024248346045328260662350506766959441224921002078134780987834208399572026751630341064076128585432804822899645121908649438079186104626223225086835026727973207053106933545605807581086173838972807582278495224764921308243782567040001796469511628492096380473864043528112882638850008499650397617471312161946511543108954047732214732604192518761590429256338997424173451213591713726558272002992188667406903054868988361837689051267067909007940396165686808795106868838064855005108089621233950746964278729824349462021919452828912238487948591999838204600089889473079741079883341554696124986558349029858769442566619743815502763844097701527868328322323759228519682883781772486449895023296886890395156702488270397371139695710695990682206272068513184347868985875355437735986097344366110220146098454669911676872422507466726262772844396687992023980393180412252620962951600694487648989414046125141359107741643310781724540269445598327469901679273898223648390876129705793738640048570537095417043836345968852339629450002567700189777577964846790537104618393106879939835732394904152494946398727231384913787249745823001738289311645727103810197887877515539458354319263104285589267937979420741297746874900374889388549539632838441649971129535778931268439396997823639463581253024656692066840729143014425564553465530544851212792187013660966565334102233995602996284005853152619476167135763137548181086648741615317739115246457134590756369220988612182423106888272870894941766951250687466714807759040424513113776919266305284693482786000012318355370203000510649927436937799090573931965327514602174194247869031450849888578824607292360533650605659626633451928857198284911861895518252866768510849271811271674466215339007494425799589266330368037494808523322072015355702363872561150182139487288949855852914835191151647214953580577992765768923589452140092207774849659839859152689402490049014646780450855312634212041907597487788211223125953582387474055558632251459696857778410866837388535466178076110887271500918602073599110008782828349828601931359617357488907006692730574882684291152057556083434843531655449794230421218496348063440957859520929997115553713223272459893888222599714500315885758226251226383907100867373946949159210357809409066172026448365305726879690385903084965357293868883959241495118228357461132946108953197235701605108487932988139914383636638205336622907254521211337215354959781371499838940575460882923252535688983896279577141026183725803250025941366898551995317381343295602432854880680913105417171403478847572564220151501945641257407256340713758082071177905993565684176799173836494504063303062708493883917533694417688017619038230941678760662664057011965910800231539615509203568524423393735430913683843408798930423887383733707108315296851537552773795174571699872575675465307644486366568477860016882096508306899350495745941389774331465508906343167441651902671583925484007761590126001875408698776037548173000537310514627103673719106723875477565011787438594421891881293376341836817387959998884620172811186225741103273771975399857174630178915209374630124944654809434747663998641540089249451044747343793912738684530410065341666811576503875144080041376635588445419247720870201032961075580784574643803834224040349515934580424790949301603070840039881116983632470720266625847077403948483754828340322244591513248736220612802960652020455915703123428617077016498405385864432461997065871572070257516341957020978947252938690548377219964733666078050655650666746125819484651038400639111767173754889123755628492910863471339403636445116498685130083472801719762230149643141939154039001471145670460001527706331210349655036747526294372592292737874421182976105207938866601322001196843303227985801688904265826419497214272904851901075941514972117647200278204061021869787050146340891622184720663076877881386068868416366064471354402886662288552667078706482937945747653944295880393480095834773836964538912825459369034140530990210203217740376593495298910340327466503746890581057555405215105328622050072523451061799089554563726997933998099494022843643356810171935875956083213474450964753871742597596389159807474860548937449082552764665987353842664468415487133534606684963715041757839071452744227881047305656416036812546802774367092104370634445506994254300152537889204524138962819925474418837039981003771901730655105777080171387232652844063902421538914319305968703705473681853824215543410439200761681559696856513162963772239637111929231912298561552761476051768494603150247970255274087139475417437666525318837699314477624013690658187766683186987513666300409400279403768729919447855449000797444258300988278145284834950271127363489449639121339411432241394689045849625215128792858854683514551608992163329455011784996662465048627222172830045774045434488939393408845069630045746752624798539349705914385735569651595250218651352938215626071449610093033419595771219480026325620655464072741468432773768618336398345601286278453384744487239861535583484271197660008924457670571819830930013272839428183593854035761689365351460247606170808699221660877760009044931352005199906138059024299151054829073698920465120899137076026463287302396273561452735416778600019151509061871192996649944544635251255611180939657319086749495750406650703523104631376966634580854744528028348644766661861687803874476503066257968690726504793327378255680115059130846355163503033322685363135207492701334721066801514843289941471860778761567160868619351770348032409538305377970923902100885197699029395920849667846810993094612702293706672531075824284278002289430678551378565535722768390905768325693280507976229296441994336567584568603213159337614179963913967140245071426406010126568343361152854104386871950971193764271019726152905089469444204032522564569281201767167121304137166303415015387588952606612266240324729892259808131866274224382895297824105879377663513842630692659814953190082184611299068317174185652027168995389944464563646470703089248951504008599661061822888980138794720478053691503842822034170421882089481794863846177177647112938213372331349756322623026418328473098549683499447951740645257409574863864331936413558808652101756422536878243221845141405047163092445446309652337506812857818627712672119535182802674487218494516703692576818047705180818990227687611349865697646700417084858530111203684655604862529459824809506713267533309643092879524182664685875269816703424729371134632023211821682751089970776093301423791969260953108431784555846379656677296878357053153268291167052511154697766977051318490201590272694053975247144309260438403419516059314927198140682502647753867140477723365586992936626937070645112420562978420494969180857648414886732128746833394793338951136536812507919552105463278901514836838501627409427945166072127389034634662362535768301094003301978784609662103275295887490983814104476638090203027305898662646192580113323630177677789960362015187337059766586565991561230871896736278580930538113756251722826451445081763396147220450655330971575966542768165629827950058414694457612107002519813422137935933949051802958001466676705218971585089638867365108649794271967367583738948555076277154455873640422096572950582246453668223767071601874641607213366687118498245978253390109444418995997750392780107313255139840820773405754956577553794289770409116409571420228693380230025899437603870877054542784281317586453244622402660818643831821462226891660685072457573117889485902077797412004329040165095691873500314917343659203512085778821361204452523870368136469226149736015867804245577929236392814077339316439392508508206520959727928573800885094153800553314652906939293080711662049154470341556580077893261447671240157808246582347200333618553406085052554826594036992297841269796755362033887196089006252179737733109866573340227539862034141198839388081519872826063624800513941456162727625429496413744959884317434642624473125087427423717708382301927586834616549672826513360395335642137272724139545454739463104449628322251582845457933561149241779746747083923968435596911958159633854978255431018163676494053663975309474153510434284245769525836241401382297602993004553331413623115085862521547918233116144303815147644592323844409534033312591272415940304330316480139161832153930361866639526736202903837917442316305999918735106124476264198373838671169480263093073548771272893575807108321917398785294490102026462289081419098788347663282918233719927616971290489000448605479370675888150101653436770635409683915707260833378384329978382632001535926336764826666274559700141470555541357684564795628678577384149810990074344641965294243434863683046415340191004534007747861984920997073615420680497487884717648552467181792474987367112740072255272190588800538041606909185987209170159028084428459572622258803142739029327530812085898685496500033506607514853966450810131101320181576443944528119520678262542516665667772589225658367409437310308191692681701747929941825717813914760436490027296973273696885927024209476169849577627774609821735288543579823873933746165418041525243007331797265820048090493269547699807443692298819880238271902411894022717673853039050368938062710759904516779434177910027164576717375460649435897796076715236692576082399547459593785193303729528593113309230666964408518409037311263890279582964898612945710648622089308353146591558950556457774096875906864050442683109026765888050991130358931886565546257805826690073777858095011429817715273547968084945928031301599621650496157700528460843729211007404968844131661389160468042607031052067033392469344850696050409683463254404597341777664046659525995421309405891229261642956098092899701872432918414448769927311569703420232182442026732210734846569104330914950450386817272059918946706401776070037996180381660685319130717417939927179286304468696802295571418406036768424662411467135246405133085173144139264971120263053400301163186239059701230700612814579728270543021200439105129200019275455983470874139432086875339191481441412783871068993677566173792832145720122415779284407893352668293527601354952001547147784537072825711418553102839482249087891314635317734225217087825443501818465384174890328302318165740210720929441444953628896415212119842220277809228686748800353964121521758377906665150446944597690876685717728001020923030903467397889577625987702929819042793216861981297432410751862996742939821733552376022465063980433476877095489742083667438482524358479920093362671073499261536442509635978938254572171219604223534931048801135996840036883391975783540196339746427425612440311802930181139735185695582618989988117958811316075569114339755674685677959928036900764099458526085419821675058081439927368204872326934818289759660177557064039914505754639463377998991196112023961998610668733702044390628198372939453058274311754372582831610070101869681001548640766315884308521449756188946074250519766336654408201104386777792941840610188681466688436569953569101322404299961245271858446078822593644050595237420612302915602324865912944449105260539387579372432085618031718618467399412498198037020264352404192035385525718637738445737173126887883204422338914094077117405266369732555437529835445785310643028012476150339747654769477837975324497094712742252873395275094111602092150983608511685837417691370283134112558707615144700568498489604300808338456202632290195011163561765905032989492777753130863693458023240541703789339666010447835045257252715361960269883076313236873064974259132787196196424346621374024293770156313978590441642103302299007768492108804953992601544447761538374942979300239461633293490859599063452938858711250267104956507397978474403844807309379849329624027364559546136877836071647522487335680653875008575928199095693277072357337411294063778987393546727728838736254506361668708240985190446787175682032349692908267897330716513026626394037531214887715636133426679029411503800223468119703782419347164191227385529597982634049525824249489381004584297761575054940323722350162865941897918756430316221819391580856455929590012212693971767718111156134625119611537902990743877562663928969788178546328357672482203921172030547393621262720602994243945777523017683019828452189053409346073165497827641541546661811362000069860698473847235166780178734237150200969007414767921753122493459975351206764098974709508735766435838682683431012232653877433197670329155560765634097241957053112237676513811748995099145494510366979493271174202393093949872617145305822631758071249330647970105326491208165384679106654911636291369219591774072389480094363848652230439114251230301727437540976950510736579784395802523606453713668228567047296997827776205239775215773191504628314406040463975203778623513875837323952073230041221445555570861593964350888214523314641085730437948149961122240579538109762211134983888224068986964032346290276113268596102816682853305085149030325552748215495078408622679995786080249473690737026626873405306052299932424760848269117260311258313117827223249824645230833114931474245579229315184388294163598032419249401202852156302446910810247963489869968887090519584584819820604231857986183254977904943759345535508926487273426498495065483968097609125809575651864291878638656751118948252709080366966631434561807538043142623374192000538765039906244673936922990047910893164884788344826422647704033143076607040449945061462449827844843237262351978322868811378865263527055247746845940099739111104279244801339852071272925084450912408790993841951019053768851872131163025230275e-7527 Imag: 9.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999991734492077330513689380886733519291674946671408282583644209407630575650409140342995770119063424352210271488566807058196756290031171701011701529100242419560679125380260046815348728362732098903518927842731522442015459899781490157079648835177627682202584659756122982363439113383649152007085924164482218692663011713033118422900677796130976701520211513973179804073856412336912323471275216141504900024638942128288839288824544726257598363037946954747139286140996713614519984757554081285137303793792145483088100948602288971711699818733272109930278128457054114530565415061361089457611033887291209920519084835617420429501562494329154552253147974514089623762774850187343594170152609219811762528104878460467096558776750623330946583444771677244364552796775144576896199277889280199535802388272476453809146195137520963756883307724481370820965867991787274442477494687690041794672045778026327225231447931749253744034412377634436209688157234219689149647284294317994921943258044571538072204447829486195849511596206224271784639410492783745179840157550028751542998761457025592735390256235705578527400511578961509524795765723846769216602534827598544122595055537226394226769431805638838745187514550034115461736580630942772164850138959148037936542189861498291521024384584123266494676747874024658839560362669842111747120028879375611008493691777571650981915670778877856645267853315285564690838429624046285747718787113149379509342042014222183841341852480309391774998684638396165511282456072444392493337309576584765407098415136785336850803964219190455543773565500777469684197751437542701373860103750789485898916431227179802485712595100542939658261974899323094584451988998658155488568870156711079534034726043188161080239900290517603670651552229213022394434669013398782592380368779183995497512609086992015378492374681966407092817332227547134424773876378525948648281148884775757472839432230631294074294857393651417886848443276428721134947690017433025325153248913554312348436870298461202014417691273225926050595655721157069521076214438687068333072136368413085150575362642380849505084969356878401879793263084179759550747507234523232524796967394646964932368873283962438332014394396863216097188067789543751355227029409828459507781736661768248386270817788738777462668450623755227855820114090319893003480371488432847626428399909880285573597465643610948439912464793498238988585822403705164724699464120093570367708183907545235654250619728658435398562469329842236925494721158612870658140940335567711385774552936065430987224064165463886334630367918901941265514184193460454762359263731196410400313787736101813913548190929105563969609759961286142020359681386163860459140037407728122634484712324766777888892902726372897663970591495119930657196797922969697210887944061663587643330053324801023593783956257456402354277193486606699249716660922975324918446780356978648595624510216079738304712240699383022296050307095326386181436021280821946461416911955375685450958489937715061871379614514046821562054907141306926601353784074703466872196180068651947552564031583073723325119005346348879674247302910116126291796234873627140147503137007974533422374524109527394238215360495340869525049828516165693842590343161532628841468410993746576931920685810272338948639553304457274026775879365952258661615059382775614181775292223518160756734117381397998310327294281664957921427253102145159956251792002855401711344356095059093958391981346288299586778567037745773029382822174296514329430014003230736855824042705301686236442231943132404408336068221178832509202932258219292166235157083978711447788498538296242489741089788080914827695688672104190838035163401216262600615305503117497776805939741146793074105899920727912025006743184132116016580499916303874955586731672402804758464494865562667287513471323587649568052936006244205760085680463430770593090850578873410241227234316091614508294678455447144047019836067353592195113713166048538989096099123016992319315487579497656388673505433776336703189690278955717239710851820593750575802469723477164557613094688400569544664083906632328980489050231322562719313784223293757762537765674860550223296757517936309514311364098955067987810622990290044610233314787617884818714076827914250567437281973993353197482098813435409755235361075969985798441337126312593960627842039335713762009624169971072389098442582559772999481646704068278833015939333008141841109592846547157670464817674715015636259928604821071600800873091706767125393539557966484701666419162259684521882415687361386486795039924743902847684455653672009117814002731651676792734130082374726001233740571183340055316038466019397098875483608598573963531197795231533784669488383851034898599484418365666793002675472937237274092221819708169896011572109417029618685719634982095269652618039495655688706701247477025085922221279463320915913152801169322907305498072372506045260384015734591714087786003188152532841620806596086485229059207306028656858812580376023757671759911398575408116826175902909012761714255020785822061641547239324759162257787133638766091640926969579726988971447262832031549969840158881281750736965695337486321518286113798675255459838184121812112863705855160188444420663952500796910005915012793169847550253471470237546273447099368433714408055787441435236639758429439479732597655485321756562005631462458721857009810050979024000068037899966702594286779178167104283058516755533365765054701540652974134109847005430449255509769055026949291803646099149879515878433840836779489981856269066079459746553702885573001041905761796108832155227593049875211636581275432414690028646002788224021719021326741635129313266431353578862049903277695541296438178355751015037527905950478175463466493661217977735637802129256184758010305201810374991157129887897715190473625904450982408862004860880041855444963715403816816697913347488258213796231921246982456767604532357185085595619225419494972418775913218700576674855723280963208713608759050096965340865736024194082384061468544443620372024356238984943165014000830348667282990315062294958017713685112256010523333894116609975573376955381051903582439853973475874565611435725163035826136598817082772994463193067795074551604854663708173556403995799133023516315205958185369423580248522691815920365525381967586987223820592252750246790439201396825488793939257364115218262946083474860523784830924246085630585814123698134711292541643275834774769358550038707208687561681985220208081109264672682079000578528676995067870581728389481309632471022423883273871806201922953495909101020921632302048824742636087442222271981512160368757719560642803804335208325783581324630948528372742921249668136599699430864754216512734431140871256723937412438287901118364475021875556572692586266456274813336044845327973318471952277585776779170860292001768956960996018482974850927059488179484838649054397646499203792764009751267114995747305588264268827308358967751669508793420818356792715912995146476118789747260518386453554450611500509015784210629209568851076896828558589881187034131236514543108275118274492680820722023356279609283225501210503057769959987375674460963328162312089681988114274578205642753188886259905616611476198660210570725766467627170954257896651227694959336299713781484347595473261014179723849074993539787736555059799178437730200032338748364417885122128183167598798597900535307364189904009133503543999817322306032952189385928613024146079985906321984855056816088225874686809752502871137380525052265954021704331783968465977328015884627311680559023242904016721675131323566011657154289720702412494820964921312304644491537963385820077869566295857565120487324475146126593516710180032457589644029021496e-01 Size: 2.395769e-15051
/*
console program in c language
based on the code by Claude Heiland-Allen
http://code.mathr.co.uk/mandelbrot-numerics
it computes center ( nucleus) of hyperbolic componnet of Mandelbrot set
using Newton method and double precision
https://wikibook.tw/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/centers
to compile
gcc m.c -Wall -std=c99 -lm
to run
./a.out
-------------
Newton method converged after 6 steps and near c= ( 0.000000 ;1.000000 ) found nucleus = center = ( -0.1225611668766536, 0.7448617666197442 ) for period = 3
Newton method converged after 5 steps and near c= ( 0.250000 ;0.500000 ) found nucleus = center = ( 0.2822713907669139, 0.5300606175785253 ) for period = 4
Newton method converged after 5 steps and near c= ( 0.250000 ;-0.500000 ) found nucleus = center = ( 0.2822713907669139, -0.5300606175785253 ) for period = 4
Newton method converged after 8 steps and near c= ( 0.300000 ;0.300000 ) found nucleus = center = ( 0.3795135880159238, 0.3349323055974976 ) for period = 5
Newton method converged after 5 steps and near c= ( -0.122561 ;0.844862 ) found nucleus = center = ( -0.1134186559494366, 0.8605694725015731 ) for period = 6
Newton method converged after 7 steps and near c= ( -1.250000 ;0.250000 ) found nucleus = center = ( -1.1380006666509646, 0.2403324012620983 ) for period = 6
Newton method converged after 35 steps and near c= ( -3.250000 ;0.250000 ) found nucleus = center = ( -1.9667732163929286, -0.0000000000000000 ) for period = 6
compare
https://gitlab.com/adammajewski/cpp-mandelbrot-center-by-knighty
https://gitlab.com/adammajewski/lawrence
https://gitlab.com/adammajewski/mandelbrot-numerics-nucleus
*/
#include <complex.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h> //
// from mandelbrot-numerics.h
enum m_newton { m_failed, m_stepped, m_converged };
typedef enum m_newton m_newton;
// from m_d_util.h
static inline bool cisfinite(double _Complex z) {
return isfinite(creal(z)) && isfinite(cimag(z));
}
static inline double cabs2(double _Complex z) {
return creal(z) * creal(z) + cimag(z) * cimag(z);
}
// last . takeWhile (\x -> 2 /= 2 + x) . iterate (/2) $ 1 :: Double
static const double epsilon = 4.440892098500626e-16;
// epsilon^2
static const double epsilon2 = 1.9721522630525295e-31;
int PrintResult(int iResult, double _Complex c_guess, double _Complex nucleus, int period, int steps)
{
static char* sResult;
switch (iResult){
case 0 : sResult = "failed"; break;
case 1 : sResult = "stepped"; break;
case 2 : sResult = "converged"; break;
default : sResult = "unknown";
}
printf("Newton method %s after %d steps and near c= ( %f ;%f ) found nucleus = center = ( %.16f, %.16f ) for period = %d \n ", sResult, steps, creal(c_guess), cimag(c_guess), creal(nucleus), cimag(nucleus), period );
return 0;
}
// ----------------------- from m_d_nucleus.c -------------------------------------------------
// nucleus = center of hyperbolic componnet of Mandelbrot set
// https://wikibook.tw/wiki/Fractals/Mathematics/Newton_method#center
//
extern m_newton m_d_nucleus_step(double _Complex *c_out, double _Complex c_guess, int period) {
double _Complex z = 0;
double _Complex dc = 0;
for (int i = 0; i < period; ++i) {
dc = 2 * z * dc + 1;
z = z * z + c_guess;
}
if (cabs2(dc) <= epsilon2) {
*c_out = c_guess;
return m_converged;
}
double _Complex c_new = c_guess - z / dc;
double _Complex d = c_new - c_guess;
if (cabs2(d) <= epsilon2) {
*c_out = c_new;
return m_converged;
}
if (cisfinite(d)) {
*c_out = c_new;
return m_stepped;
} else {
*c_out = c_guess;
return m_failed;
}
}
// compute nucleus using double precision and Newton method
extern m_newton m_d_nucleus(double _Complex *c_out, double _Complex c_guess, int period, int maxsteps) {
m_newton result = m_failed;
double _Complex c = c_guess;
int i;
for (i = 0; i < maxsteps; ++i) {
if (m_stepped != (result = m_d_nucleus_step(&c, c, period))) {
break;
}
}
*c_out = c;
PrintResult(result,c_guess, c , period, i);
return result;
}
/* ==================== main ================================================================================================*/
int main() {
double _Complex c3, c4a, c4b, c5, c3c2, c2c3;
m_d_nucleus(&c3, 0 + I * 1, 3, 64);
m_d_nucleus(&c4a, 0.25 + 0.5 * I, 4, 64);
m_d_nucleus(&c4b, 0.25 - 0.5 * I, 4, 64);
m_d_nucleus(&c5, 0.3 + 0.3 * I, 5, 64);
m_d_nucleus(&c3c2, c3 + I * 0.1, 6, 64);
m_d_nucleus(&c2c3, -1 - 0.25 + 0.25 * I, 6, 64);
m_d_nucleus(&c2c3, -3 - 0.25 + 0.25 * I, 6, 64);
return 0;
}
繪製填充有實心顏色的尾跡，使用 gmp 庫中的 mpq 型別
m-render-wakes > m-render-wakes.ppm
- 尾跡
-
週期 1 大陸附近的尾跡 -
週期 3 島附近的尾跡 -
沿主觸角的尾跡
-
逃逸路徑 1/2
cd ~/mandelbrot-numerics/c/bin
./m-interior
./m-interior: error while loading shared libraries: libmandelbrot-numerics.so: cannot open shared object file: No such file or directory
export LD_LIBRARY_PATH=${HOME}/opt/lib
export PATH=${HOME}/opt/bin:${PATH}
./m-interior
usage: ./m-interior precision z-guess-re z-guess-im c-guess-re c-guess-im interior-r interior-t period maxsteps
- ↑ askubuntu 問題：包未在 pkg-config 搜尋路徑中找到
- ↑ stackoverflow 問題：linux 載入共享庫時出錯 - 無法開啟共享物件檔案 - 沒有
- ↑ stackoverflow 問題：如何在 Linux 上永久設定路徑
- ↑ math.stackexchange 問題：曼德勃羅集中的完美圓圈
- ↑ Claude Heiland-Allen 的 atom_domain_size_estimation
- ↑ Claude Heiland-Allen 的 deriving_the_size_estimate
- ↑ Claude Heiland-Allen 的 on_the_precision_required_for_size_estimates
- ↑ fractalforums.com：deepest-e10000
- ↑ math.stackexchange 問題：測試曼德勃羅集週期為 n 的燈泡的成員資格/1151953#1151953
- ↑ math.stackexchange 問題：我渲染曼德勃羅集的方法叫什麼
- ↑ 曼德勃羅集中的內部座標
- ↑ 實用的內部距離渲染
- ↑ mathoverflow 問題：計算曼德勃羅集外部角的演算法
- ↑ fractalforums.org：re-feigenbaum-stretch