Generate mandelbrot images using many clustered computers
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

237 lines
5.1 KiB

  1. #include<stdio.h>
  2. #include<stdlib.h>
  3. #include<mpi.h>
  4. #include<math.h>
  5. #include<gmp.h>
  6. #include "types.h"
  7. void map(mpf_t out, double x, double in_min, double in_max, mpf_t out_min, mpf_t out_max, int precision)
  8. {
  9. mpf_t a, b, c;
  10. mpf_init2(a, precision);
  11. mpf_init2(b, precision);
  12. mpf_init2(c, precision);
  13. mpf_set_d(a, x);
  14. mpf_set_d(b, in_min);
  15. mpf_set_d(c, in_max);
  16. mpf_sub(a, a, b);
  17. mpf_sub(c, c, b);
  18. mpf_sub(b, out_max, out_min);
  19. mpf_mul(a, a, b);
  20. mpf_div(a, a, c);
  21. mpf_add(out, a, out_min);
  22. mpf_clear(a);
  23. mpf_clear(b);
  24. mpf_clear(c);
  25. //return (x - in_min) * (out_max - out_min) / (in_max - in_min) + out_min;
  26. }
  27. Pixel hsl_to_rgb(double h, double s, double l)
  28. {
  29. while(h >= 360) h -= 360;
  30. if(s > 1) s = 1;
  31. if(l > 1) l = 1;
  32. double c = (1 - abs(2 * l - 1)) * s;
  33. double x = c * (1 - fabs(fmod(((double)h)/60, 2) - 1));
  34. double m = l - c / 2;
  35. double r, g, b;
  36. int q = h / 60;
  37. switch(q)
  38. {
  39. case(0):
  40. r = c;
  41. g = x;
  42. b = 0;
  43. break;
  44. case(1):
  45. r = x;
  46. g = c;
  47. b = 0;
  48. break;
  49. case(2):
  50. r = 0;
  51. g = c;
  52. b = x;
  53. break;
  54. case(3):
  55. r = 0;
  56. g = x;
  57. b = c;
  58. break;
  59. case(4):
  60. r = x;
  61. g = 0;
  62. b = c;
  63. break;
  64. case(5):
  65. r = c;
  66. g = 0;
  67. b = x;
  68. }
  69. Pixel out;
  70. out.red = 255 * (r + m);
  71. out.green = 255 * (g + m);
  72. out.blue = 255 * (b + m);
  73. return out;
  74. }
  75. Pixel mandelbrot(int pixId, int w, int h,
  76. mpf_t xMin, mpf_t xMax, mpf_t yMin, mpf_t yMax,
  77. long range, long maxIter, int precision)
  78. {
  79. int x = pixId % w;
  80. int y = pixId / w;
  81. //double bounds[] = {center[0] - range, center[1] - range,
  82. // center[0] + range, center[1] + range};
  83. mpf_t Creal;
  84. mpf_init2(Creal, precision);
  85. map(Creal, x, 0, w, xMin, xMax, precision);
  86. //printf("%lf %lf %lf\r\n",mpf_get_d(rangeG), mpf_get_d(xMin), mpf_get_d(xMax));
  87. //double Creal = map(x, 0, w, bounds[0], bounds[2]);
  88. mpf_t Ci;
  89. mpf_init2(Ci, precision);
  90. map(Ci, y, 0, h, yMin, yMax, precision);
  91. //double Cimaginary = map(y, 0, h, bounds[1], bounds[3]);
  92. Pixel rtn;
  93. mpf_t real, imag, tmp, imagS;
  94. mpf_init2(real, precision);
  95. mpf_init2(imag, precision);
  96. mpf_init2(tmp, precision);
  97. mpf_init2(imagS, precision);
  98. mpf_set_d(real, 0);
  99. mpf_set_d(imag, 0);
  100. int i = 0;
  101. for(; i < maxIter; i++)
  102. {
  103. mpf_set(tmp, real);
  104. //real = real*real - imaginary*imaginary + Creal;
  105. mpf_mul(real, real, real); //real = real^2
  106. mpf_mul(imagS, imag, imag);
  107. mpf_sub(real, real, imagS);
  108. mpf_add(real, real, Creal);
  109. //imaginary = realTmp * 2 * imaginary + Cimaginary;
  110. mpf_mul_ui(imag, imag, 2);
  111. mpf_mul(imag, imag, tmp);
  112. mpf_add(imag, imag, Ci);
  113. //if(imaginary * imaginary + real * real > 4 ||
  114. // imaginary * imaginary + real * real < -4)
  115. mpf_mul(imagS, imag, imag);
  116. mpf_mul(tmp, real, real);
  117. mpf_add(tmp, imagS, tmp);
  118. if(mpf_cmp_si(tmp, 1<<16) > 0)
  119. {
  120. //We can use regular doubles here
  121. double log_zn = log2(mpf_get_d(tmp)) / 2;
  122. double nu = log2(log_zn); //this can be optimized
  123. double il = i + 1 - nu;
  124. rtn = hsl_to_rgb(il * 10, 1, 0.5);
  125. break;
  126. }
  127. }
  128. if(i == maxIter) rtn = (Pixel){.red=0, .green=0, .blue=0};
  129. mpf_clears(real, imag, tmp, imagS, Ci, Creal, NULL);
  130. return rtn;
  131. }
  132. void getConfig(Config *conf)
  133. {
  134. MPI_Recv(conf, sizeof(Config), MPI_BYTE, 0, 0,
  135. MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  136. }
  137. void slave_main(int slaveCount, int id)
  138. {
  139. Config conf;
  140. Pixel *buf = NULL;
  141. while(1)
  142. {
  143. long curPixel;
  144. MPI_Recv(&curPixel, 1, MPI_LONG, 0, 1,
  145. MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  146. if(curPixel == -1) break;
  147. if(curPixel == -2)
  148. {
  149. getConfig(&conf);
  150. if(buf == NULL)
  151. {
  152. buf = malloc(sizeof(Pixel) * conf.batch);
  153. }
  154. }
  155. else
  156. {
  157. int amtToSend = MIN(conf.batch,
  158. ((long)conf.w * conf.h) - curPixel);
  159. MPI_Recv(buf, sizeof(Pixel) * amtToSend, MPI_BYTE, 0, 0,
  160. MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  161. //Processing
  162. int precision = 4 * abs(conf.zoomExp) * MAX(10, log(conf.zoomBase));
  163. //printf("Precision: %d\r\n", precision);
  164. mpf_t cX;
  165. mpf_t cY;
  166. mpf_init2(cX, precision);
  167. mpf_init2(cY, precision);
  168. mpf_set_str(cX, conf.cX_str, 10);
  169. mpf_set_str(cY, conf.cY_str, 10);
  170. //rangeG = 10^range
  171. mpf_t rangeGx;
  172. mpf_init2(rangeGx, precision);
  173. mpf_set_d(rangeGx, conf.zoomBase);
  174. if(conf.zoomExp < 0)
  175. {
  176. mpf_pow_ui(rangeGx, rangeGx, -conf.zoomExp); //raise to 10^-range which will be huge
  177. mpf_ui_div(rangeGx, 1, rangeGx); //take reciprocol
  178. }
  179. else
  180. {
  181. mpf_pow_ui(rangeGx, rangeGx, conf.zoomExp);
  182. }
  183. mpf_t rangeGy;
  184. mpf_init2(rangeGy, precision);
  185. //Scale to screen resolution
  186. mpf_mul_ui(rangeGy, rangeGx, conf.h);
  187. mpf_mul_ui(rangeGx, rangeGx, conf.w);
  188. mpf_t xMin, xMax, yMin, yMax;
  189. mpf_init2(xMin, precision);
  190. mpf_init2(xMax, precision);
  191. mpf_init2(yMin, precision);
  192. mpf_init2(yMax, precision);
  193. mpf_sub(xMin, cX, rangeGx);
  194. mpf_add(xMax, cX, rangeGx);
  195. mpf_sub(yMin, cY, rangeGy);
  196. mpf_add(yMax, cY, rangeGy);
  197. for(int i = 0; i < amtToSend; i++)
  198. {
  199. buf[i] = mandelbrot(curPixel + i, conf.w, conf.h,
  200. xMin, xMax, yMin, yMax, conf.zoomExp,
  201. conf.maxIter, precision);
  202. }
  203. mpf_clears(cX, cY, rangeGx, rangeGy, xMin, xMax, yMin, yMax, NULL);
  204. MPI_Send(buf, sizeof(Pixel) * amtToSend, MPI_BYTE, 0, 0,
  205. MPI_COMM_WORLD);
  206. }
  207. }
  208. free(buf);
  209. }