home *** CD-ROM | disk | FTP | other *** search
/ Media Share 9 / MEDIASHARE_09.ISO / progmisc / frasrc18.zip / FRACSUBR.C < prev    next >
C/C++ Source or Header  |  1993-05-08  |  36KB  |  1,146 lines
  1. /*
  2. FRACSUBR.C contains subroutines which belong primarily to CALCFRAC.C and
  3. FRACTALS.C, i.e. which are non-fractal-specific fractal engine subroutines.
  4. */
  5.  
  6. #include <stdio.h>
  7. #ifndef XFRACT
  8. #include <stdarg.h>
  9. #else
  10. #include <varargs.h>
  11. #endif
  12. #include <float.h>
  13. #include <sys/types.h>
  14. #include <sys/timeb.h>
  15. #include <stdlib.h>
  16. #include "fractint.h"
  17. #include "fractype.h"
  18. #include "mpmath.h"
  19. #include "prototyp.h"
  20.  
  21. /* routines in this module    */
  22.  
  23. static long   _fastcall fudgetolong(double d);
  24. static double _fastcall fudgetodouble(long l);
  25. static void   _fastcall adjust_to_limits(double);
  26. static void   _fastcall smallest_add(double *);
  27. static int    _fastcall ratio_bad(double,double);
  28. static void   _fastcall plotdorbit(double,double,int);
  29. static int    _fastcall combine_worklist(void);
  30.  
  31.  
  32. extern int    calc_status;        /* status of calculations */
  33. extern char far *resume_info;        /* pointer to resume info if allocated */
  34.        int    resume_len;        /* length of resume info */
  35. static int    resume_offset;        /* offset in resume info gets */
  36. extern double plotmx1,plotmx2,plotmy1,plotmy2; /* real->screen conversion */
  37. extern int    orbit_ptr;        /* pointer into save_orbit array */
  38. extern int orbit_delay;         /* orbit delay value */
  39. extern int far *save_orbit;        /* array to save orbit values */
  40. extern int    orbit_color;        /* XOR color */
  41. extern int    num_worklist;        /* resume worklist for standard engine */
  42. extern int    fractype;         /* fractal type */
  43. extern char   stdcalcmode;        /* '1', '2', 'g', 'b'       */
  44. extern char   floatflag;        /* floating-point fractals? */
  45. extern int    integerfractal;        /* TRUE if fractal uses integer math */
  46. extern struct workliststuff worklist[MAXCALCWORK];
  47. extern int    sxdots,sydots;        /* # of dots on the physical screen    */
  48. extern int    sxoffs,syoffs;        /* physical top left of logical screen */
  49. extern int    xdots, ydots;        /* # of dots on the logical screen     */
  50. extern int    colors;            /* maximum colors available */
  51. extern long   fudge;            /* fudge factor (2**n) */
  52. extern int    bitshift;         /* bit shift for fudge */
  53. extern int    inside;            /* inside color: 1=blue     */
  54. extern int    outside;            /* outside color    */
  55. extern double xxmin,xxmax,yymin,yymax,xx3rd,yy3rd; /* corners */
  56. extern long   xmin, xmax, ymin, ymax, x3rd, y3rd;  /* integer equivs */
  57. extern int    maxit;            /* try this many iterations */
  58. extern int    attractors;        /* number of finite attractors    */
  59. extern _CMPLX  attr[];        /* finite attractor vals (f.p)    */
  60. extern _LCMPLX lattr[];     /* finite attractor vals (int)    */
  61. extern int    attrperiod[];         /* finite attractor period    */
  62. extern _CMPLX  old,new;
  63. extern _LCMPLX lold,lnew;
  64. extern double tempsqrx,tempsqry;
  65. extern long   ltempsqrx,ltempsqry;
  66. extern int    xxstart,xxstop;        /* these are same as worklist, */
  67. extern int    yystart,yystop,yybegin;    /* declared as separate items  */
  68. extern int    periodicitycheck;
  69. extern int    basin;
  70. extern int    finattract;
  71. extern int    pixelpi;            /* value of pi in pixels */
  72. extern double closenuff;
  73. extern long   lclosenuff;
  74. extern int    ixstart, ixstop, iystart, iystop;
  75. extern int    color;
  76. extern int    decomp[];
  77. extern double potparam[3];        /* three potential parameters*/
  78. extern int    distest;            /* non-zero if distance estimator */
  79. extern double param[];          /* parameters */
  80. extern int    invert;            /* non-zero if inversion active */
  81. extern int    biomorph,usr_biomorph;
  82. extern int    debugflag; /* internal use only - you didn't see this */
  83. extern long   creal, cimag;        /* for calcmand */
  84. extern long   delx, dely;        /* screen pixel increments  */
  85. extern long   delx2, dely2;        /* screen pixel increments  */
  86. extern double delxx, delyy;        /* screen pixel increments  */
  87. extern double delxx2, delyy2;        /* screen pixel increments  */
  88. extern long   delmin;            /* for calcfrac/calcmand    */
  89. extern double ddelmin;            /* same as a double        */
  90. extern int    potflag;            /* continuous potential flag */
  91. extern int    bailout;
  92. extern double rqlim;
  93. extern double  dxsize, dysize;        /* xdots-1, ydots-1        */
  94. extern int soundflag;
  95. extern int basehertz;
  96.  
  97. extern long   far *lx0, far *ly0;    /* x, y grid            */
  98. extern long   far *lx1, far *ly1;    /* adjustment for rotate    */
  99. /* note that lx1 & ly1 values can overflow into sign bit; since     */
  100. /* they're used only to add to lx0/ly0, 2s comp straightens it out  */
  101. extern double far *dx0, far *dy0;    /* floating pt equivs */
  102. extern double far *dx1, far *dy1;
  103.  
  104. extern char   usr_floatflag;
  105. extern char   usr_stdcalcmode;
  106. extern int    usr_periodicitycheck;
  107. extern int    usr_distest;
  108. extern double zwidth;
  109.  
  110. extern int neworbittype;
  111.  
  112. #define FUDGEFACTOR    29    /* fudge all values up by 2**this */
  113. #define FUDGEFACTOR2    24    /* (or maybe this)          */
  114.  
  115.  
  116. void calcfracinit() /* initialize a *pile* of stuff for fractal calculation */
  117. {
  118.    int i;
  119.    double ftemp;
  120.  
  121.    floatflag = usr_floatflag;
  122.    if (calc_status == 2) /* on resume, ensure floatflag correct */
  123.       if (curfractalspecific->isinteger)
  124.      floatflag = 0;
  125.       else
  126.      floatflag = 1;
  127.    /* if floating pt only, set floatflag for TAB screen */
  128.    if (!curfractalspecific->isinteger && curfractalspecific->tofloat == NOFRACTAL)
  129.       floatflag = 1;
  130.  
  131. init_restart:
  132.  
  133.    /* the following variables may be forced to a different setting due to
  134.       calc routine constraints;  usr_xxx is what the user last said is wanted,
  135.       xxx is what we actually do in the current situation */
  136.    stdcalcmode        = usr_stdcalcmode;
  137.    periodicitycheck = usr_periodicitycheck;
  138.    distest        = usr_distest;
  139.    biomorph        = usr_biomorph;
  140.  
  141.    potflag = 0;
  142.    if (potparam[0] != 0.0
  143.      && colors >= 64
  144.      && (curfractalspecific->calctype == StandardFractal
  145.      || curfractalspecific->calctype == calcmand
  146.      || curfractalspecific->calctype == calcmandfp)) {
  147.       potflag = 1;
  148.       distest = 0;    /* can't do distest too */
  149.       }
  150.  
  151.    if (distest)
  152.       floatflag = 1;  /* force floating point for dist est */
  153.  
  154.    if (floatflag) { /* ensure type matches floatflag */
  155.       if (curfractalspecific->isinteger != 0
  156.     && curfractalspecific->tofloat != NOFRACTAL)
  157.      fractype = curfractalspecific->tofloat;
  158.       }
  159.    else {
  160.       if (curfractalspecific->isinteger == 0
  161.     && curfractalspecific->tofloat != NOFRACTAL)
  162.      fractype = curfractalspecific->tofloat;
  163.       }
  164.    /* match Julibrot with integer mode of orbit */   
  165.    if(fractype == JULIBROTFP && fractalspecific[neworbittype].isinteger)
  166.    {
  167.       int i;
  168.       if((i=fractalspecific[neworbittype].tofloat) != NOFRACTAL)
  169.          neworbittype = i;
  170.       else
  171.          fractype = JULIBROT;   
  172.    }
  173.    else if(fractype == JULIBROT && fractalspecific[neworbittype].isinteger==0)
  174.    {
  175.       int i;
  176.       if((i=fractalspecific[neworbittype].tofloat) != NOFRACTAL)
  177.          neworbittype = i;
  178.       else
  179.          fractype = JULIBROTFP;   
  180.    }
  181.       
  182.    curfractalspecific = &fractalspecific[fractype];
  183.  
  184.    integerfractal = curfractalspecific->isinteger;
  185.  
  186. /*   if (fractype == JULIBROT)
  187.       rqlim = 4;
  188.    else */ if (potflag && potparam[2] != 0.0)
  189.       rqlim = potparam[2];
  190. /* else if (decomp[0] > 0 && decomp[1] > 0)
  191.       rqlim = (double)decomp[1]; */
  192.    else if (bailout) /* user input bailout */
  193.       rqlim = bailout;
  194.    else if (biomorph != -1) /* biomorph benefits from larger bailout */
  195.       rqlim = 100;
  196.    else
  197.       rqlim = curfractalspecific->orbit_bailout;
  198.    if (integerfractal) /* the bailout limit mustn't be too high here */
  199.       if (rqlim > 127.0) rqlim = 127.0;
  200.  
  201.    if ((curfractalspecific->flags&NOROTATE) != 0) {
  202.       /* ensure min<max and unrotated rectangle */
  203.       if (xxmin > xxmax) { ftemp = xxmax; xxmax = xxmin; xxmin = ftemp; }
  204.       if (yymin > yymax) { ftemp = yymax; yymax = yymin; yymin = ftemp; }
  205.       xx3rd = xxmin; yy3rd = yymin;
  206.       }
  207.  
  208.    /* set up bitshift for integer math */
  209.    bitshift = FUDGEFACTOR2; /* by default, the smaller shift */
  210.    if (integerfractal > 1)  /* use specific override from table */
  211.       bitshift = integerfractal;
  212.    if (integerfractal == 0) /* float? */
  213.       if ((i = curfractalspecific->tofloat) != NOFRACTAL) /* -> int? */
  214.      if (fractalspecific[i].isinteger > 1) /* specific shift? */
  215.         bitshift = fractalspecific[i].isinteger;
  216.  
  217. /* We want this code if we're using the assembler calcmand */
  218.    if (fractype == MANDEL || fractype == JULIA) { /* adust shift bits if.. */
  219.       if (potflag == 0                  /* not using potential */
  220.     && (param[0] > -2.0 && param[0] < 2.0)  /* parameters not too large */
  221.     && (param[1] > -2.0 && param[1] < 2.0)
  222.     && !invert                  /* and not inverting */
  223.     && biomorph == -1              /* and not biomorphing */
  224.     && rqlim <= 4.0               /* and bailout not too high */
  225.     && (outside > -2 || outside < -5)      /* and no funny outside stuff */
  226.     && debugflag != 1234)              /* and not debugging */
  227.      bitshift = FUDGEFACTOR;          /* use the larger bitshift */
  228.       }
  229.  
  230.    fudge = 1L << bitshift;
  231.  
  232.    l_at_rad = fudge/32768L;
  233.    f_at_rad = 1.0/32768L;
  234.  
  235.    /* now setup arrays of real coordinates corresponding to each pixel */
  236.  
  237.    adjust_to_limits(1.0); /* make sure all corners in valid range */
  238.  
  239.    delxx  = (xxmax - xx3rd) / dxsize; /* calculate stepsizes */
  240.    delyy  = (yymax - yy3rd) / dysize;
  241.    delxx2 = (xx3rd - xxmin) / dysize;
  242.    delyy2 = (yy3rd - yymin) / dxsize;
  243.  
  244.   if(fractype != CELLULAR) { /* fudgetolong fails w >10 digits in double */
  245.    creal = fudgetolong(param[0]); /* integer equivs for it all */
  246.    cimag = fudgetolong(param[1]);
  247.    xmin  = fudgetolong(xxmin);
  248.    xmax  = fudgetolong(xxmax);
  249.    x3rd  = fudgetolong(xx3rd);
  250.    ymin  = fudgetolong(yymin);
  251.    ymax  = fudgetolong(yymax);
  252.    y3rd  = fudgetolong(yy3rd);
  253.    delx  = fudgetolong(delxx);
  254.    dely  = fudgetolong(delyy);
  255.    delx2 = fudgetolong(delxx2);
  256.    dely2 = fudgetolong(delyy2);
  257.   }
  258.  
  259.    if (fractype != PLASMA) { /* skip this if plasma to avoid 3d problems */
  260.       if (integerfractal && !invert) {
  261.      if (    (delx  == 0 && delxx  != 0.0)
  262.          || (delx2 == 0 && delxx2 != 0.0)
  263.          || (dely  == 0 && delyy  != 0.0)
  264.          || (dely2 == 0 && delyy2 != 0.0) )
  265.         goto expand_retry;
  266.      lx0[0] = xmin;         /* fill up the x, y grids */
  267.      ly0[0] = ymax;
  268.      lx1[0] = ly1[0] = 0;
  269.      for (i = 1; i < xdots; i++ ) {
  270.         lx0[i] = lx0[i-1] + delx;
  271.         ly1[i] = ly1[i-1] - dely2;
  272.         }
  273.      for (i = 1; i < ydots; i++ ) {
  274.         ly0[i] = ly0[i-1] - dely;
  275.         lx1[i] = lx1[i-1] + delx2;
  276.         }
  277.      /* past max res?  check corners within 10% of expected */
  278.      if (    ratio_bad((double)lx0[xdots-1]-xmin,(double)xmax-x3rd)
  279.          || ratio_bad((double)ly0[ydots-1]-ymax,(double)y3rd-ymax)
  280.          || ratio_bad((double)lx1[(ydots>>1)-1],((double)x3rd-xmin)/2)
  281.          || ratio_bad((double)ly1[(xdots>>1)-1],((double)ymin-y3rd)/2) ) {
  282. expand_retry:
  283.         if (integerfractal        /* integer fractal type? */
  284.           && curfractalspecific->tofloat != NOFRACTAL)
  285.            floatflag = 1;        /* switch to floating pt */
  286.         else
  287.            adjust_to_limits(2.0);    /* double the size */
  288.         if (calc_status == 2)    /* due to restore of an old file? */
  289.            calc_status = 0;     /*   whatever, it isn't resumable */
  290.         goto init_restart;
  291.         }
  292.      /* re-set corners to match reality */
  293.      xmax = lx0[xdots-1] + lx1[ydots-1];
  294.      ymin = ly0[ydots-1] + ly1[xdots-1];
  295.      x3rd = xmin + lx1[ydots-1];
  296.      y3rd = ly0[ydots-1];
  297.      xxmin = fudgetodouble(xmin);
  298.      xxmax = fudgetodouble(xmax);
  299.      xx3rd = fudgetodouble(x3rd);
  300.      yymin = fudgetodouble(ymin);
  301.      yymax = fudgetodouble(ymax);
  302.      yy3rd = fudgetodouble(y3rd);
  303.      }
  304.       else {
  305.      /* set up dx0 and dy0 analogs of lx0 and ly0 */
  306.      /* put fractal parameters in doubles */
  307.      dx0[0] = xxmin;        /* fill up the x, y grids */
  308.      dy0[0] = yymax;
  309.      dx1[0] = dy1[0] = 0;
  310.      for (i = 1; i < xdots; i++ ) {
  311.         dx0[i] = dx0[i-1] + delxx;
  312.         dy1[i] = dy1[i-1] - delyy2;
  313.         }
  314.      for (i = 1; i < ydots; i++ ) {
  315.         dy0[i] = dy0[i-1] - delyy;
  316.         dx1[i] = dx1[i-1] + delxx2;
  317.         }
  318.      if (    ratio_bad(dx0[xdots-1]-xxmin,xxmax-xx3rd)
  319.          || ratio_bad(dy0[ydots-1]-yymax,yy3rd-yymax)
  320.          || ratio_bad(dx1[ydots-1],xx3rd-xxmin)
  321.          || ratio_bad(dy1[xdots-1],yymin-yy3rd))
  322.         goto expand_retry;
  323.      /* re-set corners to match reality */
  324.      xxmax = dx0[xdots-1] + dx1[ydots-1];
  325.      yymin = dy0[ydots-1] + dy1[xdots-1];
  326.      xx3rd = xxmin + dx1[ydots-1];
  327.      yy3rd = dy0[ydots-1];
  328.      }
  329.       }
  330.  
  331.    /* for periodicity close-enough, and for unity: */
  332.    /*      min(max(delx,delx2),max(dely,dely2)       */
  333.    ddelmin = fabs(delxx);
  334.    if (fabs(delxx2) > ddelmin)
  335.       ddelmin = fabs(delxx2);
  336.    if (fabs(delyy) > fabs(delyy2)) {
  337.       if (fabs(delyy) < ddelmin)
  338.      ddelmin = fabs(delyy);
  339.       }
  340.    else
  341.       if (fabs(delyy2) < ddelmin)
  342.      ddelmin = fabs(delyy2);
  343.    delmin = fudgetolong(ddelmin);
  344.  
  345.    /* calculate factors which plot real values to screen co-ords */
  346.    /* calcfrac.c plot_orbit routines have comments about this     */
  347.    ftemp = (0.0-delyy2) * delxx2 * dxsize * dysize
  348.        - (xxmax-xx3rd) * (yy3rd-yymax);
  349.    plotmx1 = delxx2 * dxsize * dysize / ftemp;
  350.    plotmx2 = (yy3rd-yymax) * dxsize / ftemp;
  351.    plotmy1 = (0.0-delyy2) * dxsize * dysize / ftemp;
  352.    plotmy2 = (xxmax-xx3rd) * dysize / ftemp;
  353.  
  354. }
  355.  
  356. static long _fastcall fudgetolong(double d)
  357. {
  358.    if ((d *= fudge) > 0) d += 0.5;
  359.    else          d -= 0.5;
  360.    return (long)d;
  361. }
  362.  
  363. static double _fastcall fudgetodouble(long l)
  364. {
  365.    char buf[30];
  366.    double d;
  367.    sprintf(buf,"%.9g",(double)l / fudge);
  368. #ifndef XFRACT
  369.    sscanf(buf,"%lg",&d);
  370. #else
  371.    sscanf(buf,"%lf",&d);
  372. #endif
  373.    return d;
  374. }
  375.  
  376.  
  377. void adjust_corner()
  378. {
  379.    /* make edges very near vert/horiz exact, to ditch rounding errs and */
  380.    /* to avoid problems when delta per axis makes too large a ratio    */
  381.    double ftemp,ftemp2;
  382.    if( (ftemp=fabs(xx3rd-xxmin)) < (ftemp2=fabs(xxmax-xx3rd)) ) {
  383.       if (ftemp*10000 < ftemp2 && yy3rd != yymax)
  384.      xx3rd = xxmin;
  385.       }
  386.    else if (ftemp2*10000 < ftemp && yy3rd != yymin)
  387.       xx3rd = xxmax;
  388.    if( (ftemp=fabs(yy3rd-yymin)) < (ftemp2=fabs(yymax-yy3rd)) ) {
  389.       if (ftemp*10000 < ftemp2 && xx3rd != xxmax)
  390.      yy3rd = yymin;
  391.       }
  392.    else if (ftemp2*10000 < ftemp && xx3rd != xxmin)
  393.       yy3rd = yymax;
  394. }
  395.  
  396. static void _fastcall adjust_to_limits(double expand)
  397. {  double cornerx[4],cornery[4];
  398.    double lowx,highx,lowy,highy,limit,ftemp;
  399.    double centerx,centery,adjx,adjy;
  400.    int i;
  401.    limit = 32767.99;
  402.    if (bitshift >= 24) limit = 31.99;
  403.    if (bitshift >= 29) limit = 3.99;
  404.    centerx = (xxmin+xxmax)/2;
  405.    centery = (yymin+yymax)/2;
  406.    if (xxmin == centerx) { /* ohoh, infinitely thin, fix it */
  407.       smallest_add(&xxmax);
  408.       xxmin -= xxmax-centerx;
  409.       }
  410.    if (yymin == centery) {
  411.       smallest_add(&yymax);
  412.       yymin -= yymax-centery;
  413.       }
  414.    if (xx3rd == centerx)
  415.       smallest_add(&xx3rd);
  416.    if (yy3rd == centery)
  417.       smallest_add(&yy3rd);
  418.    /* setup array for easier manipulation */
  419.    cornerx[0] = xxmin; cornerx[1] = xxmax;
  420.    cornerx[2] = xx3rd; cornerx[3] = xxmin+(xxmax-xx3rd);
  421.    cornery[0] = yymax; cornery[1] = yymin;
  422.    cornery[2] = yy3rd; cornery[3] = yymin+(yymax-yy3rd);
  423.    /* if caller wants image size adjusted, do that first */
  424.    if (expand != 1.0)
  425.       for (i=0; i<4; ++i) {
  426.      cornerx[i] = centerx + (cornerx[i]-centerx)*expand;
  427.      cornery[i] = centery + (cornery[i]-centery)*expand;
  428.       }
  429.    /* get min/max x/y values */
  430.    lowx = highx = cornerx[0];
  431.    lowy = highy = cornery[0];
  432.    for (i=1; i<4; ++i) {
  433.       if (cornerx[i] < lowx)  lowx  = cornerx[i];
  434.       if (cornerx[i] > highx) highx = cornerx[i];
  435.       if (cornery[i] < lowy)  lowy  = cornery[i];
  436.       if (cornery[i] > highy) highy = cornery[i];
  437.       }
  438.    /* if image is too large, downsize it maintaining center */
  439.    ftemp = highx-lowx;
  440.    if (highy-lowy > ftemp) ftemp = highy-lowy;
  441.    if ((ftemp = limit*2/ftemp) < 1.0)
  442.       for (i=0; i<4; ++i) {
  443.      cornerx[i] = centerx + (cornerx[i]-centerx)*ftemp;
  444.      cornery[i] = centery + (cornery[i]-centery)*ftemp;
  445.      }
  446.    /* if any corner has x or y past limit, move the image */
  447.    adjx = adjy = 0;
  448.    for (i=0; i<4; ++i) {
  449.       if (cornerx[i] > limit     && (ftemp = cornerx[i] - limit) > adjx)
  450.      adjx = ftemp;
  451.       if (cornerx[i] < 0.0-limit && (ftemp = cornerx[i] + limit) < adjx)
  452.      adjx = ftemp;
  453.       if (cornery[i] > limit     && (ftemp = cornery[i] - limit) > adjy)
  454.      adjy = ftemp;
  455.       if (cornery[i] < 0.0-limit && (ftemp = cornery[i] + limit) < adjy)
  456.      adjy = ftemp;
  457.       }
  458.    if (calc_status == 2 && (adjx != 0 || adjy != 0) && (zwidth == 1.0))
  459.       calc_status = 0;
  460.    xxmin = cornerx[0] - adjx;
  461.    xxmax = cornerx[1] - adjx;
  462.    xx3rd = cornerx[2] - adjx;
  463.    yymax = cornery[0] - adjy;
  464.    yymin = cornery[1] - adjy;
  465.    yy3rd = cornery[2] - adjy;
  466.    adjust_corner(); /* make 3rd corner exact if very near other co-ords */
  467. }
  468.  
  469. static void _fastcall smallest_add(double *num)
  470. {
  471.    *num += *num * 5.0e-16;
  472. }
  473.  
  474. static int _fastcall ratio_bad(double actual, double desired)
  475. {  double ftemp;
  476.    if (desired != 0)
  477.       if ((ftemp = actual / desired) < 0.95 || ftemp > 1.05)
  478.      return(1);
  479.    return(0);
  480. }
  481.  
  482.  
  483. /* Save/resume stuff:
  484.  
  485.    Engines which aren't resumable can simply ignore all this.
  486.  
  487.    Before calling the (per_image,calctype) routines (engine), calcfract sets:
  488.       "resuming" to 0 if new image, nonzero if resuming a partially done image
  489.       "calc_status" to 1
  490.    If an engine is interrupted and wants to be able to resume it must:
  491.       store whatever status info it needs to be able to resume later
  492.       set calc_status to 2 and return
  493.    If subsequently called with resuming!=0, the engine must restore status
  494.    info and continue from where it left off.
  495.  
  496.    Since the info required for resume can get rather large for some types,
  497.    it is not stored directly in save_info.  Instead, memory is dynamically
  498.    allocated as required, and stored in .fra files as a separate block.
  499.    To save info for later resume, an engine routine can use:
  500.       alloc_resume(maxsize,version)
  501.      Maxsize must be >= max bytes subsequently saved + 2; over-allocation
  502.      is harmless except for possibility of insufficient mem available;
  503.      undersize is not checked and probably causes serious misbehaviour.
  504.      Version is an arbitrary number so that subsequent revisions of the
  505.      engine can be made backward compatible.
  506.      Alloc_resume sets calc_status to 2 (resumable) if it succeeds; to 3
  507.      if it cannot allocate memory (and issues warning to user).
  508.       put_resume({bytes,&argument,} ... 0)
  509.      Can be called as often as required to store the info.
  510.      Arguments must not be far addresses.
  511.      Is not protected against calls which use more space than allocated.
  512.    To reload info when resuming, use:
  513.       start_resume()
  514.      initializes and returns version number
  515.       get_resume({bytes,&argument,} ... 0)
  516.      inverse of store_resume
  517.       end_resume()
  518.      optional, frees the memory area sooner than would happen otherwise
  519.  
  520.    Example, save info:
  521.       alloc_resume(sizeof(parmarray)+100,2);
  522.       put_resume(sizeof(int),&row, sizeof(int),&col,
  523.          sizeof(parmarray),parmarray, 0);
  524.     restore info:
  525.       vsn = start_resume();
  526.       get_resume(sizeof(int),&row, sizeof(int),&col, 0);
  527.       if (vsn >= 2)
  528.      get_resume(sizeof(parmarray),parmarray,0);
  529.       end_resume();
  530.  
  531.    Engines which allocate a large far memory chunk of their own might
  532.    directly set resume_info, resume_len, calc_status to avoid doubling
  533.    transient memory needs by using these routines.
  534.  
  535.    StandardFractal, calcmand, solidguess, and bound_trace_main are a related
  536.    set of engines for escape-time fractals.  They use a common worklist
  537.    structure for save/resume.  Fractals using these must specify calcmand
  538.    or StandardFractal as the engine in fractalspecificinfo.
  539.    Other engines don't get btm nor ssg, don't get off-axis symmetry nor
  540.    panning (the worklist stuff), and are on their own for save/resume.
  541.  
  542.    */
  543.  
  544. #ifndef XFRACT
  545. int put_resume(int len, ...)
  546. #else
  547. int put_resume(va_alist)
  548. va_dcl
  549. #endif
  550. {
  551.    va_list arg_marker;    /* variable arg list */
  552.    char *source_ptr;
  553. #ifdef XFRACT
  554.    int len;
  555. #endif
  556.  
  557.    if (resume_info == NULL)
  558.       return(-1);
  559. #ifndef XFRACT
  560.    va_start(arg_marker,len);
  561. #else
  562.    va_start(arg_marker);
  563.    len = va_arg(arg_marker,int);
  564. #endif
  565.    while (len)
  566.    {
  567.       source_ptr = va_arg(arg_marker,char *);
  568.       far_memcpy(resume_info+resume_len,source_ptr,len);
  569.       resume_len += len;
  570.       len = va_arg(arg_marker,int);
  571.    }
  572.    return(0);
  573. }
  574.  
  575. int alloc_resume(int alloclen, int version)
  576. {
  577.    if (resume_info != NULL) /* free the prior area if there is one */
  578.       farmemfree(resume_info);
  579.    if ((resume_info = farmemalloc((long)alloclen))== NULL)
  580.    {
  581.       static char msg[] = {"\
  582. Warning - insufficient free memory to save status.\n\
  583. You will not be able to resume calculating this image."};
  584.       stopmsg(0,msg);
  585.       calc_status = 3;
  586.       return(-1);
  587.    }
  588.    resume_len = 0;
  589.    put_resume(sizeof(int),&version,0);
  590.    calc_status = 2;
  591.    return(0);
  592. }
  593.  
  594. #ifndef XFRACT
  595. int get_resume(int len, ...)
  596. #else
  597. int get_resume(va_alist)
  598. va_dcl
  599. #endif
  600. {
  601.    va_list arg_marker;    /* variable arg list */
  602.    char *dest_ptr;
  603. #ifdef XFRACT
  604.    int len;
  605. #endif
  606.  
  607.    if (resume_info == NULL)
  608.       return(-1);
  609. #ifndef XFRACT
  610.    va_start(arg_marker,len);
  611. #else
  612.    va_start(arg_marker);
  613.    len = va_arg(arg_marker,int);
  614. #endif
  615.    while (len)
  616.    {
  617.       dest_ptr = va_arg(arg_marker,char *);
  618.       far_memcpy(dest_ptr,resume_info+resume_offset,len);
  619.       resume_offset += len;
  620.       len = va_arg(arg_marker,int);
  621.    }
  622.    return(0);
  623. }
  624.  
  625. int start_resume()
  626. {
  627.    int version;
  628.    if (resume_info == NULL)
  629.       return(-1);
  630.    resume_offset = 0;
  631.    get_resume(sizeof(int),&version,0);
  632.    return(version);
  633. }
  634.  
  635. void end_resume()
  636. {
  637.    if (resume_info != NULL) /* free the prior area if there is one */
  638.    {
  639.       farmemfree(resume_info);
  640.       resume_info = NULL;
  641.    }
  642. }
  643.  
  644.  
  645. /* Showing orbit requires converting real co-ords to screen co-ords.
  646.    Define:
  647.        Xs == xxmax-xx3rd           Ys == yy3rd-yymax
  648.        W  == xdots-1               D  == ydots-1
  649.    We know that:
  650.        realx == lx0[col] + lx1[row]
  651.        realy == ly0[row] + ly1[col]
  652.        lx0[col] == (col/width) * Xs + xxmin
  653.        lx1[row] == row * delxx
  654.        ly0[row] == (row/D) * Ys + yymax
  655.        ly1[col] == col * (0-delyy)
  656.   so:
  657.        realx == (col/W) * Xs + xxmin + row * delxx
  658.        realy == (row/D) * Ys + yymax + col * (0-delyy)
  659.   and therefore:
  660.        row == (realx-xxmin - (col/W)*Xs) / Xv     (1)
  661.        col == (realy-yymax - (row/D)*Ys) / Yv     (2)
  662.   substitute (2) into (1) and solve for row:
  663.        row == ((realx-xxmin)*(0-delyy2)*W*D - (realy-yymax)*Xs*D)
  664.               / ((0-delyy2)*W*delxx2*D-Ys*Xs)
  665.   */
  666.  
  667. /* sleep N * a tenth of a millisecond */
  668.  
  669. void sleepms(long ms)
  670. {
  671.     extern int tabmode;
  672.     static long scalems=0L;
  673.     int savehelpmode,savetabmode;
  674.     struct timeb t1,t2;
  675. #define SLEEPINIT 250 /* milliseconds for calibration */
  676.     savetabmode  = tabmode;
  677.     savehelpmode = helpmode;
  678.     tabmode  = 0;
  679.     helpmode = -1;
  680.     if(scalems==0L) /* calibrate */
  681.     {
  682.     /* selects a value of scalems that makes the units
  683.        10000 per sec independent of CPU speed */
  684.     char msg[80];
  685.     long millisecs;
  686.     int i,elapsed;
  687.  
  688.     scalems = 1L;
  689.     if(keypressed()) /* check at start, hope to get start of timeslice */
  690.        goto sleepexit;
  691.     /* calibrate, assume slow computer first */
  692.     do
  693.     {
  694.        scalems *= 2;
  695.        ftime(&t2);
  696.        do { /* wait for the start of a new tick */
  697.           ftime(&t1);
  698.         }
  699.         while (t2.time == t1.time && t2.millitm == t1.millitm);
  700.        sleepms(10L * SLEEPINIT); /* about 1/4 sec */
  701.        ftime(&t2);
  702.        if(keypressed()) {
  703.           scalems = 0L;
  704.           goto sleepexit;
  705.        }
  706.      }
  707.      while ((elapsed = (int)(t2.time-t1.time)*1000 + t2.millitm-t1.millitm)
  708.         < SLEEPINIT);
  709.     /* once more to see if faster (eg multi-tasking) */
  710.     do { /* wait for the start of a new tick */
  711.        ftime(&t1);
  712.        }
  713.      while (t2.time == t1.time && t2.millitm == t1.millitm);
  714.     sleepms(10L * SLEEPINIT);
  715.     ftime(&t2);
  716.     if ((i = (int)(t2.time-t1.time)*1000 + t2.millitm-t1.millitm) < elapsed)
  717.        elapsed = (i == 0) ? 1 : i;
  718.     scalems = (float)SLEEPINIT/(float)(elapsed) * scalems;
  719.     if (debugflag == 700) {
  720.        sprintf(msg,"scale factor=%ld",scalems);
  721.        stopmsg(0,msg);
  722.     }
  723.     }
  724.     if(ms > 10L * SLEEPINIT) { /* using ftime is probably more accurate */
  725.     ms /= 10;
  726.     ftime(&t1);
  727.     while(1) {
  728.        if(keypressed()) break;
  729.        ftime(&t2);
  730.        if ((t2.time-t1.time)*1000 + t2.millitm-t1.millitm >= ms) break;
  731.     }
  732.     }
  733.     else
  734.     if(!keypressed()) {
  735.        ms *= scalems;
  736.        while(ms-- >= 0);
  737.     }
  738. sleepexit:
  739.     tabmode  = savetabmode;
  740.     helpmode = savehelpmode;
  741. }
  742.  
  743.  
  744. static void _fastcall plotdorbit(double dx, double dy, int color)
  745. {
  746.    int i, j, c;
  747.    int save_sxoffs,save_syoffs;
  748.    if (orbit_ptr >= 1500) return;
  749.    i = dy * plotmx1 - dx * plotmx2; i += sxoffs;
  750.    if (i < 0 || i >= sxdots) return;
  751.    j = dx * plotmy1 - dy * plotmy2; j += syoffs;
  752.    if (j < 0 || j >= sydots) return;
  753.    save_sxoffs = sxoffs;
  754.    save_syoffs = syoffs;
  755.    sxoffs = syoffs = 0;
  756.    /* save orbit value */
  757.    if(color == -1)
  758.    {
  759.       *(save_orbit + orbit_ptr++) = i;
  760.       *(save_orbit + orbit_ptr++) = j;
  761.       *(save_orbit + orbit_ptr++) = c = getcolor(i,j);
  762.       putcolor(i,j,c^orbit_color);
  763.    }
  764.    else
  765.       putcolor(i,j,color);
  766.    sxoffs = save_sxoffs;
  767.    syoffs = save_syoffs;
  768.    if(orbit_delay > 0)
  769.       sleepms(orbit_delay);
  770.    if(soundflag==1)
  771.        snd((int)(i*1000/xdots+basehertz));
  772.    else if(soundflag > 1)
  773.        snd((int)(j*1000/ydots+basehertz));
  774.    /* placing sleepms here delays each dot */
  775. }
  776.  
  777. void iplot_orbit(ix, iy, color)
  778. long ix, iy;
  779. int color;
  780. {
  781.    plotdorbit((double)ix/fudge-xxmin,(double)iy/fudge-yymax,color);
  782. }
  783.  
  784. void plot_orbit(real,imag,color)
  785. double real,imag;
  786. int color;
  787. {
  788.    plotdorbit(real-xxmin,imag-yymax,color);
  789. }
  790.  
  791. void scrub_orbit()
  792. {
  793.    static long oldtime = 0;
  794.    int i,j,c;
  795.    int save_sxoffs,save_syoffs;
  796.    save_sxoffs = sxoffs;
  797.    save_syoffs = syoffs;
  798.    sxoffs = syoffs = 0;
  799.    while(orbit_ptr > 0)
  800.    {
  801.       c = *(save_orbit + --orbit_ptr);
  802.       j = *(save_orbit + --orbit_ptr);
  803.       i = *(save_orbit + --orbit_ptr);
  804.       putcolor(i,j,c);
  805.    }
  806.    sxoffs = save_sxoffs;
  807.    syoffs = save_syoffs;
  808.    nosnd();
  809. }
  810.  
  811.  
  812. int add_worklist(int xfrom, int xto, int yfrom, int yto, int ybegin,
  813. int pass, int sym)
  814. {
  815.    if (num_worklist >= MAXCALCWORK)
  816.       return(-1);
  817.    worklist[num_worklist].xxstart = xfrom;
  818.    worklist[num_worklist].xxstop  = xto;
  819.    worklist[num_worklist].yystart = yfrom;
  820.    worklist[num_worklist].yystop  = yto;
  821.    worklist[num_worklist].yybegin = ybegin;
  822.    worklist[num_worklist].pass      = pass;
  823.    worklist[num_worklist].sym      = sym;
  824.    ++num_worklist;
  825.    tidy_worklist();
  826.    return(0);
  827. }
  828.  
  829. static int _fastcall combine_worklist() /* look for 2 entries which can freely merge */
  830. {
  831.    int i,j;
  832.    for (i=0; i<num_worklist; ++i)
  833.       if (worklist[i].yystart == worklist[i].yybegin)
  834.      for (j=i+1; j<num_worklist; ++j)
  835.         if (worklist[j].sym == worklist[i].sym
  836.         && worklist[j].yystart == worklist[j].yybegin
  837.         && worklist[i].pass == worklist[j].pass)
  838.         {
  839.            if ( worklist[i].xxstart == worklist[j].xxstart
  840.            && worklist[i].xxstop  == worklist[j].xxstop)
  841.            {
  842.           if (worklist[i].yystop+1 == worklist[j].yystart)
  843.           {
  844.              worklist[i].yystop = worklist[j].yystop;
  845.              return(j);
  846.           }
  847.           if (worklist[j].yystop+1 == worklist[i].yystart)
  848.           {
  849.              worklist[i].yystart = worklist[j].yystart;
  850.              worklist[i].yybegin = worklist[j].yybegin;
  851.              return(j);
  852.           }
  853.            }
  854.            if ( worklist[i].yystart == worklist[j].yystart
  855.            && worklist[i].yystop  == worklist[j].yystop)
  856.            {
  857.           if (worklist[i].xxstop+1 == worklist[j].xxstart)
  858.           {
  859.              worklist[i].xxstop = worklist[j].xxstop;
  860.              return(j);
  861.           }
  862.           if (worklist[j].xxstop+1 == worklist[i].xxstart)
  863.           {
  864.              worklist[i].xxstart = worklist[j].xxstart;
  865.              return(j);
  866.           }
  867.            }
  868.         }
  869.    return(0); /* nothing combined */
  870. }
  871.  
  872. void tidy_worklist() /* combine mergeable entries, resort */
  873. {
  874.    int i,j;
  875.    struct workliststuff tempwork;
  876.    while (i=combine_worklist())
  877.    { /* merged two, delete the gone one */
  878.       while (++i < num_worklist)
  879.      worklist[i-1] = worklist[i];
  880.       --num_worklist;
  881.    }
  882.    for (i=0; i<num_worklist; ++i)
  883.       for (j=i+1; j<num_worklist; ++j)
  884.      if (worklist[j].pass < worklist[i].pass
  885.          || (worklist[j].pass == worklist[i].pass
  886.          && (worklist[j].yystart < worklist[i].yystart
  887.          || (   worklist[j].yystart == worklist[i].yystart
  888.          && worklist[j].xxstart <  worklist[i].xxstart))))
  889.      { /* dumb sort, swap 2 entries to correct order */
  890.         tempwork = worklist[i];
  891.         worklist[i] = worklist[j];
  892.         worklist[j] = tempwork;
  893.      }
  894. }
  895.  
  896.  
  897. void get_julia_attractor (double real, double imag)
  898. {
  899.    _LCMPLX lresult;
  900.    _CMPLX result;
  901.    int savper,savmaxit;
  902.    int i;
  903.  
  904.    if (attractors == 0 && finattract == 0) /* not magnet & not requested */
  905.       return;
  906.  
  907.    if (attractors >= N_ATTR)     /* space for more attractors ?  */
  908.       return;               /* Bad luck - no room left !    */
  909.  
  910.    savper = periodicitycheck;
  911.    savmaxit = maxit;
  912.    periodicitycheck = 0;
  913.    old.x = real;            /* prepare for f.p orbit calc */
  914.    old.y = imag;
  915.    tempsqrx = sqr(old.x);
  916.    tempsqry = sqr(old.y);
  917.  
  918.    lold.x = real;        /* prepare for int orbit calc */
  919.    lold.y = imag;
  920.    ltempsqrx = tempsqrx;
  921.    ltempsqry = tempsqry;
  922.  
  923.    lold.x = lold.x << bitshift;
  924.    lold.y = lold.y << bitshift;
  925.    ltempsqrx = ltempsqrx << bitshift;
  926.    ltempsqry = ltempsqry << bitshift;
  927.  
  928.    if (maxit < 500)        /* we're going to try at least this hard */
  929.       maxit = 500;
  930.    color = 0;
  931.    while (++color < maxit)
  932.       if(curfractalspecific->orbitcalc())
  933.      break;
  934.    if (color >= maxit)        /* if orbit stays in the lake */
  935.    {
  936.       if (integerfractal)   /* remember where it went to */
  937.      lresult = lnew;
  938.       else
  939.      result =  new;
  940.      for (i=0;i<10;i++) {
  941.       if(!curfractalspecific->orbitcalc()) /* if it stays in the lake */
  942.       {                /* and doen't move far, probably */
  943.      if (integerfractal)   /*   found a finite attractor    */
  944.      {
  945.         if(labs(lresult.x-lnew.x) < lclosenuff
  946.         && labs(lresult.y-lnew.y) < lclosenuff)
  947.         {
  948.            lattr[attractors] = lnew;
  949.            attrperiod[attractors] = i+1;
  950.            attractors++;   /* another attractor - coloured lakes ! */
  951.            break;
  952.         }
  953.      }
  954.      else
  955.      {
  956.         if(fabs(result.x-new.x) < closenuff
  957.         && fabs(result.y-new.y) < closenuff)
  958.         {
  959.            attr[attractors] = new;
  960.            attrperiod[attractors] = i+1;
  961.            attractors++;   /* another attractor - coloured lakes ! */
  962.            break;
  963.         }
  964.      }
  965.       } else {
  966.       break;
  967.       }
  968.      }
  969.    }
  970.    if(attractors==0)
  971.       periodicitycheck = savper;
  972.    maxit = savmaxit;
  973. }
  974.  
  975.  
  976. #define maxyblk 7    /* must match calcfrac.c */
  977. #define maxxblk 202  /* must match calcfrac.c */
  978. int ssg_blocksize() /* used by solidguessing and by zoom panning */
  979. {
  980.    int blocksize,i;
  981.    /* blocksize 4 if <300 rows, 8 if 300-599, 16 if 600-1199, 32 if >=1200 */
  982.    blocksize=4;
  983.    i=300;
  984.    while(i<=ydots)
  985.    {
  986.       blocksize+=blocksize;
  987.       i+=i;
  988.    }
  989.    /* increase blocksize if prefix array not big enough */
  990.    while(blocksize*(maxxblk-2)<xdots || blocksize*(maxyblk-2)*16<ydots)
  991.       blocksize+=blocksize;
  992.    return(blocksize);
  993. }
  994.  
  995.  
  996. /* Symmetry plot for period PI */
  997. void _fastcall symPIplot(x, y, color)
  998. int x, y, color ;
  999. {
  1000.    while(x <= xxstop)
  1001.    {
  1002.       putcolor(x, y, color) ;
  1003.       x += pixelpi;
  1004.    }
  1005. }
  1006. /* Symmetry plot for period PI plus Origin Symmetry */
  1007. void _fastcall symPIplot2J(x, y, color)
  1008. int x, y, color ;
  1009. {
  1010.    int i,j;
  1011.    while(x <= xxstop)
  1012.    {
  1013.       putcolor(x, y, color) ;
  1014.       if ((i=yystop-(y-yystart)) > iystop && i < ydots
  1015.       && (j=xxstop-(x-xxstart)) < xdots)
  1016.      putcolor(j, i, color) ;
  1017.       x += pixelpi;
  1018.    }
  1019. }
  1020. /* Symmetry plot for period PI plus Both Axis Symmetry */
  1021. void _fastcall symPIplot4J(x, y, color)
  1022. int x, y, color ;
  1023. {
  1024.    int i,j;
  1025.    while(x <= (xxstart+xxstop)/2)
  1026.    {
  1027.       j = xxstop-(x-xxstart);
  1028.       putcolor(       x , y , color) ;
  1029.       if (j < xdots)
  1030.      putcolor(j , y , color) ;
  1031.       if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1032.       {
  1033.      putcolor(x , i , color) ;
  1034.      if (j < xdots)
  1035.         putcolor(j , i , color) ;
  1036.       }
  1037.       x += pixelpi;
  1038.    }
  1039. }
  1040.  
  1041. /* Symmetry plot for X Axis Symmetry */
  1042. void _fastcall symplot2(x, y, color)
  1043. int x, y, color ;
  1044. {
  1045.    int i;
  1046.    putcolor(x, y, color) ;
  1047.    if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1048.       putcolor(x, i, color) ;
  1049. }
  1050.  
  1051. /* Symmetry plot for Y Axis Symmetry */
  1052. void _fastcall symplot2Y(x, y, color)
  1053. int x, y, color ;
  1054. {
  1055.    int i;
  1056.    putcolor(x, y, color) ;
  1057.    if ((i=xxstop-(x-xxstart)) < xdots)
  1058.       putcolor(i, y, color) ;
  1059. }
  1060.  
  1061. /* Symmetry plot for Origin Symmetry */
  1062. void _fastcall symplot2J(x, y, color)
  1063. int x, y, color ;
  1064. {
  1065.    int i,j;
  1066.    putcolor(x, y, color) ;
  1067.    if ((i=yystop-(y-yystart)) > iystop && i < ydots
  1068.        && (j=xxstop-(x-xxstart)) < xdots)
  1069.       putcolor(j, i, color) ;
  1070. }
  1071.  
  1072. /* Symmetry plot for Both Axis Symmetry */
  1073. void _fastcall symplot4(x, y, color)
  1074. int x, y, color ;
  1075. {
  1076.    int i,j;
  1077.    j = xxstop-(x-xxstart);
  1078.    putcolor(       x , y, color) ;
  1079.    if (j < xdots)
  1080.       putcolor(    j , y, color) ;
  1081.    if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1082.    {
  1083.       putcolor(    x , i, color) ;
  1084.       if (j < xdots)
  1085.      putcolor(j , i, color) ;
  1086.    }
  1087. }
  1088.  
  1089. /* Symmetry plot for X Axis Symmetry - Striped Newtbasin version */
  1090. void _fastcall symplot2basin(x, y, color)
  1091. int x, y, color ;
  1092. {
  1093.    int i,stripe;
  1094.    extern int degree;
  1095.    putcolor(x, y, color) ;
  1096.    if(basin==2 && color > 8)
  1097.       stripe=8;
  1098.    else
  1099.       stripe = 0;
  1100.    if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1101.    {
  1102.       color -= stripe;              /* reconstruct unstriped color */
  1103.       color = (degree+1-color)%degree+1;  /* symmetrical color */
  1104.       color += stripe;              /* add stripe */
  1105.       putcolor(x, i,color)  ;
  1106.    }
  1107. }
  1108.  
  1109. /* Symmetry plot for Both Axis Symmetry  - Newtbasin version */
  1110. void _fastcall symplot4basin(x, y, color)
  1111. int x, y, color ;
  1112. {
  1113.    extern int degree;
  1114.    int i,j,color1,stripe;
  1115.    if(color == 0) /* assumed to be "inside" color */
  1116.    {
  1117.       symplot4(x, y, color);
  1118.       return;
  1119.    }
  1120.    if(basin==2 && color > 8)
  1121.       stripe = 8;
  1122.    else
  1123.       stripe = 0;
  1124.    color -= stripe;          /* reconstruct unstriped color */
  1125.    color1 = degree/2+degree+2 - color;
  1126.    if(color < degree/2+2)
  1127.       color1 = degree/2+2 - color;
  1128.    else
  1129.       color1 = degree/2+degree+2 - color;
  1130.    j = xxstop-(x-xxstart);
  1131.    putcolor(       x, y, color+stripe) ;
  1132.    if (j < xdots)
  1133.       putcolor(    j, y, color1+stripe) ;
  1134.    if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1135.    {
  1136.       putcolor(    x, i, stripe + (degree+1 - color)%degree+1) ;
  1137.       if (j < xdots)
  1138.      putcolor(j, i, stripe + (degree+1 - color1)%degree+1) ;
  1139.    }
  1140. }
  1141.  
  1142. /* Do nothing plot!!! */
  1143. void _fastcall noplot(int x,int y,int color)
  1144. {
  1145. }
  1146.