home *** CD-ROM | disk | FTP | other *** search
/ Computerworld 1996 March / Computerworld_1996-03_cd.bin / idg_cd3 / grafika / fraktaly / frasr192 / fracsubr.c < prev    next >
C/C++ Source or Header  |  1995-03-26  |  45KB  |  1,453 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. static void   _fastcall adjust_to_limitsbf(double);
  32. static void   _fastcall smallest_add_bf(bf_t);
  33.        int    resume_len;        /* length of resume info */
  34. static int    resume_offset;        /* offset in resume info gets */
  35.  
  36.  
  37. #define FUDGEFACTOR    29    /* fudge all values up by 2**this */
  38. #define FUDGEFACTOR2    24    /* (or maybe this)          */
  39.  
  40. void fill_dx_array(void)
  41. {
  42.    int i;
  43.    dx0[0] = xxmin;        /* fill up the x, y grids */
  44.    dy0[0] = yymax;
  45.    dx1[0] = dy1[0] = 0;
  46.    for (i = 1; i < xdots; i++ ) {
  47.       dx0[i] = (double)(dx0[0] + i*delxx);
  48.       dy1[i] = (double)(dy1[0] - i*delyy2);
  49.    }
  50.    for (i = 1; i < ydots; i++ ) {
  51.       dy0[i] = (double)(dy0[0] - i*delyy);
  52.       dx1[i] = (double)(dx1[0] + i*delxx2);
  53.    }
  54. }
  55. void fill_lx_array(void)
  56. {
  57.    int i;
  58.    /* note that lx1 & ly1 values can overflow into sign bit; since     */
  59.    /* they're used only to add to lx0/ly0, 2s comp straightens it out  */
  60.    lx0[0] = xmin;         /* fill up the x, y grids */
  61.    ly0[0] = ymax;
  62.    lx1[0] = ly1[0] = 0;
  63.    for (i = 1; i < xdots; i++ ) {
  64.       lx0[i] = lx0[i-1] + delx;
  65.       ly1[i] = ly1[i-1] - dely2;
  66.    }
  67.    for (i = 1; i < ydots; i++ ) {
  68.       ly0[i] = ly0[i-1] - dely;
  69.       lx1[i] = lx1[i-1] + delx2;
  70.    }
  71. }
  72.  
  73. /* returns pointer to nonempty param string, NULL otherwise */
  74. char *typehasparm(int type,int parm)
  75. {
  76.    int extra;
  77.    char *ret = NULL;
  78.    if(0 <= parm && parm < 4)
  79.       ret=fractalspecific[type].param[parm];
  80.    else if(parm >= 4 && parm < MAXPARAMS)    
  81.       if((extra=find_extra_param(type)) > -1)
  82.          ret=moreparams[extra].param[parm-4];
  83.    if(ret)
  84.       if(*ret == 0)
  85.          ret = NULL;
  86.    return(ret);
  87. }
  88.  
  89. void fractal_floattobf(void)
  90. {
  91.    int i;
  92.    init_bf_dec(getprecdbl(CURRENTREZ));
  93.    floattobf(bfxmin,xxmin);
  94.    floattobf(bfxmax,xxmax);
  95.    floattobf(bfymin,yymin);
  96.    floattobf(bfymax,yymax);
  97.    floattobf(bfx3rd,xx3rd);
  98.    floattobf(bfy3rd,yy3rd);
  99.    
  100.    for (i = 0; i < MAXPARAMS; i++)
  101.       if(typehasparm(fractype,i)) 
  102.          floattobf(bfparms[i],param[i]);
  103.    calc_status = 0;
  104. }
  105.  
  106.  
  107. #ifdef _MSC_VER
  108. #if _MSC_VER == 800
  109. /* MSC8 doesn't correctly calculate the address of certain arrays here */
  110. #pragma optimize( "", off )
  111. #endif
  112. #endif
  113.  
  114. void calcfracinit(void) /* initialize a *pile* of stuff for fractal calculation */
  115. {
  116.    int tries = 0;
  117.    int i, gotprec;
  118.    double ftemp;
  119.    coloriter=oldcoloriter = 0L;
  120.    /* set up grid array compactly leaving space at end */
  121.    dx0 = MK_FP(extraseg,0);
  122.    dy1 = (dx1 = (dy0 = dx0 + xdots) + ydots) + ydots;
  123.    lx0 = (long far *) dx0;
  124.    ly1 = (lx1 = (ly0 = lx0 + xdots) + ydots) + ydots;
  125.    if(!(curfractalspecific->flags & BF_MATH))
  126.    {
  127.       int tofloat;
  128.       if((tofloat=curfractalspecific->tofloat) == NOFRACTAL)
  129.          bf_math = 0;
  130.       else if(!(fractalspecific[tofloat].flags & BF_MATH))
  131.          bf_math = 0;
  132.       else if(bf_math)
  133.       {
  134.          curfractalspecific = &fractalspecific[tofloat];
  135.          fractype = tofloat;
  136.       }
  137.    }
  138.    
  139.    /* switch back to double when zooming out if using arbitrary precision */
  140.    if(bf_math)
  141.    {
  142.       gotprec=getprecbf(CURRENTREZ);
  143.       if(gotprec <= DBL_DIG+1 && debugflag != 3200)
  144.       {
  145.          bfcornerstofloat(); 
  146.          bf_math = 0;
  147.       }
  148.       else
  149.          init_bf_dec(gotprec);
  150.    }
  151.    else if((fractype==MANDEL || fractype==MANDELFP) && debugflag==3200)
  152.    {
  153.       fractype=MANDELFP;
  154.       curfractalspecific = &fractalspecific[MANDELFP];   
  155.       fractal_floattobf();
  156.       usr_floatflag = 1;
  157.    }
  158.    else if((fractype==JULIA || fractype==JULIAFP) && debugflag==3200)
  159.    {
  160.       fractype=JULIAFP;
  161.       curfractalspecific = &fractalspecific[JULIAFP];   
  162.       fractal_floattobf();
  163.       usr_floatflag = 1;
  164.    }
  165.    else if((fractype==LMANDELZPOWER || fractype==FPMANDELZPOWER) && debugflag==3200)
  166.    {
  167.       fractype=FPMANDELZPOWER;
  168.       curfractalspecific = &fractalspecific[FPMANDELZPOWER];   
  169.       fractal_floattobf();
  170.       usr_floatflag = 1;
  171.    }   
  172.    else if((fractype==LJULIAZPOWER || fractype==FPJULIAZPOWER) && debugflag==3200)
  173.    {
  174.       fractype=FPJULIAZPOWER;
  175.       curfractalspecific = &fractalspecific[FPJULIAZPOWER];
  176.       fractal_floattobf();
  177.       usr_floatflag = 1;
  178.    }   
  179.    else
  180.       free_bf_vars();
  181.    if(bf_math)
  182.       floatflag=1;
  183.    else   
  184.       floatflag = usr_floatflag;
  185.    if (calc_status == 2) /* on resume, ensure floatflag correct */
  186.       if (curfractalspecific->isinteger)
  187.      floatflag = 0;
  188.       else
  189.      floatflag = 1;
  190.    /* if floating pt only, set floatflag for TAB screen */
  191.    if (!curfractalspecific->isinteger && curfractalspecific->tofloat == NOFRACTAL)
  192.       floatflag = 1;
  193.  
  194. init_restart:
  195.  
  196.    /* the following variables may be forced to a different setting due to
  197.       calc routine constraints;  usr_xxx is what the user last said is wanted,
  198.       xxx is what we actually do in the current situation */
  199.    stdcalcmode        = usr_stdcalcmode;
  200.    periodicitycheck = usr_periodicitycheck;
  201.    distest        = usr_distest;
  202.    biomorph        = usr_biomorph;
  203.  
  204.    potflag = 0;
  205.    if (potparam[0] != 0.0
  206.      && colors >= 64
  207.      && (curfractalspecific->calctype == StandardFractal
  208.      || curfractalspecific->calctype == calcmand
  209.      || curfractalspecific->calctype == calcmandfp)) {
  210.       potflag = 1;
  211.       distest = 0;    /* can't do distest too */
  212.       }
  213.  
  214.    if (distest)
  215.       floatflag = 1;  /* force floating point for dist est */
  216.  
  217.    if (floatflag) { /* ensure type matches floatflag */
  218.       if (curfractalspecific->isinteger != 0
  219.     && curfractalspecific->tofloat != NOFRACTAL)
  220.      fractype = curfractalspecific->tofloat;
  221.       }
  222.    else {
  223.       if (curfractalspecific->isinteger == 0
  224.     && curfractalspecific->tofloat != NOFRACTAL)
  225.      fractype = curfractalspecific->tofloat;
  226.       }
  227.    /* match Julibrot with integer mode of orbit */   
  228.    if(fractype == JULIBROTFP && fractalspecific[neworbittype].isinteger)
  229.    {
  230.       int i;
  231.       if((i=fractalspecific[neworbittype].tofloat) != NOFRACTAL)
  232.          neworbittype = i;
  233.       else
  234.          fractype = JULIBROT;   
  235.    }
  236.    else if(fractype == JULIBROT && fractalspecific[neworbittype].isinteger==0)
  237.    {
  238.       int i;
  239.       if((i=fractalspecific[neworbittype].tofloat) != NOFRACTAL)
  240.          neworbittype = i;
  241.       else
  242.          fractype = JULIBROTFP;   
  243.    }
  244.       
  245.    curfractalspecific = &fractalspecific[fractype];
  246.  
  247.    integerfractal = curfractalspecific->isinteger;
  248.  
  249. /*   if (fractype == JULIBROT)
  250.       rqlim = 4;
  251.    else */ if (potflag && potparam[2] != 0.0)
  252.       rqlim = potparam[2];
  253. /* else if (decomp[0] > 0 && decomp[1] > 0)
  254.       rqlim = (double)decomp[1]; */
  255.    else if (bailout) /* user input bailout */
  256.       rqlim = bailout;
  257.    else if (biomorph != -1) /* biomorph benefits from larger bailout */
  258.       rqlim = 100;
  259.    else
  260.       rqlim = curfractalspecific->orbit_bailout;
  261.    if (integerfractal) /* the bailout limit mustn't be too high here */
  262.       if (rqlim > 127.0) rqlim = 127.0;
  263.  
  264.    if ((curfractalspecific->flags&NOROTATE) != 0) {
  265.       /* ensure min<max and unrotated rectangle */
  266.       if (xxmin > xxmax) { ftemp = xxmax; xxmax = xxmin; xxmin = ftemp; }
  267.       if (yymin > yymax) { ftemp = yymax; yymax = yymin; yymin = ftemp; }
  268.       xx3rd = xxmin; yy3rd = yymin;
  269.       }
  270.  
  271.    /* set up bitshift for integer math */
  272.    bitshift = FUDGEFACTOR2; /* by default, the smaller shift */
  273.    if (integerfractal > 1)  /* use specific override from table */
  274.       bitshift = integerfractal;
  275.    if (integerfractal == 0) /* float? */
  276.       if ((i = curfractalspecific->tofloat) != NOFRACTAL) /* -> int? */
  277.       {
  278.      if (fractalspecific[i].isinteger > 1) /* specific shift? */
  279.         bitshift = fractalspecific[i].isinteger;
  280.       }
  281.       else
  282.          bitshift = 16;  /* to allow larger corners */
  283. /* We want this code if we're using the assembler calcmand */
  284.    if (fractype == MANDEL || fractype == JULIA) { /* adust shift bits if.. */
  285.       if (potflag == 0                  /* not using potential */
  286.     && (param[0] > -2.0 && param[0] < 2.0)  /* parameters not too large */
  287.     && (param[1] > -2.0 && param[1] < 2.0)
  288.     && !invert                  /* and not inverting */
  289.     && biomorph == -1              /* and not biomorphing */
  290.     && rqlim <= 4.0               /* and bailout not too high */
  291.     && (outside > -2 || outside < -6)      /* and no funny outside stuff */
  292.     && debugflag != 1234              /* and not debugging */
  293.     && bailoutest == Mod)              /* and bailout test = mod */
  294.      bitshift = FUDGEFACTOR;          /* use the larger bitshift */
  295.       }
  296.  
  297.    fudge = 1L << bitshift;
  298.  
  299.    l_at_rad = fudge/32768L;
  300.    f_at_rad = 1.0/32768L;
  301.  
  302.    /* now setup arrays of real coordinates corresponding to each pixel */
  303.    if(bf_math)
  304.       adjust_to_limitsbf(1.0); /* make sure all corners in valid range */
  305.    else
  306.    {
  307.       adjust_to_limits(1.0); /* make sure all corners in valid range */
  308.       delxx  = (LDBL)(xxmax - xx3rd) / (LDBL)dxsize; /* calculate stepsizes */
  309.       delyy  = (LDBL)(yymax - yy3rd) / (LDBL)dysize;
  310.       delxx2 = (LDBL)(xx3rd - xxmin) / (LDBL)dysize;
  311.       delyy2 = (LDBL)(yy3rd - yymin) / (LDBL)dxsize;
  312.       fill_dx_array();
  313.    }
  314.  
  315.    if(fractype != CELLULAR)  /* fudgetolong fails w >10 digits in double */
  316.    {
  317.       creal = fudgetolong(param[0]); /* integer equivs for it all */
  318.       cimag = fudgetolong(param[1]);
  319.       xmin  = fudgetolong(xxmin);
  320.       xmax  = fudgetolong(xxmax);
  321.       x3rd  = fudgetolong(xx3rd);
  322.       ymin  = fudgetolong(yymin);
  323.       ymax  = fudgetolong(yymax);
  324.       y3rd  = fudgetolong(yy3rd);
  325.       delx  = fudgetolong((double)delxx);
  326.       dely  = fudgetolong((double)delyy);
  327.       delx2 = fudgetolong((double)delxx2);
  328.       dely2 = fudgetolong((double)delyy2);
  329.    }
  330.  
  331.    /* skip this if plasma to avoid 3d problems */
  332.    /* skip if bf_math to avoid extraseg conflict with dx0 arrays */
  333.    if (fractype != PLASMA && bf_math == 0) 
  334.    {
  335.       if (integerfractal && !invert) 
  336.       {
  337.          if (    (delx  == 0 && delxx  != 0.0)
  338.          || (delx2 == 0 && delxx2 != 0.0)
  339.          || (dely  == 0 && delyy  != 0.0)
  340.          || (dely2 == 0 && delyy2 != 0.0) )
  341.         goto expand_retry;
  342.  
  343.          fill_lx_array();   /* fill up the x,y grids */
  344.      /* past max res?  check corners within 10% of expected */
  345.      if (    ratio_bad((double)lx0[xdots-1]-xmin,(double)xmax-x3rd)
  346.          || ratio_bad((double)ly0[ydots-1]-ymax,(double)y3rd-ymax)
  347.          || ratio_bad((double)lx1[(ydots>>1)-1],((double)x3rd-xmin)/2)
  348.          || ratio_bad((double)ly1[(xdots>>1)-1],((double)ymin-y3rd)/2) ) 
  349.          {
  350. expand_retry:
  351.         if (integerfractal        /* integer fractal type? */
  352.            && curfractalspecific->tofloat != NOFRACTAL)
  353.            floatflag = 1;        /* switch to floating pt */
  354.         else
  355.            adjust_to_limits(2.0);     /* double the size */
  356.         if (calc_status == 2)    /* due to restore of an old file? */
  357.            calc_status = 0;     /*   whatever, it isn't resumable */
  358.         goto init_restart;
  359.      } /* end if ratio bad */
  360.      
  361.      /* re-set corners to match reality */
  362.      xmax = lx0[xdots-1] + lx1[ydots-1];
  363.      ymin = ly0[ydots-1] + ly1[xdots-1];
  364.      x3rd = xmin + lx1[ydots-1];
  365.      y3rd = ly0[ydots-1];
  366.      xxmin = fudgetodouble(xmin);
  367.      xxmax = fudgetodouble(xmax);
  368.      xx3rd = fudgetodouble(x3rd);
  369.      yymin = fudgetodouble(ymin);
  370.      yymax = fudgetodouble(ymax);
  371.      yy3rd = fudgetodouble(y3rd);
  372.       } /* end if (integerfractal && !invert) */
  373.       else 
  374.       {
  375.      /* set up dx0 and dy0 analogs of lx0 and ly0 */
  376.      /* put fractal parameters in doubles */
  377.          dx0[0] = xxmin;        /* fill up the x, y grids */
  378.      dy0[0] = yymax;
  379.      dx1[0] = dy1[0] = 0;
  380.          /* this way of defining the dx and dy arrays is not the most 
  381.             accurate, but it is kept because it is used to determine
  382.             the limit of resolution */
  383.      for (i = 1; i < xdots; i++ ) {
  384.         dx0[i] = (double)(dx0[i-1] + (double)delxx);
  385.         dy1[i] = (double)(dy1[i-1] - (double)delyy2);
  386.         }
  387.      for (i = 1; i < ydots; i++ ) {
  388.         dy0[i] = (double)(dy0[i-1] - (double)delyy);
  389.         dx1[i] = (double)(dx1[i-1] + (double)delxx2);
  390.         }
  391.          if(bf_math == 0) /* redundant test, leave for now */
  392.          {
  393.             double Xctr, Yctr, Xmagfactor, Rotation, Skew;
  394.             LDBL Mag;
  395.             int dec;
  396.             cvtcentermag(&Xctr, &Yctr, &Mag, &Xmagfactor, &Rotation, &Skew);
  397.  
  398.             /* 4 digits of padding sounds good */
  399.             /* keep this aligned with CMDFILES.C */
  400.             dec = getpower10(Mag) + 4; 
  401.  
  402.             /* following is the old logic for detecting failure of double
  403.                precision. It has two advantages: it is independent of the
  404.                representation of numbers, and it is sensitive to resolution
  405.                (allows depper zooms at lower resolution. However it fails
  406.                for rotations of exactly 90 degrees, so we added a safety net
  407.                by using the magnification */
  408.             if(++tries < 2) /* for safety */
  409.             {
  410.             static FCODE err[] = {"precision-detection error"};
  411.             if(tries > 1) stopmsg(0, err);
  412.         if (   ratio_bad(dx0[xdots-1]-xxmin,xxmax-xx3rd)
  413.                 || ratio_bad(dy0[ydots-1]-yymax,yy3rd-yymax)
  414.                 || ratio_bad(dx1[ydots-1],      xx3rd-xxmin)
  415.                 || ratio_bad(dy1[xdots-1],      yymin-yy3rd)
  416.                 || (xxmax==xx3rd && yy3rd==yymax && dec > DBL_DIG+1)) 
  417.         {
  418.                if(curfractalspecific->flags & BF_MATH)
  419.                {
  420.                   fractal_floattobf();
  421.                goto init_restart;
  422.            }   
  423.            goto expand_retry;
  424.         } /* end if ratio_bad etc. */
  425.         } /* end if tries < 2 */
  426.      } /* end if bf_math == 0 */  
  427.  
  428.          /* if long double available, this is more accurate */
  429.      fill_dx_array();    /* fill up the x, y grids */
  430.  
  431.      /* re-set corners to match reality */
  432.      xxmax = dx0[xdots-1] + dx1[ydots-1];
  433.      yymin = dy0[ydots-1] + dy1[xdots-1];
  434.      xx3rd = xxmin + dx1[ydots-1];
  435.      yy3rd = dy0[ydots-1];
  436.       } /* end else */
  437.    } /* end if not plasma */
  438.  
  439.    /* for periodicity close-enough, and for unity: */
  440.    /*      min(max(delx,delx2),max(dely,dely2)       */
  441.    ddelmin = fabs((double)delxx);
  442.    if (fabs((double)delxx2) > ddelmin)
  443.       ddelmin = fabs((double)delxx2);
  444.    if (fabs((double)delyy) > fabs((double)delyy2)) 
  445.    {
  446.       if (fabs((double)delyy) < ddelmin)
  447.          ddelmin = fabs((double)delyy);
  448.    }
  449.    else if (fabs((double)delyy2) < ddelmin)
  450.       ddelmin = fabs((double)delyy2);
  451.    delmin = fudgetolong(ddelmin);
  452.  
  453.    /* calculate factors which plot real values to screen co-ords */
  454.    /* calcfrac.c plot_orbit routines have comments about this     */
  455.    ftemp = (double)((0.0-delyy2) * delxx2 * dxsize * dysize
  456.        - (xxmax-xx3rd) * (yy3rd-yymax));
  457.    if(ftemp != 0)
  458.    {
  459.       plotmx1 = (double)(delxx2 * dxsize * dysize / ftemp);
  460.       plotmx2 = (yy3rd-yymax) * dxsize / ftemp;
  461.       plotmy1 = (double)((0.0-delyy2) * dxsize * dysize / ftemp);
  462.       plotmy2 = (xxmax-xx3rd) * dysize / ftemp;
  463.    }
  464.    if(bf_math == 0)
  465.       free_bf_vars();
  466.    else
  467.    {
  468.       /* zap all of extraseg except high area to flush out bugs */
  469.       /* in production version this code can be deleted */
  470.       char far *extra;
  471.       extra = (char far *)MK_FP(extraseg,0);
  472.       far_memset(extra,0,(unsigned int)(0x10000l-(bflength+2)*22U));
  473.    }   
  474. }
  475.  
  476. #ifdef _MSC_VER
  477. #if _MSC_VER == 800
  478. #pragma optimize( "", on ) /* restore optimization options */
  479. #endif
  480. #endif
  481.  
  482. static long _fastcall fudgetolong(double d)
  483. {
  484.    if ((d *= fudge) > 0) d += 0.5;
  485.    else          d -= 0.5;
  486.    return (long)d;
  487. }
  488.  
  489. static double _fastcall fudgetodouble(long l)
  490. {
  491.    char buf[30];
  492.    double d;
  493.    sprintf(buf,"%.9g",(double)l / fudge);
  494. #ifndef XFRACT
  495.    sscanf(buf,"%lg",&d);
  496. #else
  497.    sscanf(buf,"%lf",&d);
  498. #endif
  499.    return d;
  500. }
  501.  
  502. void adjust_cornerbf(void)
  503. {
  504.    /* make edges very near vert/horiz exact, to ditch rounding errs and */
  505.    /* to avoid problems when delta per axis makes too large a ratio     */
  506.    double ftemp;
  507.    double Xmagfactor, Rotation, Skew;
  508.    LDBL Magnification;
  509.  
  510.    bf_t bftemp, bftemp2;
  511.    bf_t btmp1;
  512.    int saved; saved = save_stack();
  513.    bftemp  = alloc_stack(rbflength+2);
  514.    bftemp2 = alloc_stack(rbflength+2);
  515.    btmp1  =  alloc_stack(rbflength+2);
  516.  
  517.    /* While we're at it, let's adjust the Xmagfactor as well */
  518.    /* use bftemp, bftemp2 as bfXctr, bfYctr */
  519.    cvtcentermagbf(bftemp, bftemp2, &Magnification, &Xmagfactor, &Rotation, &Skew);
  520.    ftemp = fabs(Xmagfactor);
  521.    if (ftemp != 1 && ftemp >= (1-aspectdrift) && ftemp <= (1+aspectdrift))
  522.       {
  523.       Xmagfactor = sign(Xmagfactor);
  524.       cvtcornersbf(bftemp, bftemp2, Magnification, Xmagfactor, Rotation, Skew);
  525.       }
  526.  
  527.    /* ftemp=fabs(xx3rd-xxmin); */
  528.    abs_a_bf(sub_bf(bftemp,bfx3rd,bfxmin)); 
  529.    
  530.    /* ftemp2=fabs(xxmax-xx3rd);*/
  531.    abs_a_bf(sub_bf(bftemp2,bfxmax,bfx3rd));                               
  532.  
  533.    /* if( (ftemp=fabs(xx3rd-xxmin)) < (ftemp2=fabs(xxmax-xx3rd)) ) */
  534.    if(cmp_bf(bftemp,bftemp2) < 0) 
  535.    {
  536.       /* if (ftemp*10000 < ftemp2 && yy3rd != yymax) */
  537.       if (cmp_bf(mult_bf_int(btmp1,bftemp,10000),bftemp2) < 0 
  538.          && cmp_bf(bfy3rd,bfymax) != 0 )
  539.      /* xx3rd = xxmin; */
  540.      copy_bf(bfx3rd, bfxmin);
  541.    }
  542.  
  543.    /* else if (ftemp2*10000 < ftemp && yy3rd != yymin) */
  544.    if (cmp_bf(mult_bf_int(btmp1,bftemp2,10000),bftemp) < 0 
  545.                    && cmp_bf(bfy3rd,bfymin) != 0 )
  546.       /* xx3rd = xxmax; */
  547.       copy_bf(bfx3rd, bfxmax);
  548.  
  549.    /* ftemp=fabs(yy3rd-yymin); */     
  550.    abs_a_bf(sub_bf(bftemp,bfy3rd,bfymin)); 
  551.  
  552.    /* ftemp2=fabs(yymax-yy3rd); */
  553.    abs_a_bf(sub_bf(bftemp2,bfymax,bfy3rd));
  554.     
  555.    /* if( (ftemp=fabs(yy3rd-yymin)) < (ftemp2=fabs(yymax-yy3rd)) ) */
  556.    if(cmp_bf(bftemp,bftemp2) < 0) 
  557.    {
  558.       /* if (ftemp*10000 < ftemp2 && xx3rd != xxmax) */
  559.       if (cmp_bf(mult_bf_int(btmp1,bftemp,10000),bftemp2) < 0 
  560.                  && cmp_bf(bfx3rd,bfxmax) != 0 )
  561.      /* yy3rd = yymin; */
  562.      copy_bf(bfy3rd, bfymin);
  563.    }
  564.  
  565.    /* else if (ftemp2*10000 < ftemp && xx3rd != xxmin) */
  566.      if (cmp_bf(mult_bf_int(btmp1,bftemp2,10000),bftemp) < 0 
  567.                       && cmp_bf(bfx3rd,bfxmin) != 0 )
  568.       /* yy3rd = yymax; */
  569.       copy_bf(bfy3rd, bfymax);
  570.  
  571.  
  572.    restore_stack(saved);
  573. }
  574.  
  575. void adjust_corner(void)
  576. {
  577.    /* make edges very near vert/horiz exact, to ditch rounding errs and */
  578.    /* to avoid problems when delta per axis makes too large a ratio    */
  579.    double ftemp,ftemp2;
  580.    double Xctr, Yctr, Xmagfactor, Rotation, Skew;
  581.    LDBL Magnification;
  582.  
  583.    /* While we're at it, let's adjust the Xmagfactor as well */
  584.    cvtcentermag(&Xctr, &Yctr, &Magnification, &Xmagfactor, &Rotation, &Skew);
  585.    ftemp = fabs(Xmagfactor);
  586.    if (ftemp != 1 && ftemp >= (1-aspectdrift) && ftemp <= (1+aspectdrift))
  587.       {
  588.       Xmagfactor = sign(Xmagfactor);
  589.       cvtcorners(Xctr, Yctr, Magnification, Xmagfactor, Rotation, Skew);
  590.       }
  591.  
  592.    if( (ftemp=fabs(xx3rd-xxmin)) < (ftemp2=fabs(xxmax-xx3rd)) ) {
  593.       if (ftemp*10000 < ftemp2 && yy3rd != yymax)
  594.      xx3rd = xxmin;
  595.       }
  596.  
  597.    if (ftemp2*10000 < ftemp && yy3rd != yymin)
  598.       xx3rd = xxmax;
  599.  
  600.    if( (ftemp=fabs(yy3rd-yymin)) < (ftemp2=fabs(yymax-yy3rd)) ) {
  601.       if (ftemp*10000 < ftemp2 && xx3rd != xxmax)
  602.      yy3rd = yymin;
  603.       }
  604.  
  605.    if (ftemp2*10000 < ftemp && xx3rd != xxmin)
  606.       yy3rd = yymax;
  607.  
  608. }
  609.  
  610. static void _fastcall adjust_to_limitsbf(double expand)
  611. {  
  612.    LDBL limit;
  613.    bf_t bcornerx[4],bcornery[4];
  614.    bf_t blowx,bhighx,blowy,bhighy,blimit,bftemp;
  615.    bf_t bcenterx,bcentery,badjx,badjy,btmp1,btmp2;
  616.    bf_t bexpand;
  617.    int i;
  618.    int saved; saved = save_stack();
  619.    bcornerx[0] = alloc_stack(rbflength+2);
  620.    bcornerx[1] = alloc_stack(rbflength+2);
  621.    bcornerx[2] = alloc_stack(rbflength+2);
  622.    bcornerx[3] = alloc_stack(rbflength+2);
  623.    bcornery[0] = alloc_stack(rbflength+2);
  624.    bcornery[1] = alloc_stack(rbflength+2);
  625.    bcornery[2] = alloc_stack(rbflength+2);
  626.    bcornery[3] = alloc_stack(rbflength+2);
  627.    blowx       = alloc_stack(rbflength+2);
  628.    bhighx      = alloc_stack(rbflength+2);
  629.    blowy       = alloc_stack(rbflength+2);
  630.    bhighy      = alloc_stack(rbflength+2);
  631.    blimit      = alloc_stack(rbflength+2);
  632.    bftemp      = alloc_stack(rbflength+2);
  633.    bcenterx    = alloc_stack(rbflength+2);
  634.    bcentery    = alloc_stack(rbflength+2);
  635.    badjx       = alloc_stack(rbflength+2);
  636.    badjy       = alloc_stack(rbflength+2);
  637.    btmp1       = alloc_stack(rbflength+2);
  638.    btmp2       = alloc_stack(rbflength+2);
  639.    bexpand     = alloc_stack(rbflength+2);
  640.  
  641.    limit = 32767.99;
  642.  
  643.    if (bitshift >= 24) limit = 31.99;
  644.    if (bitshift >= 29) limit = 3.99;
  645.    floattobf(blimit,limit);
  646.    floattobf(bexpand,expand);
  647.    
  648.    add_bf(bcenterx,bfxmin,bfxmax);
  649.    half_a_bf(bcenterx);
  650.  
  651.    /* centery = (yymin+yymax)/2; */
  652.    add_bf(bcentery,bfymin,bfymax);
  653.    half_a_bf(bcentery);
  654.       
  655.    /* if (xxmin == centerx) { */
  656.    if (cmp_bf(bfxmin,bcenterx)==0) { /* ohoh, infinitely thin, fix it */
  657.       smallest_add_bf(bfxmax);
  658.       /* bfxmin -= bfxmax-centerx; */
  659.       sub_a_bf(bfxmin,sub_bf(btmp1,bfxmax,bcenterx));
  660.       }
  661.  
  662.    /* if (bfymin == centery) */
  663.    if (cmp_bf(bfymin,bcentery)==0) { 
  664.       smallest_add_bf(bfymax);
  665.       /* bfymin -= bfymax-centery; */
  666.       sub_a_bf(bfymin,sub_bf(btmp1,bfymax,bcentery));
  667.       }
  668.  
  669.    /* if (bfx3rd == centerx) */
  670.    if (cmp_bf(bfx3rd,bcenterx)==0)  
  671.       smallest_add_bf(bfx3rd);
  672.  
  673.    /* if (bfy3rd == centery) */
  674.    if (cmp_bf(bfy3rd,bcentery)==0)  
  675.       smallest_add_bf(bfy3rd);
  676.  
  677.    /* setup array for easier manipulation */
  678.    /* cornerx[0] = xxmin; */ 
  679.    copy_bf(bcornerx[0],bfxmin); 
  680.  
  681.    /* cornerx[1] = xxmax; */
  682.    copy_bf(bcornerx[1],bfxmax);
  683.  
  684.    /* cornerx[2] = xx3rd; */ 
  685.    copy_bf(bcornerx[2],bfx3rd);
  686.  
  687.    /* cornerx[3] = xxmin+(xxmax-xx3rd); */
  688.    sub_bf(bcornerx[3],bfxmax,bfx3rd);   
  689.    add_a_bf(bcornerx[3],bfxmin);   
  690.  
  691.    /* cornery[0] = yymax; */ 
  692.    copy_bf(bcornery[0],bfymax); 
  693.  
  694.    /* cornery[1] = yymin; */
  695.    copy_bf(bcornery[1],bfymin);
  696.    
  697.    /* cornery[2] = yy3rd; */ 
  698.    copy_bf(bcornery[2],bfy3rd); 
  699.    
  700.    /* cornery[3] = yymin+(yymax-yy3rd); */
  701.    sub_bf(bcornery[3],bfymax,bfy3rd);   
  702.    add_a_bf(bcornery[3],bfymin);   
  703.  
  704.    /* if caller wants image size adjusted, do that first */
  705.    if (expand != 1.0)
  706.    {
  707.       for (i=0; i<4; ++i) {
  708.      /* cornerx[i] = centerx + (cornerx[i]-centerx)*expand; */
  709.      sub_bf(btmp1,bcornerx[i],bcenterx);
  710.      mult_bf(bcornerx[i],btmp1,bexpand);
  711.      add_a_bf(bcornerx[i],bcenterx);
  712.  
  713.      /* cornery[i] = centery + (cornery[i]-centery)*expand; */
  714.      sub_bf(btmp1,bcornery[i],bcentery);
  715.      mult_bf(bcornery[i],btmp1,bexpand);
  716.      add_a_bf(bcornery[i],bcentery);
  717.       }
  718.    }
  719.  
  720.    /* get min/max x/y values */
  721.    /* lowx = highx = cornerx[0]; */
  722.    copy_bf(blowx,bcornerx[0]); copy_bf(bhighx,bcornerx[0]);   
  723.  
  724.    /* lowy = highy = cornery[0]; */
  725.    copy_bf(blowy,bcornery[0]); copy_bf(bhighy,bcornery[0]);   
  726.  
  727.    for (i=1; i<4; ++i) {
  728.       /* if (cornerx[i] < lowx)               lowx  = cornerx[i]; */
  729.       if (cmp_bf(bcornerx[i],blowx) < 0)   copy_bf(blowx,bcornerx[i]);
  730.  
  731.       /* if (cornerx[i] > highx)              highx = cornerx[i]; */
  732.       if (cmp_bf(bcornerx[i],bhighx) > 0)  copy_bf(bhighx,bcornerx[i]);
  733.  
  734.       /* if (cornery[i] < lowy)               lowy  = cornery[i]; */
  735.       if (cmp_bf(bcornery[i],blowy) < 0)   copy_bf(blowy,bcornery[i]);
  736.  
  737.       /* if (cornery[i] > highy)              highy = cornery[i]; */
  738.       if (cmp_bf(bcornery[i],bhighy) > 0)  copy_bf(bhighy,bcornery[i]);
  739.       }
  740.  
  741.    /* if image is too large, downsize it maintaining center */
  742.    /* ftemp = highx-lowx; */
  743.    sub_bf(bftemp,bhighx,blowx);
  744.    
  745.    /* if (highy-lowy > ftemp) ftemp = highy-lowy; */
  746.    if (cmp_bf(sub_bf(btmp1,bhighy,blowy),bftemp) > 0) copy_bf(bftemp,btmp1);
  747.  
  748.    /* if image is too large, downsize it maintaining center */
  749.  
  750.    floattobf(btmp1,limit*2.0);
  751.    copy_bf(btmp2,bftemp);
  752.    div_bf(bftemp,btmp1,btmp2);
  753.    floattobf(btmp1,1.0);
  754.    if (cmp_bf(bftemp,btmp1) < 0) 
  755.       for (i=0; i<4; ++i) {
  756.      /* cornerx[i] = centerx + (cornerx[i]-centerx)*ftemp; */
  757.      sub_bf(btmp1,bcornerx[i],bcenterx);
  758.      mult_bf(bcornerx[i],btmp1,bftemp);
  759.      add_a_bf(bcornerx[i],bcenterx);
  760.      
  761.      /* cornery[i] = centery + (cornery[i]-centery)*ftemp; */
  762.      sub_bf(btmp1,bcornery[i],bcentery);
  763.      mult_bf(bcornery[i],btmp1,bftemp);
  764.      add_a_bf(bcornery[i],bcentery);
  765.      }
  766.  
  767.    /* if any corner has x or y past limit, move the image */
  768.    /* adjx = adjy = 0; */
  769.    clear_bf(badjx); clear_bf(badjy);
  770.    
  771.    for (i=0; i<4; ++i) {
  772.       /* if (cornerx[i] > limit && (ftemp = cornerx[i] - limit) > adjx)
  773.      adjx = ftemp; */
  774.       if (cmp_bf(bcornerx[i],blimit) > 0 && 
  775.           cmp_bf(sub_bf(bftemp,bcornerx[i],blimit),badjx) > 0)
  776.      copy_bf(badjx,bftemp);
  777.  
  778.       /* if (cornerx[i] < 0.0-limit && (ftemp = cornerx[i] + limit) < adjx)
  779.      adjx = ftemp; */
  780.       if (cmp_bf(bcornerx[i],neg_bf(btmp1,blimit)) < 0 && 
  781.           cmp_bf(add_bf(bftemp,bcornerx[i],blimit),badjx) < 0)
  782.      copy_bf(badjx,bftemp);
  783.  
  784.       /* if (cornery[i] > limit     && (ftemp = cornery[i] - limit) > adjy)
  785.      adjy = ftemp; */
  786.       if (cmp_bf(bcornery[i],blimit) > 0 && 
  787.           cmp_bf(sub_bf(bftemp,bcornery[i],blimit),badjy) > 0)
  788.      copy_bf(badjy,bftemp);
  789.  
  790.       /* if (cornery[i] < 0.0-limit && (ftemp = cornery[i] + limit) < adjy)
  791.      adjy = ftemp; */
  792.       if (cmp_bf(bcornery[i],neg_bf(btmp1,blimit)) < 0 && 
  793.           cmp_bf(add_bf(bftemp,bcornery[i],blimit),badjy) < 0)
  794.      copy_bf(badjy,bftemp);
  795.       }
  796.  
  797.    /* if (calc_status == 2 && (adjx != 0 || adjy != 0) && (zwidth == 1.0))
  798.       calc_status = 0; */
  799.    if (calc_status == 2 && (is_bf_not_zero(badjx)|| is_bf_not_zero(badjy)) && (zwidth == 1.0))
  800.       calc_status = 0;
  801.  
  802.    /* xxmin = cornerx[0] - adjx; */
  803.    sub_bf(bfxmin,bcornerx[0],badjx);
  804.    /* xxmax = cornerx[1] - adjx; */
  805.    sub_bf(bfxmax,bcornerx[1],badjx);
  806.    /* xx3rd = cornerx[2] - adjx; */
  807.    sub_bf(bfx3rd,bcornerx[2],badjx);
  808.    /* yymax = cornery[0] - adjy; */
  809.    sub_bf(bfymax,bcornery[0],badjy);
  810.    /* yymin = cornery[1] - adjy; */
  811.    sub_bf(bfymin,bcornery[1],badjy);
  812.    /* yy3rd = cornery[2] - adjy; */
  813.    sub_bf(bfy3rd,bcornery[2],badjy);
  814.  
  815.    adjust_cornerbf(); /* make 3rd corner exact if very near other co-ords */
  816.    restore_stack(saved);
  817. }
  818.  
  819. static void _fastcall adjust_to_limits(double expand)
  820. {  
  821.    double cornerx[4],cornery[4];
  822.    double lowx,highx,lowy,highy,limit,ftemp;
  823.    double centerx,centery,adjx,adjy;
  824.    int i;
  825.  
  826.    limit = 32767.99;
  827.  
  828.    if (bitshift >= 24) limit = 31.99;
  829.    if (bitshift >= 29) limit = 3.99;
  830.  
  831.    centerx = (xxmin+xxmax)/2;
  832.    centery = (yymin+yymax)/2;
  833.       
  834.    if (xxmin == centerx) { /* ohoh, infinitely thin, fix it */
  835.       smallest_add(&xxmax);
  836.       xxmin -= xxmax-centerx;
  837.       }
  838.  
  839.    if (yymin == centery) {
  840.       smallest_add(&yymax);
  841.       yymin -= yymax-centery;
  842.       }
  843.  
  844.    if (xx3rd == centerx)
  845.       smallest_add(&xx3rd);
  846.  
  847.    if (yy3rd == centery)
  848.       smallest_add(&yy3rd);
  849.  
  850.    /* setup array for easier manipulation */
  851.    cornerx[0] = xxmin; 
  852.    cornerx[1] = xxmax;
  853.    cornerx[2] = xx3rd; 
  854.    cornerx[3] = xxmin+(xxmax-xx3rd);
  855.  
  856.    cornery[0] = yymax; 
  857.    cornery[1] = yymin;
  858.    cornery[2] = yy3rd; 
  859.    cornery[3] = yymin+(yymax-yy3rd);
  860.  
  861.    /* if caller wants image size adjusted, do that first */
  862.    if (expand != 1.0)
  863.    {
  864.       for (i=0; i<4; ++i) {
  865.      cornerx[i] = centerx + (cornerx[i]-centerx)*expand;
  866.      cornery[i] = centery + (cornery[i]-centery)*expand;
  867.       }
  868.    }
  869.    /* get min/max x/y values */
  870.    lowx = highx = cornerx[0];
  871.    lowy = highy = cornery[0];
  872.  
  873.    for (i=1; i<4; ++i) {
  874.       if (cornerx[i] < lowx)               lowx  = cornerx[i];
  875.       if (cornerx[i] > highx)              highx = cornerx[i];
  876.       if (cornery[i] < lowy)               lowy  = cornery[i];
  877.       if (cornery[i] > highy)              highy = cornery[i];
  878.       }
  879.  
  880.    /* if image is too large, downsize it maintaining center */
  881.    ftemp = highx-lowx;
  882.    
  883.    if (highy-lowy > ftemp) ftemp = highy-lowy;
  884.  
  885.    /* if image is too large, downsize it maintaining center */
  886.    if ((ftemp = limit*2/ftemp) < 1.0) {
  887.       for (i=0; i<4; ++i) {
  888.      cornerx[i] = centerx + (cornerx[i]-centerx)*ftemp;
  889.      cornery[i] = centery + (cornery[i]-centery)*ftemp;
  890.      }
  891.    }
  892.  
  893.    /* if any corner has x or y past limit, move the image */
  894.    adjx = adjy = 0;
  895.    
  896.    for (i=0; i<4; ++i) {
  897.       if (cornerx[i] > limit && (ftemp = cornerx[i] - limit) > adjx)
  898.      adjx = ftemp;
  899.       if (cornerx[i] < 0.0-limit && (ftemp = cornerx[i] + limit) < adjx)
  900.      adjx = ftemp;
  901.       if (cornery[i] > limit     && (ftemp = cornery[i] - limit) > adjy)
  902.      adjy = ftemp;
  903.       if (cornery[i] < 0.0-limit && (ftemp = cornery[i] + limit) < adjy)
  904.      adjy = ftemp;
  905.       }
  906.    if (calc_status == 2 && (adjx != 0 || adjy != 0) && (zwidth == 1.0))
  907.       calc_status = 0;
  908.    xxmin = cornerx[0] - adjx;
  909.    xxmax = cornerx[1] - adjx;
  910.    xx3rd = cornerx[2] - adjx;
  911.    yymax = cornery[0] - adjy;
  912.    yymin = cornery[1] - adjy;
  913.    yy3rd = cornery[2] - adjy;
  914.  
  915.    adjust_corner(); /* make 3rd corner exact if very near other co-ords */
  916. }
  917.  
  918. static void _fastcall smallest_add(double *num)
  919. {
  920.    *num += *num * 5.0e-16;
  921. }
  922.  
  923. static void _fastcall smallest_add_bf(bf_t num)
  924. {
  925.    bf_t btmp1;
  926.    int saved; saved = save_stack();
  927.    btmp1 = alloc_stack(bflength+2);
  928.    mult_bf(btmp1,floattobf(btmp1, 5.0e-16),num);
  929.    add_a_bf(num,btmp1);
  930.    restore_stack(saved);
  931. }
  932.  
  933. static int _fastcall ratio_bad(double actual, double desired)
  934. {  double ftemp;
  935.    ftemp = 0;
  936.    if (desired != 0 && debugflag != 3400)
  937.       ftemp = actual / desired;
  938.    if (desired != 0 && debugflag != 3400)
  939.       if ((ftemp = actual / desired) < 0.95 || ftemp > 1.05)
  940.      return(1);
  941.    return(0);
  942. }
  943.  
  944.  
  945. /* Save/resume stuff:
  946.  
  947.    Engines which aren't resumable can simply ignore all this.
  948.  
  949.    Before calling the (per_image,calctype) routines (engine), calcfract sets:
  950.       "resuming" to 0 if new image, nonzero if resuming a partially done image
  951.       "calc_status" to 1
  952.    If an engine is interrupted and wants to be able to resume it must:
  953.       store whatever status info it needs to be able to resume later
  954.       set calc_status to 2 and return
  955.    If subsequently called with resuming!=0, the engine must restore status
  956.    info and continue from where it left off.
  957.  
  958.    Since the info required for resume can get rather large for some types,
  959.    it is not stored directly in save_info.  Instead, memory is dynamically
  960.    allocated as required, and stored in .fra files as a separate block.
  961.    To save info for later resume, an engine routine can use:
  962.       alloc_resume(maxsize,version)
  963.      Maxsize must be >= max bytes subsequently saved + 2; over-allocation
  964.      is harmless except for possibility of insufficient mem available;
  965.      undersize is not checked and probably causes serious misbehaviour.
  966.      Version is an arbitrary number so that subsequent revisions of the
  967.      engine can be made backward compatible.
  968.      Alloc_resume sets calc_status to 2 (resumable) if it succeeds; to 3
  969.      if it cannot allocate memory (and issues warning to user).
  970.       put_resume({bytes,&argument,} ... 0)
  971.      Can be called as often as required to store the info.
  972.      Arguments must not be far addresses.
  973.      Is not protected against calls which use more space than allocated.
  974.    To reload info when resuming, use:
  975.       start_resume()
  976.      initializes and returns version number
  977.       get_resume({bytes,&argument,} ... 0)
  978.      inverse of store_resume
  979.       end_resume()
  980.      optional, frees the memory area sooner than would happen otherwise
  981.  
  982.    Example, save info:
  983.       alloc_resume(sizeof(parmarray)+100,2);
  984.       put_resume(sizeof(int),&row, sizeof(int),&col,
  985.          sizeof(parmarray),parmarray, 0);
  986.     restore info:
  987.       vsn = start_resume();
  988.       get_resume(sizeof(int),&row, sizeof(int),&col, 0);
  989.       if (vsn >= 2)
  990.      get_resume(sizeof(parmarray),parmarray,0);
  991.       end_resume();
  992.  
  993.    Engines which allocate a large far memory chunk of their own might
  994.    directly set resume_info, resume_len, calc_status to avoid doubling
  995.    transient memory needs by using these routines.
  996.  
  997.    StandardFractal, calcmand, solidguess, and bound_trace_main are a related
  998.    set of engines for escape-time fractals.  They use a common worklist
  999.    structure for save/resume.  Fractals using these must specify calcmand
  1000.    or StandardFractal as the engine in fractalspecificinfo.
  1001.    Other engines don't get btm nor ssg, don't get off-axis symmetry nor
  1002.    panning (the worklist stuff), and are on their own for save/resume.
  1003.  
  1004.    */
  1005.  
  1006. #ifndef XFRACT
  1007. int put_resume(int len, ...)
  1008. #else
  1009. int put_resume(va_alist)
  1010. va_dcl
  1011. #endif
  1012. {
  1013.    va_list arg_marker;    /* variable arg list */
  1014.    char *source_ptr;
  1015. #ifdef XFRACT
  1016.    int len;
  1017. #endif
  1018.  
  1019.    if (resume_info == NULL)
  1020.       return(-1);
  1021. #ifndef XFRACT
  1022.    va_start(arg_marker,len);
  1023. #else
  1024.    va_start(arg_marker);
  1025.    len = va_arg(arg_marker,int);
  1026. #endif
  1027.    while (len)
  1028.    {
  1029.       source_ptr = va_arg(arg_marker,char *);
  1030.       far_memcpy(resume_info+resume_len,source_ptr,len);
  1031.       resume_len += len;
  1032.       len = va_arg(arg_marker,int);
  1033.    }
  1034.    return(0);
  1035. }
  1036.  
  1037. int alloc_resume(int alloclen, int version)
  1038. {
  1039.    if (resume_info != NULL) /* free the prior area if there is one */
  1040.       farmemfree(resume_info);
  1041.    if ((resume_info = farmemalloc((long)alloclen))== NULL)
  1042.    {
  1043.       static FCODE msg[] = {"\
  1044. Warning - insufficient free memory to save status.\n\
  1045. You will not be able to resume calculating this image."};
  1046.       stopmsg(0,msg);
  1047.       calc_status = 3;
  1048.       return(-1);
  1049.    }
  1050.    resume_len = 0;
  1051.    put_resume(sizeof(int),&version,0);
  1052.    calc_status = 2;
  1053.    return(0);
  1054. }
  1055.  
  1056. #ifndef XFRACT
  1057. int get_resume(int len, ...)
  1058. #else
  1059. int get_resume(va_alist)
  1060. va_dcl
  1061. #endif
  1062. {
  1063.    va_list arg_marker;    /* variable arg list */
  1064.    char *dest_ptr;
  1065. #ifdef XFRACT
  1066.    int len;
  1067. #endif
  1068.  
  1069.    if (resume_info == NULL)
  1070.       return(-1);
  1071. #ifndef XFRACT
  1072.    va_start(arg_marker,len);
  1073. #else
  1074.    va_start(arg_marker);
  1075.    len = va_arg(arg_marker,int);
  1076. #endif
  1077.    while (len)
  1078.    {
  1079.       dest_ptr = va_arg(arg_marker,char *);
  1080.       far_memcpy(dest_ptr,resume_info+resume_offset,len);
  1081.       resume_offset += len;
  1082.       len = va_arg(arg_marker,int);
  1083.    }
  1084.    return(0);
  1085. }
  1086.  
  1087. int start_resume(void)
  1088. {
  1089.    int version;
  1090.    if (resume_info == NULL)
  1091.       return(-1);
  1092.    resume_offset = 0;
  1093.    get_resume(sizeof(int),&version,0);
  1094.    return(version);
  1095. }
  1096.  
  1097. void end_resume(void)
  1098. {
  1099.    if (resume_info != NULL) /* free the prior area if there is one */
  1100.    {
  1101.       farmemfree(resume_info);
  1102.       resume_info = NULL;
  1103.    }
  1104. }
  1105.  
  1106.  
  1107. /* Showing orbit requires converting real co-ords to screen co-ords.
  1108.    Define:
  1109.        Xs == xxmax-xx3rd           Ys == yy3rd-yymax
  1110.        W  == xdots-1               D  == ydots-1
  1111.    We know that:
  1112.        realx == lx0[col] + lx1[row]
  1113.        realy == ly0[row] + ly1[col]
  1114.        lx0[col] == (col/width) * Xs + xxmin
  1115.        lx1[row] == row * delxx
  1116.        ly0[row] == (row/D) * Ys + yymax
  1117.        ly1[col] == col * (0-delyy)
  1118.   so:
  1119.        realx == (col/W) * Xs + xxmin + row * delxx
  1120.        realy == (row/D) * Ys + yymax + col * (0-delyy)
  1121.   and therefore:
  1122.        row == (realx-xxmin - (col/W)*Xs) / Xv     (1)
  1123.        col == (realy-yymax - (row/D)*Ys) / Yv     (2)
  1124.   substitute (2) into (1) and solve for row:
  1125.        row == ((realx-xxmin)*(0-delyy2)*W*D - (realy-yymax)*Xs*D)
  1126.               / ((0-delyy2)*W*delxx2*D-Ys*Xs)
  1127.   */
  1128.  
  1129. /* sleep N * a tenth of a millisecond */
  1130.  
  1131. void sleepms(long ms)
  1132. {
  1133.     static long scalems = 0L;
  1134.     int savehelpmode,savetabmode;
  1135.     struct timeb t1,t2;
  1136. #define SLEEPINIT 250 /* milliseconds for calibration */
  1137.     savetabmode  = tabmode;
  1138.     savehelpmode = helpmode;
  1139.     tabmode  = 0;
  1140.     helpmode = -1;
  1141.     if(scalems==0L) /* calibrate */
  1142.     {
  1143.     /* selects a value of scalems that makes the units
  1144.        10000 per sec independent of CPU speed */
  1145.     int i,elapsed;
  1146.     scalems = 1L;
  1147.     if(keypressed()) /* check at start, hope to get start of timeslice */
  1148.        goto sleepexit;
  1149.     /* calibrate, assume slow computer first */
  1150.         showtempmsg("Calibrating timer");
  1151.     do
  1152.     {
  1153.        scalems *= 2;
  1154.        ftime(&t2);
  1155.        do { /* wait for the start of a new tick */
  1156.           ftime(&t1);
  1157.         }
  1158.         while (t2.time == t1.time && t2.millitm == t1.millitm);
  1159.        sleepms(10L * SLEEPINIT); /* about 1/4 sec */
  1160.        ftime(&t2);
  1161.        if(keypressed()) {
  1162.           scalems = 0L;
  1163.           goto sleepexit;
  1164.        }
  1165.      }
  1166.      while ((elapsed = (int)(t2.time-t1.time)*1000 + t2.millitm-t1.millitm)
  1167.         < SLEEPINIT);
  1168.     /* once more to see if faster (eg multi-tasking) */
  1169.     do { /* wait for the start of a new tick */
  1170.        ftime(&t1);
  1171.        }
  1172.      while (t2.time == t1.time && t2.millitm == t1.millitm);
  1173.     sleepms(10L * SLEEPINIT);
  1174.     ftime(&t2);
  1175.     if ((i = (int)(t2.time-t1.time)*1000 + t2.millitm-t1.millitm) < elapsed)
  1176.        elapsed = (i == 0) ? 1 : i;
  1177.     scalems = (long)((float)SLEEPINIT/(float)(elapsed) * scalems);
  1178. #if 0
  1179.     char msg[80];
  1180.     if (debugflag == 700) {
  1181.        sprintf(msg,"scale factor=%ld",scalems);
  1182.        stopmsg(0,msg);
  1183.     }
  1184. #endif
  1185.         cleartempmsg();
  1186.     }
  1187.     if(ms > 10L * SLEEPINIT) { /* using ftime is probably more accurate */
  1188.     ms /= 10;
  1189.     ftime(&t1);
  1190.     for(;;) {
  1191.        if(keypressed()) break;
  1192.        ftime(&t2);
  1193.        if ((long)((t2.time-t1.time)*1000 + t2.millitm-t1.millitm) >= ms) break;
  1194.     }
  1195.     }
  1196.     else
  1197.     if(!keypressed()) {
  1198.        ms *= scalems;
  1199.        while(ms-- >= 0);
  1200.     }
  1201. sleepexit:
  1202.     tabmode  = savetabmode;
  1203.     helpmode = savehelpmode;
  1204. }
  1205.  
  1206.  
  1207. static void _fastcall plotdorbit(double dx, double dy, int color)
  1208. {
  1209.    int i, j, c;
  1210.    int save_sxoffs,save_syoffs;
  1211.    if (orbit_ptr >= 1500) return;
  1212.    i = (int)(dy * plotmx1 - dx * plotmx2); i += sxoffs;
  1213.    if (i < 0 || i >= sxdots) return;
  1214.    j = (int)(dx * plotmy1 - dy * plotmy2); j += syoffs;
  1215.    if (j < 0 || j >= sydots) return;
  1216.    save_sxoffs = sxoffs;
  1217.    save_syoffs = syoffs;
  1218.    sxoffs = syoffs = 0;
  1219.    /* save orbit value */
  1220.    if(color == -1)
  1221.    {
  1222.       *(save_orbit + orbit_ptr++) = i;
  1223.       *(save_orbit + orbit_ptr++) = j;
  1224.       *(save_orbit + orbit_ptr++) = c = getcolor(i,j);
  1225.       putcolor(i,j,c^orbit_color);
  1226.    }
  1227.    else
  1228.       putcolor(i,j,color);
  1229.    sxoffs = save_sxoffs;
  1230.    syoffs = save_syoffs;
  1231.    if(orbit_delay > 0)
  1232.       sleepms(orbit_delay);
  1233.    if(soundflag==1)
  1234.        snd((int)(i*1000/xdots+basehertz));
  1235.    else if(soundflag > 1)
  1236.        snd((int)(j*1000/ydots+basehertz));
  1237.    /* placing sleepms here delays each dot */
  1238. }
  1239.  
  1240. void iplot_orbit(long ix, long iy, int color)
  1241. {
  1242.    plotdorbit((double)ix/fudge-xxmin,(double)iy/fudge-yymax,color);
  1243. }
  1244.  
  1245. void plot_orbit(double real,double imag,int color)
  1246. {
  1247.    plotdorbit(real-xxmin,imag-yymax,color);
  1248. }
  1249.  
  1250. void scrub_orbit(void)
  1251. {
  1252.    int i,j,c;
  1253.    int save_sxoffs,save_syoffs;
  1254.    save_sxoffs = sxoffs;
  1255.    save_syoffs = syoffs;
  1256.    sxoffs = syoffs = 0;
  1257.    while(orbit_ptr > 0)
  1258.    {
  1259.       c = *(save_orbit + --orbit_ptr);
  1260.       j = *(save_orbit + --orbit_ptr);
  1261.       i = *(save_orbit + --orbit_ptr);
  1262.       putcolor(i,j,c);
  1263.    }
  1264.    sxoffs = save_sxoffs;
  1265.    syoffs = save_syoffs;
  1266.    nosnd();
  1267. }
  1268.  
  1269.  
  1270. int add_worklist(int xfrom, int xto, int yfrom, int yto, int ybegin,
  1271. int pass, int sym)
  1272. {
  1273.    if (num_worklist >= MAXCALCWORK)
  1274.       return(-1);
  1275.    worklist[num_worklist].xxstart = xfrom;
  1276.    worklist[num_worklist].xxstop  = xto;
  1277.    worklist[num_worklist].yystart = yfrom;
  1278.    worklist[num_worklist].yystop  = yto;
  1279.    worklist[num_worklist].yybegin = ybegin;
  1280.    worklist[num_worklist].pass      = pass;
  1281.    worklist[num_worklist].sym      = sym;
  1282.    ++num_worklist;
  1283.    tidy_worklist();
  1284.    return(0);
  1285. }
  1286.  
  1287. static int _fastcall combine_worklist(void) /* look for 2 entries which can freely merge */
  1288. {
  1289.    int i,j;
  1290.    for (i=0; i<num_worklist; ++i)
  1291.       if (worklist[i].yystart == worklist[i].yybegin)
  1292.      for (j=i+1; j<num_worklist; ++j)
  1293.         if (worklist[j].sym == worklist[i].sym
  1294.         && worklist[j].yystart == worklist[j].yybegin
  1295.         && worklist[i].pass == worklist[j].pass)
  1296.         {
  1297.            if ( worklist[i].xxstart == worklist[j].xxstart
  1298.            && worklist[i].xxstop  == worklist[j].xxstop)
  1299.            {
  1300.           if (worklist[i].yystop+1 == worklist[j].yystart)
  1301.           {
  1302.              worklist[i].yystop = worklist[j].yystop;
  1303.              return(j);
  1304.           }
  1305.           if (worklist[j].yystop+1 == worklist[i].yystart)
  1306.           {
  1307.              worklist[i].yystart = worklist[j].yystart;
  1308.              worklist[i].yybegin = worklist[j].yybegin;
  1309.              return(j);
  1310.           }
  1311.            }
  1312.            if ( worklist[i].yystart == worklist[j].yystart
  1313.            && worklist[i].yystop  == worklist[j].yystop)
  1314.            {
  1315.           if (worklist[i].xxstop+1 == worklist[j].xxstart)
  1316.           {
  1317.              worklist[i].xxstop = worklist[j].xxstop;
  1318.              return(j);
  1319.           }
  1320.           if (worklist[j].xxstop+1 == worklist[i].xxstart)
  1321.           {
  1322.              worklist[i].xxstart = worklist[j].xxstart;
  1323.              return(j);
  1324.           }
  1325.            }
  1326.         }
  1327.    return(0); /* nothing combined */
  1328. }
  1329.  
  1330. void tidy_worklist(void) /* combine mergeable entries, resort */
  1331. {
  1332.    int i,j;
  1333.    WORKLIST tempwork;
  1334.    while ((i=combine_worklist()) != 0)
  1335.    { /* merged two, delete the gone one */
  1336.       while (++i < num_worklist)
  1337.      worklist[i-1] = worklist[i];
  1338.       --num_worklist;
  1339.    }
  1340.    for (i=0; i<num_worklist; ++i)
  1341.       for (j=i+1; j<num_worklist; ++j)
  1342.      if (worklist[j].pass < worklist[i].pass
  1343.          || (worklist[j].pass == worklist[i].pass
  1344.          && (worklist[j].yystart < worklist[i].yystart
  1345.          || (   worklist[j].yystart == worklist[i].yystart
  1346.          && worklist[j].xxstart <  worklist[i].xxstart))))
  1347.      { /* dumb sort, swap 2 entries to correct order */
  1348.         tempwork = worklist[i];
  1349.         worklist[i] = worklist[j];
  1350.         worklist[j] = tempwork;
  1351.      }
  1352. }
  1353.  
  1354.  
  1355. void get_julia_attractor (double real, double imag)
  1356. {
  1357.    _LCMPLX lresult;
  1358.    _CMPLX result;
  1359.    int savper;
  1360.    long savmaxit;
  1361.    int i;
  1362.  
  1363.    if (attractors == 0 && finattract == 0) /* not magnet & not requested */
  1364.       return;
  1365.  
  1366.    if (attractors >= N_ATTR)     /* space for more attractors ?  */
  1367.       return;               /* Bad luck - no room left !    */
  1368.  
  1369.    savper = periodicitycheck;
  1370.    savmaxit = maxit;
  1371.    periodicitycheck = 0;
  1372.    old.x = real;            /* prepare for f.p orbit calc */
  1373.    old.y = imag;
  1374.    tempsqrx = sqr(old.x);
  1375.    tempsqry = sqr(old.y);
  1376.  
  1377.    lold.x = (long)real;        /* prepare for int orbit calc */
  1378.    lold.y = (long)imag;
  1379.    ltempsqrx = (long)tempsqrx;
  1380.    ltempsqry = (long)tempsqry;
  1381.  
  1382.    lold.x = lold.x << bitshift;
  1383.    lold.y = lold.y << bitshift;
  1384.    ltempsqrx = ltempsqrx << bitshift;
  1385.    ltempsqry = ltempsqry << bitshift;
  1386.  
  1387.    if (maxit < 500)        /* we're going to try at least this hard */
  1388.       maxit = 500;
  1389.    color = 0;
  1390.    while (++color < maxit)
  1391.       if(curfractalspecific->orbitcalc())
  1392.      break;
  1393.    if (color >= maxit)        /* if orbit stays in the lake */
  1394.    {
  1395.       if (integerfractal)   /* remember where it went to */
  1396.      lresult = lnew;
  1397.       else
  1398.      result =  new;
  1399.      for (i=0;i<10;i++) {
  1400.       if(!curfractalspecific->orbitcalc()) /* if it stays in the lake */
  1401.       {                /* and doen't move far, probably */
  1402.      if (integerfractal)   /*   found a finite attractor    */
  1403.      {
  1404.         if(labs(lresult.x-lnew.x) < lclosenuff
  1405.         && labs(lresult.y-lnew.y) < lclosenuff)
  1406.         {
  1407.            lattr[attractors] = lnew;
  1408.            attrperiod[attractors] = i+1;
  1409.            attractors++;   /* another attractor - coloured lakes ! */
  1410.            break;
  1411.         }
  1412.      }
  1413.      else
  1414.      {
  1415.         if(fabs(result.x-new.x) < closenuff
  1416.         && fabs(result.y-new.y) < closenuff)
  1417.         {
  1418.            attr[attractors] = new;
  1419.            attrperiod[attractors] = i+1;
  1420.            attractors++;   /* another attractor - coloured lakes ! */
  1421.            break;
  1422.         }
  1423.      }
  1424.       } else {
  1425.       break;
  1426.       }
  1427.      }
  1428.    }
  1429.    if(attractors==0)
  1430.       periodicitycheck = savper;
  1431.    maxit = savmaxit;
  1432. }
  1433.  
  1434.  
  1435. #define maxyblk 7    /* must match calcfrac.c */
  1436. #define maxxblk 202  /* must match calcfrac.c */
  1437. int ssg_blocksize(void) /* used by solidguessing and by zoom panning */
  1438. {
  1439.    int blocksize,i;
  1440.    /* blocksize 4 if <300 rows, 8 if 300-599, 16 if 600-1199, 32 if >=1200 */
  1441.    blocksize=4;
  1442.    i=300;
  1443.    while(i<=ydots)
  1444.    {
  1445.       blocksize+=blocksize;
  1446.       i+=i;
  1447.    }
  1448.    /* increase blocksize if prefix array not big enough */
  1449.    while(blocksize*(maxxblk-2)<xdots || blocksize*(maxyblk-2)*16<ydots)
  1450.       blocksize+=blocksize;
  1451.    return(blocksize);
  1452. }
  1453.