store references in view and compute reference in box
[maximus:mightymandel.git] / src / fpxx_approx.c
1 #include "fpxx_approx.h"
2 #include "shader.h"
3
4 extern const char *fpxx_approx_vert;
5 const GLchar *fpxx_approx_varyings[] = {"cne", "zdz"};
6
7 void fpxx_approx_begin(struct fpxx_approx *s) {
8   s->program = 0;
9   s->iters = -1;
10   s->c0 = -1;
11   s->program = compile_program("fpxx_approx", fpxx_approx_vert, 0, 0, 2, fpxx_approx_varyings);
12   s->iters = glGetUniformLocation(s->program, "iters");D;
13   s->c0 = glGetAttribLocation(s->program, "c0");D;
14   mpfr_inits2(53
15     , s->p.x, s->p.y, s->z.x, s->z.y, s->delta4
16     , s->a.x, s->a.y, s->b.x, s->b.y, s->c.x, s->c.y
17     , s->u.x, s->u.y, s->v.x, s->v.y, s->w.x, s->w.y
18     , s->t1.x, s->t1.y, s->t2.x, s->t2.y
19     , s->t3.x, s->t3.y, s->t4.x, s->t4.y
20     , s->t5.x, s->t5.y, s->t6.x, s->t6.y
21     , (mpfr_ptr) 0);
22 }
23
24 void fpxx_approx_end(struct fpxx_approx *s) {
25   glDeleteProgram(s->program);D;
26   s->program = 0;
27   mpfr_clears
28     ( s->p.x, s->p.y, s->z.x, s->z.y, s->delta4
29     , s->a.x, s->a.y, s->b.x, s->b.y, s->c.x, s->c.y
30     , s->u.x, s->u.y, s->v.x, s->v.y, s->w.x, s->w.y
31     , s->t1.x, s->t1.y, s->t2.x, s->t2.y
32     , s->t3.x, s->t3.y, s->t4.x, s->t4.y
33     , s->t5.x, s->t5.y, s->t6.x, s->t6.y
34     , (mpfr_ptr) 0);
35 }
36
37 int fpxx_approx_do(struct fpxx_approx *s, GLuint *active_count, GLuint *vbo, GLuint query, struct view *v) {
38   mpfr_sqr(s->delta4, v->radius, MPFR_RNDN);
39   mpfr_sqr(s->delta4, s->delta4, MPFR_RNDN);
40   mpfr_prec_t p = mpfr_get_prec(v->refx);
41   c_set_prec(s->p, p);
42   c_set_prec(s->z, p);
43   c_set_prec(s->a, p);
44   c_set_prec(s->b, p);
45   c_set_prec(s->c, p);
46   c_set_prec(s->u, p);
47   c_set_prec(s->v, p);
48   c_set_prec(s->w, p);
49   c_set_prec(s->t1, p);
50   c_set_prec(s->t2, p);
51   c_set_prec(s->t3, p);
52   c_set_prec(s->t4, p);
53   c_set_prec(s->t5, p);
54   c_set_prec(s->t6, p);
55   mpfr_set(s->p.x, v->refx, MPFR_RNDN); mpfr_set(s->p.y, v->refy, MPFR_RNDN);
56   mpfr_set(s->z.x, v->refx, MPFR_RNDN); mpfr_set(s->z.y, v->refy, MPFR_RNDN);
57   mpfr_set_ui(s->a.x, 1, MPFR_RNDN); mpfr_set_ui(s->a.y, 0, MPFR_RNDN);
58   mpfr_set_ui(s->b.x, 0, MPFR_RNDN); mpfr_set_ui(s->b.y, 0, MPFR_RNDN);
59   mpfr_set_ui(s->c.x, 0, MPFR_RNDN); mpfr_set_ui(s->c.y, 0, MPFR_RNDN);
60   mpfr_set_ui(s->u.x, 0, MPFR_RNDN); mpfr_set_ui(s->u.y, 0, MPFR_RNDN);
61   mpfr_set_ui(s->v.x, 0, MPFR_RNDN); mpfr_set_ui(s->v.y, 0, MPFR_RNDN);
62   mpfr_set_ui(s->w.x, 0, MPFR_RNDN); mpfr_set_ui(s->w.y, 0, MPFR_RNDN);
63   unsigned int n = 1;
64   bool accurate = true;
65   while (accurate) {
66     s->values[0][0] = mpfr_get_d(s->a.x, MPFR_RNDN);
67     s->values[0][1] = mpfr_get_d(s->a.y, MPFR_RNDN);
68     s->values[1][0] = mpfr_get_d(s->b.x, MPFR_RNDN);
69     s->values[1][1] = mpfr_get_d(s->b.y, MPFR_RNDN);
70     s->values[2][0] = mpfr_get_d(s->c.x, MPFR_RNDN);
71     s->values[2][1] = mpfr_get_d(s->c.y, MPFR_RNDN);
72     s->values[3][0] = mpfr_get_d(s->u.x, MPFR_RNDN);
73     s->values[3][1] = mpfr_get_d(s->u.y, MPFR_RNDN);
74     s->values[4][0] = mpfr_get_d(s->v.x, MPFR_RNDN);
75     s->values[4][1] = mpfr_get_d(s->v.y, MPFR_RNDN);
76     s->values[5][0] = mpfr_get_d(s->w.x, MPFR_RNDN);
77     s->values[5][1] = mpfr_get_d(s->w.y, MPFR_RNDN);
78     n += 1;
79     // w <- 2 (a c + z w + v)
80     c_mul(s->w,  s->z, s->w, s->t3.x, s->t3.y, s->t4.x, s->t4.y);
81     c_mul(s->t1, s->a, s->c, s->t3.x, s->t3.y, s->t4.x, s->t4.y);
82     c_add(s->w,  s->w, s->t1);
83     c_add(s->w,  s->w, s->v);
84     c_mul_2ui(s->w, s->w, 1);
85     // v <- 2 (a b + z v + u)
86     c_mul(s->v,  s->z, s->v, s->t3.x, s->t3.y, s->t4.x, s->t4.y);
87     c_mul(s->t2, s->a, s->b, s->t3.x, s->t3.y, s->t4.x, s->t4.y);
88     c_add(s->v,  s->v, s->t2);
89     c_add(s->v,  s->v, s->u);
90     c_mul_2ui(s->v, s->v, 1);
91     // u <- 2 (a a + z u)
92     c_mul(s->u, s->z, s->u, s->t3.x, s->t3.y, s->t4.x, s->t4.y);
93     c_sqr(s->t1, s->a, s->t3.y, s->t4.x, s->t4.y);
94     c_add(s->u, s->u, s->t1);
95     c_mul_2ui(s->u, s->u, 1);
96     // c <- 2 (a b + z c)
97     c_mul(s->c, s->z, s->c, s->t3.x, s->t3.y, s->t4.x, s->t4.y);
98     c_add(s->c, s->t2, s->c);
99     c_mul_2ui(s->c, s->c, 1);
100     // b <- 2 z b + a a
101     c_mul(s->b, s->z, s->b, s->t3.x, s->t3.y, s->t4.x, s->t4.y);
102     c_mul_2ui(s->b, s->b, 1);
103     c_add(s->b, s->t1, s->b);
104     // a <- 2 z a + 1
105     c_mul(s->a, s->z, s->a, s->t3.x, s->t3.y, s->t4.x, s->t4.y);
106     c_mul_2ui(s->a, s->a, 1);
107     mpfr_add_ui(s->a.x, s->a.x, 1, MPFR_RNDN);
108     // z <- z z + p
109     c_sqr(s->z, s->z, s->t5.y, s->t6.x, s->t6.y);
110     c_add(s->z, s->z, s->p);
111     // accurate <- |c|^2 d^4 << |a|^2 && |w|^2 d^4 << |u|^2
112     c_mag2(s->t1.x, s->c, s->t4.x, s->t4.y);
113     c_mag2(s->t1.y, s->a, s->t4.x, s->t4.y);
114     mpfr_div(s->t1.x, s->t1.x, s->t1.y, MPFR_RNDN);
115     mpfr_mul(s->t1.x, s->t1.x, s->delta4, MPFR_RNDN);
116     c_mag2(s->t2.x, s->w, s->t4.x, s->t4.y);
117     c_mag2(s->t2.y, s->u, s->t4.x, s->t4.y);
118     mpfr_div(s->t2.x, s->t2.x, s->t2.y, MPFR_RNDN);
119     mpfr_mul(s->t2.x, s->t2.x, s->delta4, MPFR_RNDN);
120     const double eps = 1e-40;
121     accurate = mpfr_get_d(s->t1.x, MPFR_RNDN) < eps && mpfr_get_d(s->t2.x, MPFR_RNDN) < eps;
122   }
123   n -= 1;
124   fprintf(stderr, "skip: %d\n", n);
125   glEnable(GL_RASTERIZER_DISCARD);D;
126   glBindBuffer(GL_ARRAY_BUFFER, vbo[0]);D;
127   glUseProgram(s->program);D;
128   glUniform2dv(s->abcuvw, 6, &s->values[0][0]);D;
129   glUniform1i(s->iters, n);D;
130   glVertexAttribLPointer(s->c0, 2, GL_DOUBLE, 8 * sizeof(GLdouble), 0);D;
131   glEnableVertexAttribArray(s->c0);D;
132   glBindBufferBase(GL_TRANSFORM_FEEDBACK_BUFFER, 0, vbo[1]);D;
133   glBeginQuery(GL_TRANSFORM_FEEDBACK_PRIMITIVES_WRITTEN, query);D;
134   glBeginTransformFeedback(GL_POINTS);D;
135   glDrawArrays(GL_POINTS, 0, *active_count);D;
136   glEndTransformFeedback();D;
137   glEndQuery(GL_TRANSFORM_FEEDBACK_PRIMITIVES_WRITTEN);D;
138   glGetQueryObjectuiv(query, GL_QUERY_RESULT, active_count);D;
139   glBindBufferBase(GL_TRANSFORM_FEEDBACK_BUFFER, 0, 0);D;
140   glDisableVertexAttribArray(s->c0);D;
141   glUseProgram(0);D;
142   glBindBuffer(GL_ARRAY_BUFFER, 0);D;
143   glDisable(GL_RASTERIZER_DISCARD);D;
144   GLuint tmp = vbo[1]; vbo[1] = vbo[0]; vbo[0] = tmp;
145   return n;
146 }