-
Notifications
You must be signed in to change notification settings - Fork 0
/
brat-gap-hair.C
117 lines (97 loc) · 1.98 KB
/
brat-gap-hair.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
/*
* brat-gap-hair.C
*
* FUNCTION:
* The gap-tongue tries to draw the basic gaps
* of the continued fraction. Te resulting plots look like wavey hair
* or wavey seaweed.
*
* HISTORY:
* more stuff -- October 2004
*/
#include <malloc.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "brat.h"
#include "Farey.h"
#include "FareyTree.h"
void
MakeHisto (
float *glob,
int sizex,
int sizey,
double re_center,
double im_center,
double width,
double height,
int itermax,
double renorm)
{
int i,j;
int globlen = sizex*sizey;
for (i=0; i<globlen; i++) {
glob [i] = 0.0;
}
int d,n; // denom, numerator
ContinuedFraction f;
d = itermax;
for (n=1; n<d; n++)
{
int nn = n;
int dd = d;
// #define DO_TONG_RAND
#ifdef DO_TONG_RAND
nn = rand() >> 10;
dd = rand() >> 10;
#endif
nn %= dd;
if (0 == nn) continue;
if (0 == dd) continue;
int gcf = gcf64 (nn,dd);
nn /= gcf;
dd /= gcf;
f.SetRatio (nn,dd);
double x = (double)nn/(double) dd;
for (j=0; j<sizey; j++)
{
double w = (((double) (sizey-j))-0.5)/((double) sizey);
// double gap = f.ToXEven(w);
// double gap = f.ToXPlus(w);
#define REMOVE_LEADING_TERMS
#ifdef REMOVE_LEADING_TERMS
double gap = f.ToXEven(w);
gap -= x;
gap *= (double) dd;
gap *= (double) dd;
gap -= w;
gap += w*w;
gap -= 0.5*w*w*w;
gap += x;
#endif
#ifdef WHATEVER
double gap = f.ToXOdd(w);
gap -= (1.0-x);
gap *= (double) dd;
gap *= (double) dd;
gap += w;
gap += w*w;
gap -= 0.5*w*w*w;
gap += (1.0-x);
#endif
// printf ("duude x=%d/%d = %g w=%g gap=%g\n", nn, dd, x, w, gap);
i = (int) (gap * (double) sizex);
if (0>i) continue;
if (i>=sizex) continue;
glob [j*sizex +i] ++;
}
}
/* renormalize */
double r = ((double) sizex) / ((double) itermax);
for (i=0; i<sizex*sizey; i++)
{
glob [i] *= r;
}
}
/* --------------------------- END OF LIFE ------------------------- */