分形/曼德尔布罗特数值计算
< 分形
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 [email protected]: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://wikibooks.cn/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://wikibooks.cn/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