分形/数学/数值
< 分形
It will give you a clear visual perception of the gap between mathematics ... (and) your computer calculations. Lynn Wienck[1]
In the terminology of numerical analysis, computation near an attractive fixed point is stable while computation near a repelling fixed point is unstable. Generally speaking, stable computation is trustworthy while unstable computation is not. Mark McClure[2]
- 问题解对于非零但小的扰动（通常表示为 epsilon）是如何改变的
- 问题不能显式求解（使用符号方法），只能找到数值近似解
"... 收敛加速方法。假设你有一个收敛缓慢的级数，并且想了解其和（数值上）。仅仅通过对 x_1 + x_2 + ... + x_{1000} + ... x_{1000000} + ... + x_{1000000000} + ... 求和，你需要 100 年才能达到所需的精度。如果你将你的 x_n 拟合到 c_2/n^2 + c_3/n^3 + （一些更多项），你将在 1 秒内达到相同的精度。" 安德烈（理论高能物理学家）[3]
你可以在 抛物线朱莉娅集 中找到上述问题。
"你的序列被生成... 作为一个不动点迭代的输出。... 利文变换[4] ... 使用前向差分依次去除交替序列的误差项。"[5]
例子
- GSL 库[6][7]
- "检查级数逼近何时开始偏离将其运行到 Gnu 科学库中可用的利文 u 变换的结果。到目前为止，在我的测试中，这让我一直到达最终成为最低迭代的事件视界，即使在以前会导致其他方法更早地放弃级数逼近的位置，但它似乎也没有走得太远，以至于导致故障的伪影。它似乎刚刚好，全面而言，我发现这相当令人兴奋！ "（关于曼德博集的 quaz0r）[8]
- 级数计算器 [9]
- 当 2 个像素之间的距离小于：
- 当像素的大小下降到低于约 [10]
> 1. + 2 ** -53 == 1.
True
> 10. ** -324 == 0.
True
数值计算中的精度[11]
如何处理浮点精度不足的问题：[12]
- 实现比硬件原生支持更高的精度算术
- 软件模拟（使用两个浮点数模拟双精度数，定点数，...）[13]
- 使用数值上更稳定的算法
- 使用可表示的数字
区间 [0,1] 中有多少浮点数？[14]
库
- Fredrik Johansson 的 Arb[15]
查看如何使用 arb 库来检查和调整精度
/*
from arb/examples/logistic.c
public domain.
Author: Fredrik Johansson.
*/
goal = digits * 3.3219280948873623 + 3;
for (prec = 64; ; prec *= 2)
{
flint_printf("Trying prec=%wd bits...", prec);
fflush(stdout);
for (i = 0; i < n; i++)
{
// computation
if (arb_rel_accuracy_bits(x) < goal)
{
flint_printf("ran out of accuracy at step %wd\n", i);
break;
}
}
if (i == n)
{
flint_printf("success!\n");
break;
}
}
二进制数的一位需要一位：有 2 个二进制数字（0 和 1），位有 2 个状态。十进制数的一位大约需要 3.4 位。[16] 有 10 个十进制数，3 位有 8 个状态，这不够，4 位有 16 个状态，这太多了。
| 类型 | 总位数 | 精度位数 | 有效十进制数位数 |
|---|---|---|---|
| float | 32 | 24 | 6 |
| double | 64 | 53 | 15 |
| long double | 80 | 64 | 18[17] |
| mpfr_t[18] | 142 | 41 | |
| mpfr_t | 1000 | 294 |
- "计算时保留比最终答案多两位有效数字
- 只在计算的最后一步进行舍入。永远不要用舍入后的数字进行进一步的计算。
- 对于乘法和除法，找到每个因子的有效数字位数。结果将具有较少的有效数字位数。
- 对于幂和根，答案应该具有与原始数字相同的有效数字位数。" [19]
- "当你加（或减）时，保留与精度最低的数字相同的小数位数。
- 当您进行乘法（或除法）运算时，保留的 **有效数字** 数量应与精度最低的数字相同。[20]
参见
**沃尔夫·容格测试**：角度（以圈为单位）的外部参数射线
- 1/7（周期为 3）
- 321685687669320/2251799813685247（周期为 51）
- 321685687669322/2251799813685247（周期为 51）
角度大约相差 ，但对应参数射线的着陆点大约相差 0.035。[23]
对于旋转数（内部角度）为 1/34 的射线，外部角度为
落在阿尔法不动点上
它不是角度为
的射线，该射线落在以下点上：[24]
射线外部角度之差为
而射线着陆点之间的距离为
/* Maxima CAS code */ (%i1) za:0.491486549841951 +0.091874758908285*%i; (%o1) 0.091874758908285 %i + 0.491486549841951 (%i2) zb:-0.508513450158049 +0.091874758908285*%i; (%o2) 0.091874758908285 %i - 0.508513450158049 (%i3) abs(za-zb); (%o3) 1.0
Storing the angle as an exact rational (e.g. mpq_t using libgmp), allows it to be accurate for more than 53 iterations (the number of bits of accuracy of double precision floating point). Convert to double only when needing to do sin and cos for finding the target point. Double is plenty accurate for target point. ( Claude Heiland-Allen )
该测试由约翰·米尔诺提出。[25] 另请参见马克·布雷弗曼的分析 [26] 以及罗伯特·P·穆纳福的舍入误差分析[27]
马克·麦克卢尔的评论：[28] “逃逸时间算法将花费很长时间才能生成这种类型的图像，因为动态变化在那里非常缓慢。如果您想要 1/100 的分辨率，那么大约需要 2*10^8 次迭代才能通过迭代 f(z)=z+z5 将点 z0=0.01 移动到 z=2。”
这里有一个 由 gerrit 使用 1000 位小数计算的轨道，在经过大约 300000 次迭代后逃逸
让我们来看一个简单的双曲情况，其中参数 c 为
这里排斥的不动点 z_f 为
让我们来看一个简单的 抛物线情况，其中参数 c 为：[29]
这里抛物线不动点 z_f 为
让我们取 Julia 集外部但靠近不动点的点 z
其中 n 是一个正整数
检查点 z 需要多少次迭代 i 才能到达目标集（= 逃逸）
其中 ER 是 逃逸半径
显示之间的关系
- n
- 最后一次迭代
- 用于计算的数字类型（浮点数、双精度浮点数、长双精度浮点数、扩展精度浮点数、任意精度浮点数）
参见 FractalForum 获取 evaldraw 脚本[30]
| CPU 单线程 C 源代码 - 双曲线情况，浮点数。点击右侧查看 |
|---|
/*
gcc -lm -Wall p.c
./a.out
*/
#include <stdio.h>
#include <math.h> /* pow () */
#include <time.h>
float GiveLastIteration (float Zx, float Zy, float Cx, float Cy, int ER2)
{
float Zx2, Zy2; /* Zx2 = Zx * Zx; Zy2 = Zy * Zy */
float i = 0.0;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
while (Zx2 + Zy2 < ER2) /* ER2 = ER * ER */
{
Zy = 2 * Zx * Zy + Cy;
Zx = Zx2 - Zy2 + Cx;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
i += 1.0;
}
return i;
}
int main ()
{
float distance;
int n;
/* fixed point zf = zfx = zfy*I = 1.0 for c=0 */
float zfx = 1.0;
float zfy = 0.0;
/* z1 = z1x + z1y * I = point of exterior of Julia set but near parabolic fixed point zp */
float z1x;
float z1y = zfy;
float cx= 0.0;
float cy= 0.0;
/* Escape Radius ; it defines target set = { z: abs(z) > ER}
all points z in the target set are escaping to infinity */
float ER = 2.0;
float ER2;
float LastIteration;
time_t start, end;
float dif;
ER2 = ER * ER;
printf ("Using c with float and Escape Radius = %f \n", ER);
for (n = 1; n < 101; n++)
{
time (&start);
distance = pow (2.0, -n);
printf ("n= %3d distance = %e = %.20f ", n, distance, distance);
z1x = zfx + distance;
LastIteration = GiveLastIteration (z1x, z1y, cx, cy, ER2);
time (&end);
dif = difftime (end, start);
printf ("LI = %3.0f time = %2.0lf seconds\n", LastIteration, dif);
}
return 0;
}
|
| CPU 单线程 C 源代码 - 双曲线情况，双精度浮点数。点击右侧查看 |
|---|
/*
gcc -lm -Wall p.c
./a.out
*/
#include <stdio.h>
#include <math.h> /* pow () */
#include <time.h>
double GiveLastIteration(double Zx, double Zy, double Cx, double Cy, int ER2)
{
double Zx2, Zy2; /* Zx2 = Zx * Zx; Zy2 = Zy * Zy */
double i = 0.0;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
while (Zx2 + Zy2 < ER2) /* ER2 = ER * ER */
{
Zy = 2 * Zx * Zy + Cy;
Zx = Zx2 - Zy2 + Cx;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
i += 1.0;
}
return i;
}
int main ()
{
double distance;
int n;
// fixed point zf = zfx = zfy*I = 1.0 for c=0
double zfx = 1.0;
double zfy = 0.0;
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
double z1x;
double z1y = zfy;
double cx= 0.0;
double cy= 0.0;
// Escape Radius ; it defines target set = { z: abs(z)>ER}
// all points z in the target set are escaping to infinity
double ER = 2.0;
double ER2;
double LastIteration;
time_t start,end;
double dif;
ER2 = ER * ER;
printf ("Using c with double and Escape Radius = %f \n", ER);
for (n = 1; n < 101; n++)
{
time (&start);
distance = pow (2.0, -n);
printf("n= %3d distance = %e = %.20f ", n, distance, distance);
z1x = zfx + distance;
LastIteration = GiveLastIteration (z1x, z1y, cx, cy, ER2);
time (&end);
dif = difftime (end, start);
printf("LI = %3.0f time = %2.0lf seconds\n", LastIteration, dif);
}
return 0;
}
|
| CPU 单线程 C 源代码 - 双曲线情况，长双精度浮点数。点击右侧查看 |
|---|
/*
gcc -lm -Wall p.c
./a.out
*/
#include <stdio.h>
#include <math.h> /* pow () */
#include <time.h>
long double GiveLastIteration (long double Zx0, long double Zy0, long double Cx, long double Cy, int ER2)
{
long double Zx2, Zy2; /* Zx2 = Zx * Zx; Zy2 = Zy * Zy */
long double Zx = Zx0;
long double Zy = Zy0;
long double i = 0.0;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
while (Zx2 + Zy2 < ER2) /* ER2 = ER * ER */
{
Zy = 2 * Zx * Zy + Cy;
Zx = Zx2 - Zy2 + Cx;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
i += 1.0;
}
return i;
}
int main ()
{
long double distance;
int n;
/* fixed point zp = zpx = zpy * I = 1.0 for c = 0 */
long double zpx = 1.0;
long double zpy = 0.0;
/* z1 = z1x + z1y * I = point of exterior of Julia set but near parabolic fixed point zp */
long double z1x;
long double z1y = zpy;
long double cx = 0.0;
long double cy = 0.0;
/* Escape Radius ; it defines target set = { z: abs(z) > ER}
all points z in the target set are escaping to infinity */
long double ER = 2.0;
long double ER2;
long double LastIteration;
time_t start,end;
long double dif;
ER2 = ER * ER;
printf ("Using c with long double and Escape Radius = %Lf \n", ER);
for (n = 1; n < 101; n++)
{
time (&start);
distance = pow (2.0, -n);
z1x = zpx + distance;
LastIteration = GiveLastIteration (z1x, z1y, cx, cy, ER2);
time (&end);
dif = difftime (end, start);
printf("n= %3d distance = %Le = %.10Lf LI = %10.0Lf log2(LI) = %3.0Lf time = %2.0Lf seconds\a\n", n, distance, distance, LastIteration, log2l (LastIteration), dif);
}
return 0;
}
|
| CPU 单线程 C 源代码 - 抛物线情况，浮点数。点击右侧查看 |
|---|
/*
gcc -lm -Wall p.c
./a.out
*/
#include <stdio.h>
#include <math.h> /* pow () */
#include <time.h>
float GiveLastIteration (float Zx0, float Zy0, float Cx, float Cy, int ER2)
{
float Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */
float Zx = Zx0;
float Zy = Zy0;
float i = 0.0;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
while (Zx2 + Zy2 < ER2) /* ER2 = ER * ER */
{
Zy = 2 * Zx * Zy + Cy;
Zx = Zx2 - Zy2 + Cx;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
i += 1.0;
}
return i;
}
int main ()
{
float distance;
int n;
// parabolic fixed point zp = zpx = zpy*I = 0.5 for c=1/4
float zpx = 0.5;
float zpy = 0.0;
// z1 = z1x + z1y*I = point of exterior of Julia set but near parabolic fixed point zp
float z1x;
float z1y = zpy;
float cx = 0.25;
float cy = 0.0;
/* Escape Radius ; it defines target set = { z: abs(z) > ER}
all points z in the target set are escaping to infinity */
float ER = 2.0;
float ER2;
float LastIteration;
time_t start, end;
float dif;
ER2= ER*ER;
printf ("Using c with float and Escape Radius = %f \n", ER);
for (n = 1; n < 101; n++)
{
time (&start);
distance = pow (2.0, -n);
z1x = zpx + distance;
LastIteration = GiveLastIteration (z1x, z1y, cx, cy, ER2);
time (&end);
dif = difftime (end, start);
printf("n= %3d distance = %e = %.10f LI = %10.0f log2(LI) = %3.0f time = %2.0lf seconds\n", n, distance, distance, LastIteration, log2 (LastIteration), dif);
}
return 0;
}
|
| CPU 单线程 C 源代码 - 抛物线情况，双精度浮点数。点击右侧查看 |
|---|
/*
* gcc -lm -Wall p.c
* ./a.out
*/
#include <stdio.h>
#include <math.h> /* pow () */
#include <time.h>
double
GiveLastIteration(double Zx0, double Zy0, double Cx, double Cy, int ER2)
{
double Zx2 , Zy2; /* Zx2 = Zx * Zx; Zy2 = Zy * Zy */
double Zx = Zx0;
double Zy = Zy0;
double i = 0.0;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
while (Zx2 + Zy2 < ER2) { /* ER2 = ER * ER */
Zy = 2 * Zx * Zy + Cy;
Zx = Zx2 - Zy2 + Cx;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
i += 1.0;
}
return i;
}
double
max(double a, double b)
{
if (a > b)
return a;
else
return b;
}
int
main()
{
double distance;
int n;
//parabolic fixed point zp = zpx = zpy * I = 0.5 for c = 1 / 4
double zpx = 0.5;
double zpy = 0.0;
//z1 = z1x + z1y * I = point of exterior of Julia set but near parabolic fixed point zp
double z1x;
double z1y = zpy;
double cx = 0.25;
double cy = 0.0;
//Escape Radius; it defines target set = {z:abs(z) > ER}
//all points z in the target set are escaping to infinity
double ER = 2.0;
double ER2;
double LastIteration;
time_t start , end;
double dif;
ER2 = ER * ER;
for (n = 1; n < 101; n++) {
time(&start);
distance = pow(2.0, -n);
z1x = zpx + distance;
LastIteration = GiveLastIteration(z1x, z1y, cx, cy, ER2);
time(&end);
dif = difftime(end, start);
printf("n= %3d distance = %e = %.10f LI = %10.0f log2(LI) = %3.0f time = %2.0lf seconds\n",
n, distance, distance, LastIteration, log2(LastIteration), dif);
}
return 0;
}
|
| CPU 单线程 C 源代码 - 抛物线情况，长双精度浮点数。点击右侧查看 |
|---|
/*
* gcc -lm -Wall p.c
* ./a.out
*/
#include <stdio.h>
#include <math.h> /* pow () */
#include <time.h>
long double
GiveLastIteration(long double Zx0, long double Zy0, long double Cx, long double Cy, int ER2)
{
long double Zx2, Zy2; /* Zx2 = Zx * Zx; Zy2 = Zy * Zy */
long double Zx = Zx0;
long double Zy = Zy0;
long double i = 0.0;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
while (Zx2 + Zy2 < ER2) { /* ER2 = ER * ER */
Zy = 2 * Zx * Zy + Cy;
Zx = Zx2 - Zy2 + Cx;
Zx2 = Zx * Zx;
Zy2 = Zy * Zy;
i += 1.0;
}
return i;
}
int
main()
{
long double distance;
int n;
//parabolic fixed point zp = zpx = zpy * I = 0.5 for c = 1 / 4
long double zpx = 0.5;
long double zpy = 0.0;
//z1 = z1x + z1y * I = point of exterior of Julia set but near parabolic fixed point zp
long double z1x;
long double z1y = zpy;
long double cx = 0.25;
long double cy = 0.0;
//Escape Radius; it defines target set = {z:abs(z) > ER}
//all points z in the target set are escaping to infinity
long double ER = 2.0;
long double ER2;
long double LastIteration;
time_t start , end;
long double dif;
ER2 = ER * ER;
for (n = 1; n < 101; n++) {
time(&start);
distance = pow(2.0, -n);
z1x = zpx + distance;
LastIteration = GiveLastIteration(z1x, z1y, cx, cy, ER2);
time(&end);
dif = difftime(end, start);
printf("n= %3d distance = %Le = %.10Lf LI = %10.0Lf log2(LI) = %3.0Lf time = %2.0Lf seconds\a\n", n, distance, distance, LastIteration, log2l(LastIteration), dif);
}
return 0;
}
|
| CPU 单线程 C 源代码 - 抛物线情况，GSL 库。点击右侧查看 |
|---|
/*
* http://www.gnu.org/software/gsl/
*
* gcc -L/usr/lib -lgsl -lgslcblas -lm -Wall pgsl.c
* ./a.out
*
* or
* GSL_IEEE_MODE="extended-precision" ./a.out
*
*/
#include <stdio.h>
#include <math.h> /* pow () */
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <time.h>
long double
GiveLastIteration(gsl_complex z, gsl_complex c, int ER2)
{
long double i = 0.0;
while (gsl_complex_abs2(z) < ER2) { /* ER2 = ER * ER */
z = gsl_complex_add(c, gsl_complex_mul(z, z));
i += 1.0;
}
return i;
}
int
main()
{
//long double distance;
int n;
//parabolic fixed point zp = zpx = zpy * I = 0.5 for c = 1 / 4
gsl_complex zp = gsl_complex_rect(0.5, 0.0);
//long double zpx = 0.5;
//long double zpy = 0.0;
//z1 = z1x + z1y * I = point of exterior of Julia set but near parabolic fixed point zp
gsl_complex z1;
//long double z1x;
//long double z1y = zpy;
//
gsl_complex c = gsl_complex_rect(0.25, 0.0);
//long double cx = 0.25;
//long double cy = 0.0;
//Escape Radius; it defines target set = {z:abs(z) > ER}
//all points z in the target set are escaping to infinity
long double ER = 2.0;
long double ER2;
long double LastIteration;
time_t start , end;
long double dif;
ER2 = ER * ER;
for (n = 1; n < 101; n++) {
time(&start);
z1 = gsl_complex_rect(GSL_REAL(zp) + pow(2.0, -n), 0.0);
LastIteration = GiveLastIteration(z1, c, ER2);
time(&end);
dif = difftime(end, start);
printf("n= %3d LI = %10.0Lf log2(LI) = %3.0Lf time = %2.0Lf seconds\a\n", n, LastIteration, log2l(LastIteration), dif);
}
return 0;
}
|
| CPU 单线程 C 源代码 - 抛物线情况，mpfr 库。点击右侧查看 |
|---|
/*
* gcc -lm -lmpfr -lgmp -Wall -O2 pmpfr.c ./a.out
*
*/
#include <stdio.h>
#include <math.h> // log2
#include <gmp.h>
#include <mpfr.h>
#include <time.h>
// MPFR general settings
mpfr_prec_t precision = 12;
//the number of bits used to represent the significand of a floating - point number
mpfr_rnd_t RoundMode = MPFR_RNDD;
//the way to round the result of a floating - point operation, round toward minus infinity(roundTowardNegative in IEEE 754 - 2008),
double
GiveLastIteration_mpfr (mpfr_t zx,
mpfr_t zy, mpfr_t cx, mpfr_t cy, mpfr_t ER2)
{
double i = 0.0;
mpfr_t temp;
//temporary variable
mpfr_t zx2;
//zx ^ 2
mpfr_t zy2;
//zy ^ 2
// initializes MPFR variables with default_prec
mpfr_inits (temp, zx2, zy2, (mpfr_ptr) 0);
mpfr_mul (zy2, zy, zy, RoundMode); /* zy2 = zy * zy; */
mpfr_mul (zx2, zx, zx, RoundMode); /* zx2 = zx * zx; */
mpfr_add (temp, zx2, zy2, RoundMode);
//temp = zx2 + zy2
while (mpfr_greater_p (ER2, temp))
//while (Zx2 + Zy2 < ER2
{
/* zy = 2.0*zx*zy + cy; */
mpfr_mul (temp, zx, zy, RoundMode);
//temp = zx * zy
mpfr_mul_ui (temp, temp, 2.0, RoundMode);
//temp = 2 * temp = 2 * zx * zy
mpfr_add (zy, temp, cy, RoundMode);
/* zx = zx^2 - zy^2 + cx; */
mpfr_sub (temp, zx2, zy2, RoundMode);
//temp = zx ^ 2 - zy ^ 2
mpfr_add (zx, temp, cx, RoundMode);
//temp = temp + cx
mpfr_mul (zy2, zy, zy, RoundMode); /* zy2 = zy * zy; */
mpfr_mul (zx2, zx, zx, RoundMode); /* zx2 = zx * zx; */
mpfr_add (temp, zx2, zy2, RoundMode);
//temp = zx2 + zy2
i += 1.0;
}
mpfr_clears (temp, zx2, zy2, (mpfr_ptr) 0);
return i;
}
int
main (void)
{
//declares variables
double LastIteration;
time_t start, end;
double dif;
int n;
mpfr_t zpx, zpy;
//parabolic fixed point zp = zpx + zpy * I = 0.5 for c = 1 / 4 mpfr_t zx, zy;
//point z = zx + zy * I
mpfr_t distance;
//between zp and z
mpfr_t cx, cy;
//point c = cx + cy *
mpfr_t base;
mpfr_t ER;
/*
* Escape Radius; it defines target set = {z:abs(z) > ER}; all points
* z in the target set are escaping to infinity
*/
mpfr_t ER2;
/* ER2 = ER * ER */
mpfr_set_default_prec (precision);
//initializes MPFR variables with default_prec
mpfr_inits (ER, ER2, zpx, zpy, zx, zy, distance, cx, cy, base,
(mpfr_ptr) 0);
//Assignment Functions
mpfr_set_ld (zpx, 0.5, RoundMode);
mpfr_set_ld (zpy, 0.0, RoundMode);
mpfr_set_ld (zy, 0.0, RoundMode);
mpfr_set_ld (cx, 0.25, RoundMode);
mpfr_set_ld (cy, 0.0, RoundMode);
mpfr_set_ld (base, 2.0, RoundMode);
mpfr_set_ld (ER, 2.0, RoundMode);
//computations
mpfr_mul (ER2, ER, ER, RoundMode);
mpfr_printf
("Using MPFR-%s with GMP-%s with precision = %u bits and Escape Radius = %Rf \n",
mpfr_version, gmp_version, (unsigned int) precision, ER);
for (n = 1; n <= 60; n++)
{
time (&start);
mpfr_pow_si (distance, base, -n, RoundMode);
mpfr_add (zx, zpx, distance, RoundMode);
LastIteration = GiveLastIteration_mpfr (zx, zy, cx, cy, ER2);
time (&end);
dif = difftime (end, start);
mpfr_printf
("n = %3d distance = %.10Re LI = %10.0f log2(LI) = %3.0f; time = %2.0f seconds \n",
n, distance, LastIteration, log2 (LastIteration), dif);
}
//free the space used by the MPFR variables
mpfr_clears (ER, ER2, zpx, zpy, zx, zy, distance, cx, cy, base,
(mpfr_ptr) 0);
mpfr_free_cache (); /* free the cache for constants like pi */
return 0;
}
|
| float (24) | double (53) | long double (64) | MPFR (80) | MPFR (100) | |
|---|---|---|---|---|---|
| 双曲线 | 23 | 52 | 63 | 79 | 99 |
| 抛物线 | 12 | 26 | 32 |
标准 C 类型（float、double、long double）和 MPFR 精度在相同精度下结果相同
| 1 | 2 | 3 | 4 | 5 | 24 | 53 | 64 | 80 | 100 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 双曲线 | 1 | 2 | 3 | 4 | 5 | 24 | 53 | 64 | 80 | 100 |
| 抛物线 | 3 | 5 | 10 | 19 | 35 | 16 778 821 |
双曲线情况下迭代次数和计算时间之间的关系
Using MPFR-3.0.0-p8 with GMP-4.3.2 with precision = 128 bits and Escape Radius = 2.000000 n = 1 distance = 5.0000000000e-01 LI = 1 log2(LI) = 0; time = 0 seconds n = 2 distance = 2.5000000000e-01 LI = 2 log2(LI) = 1; time = 0 seconds n = 3 distance = 1.2500000000e-01 LI = 3 log2(LI) = 2; time = 0 seconds n = 4 distance = 6.2500000000e-02 LI = 4 log2(LI) = 2; time = 0 seconds n = 5 distance = 3.1250000000e-02 LI = 5 log2(LI) = 2; time = 0 seconds n = 6 distance = 1.5625000000e-02 LI = 6 log2(LI) = 3; time = 0 seconds n = 7 distance = 7.8125000000e-03 LI = 7 log2(LI) = 3; time = 0 seconds n = 8 distance = 3.9062500000e-03 LI = 8 log2(LI) = 3; time = 0 seconds n = 9 distance = 1.9531250000e-03 LI = 9 log2(LI) = 3; time = 0 seconds n = 10 distance = 9.7656250000e-04 LI = 10 log2(LI) = 3; time = 0 seconds n = 11 distance = 4.8828125000e-04 LI = 11 log2(LI) = 3; time = 0 seconds n = 12 distance = 2.4414062500e-04 LI = 12 log2(LI) = 4; time = 0 seconds n = 13 distance = 1.2207031250e-04 LI = 13 log2(LI) = 4; time = 0 seconds n = 14 distance = 6.1035156250e-05 LI = 14 log2(LI) = 4; time = 0 seconds n = 15 distance = 3.0517578125e-05 LI = 15 log2(LI) = 4; time = 0 seconds n = 16 distance = 1.5258789062e-05 LI = 16 log2(LI) = 4; time = 0 seconds n = 17 distance = 7.6293945312e-06 LI = 17 log2(LI) = 4; time = 0 seconds n = 18 distance = 3.8146972656e-06 LI = 18 log2(LI) = 4; time = 0 seconds n = 19 distance = 1.9073486328e-06 LI = 19 log2(LI) = 4; time = 0 seconds n = 20 distance = 9.5367431641e-07 LI = 20 log2(LI) = 4; time = 0 seconds n = 21 distance = 4.7683715820e-07 LI = 21 log2(LI) = 4; time = 0 seconds n = 22 distance = 2.3841857910e-07 LI = 22 log2(LI) = 4; time = 0 seconds n = 23 distance = 1.1920928955e-07 LI = 23 log2(LI) = 5; time = 0 seconds n = 24 distance = 5.9604644775e-08 LI = 24 log2(LI) = 5; time = 0 seconds n = 25 distance = 2.9802322388e-08 LI = 25 log2(LI) = 5; time = 0 seconds n = 26 distance = 1.4901161194e-08 LI = 26 log2(LI) = 5; time = 0 seconds n = 27 distance = 7.4505805969e-09 LI = 27 log2(LI) = 5; time = 0 seconds n = 28 distance = 3.7252902985e-09 LI = 28 log2(LI) = 5; time = 0 seconds n = 29 distance = 1.8626451492e-09 LI = 29 log2(LI) = 5; time = 0 seconds n = 30 distance = 9.3132257462e-10 LI = 30 log2(LI) = 5; time = 0 seconds n = 31 distance = 4.6566128731e-10 LI = 31 log2(LI) = 5; time = 0 seconds n = 32 distance = 2.3283064365e-10 LI = 32 log2(LI) = 5; time = 0 seconds n = 33 distance = 1.1641532183e-10 LI = 33 log2(LI) = 5; time = 0 seconds n = 34 distance = 5.8207660913e-11 LI = 34 log2(LI) = 5; time = 0 seconds n = 35 distance = 2.9103830457e-11 LI = 35 log2(LI) = 5; time = 0 seconds n = 36 distance = 1.4551915228e-11 LI = 36 log2(LI) = 5; time = 0 seconds n = 37 distance = 7.2759576142e-12 LI = 37 log2(LI) = 5; time = 0 seconds n = 38 distance = 3.6379788071e-12 LI = 38 log2(LI) = 5; time = 0 seconds n = 39 distance = 1.8189894035e-12 LI = 39 log2(LI) = 5; time = 0 seconds n = 40 distance = 9.0949470177e-13 LI = 40 log2(LI) = 5; time = 0 seconds n = 41 distance = 4.5474735089e-13 LI = 41 log2(LI) = 5; time = 0 seconds n = 42 distance = 2.2737367544e-13 LI = 42 log2(LI) = 5; time = 0 seconds n = 43 distance = 1.1368683772e-13 LI = 43 log2(LI) = 5; time = 0 seconds n = 44 distance = 5.6843418861e-14 LI = 44 log2(LI) = 5; time = 0 seconds n = 45 distance = 2.8421709430e-14 LI = 45 log2(LI) = 5; time = 0 seconds n = 46 distance = 1.4210854715e-14 LI = 46 log2(LI) = 6; time = 0 seconds n = 47 distance = 7.1054273576e-15 LI = 47 log2(LI) = 6; time = 0 seconds n = 48 distance = 3.5527136788e-15 LI = 48 log2(LI) = 6; time = 0 seconds n = 49 distance = 1.7763568394e-15 LI = 49 log2(LI) = 6; time = 0 seconds n = 50 distance = 8.8817841970e-16 LI = 50 log2(LI) = 6; time = 0 seconds n = 51 distance = 4.4408920985e-16 LI = 51 log2(LI) = 6; time = 0 seconds n = 52 distance = 2.2204460493e-16 LI = 52 log2(LI) = 6; time = 0 seconds n = 53 distance = 1.1102230246e-16 LI = 53 log2(LI) = 6; time = 0 seconds n = 54 distance = 5.5511151231e-17 LI = 54 log2(LI) = 6; time = 0 seconds n = 55 distance = 2.7755575616e-17 LI = 55 log2(LI) = 6; time = 0 seconds n = 56 distance = 1.3877787808e-17 LI = 56 log2(LI) = 6; time = 0 seconds n = 57 distance = 6.9388939039e-18 LI = 57 log2(LI) = 6; time = 0 seconds n = 58 distance = 3.4694469520e-18 LI = 58 log2(LI) = 6; time = 0 seconds n = 59 distance = 1.7347234760e-18 LI = 59 log2(LI) = 6; time = 0 seconds n = 60 distance = 8.6736173799e-19 LI = 60 log2(LI) = 6; time = 0 seconds n = 61 distance = 4.3368086899e-19 LI = 61 log2(LI) = 6; time = 0 seconds n = 62 distance = 2.1684043450e-19 LI = 62 log2(LI) = 6; time = 0 seconds n = 63 distance = 1.0842021725e-19 LI = 63 log2(LI) = 6; time = 0 seconds n = 64 distance = 5.4210108624e-20 LI = 64 log2(LI) = 6; time = 0 seconds n = 65 distance = 2.7105054312e-20 LI = 65 log2(LI) = 6; time = 0 seconds n = 66 distance = 1.3552527156e-20 LI = 66 log2(LI) = 6; time = 0 seconds n = 67 distance = 6.7762635780e-21 LI = 67 log2(LI) = 6; time = 0 seconds n = 68 distance = 3.3881317890e-21 LI = 68 log2(LI) = 6; time = 0 seconds n = 69 distance = 1.6940658945e-21 LI = 69 log2(LI) = 6; time = 0 seconds n = 70 distance = 8.4703294725e-22 LI = 70 log2(LI) = 6; time = 0 seconds n = 71 distance = 4.2351647363e-22 LI = 71 log2(LI) = 6; time = 0 seconds n = 72 distance = 2.1175823681e-22 LI = 72 log2(LI) = 6; time = 0 seconds n = 73 distance = 1.0587911841e-22 LI = 73 log2(LI) = 6; time = 0 seconds n = 74 distance = 5.2939559203e-23 LI = 74 log2(LI) = 6; time = 0 seconds n = 75 distance = 2.6469779602e-23 LI = 75 log2(LI) = 6; time = 0 seconds n = 76 distance = 1.3234889801e-23 LI = 76 log2(LI) = 6; time = 0 seconds n = 77 distance = 6.6174449004e-24 LI = 77 log2(LI) = 6; time = 0 seconds n = 78 distance = 3.3087224502e-24 LI = 78 log2(LI) = 6; time = 0 seconds n = 79 distance = 1.6543612251e-24 LI = 79 log2(LI) = 6; time = 0 seconds n = 80 distance = 8.2718061255e-25 LI = 80 log2(LI) = 6; time = 0 seconds n = 81 distance = 4.1359030628e-25 LI = 81 log2(LI) = 6; time = 0 seconds n = 82 distance = 2.0679515314e-25 LI = 82 log2(LI) = 6; time = 0 seconds n = 83 distance = 1.0339757657e-25 LI = 83 log2(LI) = 6; time = 0 seconds n = 84 distance = 5.1698788285e-26 LI = 84 log2(LI) = 6; time = 0 seconds n = 85 distance = 2.5849394142e-26 LI = 85 log2(LI) = 6; time = 0 seconds n = 86 distance = 1.2924697071e-26 LI = 86 log2(LI) = 6; time = 0 seconds n = 87 distance = 6.4623485356e-27 LI = 87 log2(LI) = 6; time = 0 seconds n = 88 distance = 3.2311742678e-27 LI = 88 log2(LI) = 6; time = 0 seconds n = 89 distance = 1.6155871339e-27 LI = 89 log2(LI) = 6; time = 0 seconds n = 90 distance = 8.0779356695e-28 LI = 90 log2(LI) = 6; time = 0 seconds n = 91 distance = 4.0389678347e-28 LI = 91 log2(LI) = 7; time = 0 seconds n = 92 distance = 2.0194839174e-28 LI = 92 log2(LI) = 7; time = 0 seconds n = 93 distance = 1.0097419587e-28 LI = 93 log2(LI) = 7; time = 0 seconds n = 94 distance = 5.0487097934e-29 LI = 94 log2(LI) = 7; time = 0 seconds n = 95 distance = 2.5243548967e-29 LI = 95 log2(LI) = 7; time = 0 seconds n = 96 distance = 1.2621774484e-29 LI = 96 log2(LI) = 7; time = 0 seconds n = 97 distance = 6.3108872418e-30 LI = 97 log2(LI) = 7; time = 0 seconds n = 98 distance = 3.1554436209e-30 LI = 98 log2(LI) = 7; time = 0 seconds n = 99 distance = 1.5777218104e-30 LI = 99 log2(LI) = 7; time = 0 seconds n = 100 distance = 7.8886090522e-31 LI = 100 log2(LI) = 7; time = 0 seconds n = 101 distance = 3.9443045261e-31 LI = 101 log2(LI) = 7; time = 0 seconds n = 102 distance = 1.9721522631e-31 LI = 102 log2(LI) = 7; time = 0 seconds n = 103 distance = 9.8607613153e-32 LI = 103 log2(LI) = 7; time = 0 seconds n = 104 distance = 4.9303806576e-32 LI = 104 log2(LI) = 7; time = 0 seconds n = 105 distance = 2.4651903288e-32 LI = 105 log2(LI) = 7; time = 0 seconds n = 106 distance = 1.2325951644e-32 LI = 106 log2(LI) = 7; time = 0 seconds n = 107 distance = 6.1629758220e-33 LI = 107 log2(LI) = 7; time = 0 seconds n = 108 distance = 3.0814879110e-33 LI = 108 log2(LI) = 7; time = 0 seconds n = 109 distance = 1.5407439555e-33 LI = 109 log2(LI) = 7; time = 0 seconds n = 110 distance = 7.7037197775e-34 LI = 110 log2(LI) = 7; time = 0 seconds n = 111 distance = 3.8518598888e-34 LI = 111 log2(LI) = 7; time = 0 seconds n = 112 distance = 1.9259299444e-34 LI = 112 log2(LI) = 7; time = 0 seconds n = 113 distance = 9.6296497219e-35 LI = 113 log2(LI) = 7; time = 0 seconds n = 114 distance = 4.8148248610e-35 LI = 114 log2(LI) = 7; time = 0 seconds n = 115 distance = 2.4074124305e-35 LI = 115 log2(LI) = 7; time = 0 seconds n = 116 distance = 1.2037062152e-35 LI = 116 log2(LI) = 7; time = 0 seconds n = 117 distance = 6.0185310762e-36 LI = 117 log2(LI) = 7; time = 0 seconds n = 118 distance = 3.0092655381e-36 LI = 118 log2(LI) = 7; time = 0 seconds n = 119 distance = 1.5046327691e-36 LI = 119 log2(LI) = 7; time = 0 seconds n = 120 distance = 7.5231638453e-37 LI = 120 log2(LI) = 7; time = 0 seconds n = 121 distance = 3.7615819226e-37 LI = 121 log2(LI) = 7; time = 0 seconds n = 122 distance = 1.8807909613e-37 LI = 122 log2(LI) = 7; time = 0 seconds n = 123 distance = 9.4039548066e-38 LI = 123 log2(LI) = 7; time = 0 seconds n = 124 distance = 4.7019774033e-38 LI = 124 log2(LI) = 7; time = 0 seconds n = 125 distance = 2.3509887016e-38 LI = 125 log2(LI) = 7; time = 0 seconds n = 126 distance = 1.1754943508e-38 LI = 126 log2(LI) = 7; time = 0 seconds n = 127 distance = 5.8774717541e-39 LI = 127 log2(LI) = 7; time = 0 seconds
抛物线情况
Using MPFR-3.0.0-p8 with GMP-4.3.2 with precision = 100 bits and Escape Radius = 2.000000 n = 1 distance = 5.0000000000e-01 LI = 3 log2(LI) = 2; time = 0 seconds n = 2 distance = 2.5000000000e-01 LI = 5 log2(LI) = 2; time = 0 seconds n = 3 distance = 1.2500000000e-01 LI = 10 log2(LI) = 3; time = 0 seconds n = 4 distance = 6.2500000000e-02 LI = 19 log2(LI) = 4; time = 0 seconds n = 5 distance = 3.1250000000e-02 LI = 35 log2(LI) = 5; time = 0 seconds n = 6 distance = 1.5625000000e-02 LI = 68 log2(LI) = 6; time = 0 seconds n = 7 distance = 7.8125000000e-03 LI = 133 log2(LI) = 7; time = 0 seconds n = 8 distance = 3.9062500000e-03 LI = 261 log2(LI) = 8; time = 0 seconds n = 9 distance = 1.9531250000e-03 LI = 518 log2(LI) = 9; time = 0 seconds n = 10 distance = 9.7656250000e-04 LI = 1031 log2(LI) = 10; time = 0 seconds n = 11 distance = 4.8828125000e-04 LI = 2055 log2(LI) = 11; time = 0 seconds n = 12 distance = 2.4414062500e-04 LI = 4104 log2(LI) = 12; time = 0 seconds n = 13 distance = 1.2207031250e-04 LI = 8201 log2(LI) = 13; time = 0 seconds n = 14 distance = 6.1035156250e-05 LI = 16394 log2(LI) = 14; time = 0 seconds n = 15 distance = 3.0517578125e-05 LI = 32778 log2(LI) = 15; time = 0 seconds n = 16 distance = 1.5258789062e-05 LI = 65547 log2(LI) = 16; time = 0 seconds n = 17 distance = 7.6293945312e-06 LI = 131084 log2(LI) = 17; time = 0 seconds n = 18 distance = 3.8146972656e-06 LI = 262156 log2(LI) = 18; time = 0 seconds n = 19 distance = 1.9073486328e-06 LI = 524301 log2(LI) = 19; time = 0 seconds n = 20 distance = 9.5367431641e-07 LI = 1048590 log2(LI) = 20; time = 1 seconds n = 21 distance = 4.7683715820e-07 LI = 2097166 log2(LI) = 21; time = 0 seconds n = 22 distance = 2.3841857910e-07 LI = 4194319 log2(LI) = 22; time = 2 seconds n = 23 distance = 1.1920928955e-07 LI = 8388624 log2(LI) = 23; time = 2 seconds n = 24 distance = 5.9604644775e-08 LI = 16777232 log2(LI) = 24; time = 6 seconds n = 25 distance = 2.9802322388e-08 LI = 33554449 log2(LI) = 25; time = 11 seconds n = 26 distance = 1.4901161194e-08 LI = 67108882 log2(LI) = 26; time = 21 seconds n = 27 distance = 7.4505805969e-09 LI = 134217747 log2(LI) = 27; time = 42 seconds n = 28 distance = 3.7252902985e-09 LI = 268435475 log2(LI) = 28; time = 87 seconds n = 29 distance = 1.8626451492e-09 LI = 536870932 log2(LI) = 29; time = 175 seconds n = 30 distance = 9.3132257462e-10 LI = 1073741845 log2(LI) = 30; time = 351 seconds n = 31 distance = 4.6566128731e-10 LI = 2147483669 log2(LI) = 31; time = 698 seconds n = 32 distance = 2.3283064365e-10 LI = 4294967318 log2(LI) = 32; time = 1386 seconds n = 33 distance = 1.1641532183e-10 LI = 8589934615 log2(LI) = 33; time = 2714 seconds n = 34 distance = 5.8207660913e-11 LI = 17179869207 log2(LI) = 34; time = 5595 seconds n = 35 distance = 2.9103830457e-11 LI = 34359738392 log2(LI) = 35; time = 11175 seconds n = 36 distance = 1.4551915228e-11 LI = 68719476762 log2(LI) = 36; time = 22081 seconds
双曲线情况下最大 n 几乎与有效数字的精度相同[31]
在抛物线情况下，最大 为
最后一次迭代（逃逸时间 = 逃逸迭代次数，其中 abs(zn) > ER）为：在双曲线情况下等于 n
在抛物线情况下等于 2^n
计算时间与迭代次数成正比。双曲线情况下计算时间较短。抛物线情况下随着迭代次数的增加，计算时间会迅速增长。
![{\displaystyle time[seconds]=6*2^{n-24}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/153869fa499a5447fb1bd451b8b3d890b5c9591b)
检查抛物线情况下一个点是否逃逸
- 对于 n = 34，大约需要一个小时（5595 秒）
- 对于 n = 40，大约需要一天
- 对于 n = 45，大约需要一个月
- 对于 n = 50，大约需要一年
有效数字的抵消[32]和有效数字丢失 (LOS)[33][34]
程序由于使用的数字类型的精度有限而失败。大数 (zp) 和小数 (距离) 的加法会得到一个比可保存的位数更多的位数的数字（浮点数类型只有 7 位小数）。一些最右边的位数被取消，迭代进入一个无限循环。
例如：在抛物线情况下使用浮点数时，假设
float cx = 0.25;
float Zpx = 0.5;
float Zx ;
float distance;
float Zx2;
float n = 13;
所以
distance = pow(2.0,-n); // = 1.2207031250e-04 = 0,00012207
它大于机器精度：[35]
distance > FLT_EPSILON // = pow(2, -24) = 5.96e-08 = 0,00000006
所以这个加法仍然有效
Zx = Zpx + distance; // adding big and small number gives 0,50012207
乘法之后得到
Zx2 = Zx*Zx; // = 0,250122
下一步是加法。由于浮点数格式只保存 7 位小数，所以它被截断为
Zx = Zx2 + cx; // = 0,500122 = Zp + (distance/2)
这里相对误差太大，且
d2= 0.0000000149 // distance*distance
小于 FLT_EPSILON/2.0 = 0.0000000596；
解决方案：提高精度！
#include <stdio.h>
#include <math.h> /* pow() */
#include <float.h> /* FLT_EPSILON */
#include <time.h>
#include <fenv.h> /* fegetround() */
int main()
{
float cx = 0.25;
/* Escape Radius ; it defines target set = { z: abs(z) > ER }
all points z in the target set are escaping to infinity */
float ER = 2.0;
float ER2;
time_t start, end;
float dif;
ER2= ER*ER;
float Zpx = 0.5;
float Zx; /* bad value = 0.5002; good value = 0.5004 */
float Zx2; /* Zx2=Zx*Zx */
float i = 0.0;
float d; /* distance between Zpx=1/2 and zx */
float d2; /* d2=d*d; */
int n = 13;
d = pow (2.0, -n);
Zx = Zpx + d;
d2 = d * d;
time (&start);
Zx2 = Zpx * Zpx + 2.0 * d * Zpx + d2;
printf ("Using c with float and Escape Radius = %f \n", ER);
printf ("Round mode is = %d \n", fegetround ());
printf ("i= %3.0f; Zx = %f; Zx2 = %10.8f ; d = %f ; d2 = %.10f\n", i, Zx, Zx2, d, d2);
if (d2 < (FLT_EPSILON / 2.0) )
{
printf("error : relative error to big and d2= %.10f is smaller then FLT_EPSILON/2.0 = %.10f; increase precision ! \a\n",
d2, FLT_EPSILON / 2.0);
return 1;
}
while (Zx2 < ER2) /* ER2=ER*ER */
{
Zx = Zx2 + cx;
d = Zx - Zpx;
d2 = d * d;
Zx2 = 0.25 + d + d2; /* zx2 = zx * zx = (zp + d) * (zp + d) = zp2 +2 * d * zp + d2 = 2.25 + d + d2 */
i += 1.0;
/* printf("i= %3.0f; Zx = %f; Zx2 = %10.8f ; d = %f ; d2 = %f \n", i,Zx, Zx2, d,d2); */
}
time (&end);
dif = difftime (end, start);
printf ("n = %d; distance = %3f; LI = %10.0f log2(LI) = %3.0f time = %2.0lf seconds\n", n, d, i, log2 (i), dif);
return 0;
}
波兰语解释[36]
-
不动点附近各种类型局部离散动力学的距离。参见抛物线动力学是如何复杂的 -
胖杜阿迪兔子朱利亚集不动点附近的轨道。参见分析一条轨道必须包括周期和花瓣
示例 [37]
- 应避免操作双精度格式中低于 的数字。[38]
- 像素密度[43]
- 像素间距是两个二维像素中心之间的距离
/*
precision based on pixel spacing
code by Claude Heiland-Allen
http://mathr.co.uk/
*/
static void dorender(struct view *v, struct mandelbrot_image *image) {
mpfr_div_d(v->pixel_spacing, v->radius, G.height / 2.0, GMP_RNDN);
mpfr_t pixel_spacing_log;
mpfr_init2(pixel_spacing_log, 53);
mpfr_log2(pixel_spacing_log, v->pixel_spacing, GMP_RNDN);
int pixel_spacing_bits = -mpfr_get_d(pixel_spacing_log, GMP_RNDN);
mpfr_clear(pixel_spacing_log);
int interior = 1;
int float_type = 1;
if (interior) {
if (pixel_spacing_bits > 50 / 2) {
float_type = 2;
}
if (pixel_spacing_bits > 60 / 2) {
float_type = 3;
}
} else {
if (pixel_spacing_bits > 50) {
float_type = 2;
}
if (pixel_spacing_bits > 60) {
float_type = 3;
}
}
const char *float_type_str = 0;
switch (float_type) {
case 0: float_type_str = "float"; break;
case 1: float_type_str = "double"; break;
case 2: float_type_str = "long double"; break;
case 3: float_type_str = "mpfr"; break;
default: float_type_str = "?"; break;
}
可以使用此程序检查（自动数学精度[44]）
#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>
/*
what precision of floating point numbers do I need
to draw/compute part of complex plane ( 2D rectangle ) ?
http://fraktal.republika.pl/mandel_zoom.html
https://wikibooks.cn/wiki/Fractals/Computer_graphic_techniques/2D/plane
https://wikibooks.cn/wiki/Fractals/Mathematics/Numerical
uses the code from
https://gitorious.org/~claude
by Claude Heiland-Allen
view = " cre cim radius"
view="-0.75 0 1.5"
view="0.275336142511115 6.75982538465039e-3 0.666e-5"
view="-0.16323108442468427133 1.03436384057374316916 1e-5"
gcc p.c -lmpfr -lgmp -Wall
*/
int main()
{
//declare variables
int height = 720;
int float_type ;
int pixel_spacing_bits;
mpfr_t radius;
mpfr_init2(radius, 53); //
mpfr_set_d(radius, 1.0e-15, GMP_RNDN);
mpfr_t pixel_spacing;
mpfr_init2(pixel_spacing, 53); //
mpfr_t pixel_spacing_log;
mpfr_init2(pixel_spacing_log, 53);
printf ("radius = "); mpfr_out_str (stdout, 10, 0, radius, MPFR_RNDD); putchar ('\n');
// compute
// int mpfr_div_d (mpfr_t rop, mpfr_t op1, double op2, mpfr_rnd_t rnd)
mpfr_div_d(pixel_spacing, radius, height / 2.0, GMP_RNDN);
printf ("pixel_spacing = "); mpfr_out_str (stdout, 10, 0, pixel_spacing, MPFR_RNDD); putchar ('\n');
mpfr_log2(pixel_spacing_log, pixel_spacing, MPFR_RNDN);
printf ("pixel_spacing_log = "); mpfr_out_str (stdout, 10, 0, pixel_spacing_log, MPFR_RNDD); putchar ('\n');
pixel_spacing_bits = -mpfr_get_d(pixel_spacing_log, GMP_RNDN);
printf ("pixel_spacing_bits = %d \n", pixel_spacing_bits);
float_type = 0;
if (pixel_spacing_bits > 40) float_type = 1;
if (pixel_spacing_bits > 50) float_type = 2;
if (pixel_spacing_bits > 60) float_type = 3;
switch (float_type) {
case 0: fprintf(stderr, "render using float \n"); break;
case 1: fprintf(stderr, "render using double \n"); break;
case 2: fprintf(stderr, "render using long double \n"); break;
case 3: fprintf(stderr, "render using MPFR - arbitrary precision \n");
}
return 0;
}
使用参数 m = 2 的帐篷映射计算轨道时出现的数值误差。经过 50 次双精度数迭代后，所有轨道都归零。
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
/*
https://math.stackexchange.com/questions/2453939/is-this-characteristic-of-tent-map-usually-observed
gcc t.c -Wall
a@zelman:~/c/varia/tent$ ./a.out
*/
/* ------------ constants ---------------------------- */
double m = 2.0; /* parameter of tent map */
double a = 1.0; /* upper bound for randum number generator */
int iMax = 100;
/* ------------------- functions --------------------------- */
/*
tent map
*/
double f(double x0, double m){
double x1;
if (x0 < 0.5)
x1 = m*x0;
else x1 = m*(1.0 - x0);
return x1;
}
/* random double from 0.0 to a
https://stackoverflow.com/questions/13408990/how-to-generate-random-float-number-in-c
*/
double GiveRandom(double a){
srand((unsigned int)time(NULL));
return (((double)rand()/(double)(RAND_MAX)) * a);
}
int main(void){
int i = 0;
double x = GiveRandom(a); /* x0 = random */
for (i = 0; i<iMax; i++){
printf("i = %3d \t x = %.16f\n",i, x);
x = f(x,m); /* iteration of the tent map */
}
return 0;
}
结果
i = 0 x = 0.1720333817284710 i = 1 x = 0.3440667634569419 i = 2 x = 0.6881335269138839 i = 3 x = 0.6237329461722323 i = 4 x = 0.7525341076555354 i = 5 x = 0.4949317846889292 i = 6 x = 0.9898635693778584 i = 7 x = 0.0202728612442833 i = 8 x = 0.0405457224885666 i = 9 x = 0.0810914449771332 i = 10 x = 0.1621828899542663 i = 11 x = 0.3243657799085327 i = 12 x = 0.6487315598170653 i = 13 x = 0.7025368803658694 i = 14 x = 0.5949262392682613 i = 15 x = 0.8101475214634775 i = 16 x = 0.3797049570730451 i = 17 x = 0.7594099141460902 i = 18 x = 0.4811801717078197 i = 19 x = 0.9623603434156394 i = 20 x = 0.0752793131687213 i = 21 x = 0.1505586263374425 i = 22 x = 0.3011172526748851 i = 23 x = 0.6022345053497702 i = 24 x = 0.7955309893004596 i = 25 x = 0.4089380213990808 i = 26 x = 0.8178760427981615 i = 27 x = 0.3642479144036770 i = 28 x = 0.7284958288073540 i = 29 x = 0.5430083423852921 i = 30 x = 0.9139833152294159 i = 31 x = 0.1720333695411682 i = 32 x = 0.3440667390823364 i = 33 x = 0.6881334781646729 i = 34 x = 0.6237330436706543 i = 35 x = 0.7525339126586914 i = 36 x = 0.4949321746826172 i = 37 x = 0.9898643493652344 i = 38 x = 0.0202713012695312 i = 39 x = 0.0405426025390625 i = 40 x = 0.0810852050781250 i = 41 x = 0.1621704101562500 i = 42 x = 0.3243408203125000 i = 43 x = 0.6486816406250000 i = 44 x = 0.7026367187500000 i = 45 x = 0.5947265625000000 i = 46 x = 0.8105468750000000 i = 47 x = 0.3789062500000000 i = 48 x = 0.7578125000000000 i = 49 x = 0.4843750000000000 i = 50 x = 0.9687500000000000 i = 51 x = 0.0625000000000000 i = 52 x = 0.1250000000000000 i = 53 x = 0.2500000000000000 i = 54 x = 0.5000000000000000 i = 55 x = 1.0000000000000000 i = 56 x = 0.0000000000000000 i = 57 x = 0.0000000000000000 i = 58 x = 0.0000000000000000 i = 59 x = 0.0000000000000000 i = 60 x = 0.0000000000000000 i = 61 x = 0.0000000000000000 i = 62 x = 0.0000000000000000 i = 63 x = 0.0000000000000000 i = 64 x = 0.0000000000000000 i = 65 x = 0.0000000000000000 i = 66 x = 0.0000000000000000 i = 67 x = 0.0000000000000000 i = 68 x = 0.0000000000000000 i = 69 x = 0.0000000000000000 i = 70 x = 0.0000000000000000 i = 71 x = 0.0000000000000000 i = 72 x = 0.0000000000000000 i = 73 x = 0.0000000000000000 i = 74 x = 0.0000000000000000 i = 75 x = 0.0000000000000000 i = 76 x = 0.0000000000000000 i = 77 x = 0.0000000000000000 i = 78 x = 0.0000000000000000 i = 79 x = 0.0000000000000000 i = 80 x = 0.0000000000000000 i = 81 x = 0.0000000000000000 i = 82 x = 0.0000000000000000 i = 83 x = 0.0000000000000000 i = 84 x = 0.0000000000000000 i = 85 x = 0.0000000000000000 i = 86 x = 0.0000000000000000 i = 87 x = 0.0000000000000000 i = 88 x = 0.0000000000000000 i = 89 x = 0.0000000000000000 i = 90 x = 0.0000000000000000 i = 91 x = 0.0000000000000000 i = 92 x = 0.0000000000000000 i = 93 x = 0.0000000000000000 i = 94 x = 0.0000000000000000 i = 95 x = 0.0000000000000000 i = 96 x = 0.0000000000000000 i = 97 x = 0.0000000000000000 i = 98 x = 0.0000000000000000 i = 99 x = 0.0000000000000000
- Mark McClure 关于计算的一些说明
- Gabriel Peyré 的《数据科学的数值之旅》汇集了 Matlab、Python、Julia 和 R 实验，探索现代数学数据科学。
- 逻辑映射计算中的数值误差
- Keisan 在线计算器
- Gonzalo Galiano Casas 和 Esperanza García Gonzalo 撰写的《Python 中的有限算术和误差分析》
- math.stackexchange 问题：is-there-an-equation-for-the-number-of-iterations-required-to-escape-the-mandelb ?
- ↑ mathoverflow 问题 : /rounding-errors-in-images-of-julia-sets
- ↑ Mark McClure 的 ChaosAndFractals。2.11 关于计算的一些说明
- ↑ gmane.org 讨论中的级数极限
- ↑ math.stackexchange 问题 : levins-u-transformation
- ↑ math.stackexchange 问题 : accelerating-convergence-of-a-sequence
- ↑ gsl-1.0 : 级数加速
- ↑ gsl 手册 : 级数加速
- ↑ 分形论坛 : series-acceleration
- ↑ wolframalpha 级数计算器
- ↑ Fractalshades 文档
- ↑ 你永远不想知道的关于浮点数的信息，但你将被迫从 volkerschatz 那里了解它
- ↑ stackoverflow : 处理 OpenCL 粒子系统中缺乏浮点数精度的问题
- ↑ 使用 GLSL 进行密集计算 - 第 5 部分：由 Henry Thasler 模拟的四精度
- ↑ 区间 [0,1] 中有多少个浮点数？由 Daniel Lemire
- ↑ Fredrik Johansson 的 Arb
- ↑ 使用 Chudnovsky 和 GMP 计算 π 由 Beej Jorgensen
- ↑ [:w:Extended precision: Wikipedia 中的扩展精度]
- ↑ [:w:GNU MPFR|Wikipedia 中的 GNU MPFR]
- ↑ 有效数字和舍入 版权所有 © 2003–2014 Stan Brown
- ↑ mathforum : 有效数字和十进制位数的规则
- ↑ 舍入误差 Robert P. Munafo，1996 年 12 月 3 日。
- ↑ wikipedia : Shadowing_lemma
- ↑ Wolf Jung 的参数射线绘制精度测试
- ↑ 可以使用 Wolf Jung 的程序 Mandel 通过“Ray to point menu position”（或 Y 键）找到它
- ↑ 一个复变量中的动力学：9-5-91 版的导论讲义 附录 G 由 John W. Milnor
- ↑ 抛物线 Julia 集是多项式时间可计算的 由 Mark Braverman
- ↑ 舍入误差 由 Robert P. Munafo，1996 年 12 月 3 日。
- ↑ math.stackexchange 问题 : what-is-the-shape-of-parabolic-critical-orbit
- ↑ 抛物线 Julia 集是多项式时间可计算的 由 Mark Braverman
- ↑ fractalforums : numerical-problem-with-bailout-test
- ↑ wikipedia : 浮点数，IEEE_754
- ↑ wikipedia : 有效数字
- ↑ wikipedia : 精度损失
- ↑ Oracle 的《数值计算指南》中的 IEEE 算术
- ↑ 机器 epsilon
- ↑ pl.comp.os.linux.programowanie 上的波兰语讨论
- ↑ 使用 DEM/M 的精度和曼德尔缩放
- ↑ Fractalshades-doc : math
- ↑ reenigne 博客 : arbitrary precision mandelbrot sets
- ↑ hpdz : Bignum 由 Michael Condron
- ↑ fractint : 高精度和深度缩放
- ↑ chaospro 文档 : parmparm
- ↑ 像素密度
- ↑ 自动数学精度 Robert P. Munafo，1993 年 4 月 15 日。