home *** CD-ROM | disk | FTP | other *** search
/ Shareware Overload / ShartewareOverload.cdr / windows / winsrc.zip / MPMATH_C.C < prev    next >
Text File  |  1990-11-01  |  11KB  |  445 lines
  1. /* MPMath_c.c (C) 1989, Mark C. Peterson, CompuServe [70441,3353]
  2.      All rights reserved.
  3.  
  4.    Code may be used in any program provided the author is credited
  5.      either during program execution or in the documentation.  Source
  6.      code may be distributed only in combination with public domain or
  7.      shareware source code.  Source code may be modified provided the
  8.      copyright notice and this message is left unchanged and all
  9.      modifications are clearly documented.
  10.  
  11.      I would appreciate a copy of any work which incorporates this code,
  12.      however this is optional.
  13.  
  14.      Mark C. Peterson
  15.      253 West St., H
  16.      Planstville, CT -6479
  17.      (203) 276-9474
  18. */
  19.  
  20. #include "mpmath.h"
  21.  
  22. int MPaccuracy = 32;
  23.  
  24. struct MP MPsub(struct MP x, struct MP y) {
  25.    y.Exp ^= 0x8000;
  26.    return(MPadd(x, y));
  27. }
  28.  
  29. /* added by TW */
  30. struct MP MPsub086(struct MP x, struct MP y) {
  31.    y.Exp ^= 0x8000;
  32.    return(MPadd086(x, y));
  33. }
  34.  
  35. /* added by TW */
  36. struct MP MPsub386(struct MP x, struct MP y) {
  37.    y.Exp ^= 0x8000;
  38.    return(MPadd386(x, y));
  39. }
  40.  
  41. struct MP MPabs(struct MP x) {
  42.    x.Exp &= 0x7fff;
  43.    return(x);
  44. }
  45.  
  46. struct MPC MPCsqr(struct MPC x) {
  47.    struct MPC z;
  48.  
  49.     z.x = pMPsub(pMPmul(x.x, x.x), pMPmul(x.y, x.y));
  50.     z.y = pMPmul(x.x, x.y);
  51.     z.y.Exp++;
  52.    return(z);
  53. }
  54.  
  55. struct MP MPCmod(struct MPC x) {
  56.     return(pMPadd(pMPmul(x.x, x.x), pMPmul(x.y, x.y)));
  57. }
  58.  
  59. struct MPC MPCmul(struct MPC x, struct MPC y) {
  60.    struct MPC z, t;
  61.  
  62.     z.x = pMPsub(pMPmul(x.x, y.x), pMPmul(x.y, y.y));
  63.     z.y = pMPadd(pMPmul(x.x, y.y), pMPmul(x.y, y.x));
  64.    return(z);
  65. }
  66.  
  67. struct MPC MPCdiv(struct MPC x, struct MPC y) {
  68.    struct MP mod;
  69.  
  70.    mod = MPCmod(y);
  71.     y.y.Exp ^= 0x8000;
  72.     y.x = pMPdiv(y.x, mod);
  73.     y.y = pMPdiv(y.y, mod);
  74.    return(MPCmul(x, y));
  75. }
  76.  
  77. struct MPC MPCadd(struct MPC x, struct MPC y) {
  78.    struct MPC z;
  79.  
  80.     z.x = pMPadd(x.x, y.x);
  81.     z.y = pMPadd(x.y, y.y);
  82.    return(z);
  83. }
  84.  
  85. struct MPC MPCsub(struct MPC x, struct MPC y) {
  86.    struct MPC z;
  87.  
  88.     z.x = pMPsub(x.x, y.x);
  89.     z.y = pMPsub(x.y, y.y);
  90.    return(z);
  91. }
  92.  
  93. struct MPC MPCone = { 0x3fff, 0x80000000l, 0, 0l };
  94.  
  95. struct MPC MPCpow(struct MPC x, int exp) {
  96.    struct MPC z;
  97.    struct MPC zz;
  98.  
  99.    if(exp & 1)
  100.       z = x;
  101.    else
  102.       z = MPCone;
  103.    exp >>= 1;
  104.    while(exp) {
  105.         zz.x = pMPsub(MPmul(x.x, x.x), pMPmul(x.y, x.y));
  106.         zz.y = pMPmul(x.x, x.y);
  107.         zz.y.Exp++;
  108.       x = zz;
  109.       if(exp & 1) {
  110.             zz.x = pMPsub(pMPmul(z.x, x.x), pMPmul(z.y, x.y));
  111.             zz.y = pMPadd(pMPmul(z.x, x.y), pMPmul(z.y, x.x));
  112.      z = zz;
  113.       }
  114.       exp >>= 1;
  115.    }
  116.    return(z);
  117. }
  118.  
  119. int MPCcmp(struct MPC x, struct MPC y) {
  120.    struct MPC z;
  121.  
  122.     if(pMPcmp(x.x, y.x) || pMPcmp(x.y, y.y)) {
  123.         z.x = MPCmod(x);
  124.         z.y = MPCmod(y);
  125.         return(pMPcmp(z.x, z.y));
  126.    }
  127.    else
  128.       return(0);
  129. }
  130.  
  131. struct complex MPC2cmplx(struct MPC x) {
  132.    struct complex z;
  133.  
  134.     z.x = *pMP2d(x.x);
  135.     z.y = *pMP2d(x.y);
  136.    return(z);
  137. }
  138.  
  139. struct MPC cmplx2MPC(struct complex z) {
  140.    struct MPC x;
  141.  
  142.     x.x = pd2MP(z.x);
  143.     x.y = pd2MP(z.y);
  144.    return(x);
  145. }
  146.  
  147. /* function pointer versions added by Tim Wegner 12/07/89 */
  148. int       (*ppMPcmp)() = MPcmp086;
  149. int       (*pMPcmp)(struct MP x, struct MP y) = MPcmp086;
  150. struct MP  (*pMPmul)(struct MP x, struct MP y) = MPmul086;
  151. struct MP  (*pMPdiv)(struct MP x, struct MP y) = MPdiv086;
  152. struct MP  (*pMPadd)(struct MP x, struct MP y) = MPadd086;
  153. struct MP  (*pMPsub)(struct MP x, struct MP y) = MPsub086;
  154. struct MP  (*pd2MP)(double x)               = d2MP086 ;
  155. double *   (*pMP2d)(struct MP m)           = MP2d086 ;
  156. struct MP  (*pfg2MP)(long x, int fg)           = fg2MP086;
  157.  
  158. void setMPfunctions(void) {
  159.    if(cpu == 386)
  160.    {
  161.       pMPmul = MPmul386;
  162.       pMPdiv = MPdiv386;
  163.       pMPadd = MPadd386;
  164.       pMPsub = MPsub386;
  165.       pMPcmp = MPcmp386;
  166.       pd2MP  = d2MP386 ;
  167.       pMP2d  = MP2d386 ;
  168.       pfg2MP = fg2MP386;
  169.    }
  170.    else
  171.    {
  172.       pMPmul = MPmul086;
  173.       pMPdiv = MPdiv086;
  174.       pMPadd = MPadd086;
  175.       pMPsub = MPsub086;
  176.       pMPcmp = MPcmp086;
  177.       pd2MP  = d2MP086 ;
  178.       pMP2d  = MP2d086 ;
  179.       pfg2MP = fg2MP086;
  180.    }
  181. }
  182.  
  183. #define sqr(x) ((x) * (x))
  184.  
  185. extern int debugflag, fpu;
  186.  
  187. struct complex ComplexPower(struct complex x, struct complex y) {
  188.    struct complex z, cLog, t;
  189.    double dist, e2x, siny, cosy;
  190.  
  191.    FPUcplxlog(&x, &cLog);
  192.    FPUcplxmul(&cLog, &y, &t);
  193.  
  194.    if(fpu == 387)
  195.       FPUcplxexp387(&t, &z);
  196.    else {
  197.       e2x = exp(t.x);
  198.       FPUsincos(&t.y, &siny, &cosy);
  199.       z.x = e2x * cosy;
  200.       z.y = e2x * siny;
  201.    }
  202.    return(z);
  203. }
  204.  
  205.  
  206. /***** FRACTINT specific routines and variables *****/
  207.  
  208. #ifndef TESTING_MATH
  209.  
  210. #include "fractint.h"
  211.  
  212. extern double param[];
  213. extern struct complex old, new, init;
  214. extern double threshold, roverd, d1overd, dx0[], dy0[];
  215. extern int periodicitycheck, row, col, debugflag;
  216.  
  217.  
  218. #include <float.h>
  219. #include <stdlib.h>
  220.  
  221. extern int  xdots, ydots;     /* coordinates of dots on the screen  */
  222. extern int  colors;          /* maximum colors available */
  223. extern int  maxit;
  224.  
  225. unsigned char far *LogTable = (unsigned char far *)0;
  226. extern int LogFlag;
  227.    /* LogFlag == 1  -- standard log palettes
  228.       LogFlag == -1 -- 'old' log palettes
  229.       LogFlag >  1  -- compress counts < LogFlag into color #1
  230.       LogFlag < -1  -- use quadratic palettes based on square roots && compress
  231.    */
  232.  
  233. void SetupLogTable(void) {
  234.    float l, f, c, m;
  235.    unsigned n, x, prev, limit, lf;
  236.  
  237.    if(colors < 16) {
  238.       stopmsg(0,"Need at least 16 colors to use logaritmic/quadratic palettes.");
  239.       LogFlag = 0;
  240.       return;
  241.    }
  242.    if(!(LogTable = farmemalloc((long)maxit + 1))) {
  243.       stopmsg(0,"Insufficient memory for logarithmic palette with this maxiter.");
  244.       LogFlag = 0;
  245.       return;
  246.    }
  247.    if (LogFlag > -2) {
  248.       lf = (LogFlag > 1) ? LogFlag : 0;
  249.       if (lf >= maxit)
  250.      lf = maxit - 1;
  251.       Fg2Float((long)(maxit-lf), 0, m);
  252.       fLog14(m, m);
  253.       Fg2Float((long)(colors-(lf?2:1)), 0, c);
  254.       fDiv(m, c, m);
  255.       for (prev = 1; prev <= lf; prev++)
  256.      LogTable[prev] = 1;
  257.       for (n = (lf?2:1); n < colors; n++) {
  258.      Fg2Float((long)n, 0, f);
  259.      fMul16(f, m, f);
  260.      fExp14(f, l);
  261.      limit = Float2Fg(l, 0) + lf;
  262.      if (limit > maxit || n == colors-1)
  263.         limit = maxit;
  264.      while (prev <= limit)
  265.         LogTable[prev++] = n;
  266.       }
  267.    } else {
  268.       if ((lf = 0 - LogFlag) >= maxit)
  269.      lf = maxit - 1;
  270.       Fg2Float((long)(maxit-lf), 0, m);
  271.       fSqrt14(m, m);
  272.       Fg2Float((long)(colors-2), 0, c);
  273.       fDiv(m, c, m);
  274.       for (prev = 1; prev <= lf; prev++)
  275.      LogTable[prev] = 1;
  276.       for (n = 2; n < colors; n++) {
  277.      Fg2Float((long)n, 0, f);
  278.      fMul16(f, m, f);
  279.      fMul16(f, f, l);
  280.      limit = Float2Fg(l, 0) + lf;
  281.      if (limit > maxit || n == colors-1)
  282.         limit = maxit;
  283.      while (prev <= limit)
  284.         LogTable[prev++] = n;
  285.       }
  286.    }
  287.    LogTable[0] = 0;
  288.    if (LogFlag != -1)
  289.       for (n = 1; n < maxit; n++) /* spread top to incl unused colors */
  290.      if (LogTable[n] > LogTable[n-1])
  291.         LogTable[n] = LogTable[n-1]+1;
  292. }
  293.  
  294. long far ExpFloat14(long x) {
  295.    static float fLogTwo = (float)0.6931472;
  296.    int f;
  297.    long Ans;
  298.  
  299.    f = 23 - (int)RegFloat2Fg(RegDivFloat(x, *(long*)&fLogTwo), 0);
  300.    Ans = ExpFudged(RegFloat2Fg(x, 16), f);
  301.    return(RegFg2Float(Ans, (char)f));
  302. }
  303.  
  304. extern struct complex tmp;
  305. extern int color, colors;
  306. double TwoPi;
  307. struct complex temp, t2, BaseLog;
  308. struct complex cdegree = { 3.0, 0.0 },
  309.            croot   = { 1.0, 0.0 };
  310.  
  311. int ComplexNewtonSetup(void) {
  312.    threshold = .001;
  313.    periodicitycheck = 0;
  314.    if(param[0] != 0.0 || param[1] != 0.0 || param[2] != 0.0 ||
  315.       param[3] != 0.0) {
  316.       croot.x = param[2];
  317.       croot.y = param[3];
  318.       cdegree.x = param[0];
  319.       cdegree.y = param[1];
  320.       FPUcplxlog(&croot, &BaseLog);
  321.       TwoPi = asin(1.0) * 4;
  322.    }
  323.    return(1);
  324. }
  325.  
  326. int ComplexNewton(void) {
  327.    struct complex cd1;
  328.    double dist, mod, Angle, e2x;
  329.  
  330.    /* new = ((cdegree-1) * old**cdegree) + croot
  331.         ----------------------------------
  332.          cdegree * old**(cdegree-1)        */
  333.  
  334.    cd1.x = cdegree.x - 1.0;
  335.    cd1.y = cdegree.y;
  336.  
  337.    temp = ComplexPower(old, cd1);
  338.    FPUcplxmul(&temp, &old, &new);
  339.  
  340.    tmp.x = new.x - croot.x;
  341.    tmp.y = new.y - croot.y;
  342.    if((sqr(tmp.x) + sqr(tmp.y)) < threshold)
  343.       return(1);
  344.  
  345.    FPUcplxmul(&new, &cd1, &tmp);
  346.    tmp.x += croot.x;
  347.    tmp.y += croot.y;
  348.  
  349.    FPUcplxmul(&temp, &cdegree, &t2);
  350.    FPUcplxdiv(&tmp, &t2, &old);
  351.    new = old;
  352.    return(0);
  353. }
  354.  
  355. int ComplexBasin(void) {
  356.    struct complex cd1;
  357.    double dist, mod, Angle, e2x;
  358.  
  359.    /* new = ((cdegree-1) * old**cdegree) + croot
  360.         ----------------------------------
  361.          cdegree * old**(cdegree-1)        */
  362.  
  363.    cd1.x = cdegree.x - 1.0;
  364.    cd1.y = cdegree.y;
  365.  
  366.    temp = ComplexPower(old, cd1);
  367.    FPUcplxmul(&temp, &old, &new);
  368.  
  369.    tmp.x = new.x - croot.x;
  370.    tmp.y = new.y - croot.y;
  371.    if((sqr(tmp.x) + sqr(tmp.y)) < threshold) {
  372.       if(fabs(old.y) < .01)
  373.      old.y = 0.0;
  374.       FPUcplxlog(&old, &temp);
  375.       FPUcplxmul(&temp, &cdegree, &tmp);
  376.       mod = tmp.y/TwoPi;
  377.       color = (int)mod;
  378.       if(fabs(mod - color) > 0.5) {
  379.      if(mod < 0.0)
  380.         color--;
  381.      else
  382.         color++;
  383.       }
  384.       color += 2;
  385.       if(color < 0)
  386.      color += 128;
  387.       return(1);
  388.    }
  389.  
  390.    FPUcplxmul(&new, &cd1, &tmp);
  391.    tmp.x += croot.x;
  392.    tmp.y += croot.y;
  393.  
  394.    FPUcplxmul(&temp, &cdegree, &t2);
  395.    FPUcplxdiv(&tmp, &t2, &old);
  396.    new = old;
  397.    return(0);
  398. }
  399.  
  400. extern int Distribution, Offset, Slope;
  401. extern long con;
  402.  
  403. /*** PB, commented this out, it was unused, actual work is in prompts.c
  404. int Starfield(void) {
  405.    int c;
  406.  
  407.    plasma();
  408.    Distribution = (int)param[1];
  409.    con = (long)(param[2] / 100 * (1L << 16));
  410.    Slope = (int)param[3];
  411.    for(row = 0; row < ydots; row++) {
  412.       for(col = 0; col < xdots; col++) {
  413.      if(check_key())
  414.         return(-1);
  415.      c = getcolor(col, row);
  416.      putcolor(col, row, GausianNumber(c, colors));
  417.       }
  418.    }
  419.    return(0);
  420. }
  421.   ***/
  422.  
  423. int GausianNumber(int Probability, int Range) {
  424.    int n, r;
  425.    long Accum = 0, p;
  426.  
  427.    p = divide((long)Probability << 16, (long)Range << 16, 16);
  428.    p = multiply(p, con, 16);
  429.    p = multiply((long)Distribution << 16, p, 16);
  430.    if(!(rand() % (Distribution - (int)(p >> 16) + 1))) {
  431.       for(n = 0; n < Slope; n++)
  432.      Accum += rand();
  433.       Accum /= Slope;
  434.       r = (int)(multiply((long)Range << 15, Accum, 15) >> 14);
  435.       r = r - Range;
  436.       if(r < 0)
  437.      r = -r;
  438.       return(Range - r + Offset);
  439.    }
  440.    return(Offset);
  441. }
  442.  
  443. #endif
  444.  
  445.