fork download
copy 
  1. #include <string.h>
  2. #include <stdio.h>
  3.  
  4. /* triple-single functions adapted from ENS of Lyon's CRlibm, 
  5.    http://l...content-available-to-author-only...n.fr/projects/crlibm */
  6.  
  7. void Add12(float *s, float *r, float a, float b)          
  8. {
  9.   float _z, _a=a, _b=b;    
  10.   *s = _a + _b;              
  11.   _z = *s - _a;              
  12.   *r = _b - _z;   
  13. }
  14.  
  15. void Add22(float *zh, float *zl, float xh, float xl, float yh, float yl)  
  16. {
  17.   float _r,_s;                            
  18.   _r = (xh)+(yh);                          
  19.   _s = ((((xh)-_r) +(yh)) + (yl)) + (xl);  
  20.   *zh = _r+_s;                             
  21.   *zl = (_r - (*zh)) + _s;    
  22. }
  23.  
  24. #define ABS(x) (((x)>0) ? (x) : (-(x)))
  25.  
  26. void Sub22Cond(float *zh, float *zl, float xh, float xl, float yh, float yl)  
  27. {
  28.   float _r,_s;
  29.   _r = (xh)+(-yh);
  30.   _s = ((ABS(xh)) > (ABS(yh)))? 
  31.        ((xh)-_r+(-yh)+(-yl)+(xl)) : ((-yh)-_r+(xh)+(xl)+(-yl));
  32.   *zh = _r+_s;
  33.   *zl = (_r - (*zh)) + _s;
  34. }
  35.  
  36. void Add12Cond(float *s, float *r, float a, float b)
  37. {
  38.   float _z, _a=a, _b=b;
  39.   *s = _a + _b;
  40.   if (ABS(a) > ABS(b)){            
  41.     _z = *s - _a;                   
  42.     *r = _b - _z;                   
  43.   }else {                          
  44.     _z = *s - _b;                   
  45.     *r = _a - _z;}
  46. }   
  47.  
  48. void Mul12(float *rh, float *rl, float u, float v)                        
  49. {                                               
  50.   const float c  = 4097.0f; /* 1<<12 + 1 */   
  51.   float up, u1, u2, vp, v1, v2;                
  52.   float _u=u, _v=v;                            
  53.   up = _u*c;        vp = _v*c;                  
  54.   u1 = (_u-up)+up;  v1 = (_v-vp)+vp;            
  55.   u2 = _u-u1;       v2 = _v-v1;                 
  56.  
  57.   *rh = _u*_v;                                  
  58.   *rl = (((u1*v1-*rh)+(u1*v2))+(u2*v1))+(u2*v2);
  59. }
  60.  
  61.  
  62. void Mul23(float *resh, float *resm, float *resl, 
  63.            float ah, float al, float bh, float bl)                
  64. {                                                              
  65.     float _t1, _t2, _t3, _t4, _t5, _t6, _t7, _t8, _t9, _t10;  
  66.  
  67.     Mul12((resh),&_t1,(ah),(bh));                              
  68.     Mul12(&_t2,&_t3,(ah),(bl));                                
  69.     Mul12(&_t4,&_t5,(al),(bh));                                
  70.     _t6 = (al) * (bl);                                         
  71.     Add22(&_t7,&_t8,_t2,_t3,_t4,_t5);                      
  72.     Add12(&_t9,&_t10,_t1,_t6);                                   
  73.     Add22((resm),(resl),_t7,_t8,_t9,_t10);                 
  74. }
  75.  
  76. void Add33Cond(float *resh, float *resm, float *resl,
  77.                float ah, float am, float al, float bh, float bm, float bl)      
  78. {                                                            
  79.     float _t1, _t2, _t3, _t4, _t5, _t6, _t7, _t8;           
  80.  
  81.     Add12Cond(resh,&_t1,(ah),(bh));                          
  82.     Add12Cond(&_t2,&_t3,(am),(bm));                            
  83.     _t6 = (al) + (bl);                                       
  84.     Add12Cond(&_t7,&_t4,_t1,_t2);                              
  85.     _t5 = _t3 + _t4;                                         
  86.     _t8 = _t5 + _t6;                                         
  87.     Add12Cond(resm,resl,_t7,_t8);                  
  88. }
  89.  
  90. void DoRenormalize3(float *resh, float *resm, float *resl, 
  91.                     float ah, float am, float al)     
  92. {                                                      
  93.     float _t1h, _t1l, _t2l;                           
  94.  
  95.     Add12(&_t1h, &_t1l, (am), (al));                     
  96.     Add12(resh, &_t2l, (ah), (_t1h));              
  97.     Add12(resm, resl, _t2l, _t1l);           
  98. }
  99.  
  100. /* Parse an integer into a double-single. */
  101. void parse(float *h, float *l, const char *s)
  102. {
  103.   float acc = 0;
  104.  
  105.   for ( ; ; )
  106.     {
  107.       acc = 10.0f * acc + (*s - '0');
  108.  
  109.       s++;
  110.  
  111.       if (*s == 0) 
  112.         {
  113.           *h = acc;
  114.           *l = 0;
  115.           break;
  116.         }
  117.  
  118.       /* test acc rather than s; Regehr is the sort of
  119.          person to test long suites of leading zeroes: */
  120.       if (acc >= 100000.0f)
  121.         {
  122.           float acc2 = 0;
  123.           float factor = 1;
  124.           for ( ; ; )
  125.             {
  126.               /* XXX I'll bet this can be simplified by computing
  127.                  acc *= 10.0f; instead of factor *= 10.0f; and a simple
  128.                  Add12 of acc and acc2 at the end, for the right threshold
  129.                  instead of 100000. */
  130.               factor *= 10.0f;
  131.               acc2 = 10.0f * acc2 + (*s - '0');
  132.               s++;
  133.               if (*s == 0) break;
  134.             }
  135.           Mul12(h, l, acc, factor);
  136.           Add22(h, l, *h, *l, acc2, 0);
  137.           break;
  138.         } 
  139.     }
  140. }
  141.  
  142. int main(int c, char **v)
  143. {
  144.   if (c != 7)
  145.     {
  146.       printf("Usage!\n");
  147.       return 1;
  148.     }
  149.   float xh1, xl1, yh1, yl1, xh2, xl2, yh2, yl2, xh3, xl3, yh3, yl3;
  150.   parse(&xh1, &xl1, v[1]);
  151.   parse(&yh1, &yl1, v[2]);
  152.   parse(&xh2, &xl2, v[3]);
  153.   parse(&yh2, &yl2, v[4]);
  154.   parse(&xh3, &xl3, v[5]);
  155.   parse(&yh3, &yl3, v[6]);
  156.  
  157.   float dxh1, dxl1, dyh1, dyl1, dxh2, dxl2, dyh2, dyl2, dxh3, dxl3, dyh3, dyl3;
  158.  
  159.   Sub22Cond(&dxh1, &dxl1, xh1, xl1, xh2, xl2);
  160.   Sub22Cond(&dyh1, &dyl1, yh1, yl1, yh2, yl2);
  161.   Sub22Cond(&dxh2, &dxl2, xh2, xl2, xh3, xl3);
  162.   Sub22Cond(&dyh2, &dyl2, yh2, yl2, yh3, yl3);
  163.   Sub22Cond(&dxh3, &dxl3, xh3, xl3, xh1, xl1);
  164.   Sub22Cond(&dyh3, &dyl3, yh3, yl3, yh1, yl1);
  165.  
  166.   float sdxh1, sdxm1, sdxl1, sdyh1, sdym1, sdyl1;
  167.   float sdxh2, sdxm2, sdxl2, sdyh2, sdym2, sdyl2;
  168.   float sdxh3, sdxm3, sdxl3, sdyh3, sdym3, sdyl3;
  169.  
  170.   // Are vectors colinear? Use cross-product v1x*v2y == v2x*v1y 
  171.  
  172.   Mul23(&sdxh1, &sdxm1, &sdxl1, dxh1, dxl1, dyh2, dyl2);
  173.   Mul23(&sdxh2, &sdxm2, &sdxl2, dxh2, dxl2, dyh1, dyl1);
  174.   DoRenormalize3(&sdxh1, &sdxm1, &sdxl1, sdxh1, sdxm1, sdxl1);
  175.   DoRenormalize3(&sdxh2, &sdxm2, &sdxl2, sdxh2, sdxm2, sdxl2);
  176.  
  177.   /* XXX Compare. Products may be acurate to 2^-68 or so
  178.      if we interpolate from CRlibm's error computations for
  179.      double. We need 2^62. Can we compare for exact equality
  180.      or do we need to COMPARE UP TO EPSILON? */
  181.  
  182.   // Now reuse the six triple-singles to compute squares of lengths of sides
  183.  
  184.   Mul23(&sdxh1, &sdxm1, &sdxl1, dxh1, dxl1, dxh1, dxl1);
  185.   Mul23(&sdyh1, &sdym1, &sdyl1, dyh1, dyl1, dyh1, dyl1);
  186.   Mul23(&sdxh2, &sdxm2, &sdxl2, dxh2, dxl2, dxh2, dxl2);
  187.   Mul23(&sdyh2, &sdym2, &sdyl2, dyh2, dyl2, dyh2, dyl2);
  188.   Mul23(&sdxh3, &sdxm3, &sdxl3, dxh3, dxl3, dxh3, dxl3);
  189.   Mul23(&sdyh3, &sdym3, &sdyl3, dyh3, dyl3, dyh3, dyl3);
  190.  
  191.   DoRenormalize3(&sdxh1, &sdxm1, &sdxl1, sdxh1, sdxm1, sdxl1);
  192.   DoRenormalize3(&sdxh2, &sdxm2, &sdxl2, sdxh2, sdxm2, sdxl2);
  193.   DoRenormalize3(&sdxh3, &sdxm3, &sdxl3, sdxh3, sdxm3, sdxl3);
  194.   DoRenormalize3(&sdyh1, &sdym1, &sdyl1, sdyh1, sdym1, sdyl1);
  195.   DoRenormalize3(&sdyh2, &sdym2, &sdyl2, sdyh2, sdym2, sdyl2);
  196.   DoRenormalize3(&sdyh3, &sdym3, &sdyl3, sdyh3, sdym3, sdyl3);
  197.  
  198.   Add33Cond(&sdxh1, &sdxm1, &sdxl1, sdxh1, sdxm1, sdxl1, sdyh1, sdym1, sdyl1);
  199.   Add33Cond(&sdxh2, &sdxm2, &sdxl2, sdxh2, sdxm2, sdxl2, sdyh2, sdym2, sdyl2);
  200.   Add33Cond(&sdxh3, &sdxm3, &sdxl3, sdxh3, sdxm3, sdxl3, sdyh3, sdym3, sdyl3);
  201.  
  202.   DoRenormalize3(&sdxh1, &sdxm1, &sdxl1, sdxh1, sdxm1, sdxl1);
  203.   DoRenormalize3(&sdxh2, &sdxm2, &sdxl2, sdxh2, sdxm2, sdxl2);
  204.   DoRenormalize3(&sdxh3, &sdxm3, &sdxl3, sdxh3, sdxm3, sdxl3);
  205.  
  206.   // We have the squares of the lengths of each side, test them for equality:
  207.  
  208.   float n = 0.0f;
  209.  
  210.   if (sdxh1 == sdxh2 && sdxm1 == sdxm2 && sdxl1 == sdxl2)
  211.     // See XXX above.
  212.     n = 1.0f;
  213.  
  214.   if (sdxh1 == sdxh3 && sdxm1 == sdxm3 && sdxl1 == sdxl3)
  215.     n++;
  216.  
  217.   if (sdxh2 == sdxh3 && sdxm2 == sdxm3 && sdxl2 == sdxl3)
  218.     n++;
  219.  
  220.   if (n == 1.0f) printf("isosceles");
  221.   else if (n == 3.0f) printf("equilateral");
  222.   else printf("scalene");
  223.  
  224.   printf(" ");
  225.  
  226. #define GT(n, p) \
  227.   (sdxh##n > sdxh##p ||  \
  228.       (sdxh##n == sdxh##p && (sdxm##n > sdxm##p ||  \
  229.                           (sdxm##n == sdxm##p && sdxl##n > sdxl##p))))
  230.  
  231.   float tmp;
  232.  
  233. #define SWAP(n, p) \
  234.   tmp = sdxh##n; sdxh##n = sdxh##p; sdxh##p = tmp; \
  235.   tmp = sdxm##n; sdxm##n = sdxm##p; sdxm##p = tmp; \
  236.   tmp = sdxl##n; sdxl##n = sdxl##p; sdxl##p = tmp;
  237.  
  238.   if (GT(2, 1))
  239.     {
  240.       SWAP(1, 2)
  241.     }
  242.  
  243.   if (GT(3, 1))
  244.     {
  245.       SWAP(1, 3)
  246.     }
  247.  
  248.   /* (sdxh1, sdxm1, sdxl1) contains the longest square of side length,
  249.      the other two triple-singles contain the other two sides. */
  250.  
  251.   Add33Cond(&sdxh2, &sdxm2, &sdxl2, sdxh2, sdxm2, sdxl2, sdxh3, sdxm3, sdxl3);
  252.   DoRenormalize3(&sdxh2, &sdxm2, &sdxl2, sdxh2, sdxm2, sdxl2);
  253.  
  254.   if (GT(2, 1))
  255.     printf("acute");
  256.   else if (GT(1, 2))
  257.     printf("obtuse");
  258.   else printf("right");
  259.  
  260.   printf("\n");
  261.   return 0;
  262. }
  263.  
Runtime error #stdin #stdout 0s 2252KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Usage!
https://ideone.com/1uBk3Z
language:
C99 (gcc 8.3)
created:
10 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!