分形/數學/數值
< 分形
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 |
另請參閱：Rick Regan 的二進位制浮點數的十進位制精度
- "計算時保留比最終答案多兩位有效數字
- 只在計算的最後一步進行舍入。永遠不要用舍入後的數字進行進一步的計算。
- 對於乘法和除法，找到每個因子的有效數字位數。結果將具有較少的有效數字位數。
- 對於冪和根，答案應該具有與原始數字相同的有效數字位數。" [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://wikibook.tw/wiki/Fractals/Computer_graphic_techniques/2D/plane
https://wikibook.tw/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 日。