home *** CD-ROM | disk | FTP | other *** search
/ APDL Public Domain 1 / APDL_PD1A.iso / fractal / tesseract / !Tesseract / source / c / tessera < prev   
Encoding:
Text File  |  1991-12-18  |  18.3 KB  |  864 lines
  1. /*  > Tessera
  2.  *  David McQuillan  Aug 83
  3.  */
  4.  
  5. /* #define TRACE */
  6. /* #define TIMING */
  7.  
  8. #include <setjmp.h>
  9. #include <signal.h>
  10. #include <math.h>
  11. #include <stdio.h>
  12. #include <string.h>
  13. #include <stdlib.h>
  14.  
  15. #ifdef TIMING
  16. #include <time.h>
  17. #endif
  18.  
  19. #include "swis.h"
  20. #include "kernel.h"
  21. #include "ploth.h"
  22.  
  23. typedef enum{FALSE, TRUE} bool;
  24. typedef unsigned char small;
  25.  
  26. static small v1[48] = {
  27.      0,  1,  2,  3,  4,  5,  6,  7,
  28.      8,  9, 10, 11, 12, 13, 14, 15,
  29.      0,  0,  0,  0,  3,  3,  3,  3,
  30.      5,  5,  5,  5,  6,  6,  6,  6,
  31.      9,  9,  9,  9, 10, 10, 10, 10,
  32.     12, 12, 12, 12, 15, 15, 15, 15};
  33.  
  34. static small v2[48] = {
  35.      0,  1,  2,  3,  4,  5,  6,  7,
  36.      8,  9, 10, 11, 12, 13, 14, 15,
  37.      1,  2,  4,  8,  1,  2,  7, 11,
  38.      4,  7,  1, 13,  7,  4,  2, 14,
  39.      8, 11, 13,  1, 11,  8, 14,  2,
  40.     13, 14,  8,  4, 14, 13, 11,  7};
  41.  
  42. static small p2[4] = {1, 2, 4, 8};
  43.  
  44. static small dimcol[4][2] = {{8,6},{9,14},{10,13},{11,12}};
  45.  
  46. static unsigned colour[16];
  47.  
  48. typedef int ifloat;
  49.  
  50. #define ISCALE 30
  51. #define IMULT (1<<ISCALE)
  52.  
  53. #define ICONV(x) ((ifloat)(IMULT * (x)))
  54. #define ITIMES(x,y) MulScale(x, y, ISCALE)
  55. #define IDIV(x,y) DivScale(x, y, ISCALE)
  56. #define ITIMES0(x,y) ((x)>>(ISCALE/2))*((y)>>(ISCALE/2))
  57.  
  58. typedef int ffloat;
  59.  
  60. #define FSCALE 20
  61. #define FMULT (1<<FSCALE)
  62.  
  63. #define FCONV(x) ((ffloat)(FMULT * (x)))
  64. #define FTIMES(x,y) MulScale(x, y, FSCALE)
  65. #define FTIMES0(x,y) (((x)*(y)) >> FSCALE)
  66.  
  67. #define ITOF(x) ((x) >> (ISCALE-FSCALE))
  68.  
  69. #define R2 1.414213962373095
  70. #define PI 3.141592653589793
  71.  
  72. #define FPI2 (ffloat)(2*PI*FMULT + 0.5)
  73. #define FPI (ffloat)(PI*FMULT + 0.5)
  74.  
  75. #define QSCALE 8
  76. #define QMULT (1 << QSCALE)
  77. #define QSIZE (int)(PI/2*QMULT + 2)
  78.  
  79. static ifloat icos1[QSIZE], icos2[QMULT];
  80. static ifloat isin1[QSIZE], isin2[QMULT];
  81.  
  82. static ffloat recip[1000];
  83.  
  84. static ifloat t[4][4] = {
  85.     IMULT, 0, 0, 0,
  86.     0, IMULT, 0, 0,
  87.     0, 0, IMULT, 0,
  88.     0, 0, 0, IMULT};
  89.  
  90. #define SIZE 512
  91. #define EPS 0.1
  92.  
  93. static void isincos(ffloat t, ifloat *s, ifloat *c);
  94. static ifloat iasin(ifloat x);
  95. static void ptor(int x, int y, ffloat *t, int *r);
  96. static double max_width(double hd, double rd);
  97. static void rot(int a, int b, ifloat s, ifloat c);
  98.  
  99.  
  100. /*
  101.  *  interrupt
  102.  */
  103.  
  104. typedef struct
  105. {
  106.     jmp_buf js;
  107. } jmp_struct;
  108.  
  109. static jmp_struct err_jmp;
  110.  
  111. static void interrupt(int type)
  112. {
  113.     switch (type)
  114.     {
  115.         case SIGINT:
  116.         fprintf(stderr, "\nEscape pressed\n");
  117.         break;
  118.     default:
  119.         fprintf(stderr, "\nSignal type %d\n", type);
  120.     }
  121.     longjmp(err_jmp.js, 1);
  122. }
  123.  
  124. static void vdu(int c)
  125. {
  126.     _kernel_oswrch(c);
  127. }
  128.  
  129. static void palette(int l, int p, int r, int g, int b)
  130. {
  131.     vdu(19);
  132.     vdu(l);
  133.     vdu(p);
  134.     vdu(r);
  135.     vdu(g);
  136.     vdu(b);
  137. }
  138.  
  139. #define SolidBoth       0x00
  140. #define TriangleFill    0x50
  141. #define Circle          0x90
  142. #define CircleFill      0x98
  143. #define DrawRelFore     1
  144. #define MoveCursorAbs   4
  145. #define DrawAbsFore     5
  146.  
  147. static void mouse(int *x, int *y, int *b)
  148. {
  149.     _kernel_swi_regs regs;
  150.  
  151.     _kernel_swi(OS_Mouse, ®s, ®s);
  152.     *x = regs.r[0];
  153.     *y = regs.r[1];
  154.     *b = regs.r[2];
  155. }
  156.  
  157.  
  158. int main(void)
  159. {
  160.     int i, j, k;
  161.     bool reveal = FALSE;
  162.     double fact = 1;
  163.     double hd = 4;
  164.     double rd = 6;
  165.     double scmax;
  166.     double scmult;
  167.     ffloat ihd;
  168.     ffloat ird;
  169.     ffloat iscmult;
  170.     int old_x, old_y;
  171.     ffloat r1, t1;
  172.     ifloat st1, ct1;
  173.     ifloat sr1, cr1;
  174.     bool bp[16];
  175.     ffloat p[16][4];
  176.     int ff[16], px[16], py[16];
  177.     int pz[48];
  178.     small seq[48];
  179.  
  180. #ifdef TIMING
  181.     clock_t clockval;
  182.     int iter = 0;
  183. #endif
  184.  
  185.     if (setjmp(err_jmp.js))
  186.     {
  187.         plot_end();
  188.         exit(0);
  189.     }
  190.     signal(SIGINT, *interrupt);
  191.  
  192.     /* set up trig values */
  193.  
  194.     for (i=0; i < QSIZE; ++i)
  195.     {
  196.         double qi = i*1.0 / (1 << QSCALE);
  197.  
  198.         isin1[i] = ICONV(sin(qi));
  199.         icos1[i] = ICONV(cos(qi));
  200.     }
  201.  
  202.     for (i=0; i < QMULT; ++i)
  203.     {
  204.         double qi = i*1.0 / (1 << (2*QSCALE));
  205.  
  206.         isin2[i] = ICONV(sin(qi));
  207.         icos2[i] = ICONV(cos(qi));
  208.     }
  209.  
  210.     recip[0] = FMULT;
  211.     for (i=1; i < 999; i++)
  212.     {
  213.         recip[i] = FMULT/i;
  214.     }
  215.  
  216.     /* Set up Screen */
  217.     /* two banks */
  218.     /* origin at centre */
  219.     /* mouse on */
  220.  
  221.     if (plot_init() != 0)
  222.     {
  223.         plot_end();
  224.         return 0;
  225.     }
  226.  
  227.     /* Set up initial conditions */
  228.  
  229.     scmax = fact * max_width(hd, rd);
  230.     scmult = 1 / scmax;
  231.     ihd = FCONV(hd);
  232.     ird = FCONV(rd);
  233.     iscmult = FCONV(scmult);
  234.  
  235.     /* Set up colours */
  236.  
  237.     palette( 0, 16, 180, 221, 221);   /* centre of lines */
  238.     palette( 1, 16, 221, 221, 221);   /* light background circle */
  239.     palette( 2, 16, 255, 255, 255);   /* background colour */
  240.     palette( 7, 16,   0,   0,   0);
  241.  
  242.     palette( 6, 16,  80,  80,  80);
  243.     palette( 8, 16, 255,   0,   0);
  244.     palette( 9, 16,   0, 255,   0);
  245.     palette(10, 16,   0,   0, 255);
  246.     palette(11, 16,   0,   0,   0);
  247.     palette(12, 16, 255, 255, 255);
  248.     palette(13, 16,  80,  80,  80);
  249.     palette(14, 16,  80,  80,  80);
  250.  
  251.     palette(15, 16,  80,   0,   0);   /* outside of lines */
  252.  
  253.     for (i = 0; i <= 15; i++)
  254.     {
  255.         unsigned col = 0;
  256.  
  257.         for (j = 0; j <= 3; j++)
  258.         {
  259.             col |= dimcol[j][(i & (1<<j)) ? 1 : 0] << (4*j);
  260.         }
  261.         col = ((col & 0xFF00) << 8) | (col & 0xFF);
  262.         colour[i] = col | (col << 8);
  263.     }
  264.  
  265.     for (i = 0; i <= 47; i++)
  266.     {
  267.         seq[i] = i;
  268.     }
  269.  
  270.     /* Set up Initial Position */
  271.  
  272.     for (i = 0; i <= 5; i++)
  273.     {
  274.         ifloat s, c;
  275.  
  276.         s = ICONV(sin(0.1));
  277.         c = ICONV(cos(0.1));
  278.  
  279.         for (j = 1; j <= 3; j++)
  280.         {
  281.             rot(j, 0, s, c);
  282.         }
  283.     }
  284.  
  285.     /* button is pressed on entry ! */
  286.  
  287.     {
  288.         int b;
  289.  
  290.         mouse(&old_x, &old_y, &b);
  291.  
  292.         ptor(old_x, old_y, &r1, &t1);
  293.         r1 = r1 * FCONV(PI/800);
  294.  
  295.         isincos(r1, &sr1, &cr1);
  296.         isincos(t1, &st1, &ct1);
  297.     }
  298.  
  299. #ifdef TIMING
  300.     /*  init timer */
  301.  
  302.     clockval = clock();
  303. #endif
  304.  
  305.     /* Loop Displaying Tesseract */
  306.  
  307.     do
  308.     {
  309.         /* Get mouse input */
  310.  
  311.         {
  312.             int b, x, y;
  313.             ffloat r2, t2;
  314.             ifloat st2, ct2;
  315.             ifloat sr2, cr2;
  316.  
  317.             mouse(&x, &y, &b);
  318.             b &= 7;
  319.  
  320.             ptor(x, y, &r2, &t2);
  321.             r2 = r2 * FCONV(PI/800);
  322.  
  323.             isincos(t2, &st2, &ct2);
  324.             isincos(r2, &sr2, &cr2);
  325.  
  326.             if (b == 6)
  327.             {
  328.                 reveal = TRUE;
  329.             }
  330.             else if (b == 3)
  331.             {
  332.                 reveal = FALSE;
  333.             }
  334.             else if ((b == 1 || b == 2 || b == 4)
  335.                  && (x != old_x || y != old_y))
  336.             {
  337.                 ffloat r3, t3;
  338.                 ifloat x2, y2, z2;
  339.                 ifloat x3, y3, z3;
  340.                 bool in_out = ((b & 2) != 0);
  341.                 ifloat st3, ct3;
  342.                 ifloat sr3, cr3;
  343.  
  344.                 /* turn 1 on z axis till on xz plane */
  345.  
  346.                 x2 = ITIMES0(sr2, ITIMES0(ct2, ct1) + ITIMES0(st2, st1));
  347.                 y2 = ITIMES0(sr2, ITIMES0(st2, ct1) - ITIMES0(ct2, st1));
  348.                 z2 = cr2;
  349.  
  350.                 /* turn 1 on y axis till z == 1 */
  351.  
  352.                 x3 = ITIMES0(x2, cr1) - ITIMES0(z2, sr1);
  353.                 y3 = y2;
  354.                 z3 = ITIMES0(x2, sr1) + ITIMES0(z2, cr1);
  355.  
  356.                 /* turn 2 on z axis till y == 0 */
  357.  
  358.                 ptor(x3, y3, &r3, &t3);
  359.  
  360.                 /* turn till 1 goes to 2 on y axis */
  361.  
  362.                 r3 = ITOF(iasin(r3));
  363.  
  364.                 isincos(t3, &st3, &ct3);
  365.                 isincos(r3, &sr3, &cr3);
  366.  
  367.                 if (in_out)
  368.                 {
  369.                     double nfact, nhd, nrd;
  370.                     double lmax;
  371.  
  372.                     nfact  = fact * (1 + EPS*ct1/IMULT*ct3/IMULT*r3/FMULT);
  373.                     nhd = hd * (1 - EPS*st1/IMULT*ct3/IMULT*r3/FMULT);
  374.                     nrd = rd * (1 - EPS*r1/FMULT*st3/IMULT*r3/FMULT);
  375.  
  376.                     lmax = nfact * max_width(nhd, nrd);
  377.  
  378.                     if (lmax < scmax && lmax > scmax/100)
  379.                     {
  380.                         fact  = nfact;
  381.                         hd = nhd;
  382.                         rd = nrd;
  383.  
  384.                         scmult = fact / scmax;
  385.  
  386.                         ihd = FCONV(hd);
  387.                         ird = FCONV(rd);
  388.                         iscmult = FCONV(scmult);
  389.                     }
  390.                 }
  391.                 else
  392.                 {
  393.                     bool hyper = ((b & 4) != 0);
  394.  
  395.                     /* position as described above */
  396.  
  397.                     rot(1, 2, -st1, ct1);
  398.  
  399.                     rot(3, 1, -sr1, cr1);
  400.  
  401.                     rot(1, 2, -st3, ct3);
  402.  
  403.                     /* if hyper then cross-product to x axis */
  404.  
  405.                     if (hyper)
  406.                         rot(1, 0, -sr3, cr3);
  407.                     else
  408.                         rot(3, 1, sr3, cr3);
  409.  
  410.                     /* reverse manipulations */
  411.  
  412.                     rot(1, 2, st3, ct3);
  413.  
  414.                     rot(3, 1, sr1, cr1);
  415.  
  416.                     rot(1, 2, st1, ct1);
  417.                 }
  418.             }
  419.  
  420.             old_x = x;
  421.             old_y = y;
  422.             t1 = t2;
  423.             r1 = r2;
  424.             ct1 = ct2;
  425.             st1 = st2;
  426.             cr1 = cr2;
  427.             sr1 = sr2;
  428.         }
  429.  
  430.         /* Generate Vertices */
  431.  
  432.         for (j = 0; j <= 3; j++)
  433.         {
  434.             p[0][j] = ITOF(t[0][j]) + ITOF(t[1][j])
  435.                     + ITOF(t[2][j]) + ITOF(t[3][j]);
  436.  
  437.             for (k = 0; k <= 3; k++)
  438.             {
  439.                 ffloat x = p2[k];
  440.                 ffloat r = 2 * ITOF(t[k][j]);
  441.  
  442.                 for (i=0; i < x; i++)
  443.                     p[x+i][j] = p[i][j] - r;
  444.             }
  445.         }
  446.  
  447.         /* See which cubes and so which vertices are visible */
  448.  
  449.         {
  450.             int vism = 0;
  451.             int posm = 0;
  452.  
  453.             if (reveal)
  454.             {
  455.                 for (j = 0; j <= 15; j++)
  456.                 {
  457.                     bp[j] = TRUE;
  458.                 }
  459.             }
  460.             else
  461.             {
  462.                 for (k = 0; k <= 3; k++)
  463.                 {
  464.                     ffloat tk0 = ITOF(t[k][0]);
  465.  
  466.                     if (FTIMES((tk0 > 0 ? tk0 : -tk0), ihd) > FMULT)
  467.                         vism |= p2[k];
  468.  
  469.                     if (tk0 > 0)
  470.                         posm |= p2[k];
  471.                 }
  472.  
  473.                 for (j = 0; j <= 15; j++)
  474.                 {
  475.                     bp[j] = (((j ^ posm) & vism) != 0);
  476.                 }
  477.             }
  478.         }
  479.  
  480.         /* Generate X Perspective */
  481.  
  482.         {
  483.             for (k = 0; k <= 15; k++)
  484.             {
  485.                 int g = iscmult / 
  486.                     ((FTIMES(ihd-p[k][0], ird) - p[k][3]) / SIZE);
  487.  
  488.                 ff[k] = g / 4;
  489.                 px[k] = FTIMES0(g, p[k][1]);
  490.                 py[k] = FTIMES0(g, p[k][2]);
  491.                 pz[k] = FTIMES0(g, p[k][3]);
  492.             }
  493.         }
  494.  
  495.         /*
  496.          * Order vertices and lines by distance
  497.          * Not perfect 3D but seems to work okay mostly
  498.          * Should be almost in order already
  499.          */
  500.  
  501.         {
  502.             small *pv1 = v1;
  503.             small *pv2 = v2;
  504.             int *ppz = pz;
  505.  
  506.             for (i = 16; i <= 47; i++ )
  507.                 ppz[i] = (ppz[pv1[i]] + ppz[pv2[i]]) >> 1;
  508.         }
  509.  
  510.         {
  511.             int *ppz = pz;
  512.             small *pseq = seq;
  513.             int sk = pseq[0];
  514.             int pzk = ppz[sk];
  515.  
  516.             for (k = 1; k <= 47; k++)
  517.             {
  518.                 int sn;
  519.                 int pzn;
  520.  
  521.                 if (pzk > (pzn = ppz[sn = pseq[k]]))
  522.                 {
  523.                     int n = k;
  524.  
  525.                     do
  526.                     {
  527.                         pseq[n] = sk;
  528.                         n --;
  529.                     } while (n >= 1 && ppz[sk = pseq[n-1]] > pzn);
  530.  
  531.                     pseq[n] = sn;
  532.  
  533.                     sn = pseq[k];
  534.                     pzn = ppz[sn];
  535.                 }
  536.  
  537.                 sk = sn;
  538.                 pzk = pzn;
  539.             }
  540.         }
  541.  
  542.         /* Start New Screen */
  543.  
  544.         plot_begin();
  545.  
  546. #ifdef TIMING
  547.         if (iter == 100)
  548.         {      
  549.             clock_t oldclock = clockval;
  550.  
  551.             vdu(31);
  552.             vdu(0);
  553.             vdu(0);
  554.             clockval = clock();
  555.  
  556.             printf("%4d\n", (100*CLOCKS_PER_SEC)/(clockval-oldclock));
  557.             iter=0;
  558.         }
  559.         iter++;
  560. #endif
  561.  
  562.         /* Draw Vertex Spheres and lines */
  563.  
  564.         {
  565.             for (k = 0; k <= 47; k++)
  566.             {
  567.                 int i = seq[k];
  568.  
  569.                 if (i <= 15 && bp[i])
  570.                 {
  571.                     plot_circlefill(px[i], py[i], ff[i], colour[i]);
  572.                 }
  573.                 else if (i > 15 && bp[v1[i]] && bp[v2[i]])
  574.                 {
  575.                     int a = v1[i];
  576.                     int b = v2[i];
  577.                     int xa = px[a];
  578.                     int xb = px[b];
  579.                     int ya = py[a];
  580.                     int yb = py[b];
  581.                     int za = pz[a];
  582.                     int zb = pz[b];
  583.                     int fa = ff[a];
  584.                     int fb = ff[b];
  585.                     int dx = xb - xa;
  586.                     int dy = yb - ya;
  587.                     int dz = zb - za;
  588.  
  589.                     int dr2 = dx * dx + dy * dy;
  590.                     int ds2 = dr2 + dz * dz;
  591.                     int ds;
  592.                     int w1, w2;
  593.  
  594.                     ffloat ids;
  595.  
  596.                     ds = 4 * (fa + fb);
  597.                     ds = ((ds + FTIMES0(ds2, recip[ds])) >> 1);
  598.                     ds = ((ds + FTIMES0(ds2, recip[ds])) >> 1);
  599.  
  600.                     ids = recip[ds+7];
  601.                     w1 = fa*ids;
  602.                     w2 = fb*ids;
  603.                     xa += FTIMES0(dx, w1);
  604.                     ya += FTIMES0(dy, w1);
  605.                     xb -= FTIMES0(dx, w2);
  606.                     yb -= FTIMES0(dy, w2);
  607.  
  608.                     plot_line(xa, ya, xb, yb);
  609. /*
  610.                     {
  611.                     int dr;
  612.                     int qxa, qya, qxb, qyb;
  613.                     ffloat idr;
  614.  
  615.                     dr = dx * dx + dy * dy;
  616.                     dr = (dx > 0 ? dx : -dx) + (dy > 0 ? dy : -dy) + 1;
  617.                     dr = (dr + dr2*recip[dr]/FMULT) >> 1;
  618.                     idr = recip[dr];
  619.  
  620.                     qxa =  fa*dy*idr/FMULT >> 3;
  621.                     qya = -fa*dx*idr/FMULT >> 3;
  622.                     qxb =  fb*dy*idr/FMULT >> 3;
  623.                     qyb = -fb*dx*idr/FMULT >> 3;
  624.  
  625.                     gcol(0, 5);
  626.                     plot(MoveCursorAbs, xa+qxa, ya+qya);
  627.                     plot(MoveCursorAbs, xa-qxa, ya-qya);
  628.                     plot(DrawAbsFore+TriangleFill, xb+qxb, yb+qyb);
  629.                     plot(DrawAbsFore+TriangleFill, xb-qxb, yb-qyb);
  630.  
  631.                     gcol(0, 2);
  632.                     plot(MoveCursorAbs, xa, ya);
  633.                     plot(DrawAbsFore+SolidBoth, xb, yb);
  634.                     }
  635. */
  636.                 }
  637.             }
  638.         }
  639.  
  640.         /* screen done - display it */
  641.         /* escape pressed if 1 returned */
  642.     }
  643.     while (plot_show() == 0);
  644.  
  645.     plot_end();
  646.  
  647.     return 0;
  648. }
  649.  
  650. /*
  651.  * Max width of Tesseract
  652.  * so all checks can be removed
  653.  * that points are on screen if
  654.  * desired.
  655.  */
  656.  
  657. static double max_width(double hd, double rd)
  658. {
  659.     double r1 = 2.25;    /* diagonal plus ball */
  660.     double r2 = r1 / sqrt(hd*hd - r1*r1);
  661.     double r3 = r2 / sqrt(rd*rd - r2*r2);
  662.  
  663.     return r3;
  664. }
  665.  
  666. /*
  667.  * Move on hypersphere
  668.  * Parallel displacement X, Y or Z
  669.  * or Rotation X, Y or Z
  670.  */
  671.  
  672. static void rot(int a, int b, ifloat s, ifloat c)
  673. {
  674.     int i;
  675.  
  676.     for (i = 0; i <= 3; i++)
  677.     {
  678. /*
  679.         ifloat w0 = t[i][b];
  680.         ifloat w1 = t[i][a];
  681.  
  682.         t[i][b] =  ITIMES(c, w0) + ITIMES(s, w1);
  683.         t[i][a] = -ITIMES(s, w0) + ITIMES(c, w1);
  684. */
  685.         CrossMul(s, c, &t[i][b], &t[i][a]);
  686.     }
  687. }
  688.  
  689. /* sin and cos of angle */
  690.  
  691. static void isincos(ffloat t, ifloat *s, ifloat *c)
  692. {
  693.     ffloat tp, tq, t1, t2;
  694.     ifloat s1, c1, s2, c2, ss, cc;
  695.  
  696.     while (t <= -FPI) t += FPI2;
  697.     while (t > FPI) t -= FPI2;
  698.  
  699.     tp = (t >= 0 ? t : -t);
  700.     tq = (tp >= (FPI/2) ? FPI-tp : tp);
  701.  
  702. #if (FSCALE > 2*QSCALE)
  703.     t1 = tq >> (FSCALE - 2*QSCALE);
  704. #else
  705.     t1 = tq << (2*QSCALE - FSCALE);
  706. #endif
  707.  
  708.     t2 = t1 & (QMULT-1);
  709.     t1 = t1 >> QSCALE;
  710.  
  711.     s1 = isin1[t1];
  712.     c1 = icos1[t1];
  713.     s2 = isin2[t2];
  714.     c2 = icos2[t2];
  715.     ss = ITIMES(s1, c2) + ITIMES(c1, s2);
  716.     cc = ITIMES(c1, c2) - ITIMES(s1, s2);
  717.     if (t < 0) ss = -ss;
  718.     if (tp >= (FPI/2)) cc = -cc;
  719.     *s = ss;
  720.     *c = cc;
  721. }
  722.  
  723.  
  724. static ifloat sqrt2[] = {
  725.     ICONV( 0.00098843065718),
  726.     ICONV( 0.00079479950957),
  727.     ICONV(-0.0035890535377),
  728.     ICONV( 0.011028809744),
  729.     ICONV(-0.044195203560),
  730.     ICONV( 0.35355338194),
  731.     ICONV( 1.41421356237)
  732.     };
  733.  
  734. static ifloat sqrt1[] = {
  735.     ICONV( 0.0135199291026),
  736.     ICONV(-0.0226657767832),
  737.     ICONV( 0.0278720776889),
  738.     ICONV(-0.0389582788321),
  739.     ICONV( 0.0624811144548),
  740.     ICONV(-0.125001503933),
  741.     ICONV( 0.5),
  742.     ICONV( 1.0)
  743.     };
  744.  
  745. static ifloat atn[] = {
  746.     ICONV( 0.0805374449538),
  747.     ICONV(-0.138776856032),
  748.     ICONV( 0.199777106478),
  749.     ICONV(-0.333329491539),
  750.     ICONV( 1.0)
  751.     };
  752.  
  753. static ifloat asn[] = {
  754.     ICONV( 0.042163199048),
  755.     ICONV( 0.024181311049),
  756.     ICONV( 0.045470025998),
  757.     ICONV( 0.074953002686),
  758.     ICONV( 0.16666752422),
  759.     ICONV( 1.0)
  760.     };
  761.  
  762. /* Low precision polynomial evaluate */
  763.  
  764. static ifloat i0polevl(ifloat x, ifloat *coef, int N)
  765. {
  766.     ifloat ans = *coef++;
  767.  
  768.     do
  769.         ans = ITIMES0(ans, x) + *coef++;
  770.     while (--N);
  771.  
  772.     return ans;
  773. }
  774.  
  775. /* low precision asin */
  776.  
  777. static ifloat iasin(ifloat x)
  778. {
  779.     ifloat x2 = ITIMES0(x, x);
  780.  
  781.     return ITIMES(x, i0polevl(x2, asn, 5));
  782. }
  783.  
  784. /* low precision polar to rect */
  785.  
  786. static void ptor(int x, int y, int *r, ffloat *t)
  787. {
  788.     int w;
  789.     ifloat d, r2, t1;
  790.     int code = 0;
  791.  
  792.     if (x < 0)
  793.     {
  794.         code = 2;
  795.         x = -x;
  796.     }
  797.  
  798.     if (y < 0)
  799.     {
  800.        code |= 1;
  801.        y = -y;
  802.     }
  803.  
  804.     if (x == 0)
  805.     {
  806.         *r = y;
  807.  
  808.         if (code & 1)
  809.             *t = FCONV(-PI/2);
  810.         else if (y == 0)
  811.             *t = 0;
  812.         else
  813.             *t = FCONV(PI/2);
  814.  
  815.         return;
  816.     }
  817.  
  818.     if (y == 0)
  819.     {
  820.         *r = x;
  821.  
  822.         if (code & 2)
  823.             *t = FCONV(PI);
  824.         else
  825.             *t = 0;
  826.  
  827.         return;
  828.     }
  829.  
  830.     if (x < y)
  831.     {
  832.         code |= 4;
  833.         w = y;
  834.         y = x;
  835.         x = w;
  836.     }
  837.  
  838.     d = IDIV(y, x);
  839.     r2 = ITIMES0(d, d);
  840.  
  841.     if (r2 > ICONV(R2-1))
  842.     {
  843.         r2 = r2 - IMULT;
  844.         *r = ITIMES(x, i0polevl(r2, sqrt2, 6));
  845.  
  846.         code |= 8;
  847.         d = IDIV(d-IMULT, d+IMULT);
  848.         r2 = ITIMES0(d, d);
  849.     }
  850.     else
  851.     {
  852.         *r = ITIMES(x, i0polevl(r2, sqrt1, 7));
  853.     }
  854.  
  855.     t1 = ITOF(ITIMES(d, i0polevl(r2, atn, 4)));
  856.  
  857.     if (code & 8) t1 += FCONV(PI/4);
  858.     if (code & 4) t1 = FCONV(PI/2) - t1;
  859.     if (code & 2) t1 = FCONV(PI) - t1;
  860.     if (code & 1) t1 = -t1;
  861.  
  862.     *t = t1;
  863. }
  864.