Mail Archives: djgpp/2001/06/16/07:45:21
From: | invalid AT erehwon DOT invalid (Graaagh the Mighty)
|
Newsgroups: | comp.os.msdos.djgpp
|
Subject: | Mysterious NaNs...
|
Organization: | Low Charisma Anonymous
|
Message-ID: | <3b2b4290.163568911@news.primus.ca>
|
X-Newsreader: | Forte Free Agent 1.11/32.235
|
Lines: | 655
|
Date: | Sat, 16 Jun 2001 11:36:51 GMT
|
NNTP-Posting-Host: | 207.176.153.100
|
X-Complaints-To: | news AT primus DOT ca
|
X-Trace: | news2.tor.primus.ca 992691528 207.176.153.100 (Sat, 16 Jun 2001 07:38:48 EDT)
|
NNTP-Posting-Date: | Sat, 16 Jun 2001 07:38:48 EDT
|
To: | djgpp AT delorie DOT com
|
DJ-Gateway: | from newsgroup comp.os.msdos.djgpp
|
Reply-To: | djgpp AT delorie DOT com
|
This is seriously weird.
This is a hacked-up-for-debugging version of a Mandelbrot set related
program I am working on. This altered version just runs
"f_find_components" and prints loads of debugging info. The odd stuff
is with a certain loop below. Compile with:
gcc periods.c -o per.exe -lalleg -g -O6 -ffast-math
-fomit-frame-pointer -fforce-addr -funroll-loops -march=k6
-malign-double
(gcc version 2.95.2; gas version 2.8.1.)
I suspect it may be an optimizer bug.
There are comments in the trouble area to help explain better; these
start "// DEBUG".
/*
* periods.c
* (c) 2001 Paul Derbyshire
* Lesser Gnu Public License. See http://www.fsf.org.
* Quick b&w Mandelbrot renderer. Finds components of period n and
writes images of them.
* Uses "MDisk" algorithm from _The_Science_of_Fractal_Images_, ISBN
0-387-96608-0
* Requires Allegro and a gcc port to compile. Allegro can be had at
http://www.allegro.cc.
*/
#include <stdio.h>
#include <allegro.h>
#include <math.h>
#include <unistd.h>
#include <stdlib.h>
#include <values.h> // MAXINT
#include <bios.h> // biostime
typedef double ldouble; // Change this when porting to
machines/compilers with no long double.
typedef unsigned int word;
typedef int bool;
int word_bits = 8*sizeof(word);
int ldouble_bits = 64; // Change this when porting to machines with
different long double mantissa size.
int width = 1024;
int height = 768;
double max_ldouble_amag = 0.0;
PALETTE palette;
int p_greys = 0;
int p_void = 255;
int max_iters = 0;
int max_components = 0;
ldouble max_radius_squared = 9.0;
//double edge_width = 0.1;
double edge_width = 0.02;
double min_disk_radius = 1.0;
ldouble per_width = 1.0;
int recursive_disk_radius = 1;
int aa_degree = 1;
int num_components = 0;
ldouble f_cx, f_cy, f_mag, f_x1, f_y1, f_x2, f_y2, f_pix_size,
f_aa_size, *f_x_orbit, *f_y_orbit, *f_x_orbit2, *f_y_orbit2;
int precision_bits = 0;
int num_words = 0;
bool do_per;
bool resume;
int symmetry_line;
BITMAP *main_screen;
int preview_grid_size = 64;
void set_max_iters (int);
void f_find_components (int);
int num_components_of_period (int);
void f_set_cmag (ldouble, ldouble, ldouble);
void image_draw_loop (int);
void draw_image (void);
void q_wait_key (void);
void f_draw_image (void);
void f_do_pixel (int, int, ldouble, ldouble);
void f_do_disk_1 (int, int, ldouble, ldouble);
void f_do_disk (int, int, ldouble, ldouble);
double f_calc_dist (ldouble, ldouble);
double f_calc_dist_per (ldouble, ldouble);
ldouble f_jitter (void);
bool f_check_known_components (ldouble, ldouble);
int main (void) {
int i, per, image_max_iters;
printf("\nInitializing...\n");
srandom(biostime(0, 0));
max_ldouble_amag = pow(2.0, (double)(ldouble_bits - 3));
// ---------------------------------------
// TUNABLE PARAMETERS
// DON'T CHANGE THESE to debug the NaN mystery.
aa_degree = 5;
set_max_iters(16384);
image_max_iters = 2048; // Must be less than arg to set_max_iters!!!
max_components = 10;
edge_width = 0.1;
per = 13;
resume = 0;
// END TUNABLE PARAMETERS
// ---------------------------------------
i = num_components_of_period(per);
if (i <= max_components)
printf("Locating all %d components of period %d.\n", i, per);
else
printf("Locating %d random components of period %d.\n",
max_components, per);
f_find_components(per);
exit(0);
max_iters = image_max_iters;
if (allegro_init()) {
fprintf(stderr, "Initialization error!\n(non-zero return from
allegro_init())\n\n");
return 1;
}
set_color_depth(8);
if (set_gfx_mode(GFX_AUTODETECT, width, height, 0, 0)) {
fprintf(stderr, "Initialization error!\nGraphics mofe %dx%dx256
unavailable.\n(non-zero return from set_gfx_mode())\n\n", width,
height);
return 1;
}
palette[p_void].r = 31;
palette[p_void].g = 0;
palette[p_void].b = 15;
for (i = 0; i < 64; ++i) {
palette[p_greys+i].r = i;
palette[p_greys+i].g = i;
palette[p_greys+i].b = i;
}
set_palette(palette);
main_screen = create_sub_bitmap(screen, 0, 0, width, height);
image_draw_loop(per);
set_gfx_mode(GFX_TEXT, 80, 25, 0, 0);
if (max_iters != 0) {
free(f_x_orbit);
free(f_y_orbit);
free(f_x_orbit2);
free(f_y_orbit2);
}
printf("Done.\n\n");
return 0;
}
void set_max_iters (int max) {
if (max_iters != 0) {
free(f_x_orbit);
free(f_y_orbit);
free(f_x_orbit2);
free(f_y_orbit2);
}
max_iters = max;
f_x_orbit = malloc(max_iters*sizeof(ldouble));
f_y_orbit = malloc(max_iters*sizeof(ldouble));
f_x_orbit2 = malloc(max_iters*sizeof(ldouble));
f_y_orbit2 = malloc(max_iters*sizeof(ldouble));
}
void f_find_components (int per) {
ldouble mag = ((ldouble)(per*per))*0.25;
ldouble cx, cy, ocx, ocy, x, y, xp, yp, t, *x1, *y1;
int i, iters, i2;
int m = 0;
int n = 0;
int q = num_components_of_period(per);
int expected = (q > max_components)?max_components:q;
int seed = random();
int diverged = 0;
int smaller = 0;
int duplicates = 0;
int maxi = 0;
f_pix_size = 1/(width*mag);
f_aa_size = f_pix_size/aa_degree;
if (per == 1) {
(*f_x_orbit2) = 0.0;
(*f_y_orbit2) = 0.0;
num_components = 1;
srandom(seed);
return;
}
if (per == 2) {
(*f_x_orbit2) = -1.0;
(*f_y_orbit2) = 0.0;
num_components = 1;
srandom(seed);
return;
}
srandom(1746573498L);
for (i = 0; i < max_iters; ++i) {
cx = (((ldouble)(random()))/((ldouble)(MAXINT)))*4.0;
cy = (((ldouble)(random()))/((ldouble)(MAXINT)))*3.0;
// This following loop is the oddly behaving one.
for (iters = 0; iters < max_iters; ++iters) {
ocx = cx;
ocy = cy;
x = 0.0;
y = 0.0;
xp = 0.0;
yp = 0.0;
for (i2 = 0; i2 < per; ++i2) {
t = 2.0*(x*xp - y*yp) + 1.0;
yp = 2.0*(x*yp + y*xp);
xp = t;
t = x*x - y*y + cx;
y = 2.0*x*y + cy;
x = t;
}
t = xp*xp + yp*yp;
cx -= (xp*x + yp*y)/t;
cy -= (xp*y - yp*x)/t;
if (!((cx * cx) >= 0.0)) {
// By experiment I've found that this detects NaNs.
// DEBUG: The following will never be executed.
printf("NaNs detected!\ncx = %1.16f, cy = %1.16f, t = %1.16f,
iters = %d.\n", cx, cy, t, iters);
iters = max_iters;
break;
}
if ((fabs(cx - ocx) < f_aa_size) && (fabs(cy - ocy) <
f_aa_size)) {
--iters;
break;
}
}
if (iters > maxi) {
if (iters < max_iters)
maxi = iters + 1;
else
maxi = iters;
}
if (!((cx * cx) >= 0.0)) printf("NaNs detected!\n"); // DEBUG:
This will!
// DEBUG: The weirdness is that between executing the first NaN
// detector in the final loop iteration and the second one here,
// cx is not assigned to. Also, as cx is a local, it's on the
// stack and not the heap so if there is a bug in the pointer
// arithmetic for the f_x_orbit stuff, which I doubt, it's not
// clobbering cx. Also, no pointer arithmetic is invoked between
// the two NaN detectors.
// If the optimizer is doing something weird here, it constitutes
// a bug, since it is clearly altering the program semantics!
// I am aware that f_fast_math may affect the exact values of
// floating point operations, but it shouldn't make floats
// mysteriously change values when they aren't assigned to.
printf("(%1.16f, %1.16f)\n", cx, cy);
printf("(%1.16f, %1.16f)\n", ocx, ocy);
printf("(%1.16f, %1.16f)\n", fabs(cx-ocx), fabs(cy-ocy));
if (iters < max_iters) {
// if ((fabs(cx - ocx) < f_aa_size) && (fabs(cy - ocy) <
f_aa_size)) {
if ((fabs(cx) > 0.2) || (fabs(cy) > 0.2)) {
if (f_calc_dist_per(cx, cy) < (1/(mag*mag) - f_aa_size)) {
// Check for duplication, including symmetric.
x1 = f_x_orbit2;
y1 = f_y_orbit2;
for (i2 = 0; i2 < n; ++i2) {
if ((fabs(cx - (*x1)) < f_aa_size) && (fabs(cy - (*y1)) <
f_aa_size)) break;
if ((fabs(cx - (*x1)) < f_aa_size) && (fabs(cy + (*y1)) <
f_aa_size)) break;
++x1;
++y1;
}
if (i2 == n) {
// It's new!
(*x1) = cx;
(*y1) = fabs(cy);
++n;
++m;
if (fabs(cy) > f_aa_size) ++m;
printf("Got %d of %d.\n", m, expected);
printf("Mag: %1.16f.\n", 1/(64*f_calc_dist_per(*x1,
*y1)));
if (m == q) break; // We won't find any more.
if (m == max_components) break; // We have enough.
} else {
++duplicates;
}
} else {
++smaller;
}
} else {
++smaller;
}
} else {
++diverged;
}
if (!(i & 0x03ff)) printf("Iteration %d: found %d (%d), %d
diverged, %d smaller, %d duplicates; %d maxi.\n", i, n, m, diverged,
smaller, duplicates, maxi);
}
printf("Iteration %d: found %d (%d), %d diverged, %d smaller, %d
duplicates; %d maxi.\n", i, n, m, diverged, smaller, duplicates,
maxi);
if (m < expected) {
fprintf(stderr, "Failed to find %d components!\n\n", expected -
m);
exit(4);
}
num_components = n;
srandom(seed);
}
int num_components_of_period (int per) {
int i;
int n;
if (per == 1) return 1;
if (per > 32) return max_components + 1;
n = (1 << (per - 1)) - 1;
for (i = 2; i < per; ++i) {
if ((per/i)*i == per)
n -= num_components_of_period(i);
}
return n;
}
void f_set_cmag (ldouble cx, ldouble cy, ldouble mag) {
f_cx = cx;
f_cy = cy;
f_mag = mag;
f_pix_size = 1/(width*mag);
f_x1 = f_cx - (width/2)*f_pix_size;
f_x2 = f_cx + (width/2)*f_pix_size;
f_y1 = f_cy - (height/2)*f_pix_size;
f_y2 = f_cy + (height/2)*f_pix_size;
f_aa_size = f_pix_size/aa_degree;
}
void image_draw_loop (int per) {
int i;
ldouble *x1 = f_x_orbit2;
ldouble *y1 = f_y_orbit2;
ldouble mag;
FILE *foo = NULL;
char filename[12];
for (i = 0; i < num_components; ++i) {
sprintf(filename, "p%07d.bmp", i);
if (resume) foo = fopen(filename, "rb");
if (foo != NULL) {
fclose(foo);
} else {
f_set_cmag(*x1, *y1, 1/(64*f_calc_dist_per(*x1, *y1)));
draw_image();
filename[12] = 0;
save_bmp(filename, main_screen, palette); // Main screen turn
on.
mag = ((ldouble)(per*per))*0.25;
f_pix_size = 1/(width*mag);
f_aa_size = f_pix_size/aa_degree;
}
++x1;
++y1;
}
}
void draw_image (void) {
// if (use_float)
f_draw_image();
// else
// i_draw_image();
}
void q_wait_key (void) {
char foo[256];
gets(foo);
}
void f_draw_image (void) {
ldouble cx, cy;
int x, y;
int p = preview_grid_size;
do_per = 0;
rectfill(screen, 0, 0, width-1, height-1, p_void);
cy = f_y2;
symmetry_line = -1;
if ((f_y1 < 0.0) && (f_y2 > 0.0)) {
for (y = 0; y < height; ++y) {
cy -= f_pix_size;
if (cy < 0) {
symmetry_line = 2*y;
break;
}
}
cy = f_y2;
}
// Build "previews".
while (p > 0) {
for (y = 0; y < height; y+=p) {
cx = f_x1;
for (x = 0; x < width; x+=p) {
f_do_pixel(x, y, cx, cy);
cx += p*f_pix_size;
}
cy -= p*f_pix_size;
}
cy = f_y2;
if (p > 1) {
for (y = 0; y < height; y+=p) {
cx = f_x1;
for (x = 0; x < width; x+=(p/2)) {
f_do_pixel(x, y, cx, cy);
cx += p*f_pix_size/2;
}
cy -= p*f_pix_size;
}
cy = f_y2;
}
p/=2;
}
}
void f_do_pixel (int x, int y, ldouble cx, ldouble cy) {
if ((cx < -2.0) || (cx > 0.5) || (cy < -1.2) || (cy > 1.2)) {
putpixel(screen, x, y, p_greys);
if (symmetry_line != -1) putpixel(screen, x, symmetry_line - y,
p_greys);
} else if (getpixel(screen, x, y) == p_void) {
f_do_disk(x, y, cx, cy);
}
}
void f_do_disk_1 (int x, int y, ldouble cx, ldouble cy) {
if ((x < 0) || (x >= width) || (y < 0) || (y >= height)) return;
if ((cx < -2.0) || (cx > 0.5) || (cy < -1.2) || (cy > 1.2)) return;
if (getpixel(screen, x, y) != p_void) return;
f_do_disk(x, y, cx, cy);
}
void f_do_disk (int x, int y, ldouble cx, ldouble cy) {
int radius, aa_x, aa_y, aa_count;
double f_radius, dist, aa_cx, aa_cy;
dist = f_calc_dist(cx, cy);
if (dist < min_disk_radius*f_pix_size) {
aa_count = 0;
aa_cx = cx;
for (aa_x = 0; aa_x < aa_degree; ++aa_x) {
aa_cy = cy;
for (aa_y = 0; aa_y < aa_degree; ++aa_y) {
if ((aa_x != 0) || (aa_y != 0)) dist = f_calc_dist(aa_cx,
aa_cy);
if (dist < edge_width*f_pix_size) ++aa_count;
aa_cy += f_aa_size;
}
aa_cx += f_aa_size;
}
aa_count*=63;
aa_count/=(aa_degree*aa_degree);
putpixel(screen, x, y, p_greys+aa_count);
if (symmetry_line != -1) putpixel(screen, x, symmetry_line - y,
p_greys+aa_count);
} else {
f_radius = dist/f_pix_size;
if (f_radius > (double)width) f_radius = width;
radius = f_radius;
circlefill(screen, x, y, radius, p_greys);
if (symmetry_line != -1) circlefill(screen, x, symmetry_line - y,
radius, p_greys);
if (radius >= recursive_disk_radius) {
f_do_disk_1(x-radius-1, y, cx-(radius+1)*f_pix_size, cy);
f_do_disk_1(x+radius+1, y, cx+(radius+1)*f_pix_size, cy);
f_do_disk_1(x, y-radius-1, cx, cy+(radius+1)*f_pix_size);
f_do_disk_1(x, y+radius+1, cx, cy-(radius+1)*f_pix_size);
}
}
}
double f_calc_dist (ldouble pcx, ldouble pcy) {
ldouble *x1, *y1;
ldouble x = 0.0;
ldouble y = 0.0;
ldouble x2 = 0.0;
ldouble y2 = 0.0;
ldouble xp, yp, t;
int iters, i2;
bool escaped = 0;
ldouble cx = pcx + f_jitter();
ldouble cy = pcy + f_jitter();
if (f_check_known_components(cx, cy)) return 0.5*(edge_width +
min_disk_radius)*f_pix_size;
if (do_per) return f_calc_dist_per(cx, cy);
x1 = f_x_orbit;
y1 = f_y_orbit;
for (iters = 0; iters < max_iters; ++iters) {
y = 2.0*x*y + cy;
x = x2 - y2 + cx;
x2 = x*x;
y2 = y*y;
(*x1) = x;
(*y1) = y;
++x1;
++y1;
if ((x2 + y2) > max_radius_squared) {
escaped = 1;
break;
}
}
if (escaped) {
xp = 1.0;
yp = 0.0;
x1 = f_x_orbit;
y1 = f_y_orbit;
for (i2 = 0; i2 < iters; ++i2) {
t = 2.0*((*x1)*xp - (*y1)*yp) + 1.0;
yp = 2.0*((*x1)*yp + (*y1)*xp);
xp = t;
++x1;
++y1;
}
t = sqrt(x2 + y2);
return t*log(t)*0.5/sqrt(xp*xp + yp*yp);
}
do_per = 1;
return 0.5*(edge_width + min_disk_radius)*f_pix_size;
}
double f_calc_dist_per (ldouble cx, ldouble cy) {
ldouble *x1, *y1;
ldouble x = 0.0;
ldouble y = 0.0;
ldouble x2 = 0.0;
ldouble y2 = 0.0;
ldouble pcx = 0.0;
ldouble pcy = 0.0;
ldouble pcx2 = 0.0;
ldouble pcy2 = 0.0;
ldouble xz, yz, xc, yc, xzz, yzz, xcz, ycz, xp, yp, t, u, z;
int iters, i2, pi;
bool even = 1;
bool escaped = 0;
bool periodic = 0;
x1 = f_x_orbit;
y1 = f_y_orbit;
for (iters = 0; iters < max_iters; ++iters) {
y = 2.0*x*y + cy;
x = x2 - y2 + cx;
x2 = x*x;
y2 = y*y;
(*x1) = x;
(*y1) = y;
++x1;
++y1;
even = 1 - even;
if (even) {
pcy = 2.0*pcx*pcy + cy;
pcx = pcx2 - pcy2 + cx;
pcx2 = pcx*pcx;
pcy2 = pcy*pcy;
}
if ((fabs(x - pcx) < per_width*f_aa_size) && (fabs(y - pcy) <
per_width*f_aa_size)) {
periodic = 1;
break;
}
if ((x2 + y2) > max_radius_squared) {
escaped = 1;
break;
}
}
if (periodic) {
// Find start of period.
--x1;
--y1; // Point back to x and y.
for (pi = iters - 1; pi > 0; --pi) {
--x1;
--y1;
if ((fabs((*x1) - x) < per_width*f_aa_size) && (fabs((*y1) - y)
< per_width*f_aa_size)) break;
}
// Now pi is such that pi...iters-1 is one repetition of the
cycle, and iters-pi is its period. Also, x1 and y1 now
// point at z_{pi}.
xz = 1.0;
yz = 0.0;
xc = 0.0;
yc = 0.0;
xzz = 0.0;
yzz = 0.0;
xcz = 0.0;
ycz = 0.0;
for (i2 = pi; i2 < iters; ++i2) {
t = 2.0*(xz*xz - yz*yz + (*x1)*xzz - (*y1)*yzz);
yzz = 2.0*(2.0*xz*yz + (*x1)*yzz + (*y1)*xzz);
xzz = t;
t = 2.0*(xz*xc - yz*yc + (*x1)*xcz - (*y1)*ycz);
ycz = 2.0*(xz*yc + yz*xc + (*x1)*ycz + (*y1)*xcz);
xcz = t;
t = 2.0*((*x1)*xz - (*y1)*yz);
yz = 2.0*((*x1)*yz + (*y1)*xz);
xz = t;
t = 2.0*((*x1)*xc - (*y1)*yc) + 1.0;
yc = 2.0*((*x1)*yc + (*y1)*xc);
xc = t;
++x1;
++y1;
}
t = 1.0 - xz;
u = -yz;
z = t*t + u*u;
x = (xc*t + yc*u)/z;
y = (yc*t - xc*u)/z;
xp = xcz + xzz*x - yzz*y;
yp = ycz + xzz*y + yzz*x;
return 0.25*(1.0 - xz*xz - yz*yz)/sqrt(xp*xp + yp*yp);
} else if (escaped) {
do_per = 0;
xp = 1.0;
yp = 0.0;
x1 = f_x_orbit;
y1 = f_y_orbit;
for (i2 = 0; i2 < iters; ++i2) {
t = 2.0*((*x1)*xp - (*y1)*yp) + 1.0;
yp = 2.0*((*x1)*yp + (*y1)*xp);
xp = t;
++x1;
++y1;
}
t = sqrt(x2 + y2);
return t*log(t)*0.5/sqrt(xp*xp + yp*yp);
}
return 0.5*(edge_width + min_disk_radius)*f_pix_size;
}
ldouble f_jitter (void) {
return ((((ldouble)(random()))/((ldouble)(MAXINT))) -
0.5)*f_aa_size;
}
bool f_check_known_components (ldouble cx, ldouble cy) {
ldouble t = cx + 1;
ldouble u, v, w, a, m;
if ((t*t + cy*cy) <= 0.0625) return 1;
t = 1 - 4*cx;
u = -4*cy;
a = atan2(u, t)/2;
m = sqrt(sqrt(t*t + u*u));
v = 1 + m*cos(a);
w = m*sin(a);
if ((v*v + w*w) <= 1) return 1;
v = 2 - v;
if ((v*v + w*w) <= 1) return 1;
return 0;
}
--
Bill Gates: "No computer will ever need more than 640K of RAM." -- 1980
"There's nobody getting rich writing software that I know of." -- 1980
"This antitrust thing will blow over." -- 1998
Combine neo, an underscore, and one thousand sixty-one to make my hotmail addy.
- Raw text -