home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_01_01 / 1n01031a < prev    next >
Text File  |  1990-05-17  |  6KB  |  249 lines
  1. /**********
  2. **
  3. ** Listing one: NEWT.C
  4. ** Lee Daniel Crocker
  5. ** 3/15/90
  6. */
  7.  
  8. #include <stdlib.h>
  9. #include <conio.h>
  10. #include <graph.h>
  11. #include <float.h>
  12. #include <math.h>
  13.  
  14. typedef struct { double r, i; } COMPLEX;
  15. typedef double REAL;
  16.  
  17. /* #define ASSEMBLY 1 */
  18.  
  19. #define DEGREE  8
  20. #define ILIMIT  1000    /* Iteration count limit.   */
  21. #define EPSILON 0.000001
  22.  
  23. #define PI      3.14159265358979323846
  24.  
  25. #define WIDTH   320     /* Assume CGA.              */
  26. #define HEIGHT  200
  27. #define MODE    _MRES4COLOR
  28.  
  29. #define YMIN    -3.0    /* Rectangular section      */
  30. #define YMAX    3.0     /* of complex plane.        */
  31. #define XMIN    -4.0
  32. #define XMAX    4.0
  33.  
  34. COMPLEX roots[DEGREE];  /* Precalculated values.    */
  35. REAL eps2;
  36.  
  37. #ifdef ASSEMBLY
  38.  
  39. extern void cmult(COMPLEX *, COMPLEX *, COMPLEX *);
  40. extern int cdiv(COMPLEX *, COMPLEX *, COMPLEX *);
  41. extern void cpower(COMPLEX *, unsigned, COMPLEX *);
  42. extern int calcnewton(COMPLEX *);
  43.  
  44. #else
  45.  
  46. /* In all of the complex math procedures below, complex
  47. ** arguments are passed through pointers for efficiency.
  48. ** Any or all of the arguments passed may point to the
  49. ** same structure without affecting the result, and only
  50. ** the structure pointed to by the result argument is
  51. ** modified (unless, of course, the same pointer is also
  52. ** a non-result argument).
  53. */
  54. /* Multiply two complex numbers producing a     complex
  55. ** result.  Knuth's trick of doing it with three
  56. ** multiplies and seven adds doesn't help here because
  57. ** the 87's multiply is not much slower that its add.
  58. */
  59. void cmult(COMPLEX *z1, COMPLEX *z2, COMPLEX *result)
  60. {
  61.     REAL tr;
  62.  
  63.     tr = z1->r * z2->r - z1->i * z2->i;
  64.     result->i = z1->r * z2->i + z1->i * z2->r;
  65.     result->r = tr;
  66. }
  67.  
  68. /* Divide two complex numbers producing a complex result.
  69. ** returns 0 if all is well, -1 if overflow would occur.
  70. */
  71. int cdiv(COMPLEX *z1, COMPLEX *z2, COMPLEX *result)
  72. {
  73.     COMPLEX tc;
  74.     REAL m, tr;
  75.  
  76.     m = (z2->r * z2->r) + (z2->i * z2->i);
  77.     if (m < FLT_MIN) return -1;
  78.     else m = 1.0 / m;
  79.  
  80.     tc = *z2;
  81.     tr = m * (z1->r * tc.r + z1->i * tc.i);
  82.     result->i = m * (z1->i * tc.r - z1->r * tc.i);
  83.     result->r = tr;
  84.  
  85.     return 0;
  86. }
  87.  
  88. /* Raise complex number to integer power.  Operates by
  89. ** repeatedly squaring the base, multiplying that by the
  90. ** result for every bit set in the exponent, therefore
  91. ** running in O(log2(exp)) time.
  92. */
  93. void cpower(COMPLEX *base, unsigned exp, COMPLEX *result)
  94. {
  95.     REAL xt, yt, t2;
  96.  
  97.     xt = base->r;   yt = base->i;
  98.  
  99.     if (exp & 1) {
  100.         result->r = xt;
  101.         result->i = yt;
  102.     } else {
  103.         result->r = 1.0;
  104.         result->i = 0.0;
  105.     }
  106.     exp >>= 1;
  107.  
  108.     while (exp) {
  109.         t2 = (xt + yt) * (xt - yt);
  110.         yt = 2 * xt * yt;
  111.         xt = t2;
  112.  
  113.         if (exp & 1) {
  114.             t2 = xt * result->r - yt * result->i;
  115.             result->i = result->i * xt + yt * result->r;
  116.             result->r = t2;
  117.         }
  118.         exp >>= 1;
  119.     }
  120. }
  121.  
  122. /* Calculate one iteration.  The formula is z(n+1)=
  123. ** z(n) - f(z(n)) / f'(z(n), where z(n)=old, z(n+1)=new,
  124. ** f(z) = z^DEGREE, and f'(z) is the first derivative of
  125. ** f(z), in this case DEGREE*z^(DEGREE-1).  Both old and
  126. ** new are modified, and must point to different
  127. ** structures.  Returns 1 if new point is sufficently
  128. ** close to old, or if a divide overflow occurred.
  129. */
  130. int onenewton(COMPLEX *old, COMPLEX *new)
  131. {
  132.     COMPLEX tc;             /* Temporary complex */
  133.  
  134.     /* Set new=old^DEGREE, tc=old^(DEGREE-1).
  135.     */
  136.     cpower(old, DEGREE-1, &tc);
  137.     cmult(&tc, old, new);
  138.  
  139.     new->r = new->r * (DEGREE-1) + 1.0;
  140.     new->i *= (DEGREE-1);
  141.  
  142.     tc.r *= DEGREE;
  143.     tc.i *= DEGREE;
  144.  
  145.     cdiv(new, &tc, new);
  146.  
  147.     if (fabs(new->r - old->r) < EPSILON &&
  148.         fabs(new->i - old->i) < EPSILON) return 1;
  149.     else return 0;
  150. }
  151.  
  152. /* Calculate the color of the point given.  Argument is
  153. ** not modified.
  154. */
  155. int calcnewton(COMPLEX *point)
  156. {
  157.     int i, root, iterations;
  158.     COMPLEX old, new;
  159.  
  160.     old = *point;
  161.     iterations = 0;
  162.  
  163.     while (++iterations < ILIMIT) {
  164.         if (onenewton(&old, &new)) break;
  165.         else old = new;
  166.     }
  167.  
  168.     /* To which root did we converge?  Begin by assuming
  169.     ** root 0, then test each of the other roots exiting
  170.     ** if one of them works.  Thus, we don't ever have
  171.     ** to test root 0 explicitly.
  172.     */
  173.     root = 0;
  174.     for (i=DEGREE-1; i>0; --i) {
  175.         if (fabs(new.r - roots[i].r) < eps2 &&
  176.             fabs(new.i - roots[i].i) < eps2) {
  177.             root = i;
  178.             break;
  179.         }
  180.     }
  181.  
  182.     /* At this point, <iterations> has the iteration count and
  183.     ** <root> has the root number.  These can be combined in
  184.     ** various ways to produce a color for interesting effects.
  185.     ** in this case, we will simply assign a color based on
  186.     ** the root, ignoring the iteration count.
  187.     */
  188.     return root;
  189. }
  190.  
  191. #endif  /* ASSEMBLY */
  192.  
  193. /* Main program:
  194. */
  195. int newtonplot(void)
  196. {
  197.     int i, px, py, color;
  198.     COMPLEX point;
  199.     REAL xinc, yinc, theta, dtheta;
  200.  
  201.     point.r = XMIN;
  202.     point.i = YMIN;
  203.  
  204.     /* X and Y increments */
  205.  
  206.     xinc = (XMAX - XMIN) / (WIDTH - 1);             
  207.     yinc = (YMAX - YMIN) / (HEIGHT - 1);
  208.  
  209.     eps2 = 0.5 / DEGREE;
  210.     theta = dtheta = (2.0 * PI) / DEGREE;
  211.     roots[0].r = 1.0; roots[0].i = 0.0;
  212.  
  213.     for (i = DEGREE-1; i > 0; --i) {
  214.         roots[i].r = cos(theta);
  215.         roots[i].i = sin(theta);
  216.         theta += dtheta;
  217.     }
  218.  
  219.     _setvideomode(MODE);
  220.  
  221.     for (py = 0; py < HEIGHT; ++py) {
  222.         for (px = 0; px < WIDTH; ++px) {
  223.  
  224.             color = calcnewton(&point);
  225.  
  226.             _setcolor(color);
  227.             _setpixel(px, py);
  228.  
  229.             point.r += xinc;
  230.  
  231.             if (kbhit()) break;
  232.         }
  233.         if (kbhit()) break;
  234.  
  235.         point.i += yinc;
  236.         point.r = XMIN;
  237.     }
  238.  
  239.     getch();
  240.     _setvideomode(_DEFAULTMODE);
  241.  
  242.     return 0;
  243. }
  244.  
  245. int main(int argc, char *argv[])
  246. {
  247.     return newtonplot();
  248. }
  249.