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