home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #27 / NN_1992_27.iso / spool / comp / theory / cellaut / 544 < prev    next >
Encoding:
Text File  |  1992-11-21  |  14.2 KB  |  540 lines
  1. Newsgroups: comp.theory.cell-automata
  2. Path: sparky!uunet!zaphod.mps.ohio-state.edu!cis.ohio-state.edu!news.sei.cmu.edu!fs7.ece.cmu.edu!crabapple.srv.cs.cmu.edu!hopkins
  3. From: hopkins@cs.cmu.edu (Don Hopkins)
  4. Subject: fast cellular automata implementation
  5. In-Reply-To: eeklund@nyx.cs.du.edu's message of Wed, 18 Nov 92 19:17:44 GMT
  6. Message-ID: <By30pu.E9s.1@cs.cmu.edu>
  7. Originator: hopkins@ECP.GARNET.CS.CMU.EDU
  8. Sender: news@cs.cmu.edu (Usenet News System)
  9. Nntp-Posting-Host: ecp.garnet.cs.cmu.edu
  10. Organization: School of Computer Science, Carnegie Mellon University
  11. References: <1992Nov17.122445.14209@odin.diku.dk>
  12.     <1992Nov18.191744.16470@mnemosyne.cs.du.edu>
  13. Date: Sat, 21 Nov 1992 19:56:17 GMT
  14. Lines: 524
  15.  
  16. eeklund@nyx.cs.du.edu (Eivind Eklund) writes:
  17.  
  18.    [...] And this is also nice for doing FAST implementations - see if
  19.    you can find out how... (This might be common knowledge allready,
  20.    but I have a hack which I believe to be VERY fast - somebody tell
  21.    me how many cells/second they are able to do on what kinds of
  22.    computers 8-) I am able to do approx 1.3 million cells/second on a
  23.    1 MIPS Amiga. (7.14 Mhz 68000) 
  24.  
  25. It's important to reduce the number of times you access memory, and
  26. keep as much as you can in registers.  You could implement a pure
  27. counting rule pretty efficiently by moving a 3x3 window over the
  28. cells, and keeping a count of the 1's under it, incrementing your
  29. count by 0..3 depending on the leading edge and decrementing it by
  30. that much three steps later (when those cells fall off the trailing
  31. edge). 
  32.  
  33. Another way to make it go faster (that would work well with the above
  34. technique) is to eliminate the edge conditions in the loop, by making
  35. the temporary array 2x2 cells bigger than the destination array.
  36. Before applying the rule, copy the cells into the temporary array
  37. inset by 1,1 (so they're centered with a 1 cell perimeter all
  38. around), then set the cells around the perimeter to the corresponding
  39. wrapping edge (to form a topological torus), and loop over the inset
  40. cells to compute the output array.  Since the loop's inset by one, you
  41. can address all neighbors uniformly like "west = src[-1], south =
  42. src[srcline]", instead of "west = src[x==0 ? srcline-1 : -1]", etc. 
  43.  
  44. I used this technique in my cellular automata machine implementation,
  45. that runs on HyperLook under OpenWindows 3.0.  The SPARC binary (with
  46. the user interface that lets you draw with the mouse in the cells
  47. while it's running) is available via anonymous ftp from ftp.uu.net, in
  48. the "graphics/NeWS" directory (be sure to get the HyperLook runtime
  49. too).  On a Sparc 2 it's quite fluid, and makes a pretty effective
  50. psychodellic doodle pad.  (It also includes a free lava lamp.) 
  51.  
  52. Here's the source code for the guts of my CAM simulator.  
  53. I've removed the cryptic user interface stuff and distilled it down,
  54. so the code can be easily reused.  I generate the rule lookup tables
  55. using SPARC Forth, in a language mostly compatible with the one used
  56. in Margolis & Toffoli's Cellular Automata Machines.  The simulator has
  57. a few built-in neighborhoods that index into the lookup table, and
  58. also several built-in rules implemented directory.  You have to
  59. initialize the cells, rule, and lookup table yourself (although the
  60. build-in rules don't use the lookup table), and roll your own display
  61. and user interface (unless you just enjoy the intellectual
  62. satisfaction of knowing it's applying all those rules).
  63.  
  64. Some day when I have the time I will port this to X11, probably using
  65. the Tk toolkit and TCL as an extension language, so you can define
  66. rules in TCL instead of Forth.  (It doesn't matter how efficient the
  67. high level rule language is because it's just used to generate a
  68. lookup table which is used in the inner loop).  
  69.  
  70. My apologies for the gnarly cpp macros, but C is such a cheezy
  71. language, with a brain dead macro facility tacked on as an
  72. afterthought, that once you get past anything more complicated than
  73. "#define MAX_BUFFER_SIZE 32", you end up writing in two languages at
  74. once, interleaving them with backslashes, and omitting semicolons in
  75. strategic places.  Note that "#define foo (bar)" is fundamentally
  76. different than "#define foo(bar)".  The idea behind the design of the
  77. macros is that it should be possible to re-use and combine rule
  78. definitions in useful ways (like n_margolis_ph or n_eco).  
  79.  
  80.     -Don
  81.  
  82. /*
  83.  * Cellular Automata Machine Simulator.
  84.  * Copyright (C) 1992 by Don Hopkins.  All rights reserved. 
  85.  * This program is provided for unrestricted use, provided that this 
  86.  * copyright message is preserved. There is no warranty, and no author 
  87.  * or distributer accepts responsibility for any damage caused by this 
  88.  * program. 
  89.  */
  90.  
  91.  
  92. #include <sys/types.h>
  93.  
  94. u_char *ruletable;
  95. int ruletablesize;
  96. u_char *src, *dst;
  97. u_char srcline, dstline;
  98. int width, height;
  99. u_char phase = 0;
  100. u_char wrap = 3;
  101. int steps = 1;
  102. int flow;
  103. void (*rule)();
  104.  
  105.  
  106. /* Generate a fast Cellular Automata Machine inner loop body */
  107. #define CAM_LOOP_BODY(BODY)                        \
  108.       { int y;                                \
  109.     for (y=0; y<height; y++) {                    \
  110.       long l0 = (src[0]<<8) + src[1],                \
  111.            l1 = (src[srcline]<<8) + src[srcline+1],            \
  112.            l2 = (src[srcline+srcline]<<8) + src[srcline+srcline+1];    \
  113.       int x;                            \
  114.       for (x=0; x<width; x++) {                    \
  115.         l0 = (l0<<8) + src[2];                    \
  116.         l1 = (l1<<8) + src[srcline+2];                \
  117.         l2 = (l2<<8) + src[srcline+srcline+2];            \
  118.         BODY;                            \
  119.         src++; dst++;                        \
  120.       }                                \
  121.       src += srcline - width; dst += dstline - width;        \
  122.     }                                \
  123.       }
  124.  
  125. /* Apply the rule and output the result */
  126. #define CAM_LOOP(RULE)                            \
  127.     CAM_LOOP_BODY(*dst = (RULE))
  128.  
  129. /* Rule defined by lookup table indexed by some neighborhood */
  130. #define CAM_TABLE_LOOP(NEIGHBORHOOD)                    \
  131.       CAM_LOOP(ruletable[(NEIGHBORHOOD)])
  132.  
  133. /* Extract neighboring cells from registers */
  134. #define NORTHWEST    ((u_char)((l0>>16) & 0xff))
  135. #define NORTH        ((u_char)((l0>>8) & 0xff))
  136. #define NORTHEAST    ((u_char)(l0 & 0xff))
  137. #define WEST        ((u_char)((l1>>16) & 0xff))
  138. #define CENTER        ((u_char)((l1>>8) & 0xff))
  139. #define EAST        ((u_char)(l1 & 0xff))
  140. #define SOUTHWEST    ((u_char)((l2>>16) & 0xff))
  141. #define SOUTH        ((u_char)((l2>>8) & 0xff))
  142. #define SOUTHEAST    ((u_char)(l2 & 0xff))
  143.  
  144. /* Eight neighbor sum of plane p (0..7) */
  145. #define SUM8p(p)    (((l0>>p)&1) + ((l0>>(p+8))&1) + ((l0>>(p+16))&1) + \
  146.              ((l1>>p)&1) +              ((l1>>(p+16))&1) + \
  147.              ((l2>>p)&1) + ((l2>>(p+8))&1) + ((l2>>(p+16))&1))
  148.  
  149. /* Nine neighbor sum of plane p (0..7) */
  150. #define SUM9p(p)    (SUM8p(p) + ((l1>>(p+8))&1))
  151.  
  152. /* Eight and nine neighbor sum of plane 0 */
  153. #define SUM8        SUM8p(0)
  154. #define SUM9        SUM9p(0)
  155.  
  156.  
  157. init_cam(int w, int h)
  158. {
  159.   void n_dheat();
  160.  
  161.   width = w; height = h;
  162.  
  163.   /* allocate output matrix of size width x height */
  164.   dstline = width;
  165.   dst = (u_char *)malloc(dstline * height);
  166.  
  167.   /* allocate temporary matrix of size width+2 x height+2 */
  168.   srcline = width + 2;
  169.   src = (u_char *)malloc(srcline * (height + 2) + (width + 2));
  170.  
  171.   rule = n_dheat;
  172. }
  173.  
  174.  
  175. do_rule()
  176. {
  177.     int step;
  178.  
  179.     for (step=0; step < steps; step++) {
  180.     int x, y;
  181.     u_char *s = src + srcline + 1, /* offset by 1,1 */
  182.            *d = dst;
  183.  
  184. /*
  185.  * Before applying the rule, we copy from the cells back from the dst 
  186.  * matrix to the temporary src matrix.  The temporary matrix is 2x2 larger 
  187.  * than the actual matrix we're computing, which is copied into the center, 
  188.  * and the wrapping edges copied around the perimeter, to eliminate edge 
  189.  * condition tests from the inner loop.
  190.  *
  191.  *    ff    f0 f1 ... fe ff        f0
  192.  *
  193.  *    0f    00 01 ... 0e 0f        00
  194.  *    1f    10 11 ... 1e 1f        10
  195.  *    ..    .. ..     .. ..        ..
  196.  *    ef    e0 e1 ... ee ef        e0
  197.  *    ff    f0 f1 ... fe ff        f0
  198.  *
  199.  *    0f    00 01 ... 0e 0f        00
  200.  *
  201.  * wrap value:    effect:
  202.  *    0    no effect
  203.  *    1    copy cells, no wrap
  204.  *    2    no copy, wrap edges
  205.  *    3    copy cells, wrap edges
  206.  *    4    copy cells, same edges
  207.  */
  208.  
  209.     switch (wrap) {
  210.     case 0:
  211.       break;
  212.     case 1:
  213.       for (y=0; y<height; y++) {
  214.         bcopy(d, s, width);
  215.         s += srcline;
  216.         d += dstline;
  217.       }
  218.       break;
  219.     case 2:
  220.       for (y=0; y<height; y++) {
  221.         s[-1] = s[width-1];
  222.         s[width] = s[0];
  223.         s += srcline;
  224.         d += dstline;
  225.       }
  226.       bcopy(src + srcline*height, src, srcline);
  227.       bcopy(src + srcline, src + srcline*(height+1), srcline);
  228.       break;
  229.     case 3:
  230.       for (y=0; y<height; y++) {
  231.         bcopy(d, s, width);
  232.         s[-1] = s[width-1];
  233.         s[width] = s[0];
  234.         s += srcline;
  235.         d += dstline;
  236.       }
  237.       bcopy(src + srcline*height, src, srcline);
  238.       bcopy(src + srcline, src + srcline*(height+1), srcline);
  239.       break;
  240.     case 4:
  241.       for (y=0; y<height; y++) {
  242.         bcopy(d, s, width);
  243.         s[-1] = s[0];
  244.         s[width] =  s[width-1];
  245.         s += srcline;
  246.         d += dstline;
  247.       }
  248.       bcopy(src + srcline*height, src + (srcline * (height + 1)), srcline);
  249.       bcopy(src + srcline, src, srcline);
  250.       break;
  251.     }
  252.  
  253.     (*rule)();
  254.  
  255.     phase = !phase;
  256.  
  257.     } /* for step */
  258. }
  259.  
  260.  
  261. void
  262. n_moore_a()
  263. {
  264.     /* 0    1    2    3    4    5    6    7    8     9     */
  265.     /* c    c'   se   sw   ne   nw   e    w    s     n     */
  266.     /* 0x1  0x2  0x4  0x8  0x10 0x20 0x40 0x80 0x100 0x200 */
  267.  
  268. #define MOORE_A (                            \
  269.     ((NORTHWEST&1)<<5) |    ((NORTH&1)<<9) |((NORTHEAST&1)<<4) |    \
  270.     ((WEST&1)<<7) |        (CENTER&3) |    ((EAST&1)<<6) |        \
  271.     ((SOUTHWEST&1)<<3) |    ((SOUTH&1)<<8) |((SOUTHEAST&1)<<2)    \
  272.     )
  273.  
  274.     CAM_TABLE_LOOP(MOORE_A)
  275. }
  276.  
  277. void
  278. n_moore_ab()
  279. {
  280.     /* 0    1    2    3    4    5    6    7    8     9     10    11    */
  281.     /* c    c'   se   sw   ne   nw   e    w    s     n     &c    &c'   */
  282.     /* 0x1  0x2  0x4  0x8  0x10 0x20 0x40 0x80 0x100 0x200 0x400 0x800 */
  283.  
  284. #define MOORE_AB (MOORE_A | ((CENTER&12)<<8))
  285.  
  286.     CAM_TABLE_LOOP(MOORE_AB)
  287. }
  288.  
  289. void
  290. n_vonn_neumann()
  291. {
  292.     /* 0    1    2    3    4    5    6    7    8     9     */
  293.     /* c    c'   e'   w'   s'   n'   e    w    s     n     */
  294.     /* 0x1  0x2  0x4  0x8  0x10 0x20 0x40 0x80 0x100 0x200 */
  295.  
  296. #define VON_NEUMANN (                            \
  297.     (CENTER&3) |                            \
  298.      ((EAST&1)<<6) | ((EAST&2)<<1) |                    \
  299.     ((WEST&1)<<7) | ((WEST&2)<<2) |                    \
  300.     ((SOUTH&1)<<8) | ((SOUTH&2)<<3) |                \
  301.     ((NORTH&1)<<9) | ((NORTH&2)<<4)                    \
  302.     )
  303.  
  304.     CAM_TABLE_LOOP(VON_NEUMANN)
  305. }
  306.  
  307. void
  308. n_margolis()
  309. {
  310.     register u_char i;
  311.  
  312.     /* 0    1    2    3    4    5    6    7    8     9     */
  313.     /* c    c'   cw   ccw  opp  cw'  ccw' opp'             */
  314.     /* 0x1  0x2  0x4  0x8  0x10 0x20 0x40 0x80 0x100 0x200 */
  315.  
  316. #define MARGOLIS_ODD (                            \
  317.     (CENTER & 3) |                            \
  318.     (i=(x&1 ? (y&1 ? (EAST) : (NORTH))                \
  319.         : (y&1 ? (SOUTH) : (WEST))),                \
  320.      (((i&1)<<2) | ((i&2)<<4))) |                    \
  321.     (i=(x&1 ? (y&1 ? (SOUTH) : (EAST))                \
  322.         : (y&1 ? (WEST) : (NORTH))),                \
  323.      (((i&1)<<3) | ((i&2)<<5))) |                    \
  324.     (i=(x&1 ? (y&1 ? (SOUTHEAST):(NORTHEAST))            \
  325.         : (y&1 ? (SOUTHWEST):(NORTHWEST))),            \
  326.      (((i&1)<<4) | ((i&2)<<6)))                    \
  327.     )
  328.  
  329. #define MARGOLIS_EVEN (                            \
  330.     (CENTER & 3) |                            \
  331.     (i=(x&1 ? (y&1 ? (WEST) : (SOUTH))                \
  332.         : (y&1 ? (NORTH) : (EAST))),                \
  333.      (((i&1)<<2) | ((i&2)<<4))) |                    \
  334.     (i=(x&1 ? (y&1 ? (NORTH) : (WEST))                \
  335.         : (y&1 ? (EAST) : (SOUTH))),                \
  336.      (((i&1)<<3) | ((i&2)<<5))) |                    \
  337.     (i=(x&1 ? (y&1 ? (NORTHWEST) : (SOUTHWEST))            \
  338.         : (y&1 ? (NORTHEAST) : (SOUTHEAST))),            \
  339.      (((i&1)<<4) | ((i&2)<<6)))                    \
  340.     )
  341.  
  342.     if (phase) {
  343.         CAM_TABLE_LOOP(MARGOLIS_ODD)
  344.     } else {
  345.         CAM_TABLE_LOOP(MARGOLIS_EVEN)
  346.     }
  347. }
  348.  
  349. void
  350. n_margolis_ph()
  351. {
  352.     register u_char i;
  353.  
  354.     /* 0    1    2    3    4    5    6    7    8    9      */
  355.     /* c    c'   cw   ccw  opp  cw'  ccw' opp' pha   pha'  */
  356.     /* 0x1  0x2  0x4  0x8  0x10 0x20 0x40 0x80 0x100 0x200 */
  357.  
  358. #define MARGOLIS_ODD_PH (MARGOLIS_ODD | 0x100)
  359. #define MARGOLIS_EVEN_PH (MARGOLIS_EVEN | 0x200)
  360.  
  361.     if (phase) {
  362.     CAM_TABLE_LOOP(MARGOLIS_ODD_PH)
  363.     } else {
  364.     CAM_TABLE_LOOP(MARGOLIS_EVEN_PH)
  365.     }
  366. }
  367.  
  368. void
  369. n_margolis_hv()
  370. {
  371.     register u_char i;
  372.  
  373.     /* 0    1    2    3    4    5    6    7    8    9      */
  374.     /* c    c'   cw   ccw  opp  cw'  ccw' opp' horz  vert  */
  375.     /* 0x1  0x2  0x4  0x8  0x10 0x20 0x40 0x80 0x100 0x200 */
  376.  
  377. #define MARGOLIS_ODD_HV (MARGOLIS_ODD | ((x&1)<<8) | ((y&1)<<9))
  378. #define MARGOLIS_EVEN_HV (MARGOLIS_EVEN | ((x&1)<<8) | ((y&1)<<9))
  379.  
  380.     if (phase) {
  381.     CAM_TABLE_LOOP(MARGOLIS_ODD_HV)
  382.     } else {
  383.     CAM_TABLE_LOOP(MARGOLIS_EVEN_HV)
  384.     }
  385. }
  386.  
  387. void
  388. n_life()
  389. {
  390.   int s;
  391.  
  392. #define LIFE (                            \
  393.      ((CENTER&1) ? (((s = SUM8) == 2) || (s == 3))        \
  394.               : (SUM8 == 3)) |                \
  395.      (CENTER<<1)                        \
  396.     )
  397.  
  398.   CAM_LOOP(LIFE)
  399. }
  400.  
  401. void
  402. n_brain()
  403. {
  404.   int s;
  405.  
  406. #define BRAIN (                            \
  407.      (((((s = CENTER)&3) == 0) && (SUM8 == 2)) ? 1 : 0) |    \
  408.      (s<<1)                            \
  409.     )
  410.  
  411.   CAM_LOOP(BRAIN)
  412. }
  413.  
  414. void
  415. n_heat()
  416. {
  417.  
  418. #define HEAT (                                \
  419.     ((long)(NORTHWEST + NORTH + NORTHEAST +                \
  420.         WEST +             EAST +                \
  421.         SOUTHWEST + SOUTH + SOUTHEAST + flow)) >> 3        \
  422.     )
  423.  
  424.     CAM_LOOP(HEAT)
  425. }
  426.  
  427. void
  428. n_dheat()
  429. {
  430.   int last = 0;
  431.  
  432. #define DHEAT (                                \
  433.     ((last = (long)(NORTHWEST + NORTH + NORTHEAST +            \
  434.             WEST +             EAST +            \
  435.             SOUTHWEST + SOUTH + SOUTHEAST + flow        \
  436.             + (last&0x07))), last >> 3)            \
  437.     )
  438.  
  439.     CAM_LOOP(DHEAT)
  440. }
  441.  
  442. void
  443. n_lheat()
  444. {
  445.  
  446. #define LHEAT (                                \
  447.     ((long)(NORTH + WEST + EAST + SOUTH + flow)) >> 2        \
  448.     )
  449.  
  450.     CAM_LOOP(LHEAT)
  451. }
  452.  
  453. void
  454. n_ldheat()
  455. {
  456.   int last = 0;
  457.  
  458. #define LDHEAT (                            \
  459.     ((last = (long)(NORTH + WEST + EAST + SOUTH + flow        \
  460.             + (last&0x03))), last >> 2)            \
  461.   )
  462.  
  463.   CAM_LOOP(LDHEAT)
  464. }
  465.  
  466. void
  467. n_anneal()
  468. {
  469.   int s;
  470.  
  471. #define ANNEAL (                            \
  472.       ((s = SUM9) > 5) || (s == 4)                    \
  473.     )
  474.   CAM_LOOP(ANNEAL)
  475. }
  476.  
  477. void
  478. n_anneal4()
  479. {
  480.   int s;
  481.  
  482. #define ANNEAL4 (                            \
  483.       ((((s = SUM9p(0)) > 5) || (s == 4)) ? 1 : 0) |        \
  484.       ((((s = SUM9p(1)) > 5) || (s == 4)) ? 2 : 0) |        \
  485.       ((((s = SUM9p(2)) > 5) || (s == 4)) ? 4 : 0) |        \
  486.       ((((s = SUM9p(3)) > 5) || (s == 4)) ? 8 : 0) |        \
  487.       (CENTER << 4)                            \
  488.     )
  489.   CAM_LOOP(ANNEAL4)
  490. }
  491.  
  492. void
  493. n_anneal8()
  494. {
  495.   int s;
  496.  
  497. #define ANNEAL8 (                            \
  498.       ((((s = SUM9p(0)) > 5) || (s == 4)) ? 1 : 0) |        \
  499.       ((((s = SUM9p(1)) > 5) || (s == 4)) ? 2 : 0) |        \
  500.       ((((s = SUM9p(2)) > 5) || (s == 4)) ? 4 : 0) |        \
  501.       ((((s = SUM9p(3)) > 5) || (s == 4)) ? 8 : 0) |        \
  502.       ((((s = SUM9p(4)) > 5) || (s == 4)) ? 16 : 0) |        \
  503.       ((((s = SUM9p(5)) > 5) || (s == 4)) ? 32 : 0) |        \
  504.       ((((s = SUM9p(6)) > 5) || (s == 4)) ? 64 : 0) |        \
  505.       ((((s = SUM9p(7)) > 5) || (s == 4)) ? 128 : 0)        \
  506.     )
  507.   CAM_LOOP(ANNEAL8)
  508. }
  509.  
  510. void
  511. n_eco()
  512. {
  513.   int s;
  514.  
  515. #define ANTILIFE (                            \
  516.      ((CENTER&1) ? (SUM8 != 5)                    \
  517.               : (((s = SUM8) != 5) && (s != 6))) |        \
  518.      (CENTER<<1)                            \
  519.     )
  520.  
  521. #define ECO (                                \
  522.       (((s = SUM9p(7)) > 5) || (s == 4) ? 128 : 0) |        \
  523.       ((CENTER&128) ? ((ANTILIFE)&127) : ((BRAIN)&127))        \
  524.     )
  525.   CAM_LOOP(ECO)
  526. }
  527.  
  528. void
  529. n_torben()
  530. {
  531.   int s;
  532.  
  533. /* 0 0 0 1 0 1 0 1 1 1 */
  534.  
  535. #define TORBEN (                            \
  536.       ((s = SUM9) > 6) || (s == 5) || (s == 3)            \
  537.     )
  538.   CAM_LOOP(TORBEN)
  539. }
  540.