home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / compsrcs / misc / volume08 / interval.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1991-08-27  |  11.5 KB  |  429 lines

  1. From decwrl!sun-barr!cs.utexas.edu!uunet!allbery Sat Oct 28 16:26:05 PDT 1989
  2. Article 1145 of comp.sources.misc:
  3. Path: decwrl!sun-barr!cs.utexas.edu!uunet!allbery
  4. From: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  5. Newsgroups: comp.sources.misc
  6. Subject: v08i098: Interval Arithmetic in C++
  7. Message-ID: <70936@uunet.UU.NET>
  8. Date: 28 Oct 89 18:14:26 GMT
  9. Sender: allbery@uunet.UU.NET
  10. Reply-To: @uunet.uu.net:Dan.McCue%newcastle.ac.uk@nsfnet-relay.ac.uk
  11. Lines: 413
  12. Approved: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  13.  
  14. Posting-number: Volume 8, Issue 98
  15. Submitted-by: @uunet.uu.net:Dan.McCue%newcastle.ac.uk@nsfnet-relay.ac.uk
  16. Archive-name: interval.c++
  17.  
  18.     Enclosed is a shar archive for a simple interval
  19. arithmetic package written in C++.  There is no
  20. restriction on use or re-distribution.  Comments and
  21. improvements may be addressed to:
  22.   From U.K.: Dan.McCue@uk.ac.newcastle
  23.   Elsewhere: Dan.McCue@newcastle.ac.uk
  24. ------------------------ cut here ----------------------------
  25. #! /bin/sh
  26. # This is a shell archive, meaning:
  27. # 1. Remove everything above the #! /bin/sh line.
  28. # 2. Save the resulting text in a file.
  29. # 3. Execute the file with /bin/sh (not csh) to create:
  30. #    Makefile
  31. #    README
  32. #    Tinterval.c
  33. #    interval.c
  34. #    interval.h
  35. # This archive created: Thu Aug 24 13:13:29 1989
  36. export PATH; PATH=/bin:/usr/bin:$PATH
  37. echo shar: "extracting 'Makefile'" '(715 characters)'
  38. if test -f 'Makefile'
  39. then
  40.     echo shar: "will not over-write existing file 'Makefile'"
  41. else
  42. cat << \SHAR_EOF > 'Makefile'
  43. HDRS          = interval.h
  44. CC          = CC
  45. LINKER          = CC
  46. MAKEFILE      = Makefile
  47. OBJS          = Tinterval.o \
  48.         interval.o
  49. PRINT          = pr
  50. PROGRAM          = Tinterval
  51. SRCS          = Tinterval.c \
  52.         interval.c
  53.  
  54. all:        $(PROGRAM)
  55.  
  56. $(PROGRAM):     $(OBJS) $(LIBS)
  57.         @echo -n "Linking $(PROGRAM) ... "
  58.         @$(LINKER) $(LDFLAGS) $(OBJS) $(LIBS) -o $(PROGRAM)
  59.         @echo "done"
  60.  
  61. clean:;        @rm -f $(OBJS) core a.out
  62.  
  63. depend:;    @mkmf -f $(MAKEFILE) PROGRAM=$(PROGRAM) DEST=$(DEST)
  64.  
  65. index:;        @ctags -wx $(HDRS) $(SRCS)
  66.  
  67. print:;        @$(PRINT) $(HDRS) $(SRCS)
  68.  
  69. tags:           $(HDRS) $(SRCS); @ctags $(HDRS) $(SRCS)
  70.  
  71. ###
  72. Tinterval.o: /usr/include/CC/stdio.h /usr/include/CC/string.h interval.h
  73. interval.o: /usr/include/CC/stdio.h interval.h
  74. SHAR_EOF
  75. if test 715 -ne "`wc -c < 'Makefile'`"
  76. then
  77.     echo shar: "error transmitting 'Makefile'" '(should have been 715 characters)'
  78. fi
  79. fi
  80. echo shar: "extracting 'README'" '(498 characters)'
  81. if test -f 'README'
  82. then
  83.     echo shar: "will not over-write existing file 'README'"
  84. else
  85. cat << \SHAR_EOF > 'README'
  86.     This is a simple implementation of an interval
  87. arithmetic package I wrote while trying to learn C++.
  88. The accuracy is dubious (since proper attention has not
  89. been paid to rounding) and there is no attempt to recover
  90. from errors.  The result is not too robust, but fun, useful
  91. and (in my case), instructive.
  92.     Be sure to modify the makefile to use your c++ compiler.
  93. You may also need to change the pathnames of include files
  94. (e.g., <CC/stdio.h>) in interval.c and Tinterval.c for your
  95. system.
  96. SHAR_EOF
  97. if test 498 -ne "`wc -c < 'README'`"
  98. then
  99.     echo shar: "error transmitting 'README'" '(should have been 498 characters)'
  100. fi
  101. fi
  102. echo shar: "extracting 'Tinterval.c'" '(1020 characters)'
  103. if test -f 'Tinterval.c'
  104. then
  105.     echo shar: "will not over-write existing file 'Tinterval.c'"
  106. else
  107. cat << \SHAR_EOF > 'Tinterval.c'
  108. // A simple test program for the interval class
  109.  
  110. #include <CC/stdio.h>
  111. #include <CC/string.h>
  112. #include "interval.h"
  113.  
  114. void Usage(char *name)
  115. {
  116.     fprintf(stderr, "usage: %s real real +|-|*|/ real real\n", name);
  117. }
  118.  
  119. main(int argc, char **argv)
  120. {
  121.     float a,b,c,d;
  122.  
  123.     if (argc != 6) {
  124.     Usage(argv[0]);
  125.     exit(1);
  126.     }
  127.  
  128.     a = atof(argv[1]);
  129.     b = atof(argv[2]);
  130.     c = atof(argv[4]);
  131.     d = atof(argv[5]);
  132.     {
  133.     interval I1(a,b);
  134.     interval I2(c,d);
  135.     interval Result, Result2;
  136.  
  137.     switch (argv[3][0]) {
  138.     case '+': Result = I1 + I2; break;
  139.     case '-': Result = I1 - I2; break;
  140.     case '*': Result = I1 * I2; break;
  141.     case '/': Result = I1 / I2; break;
  142.     default: 
  143.         Usage(argv[0]);
  144.         exit(2);
  145.     }
  146.  
  147.     printf("[%f, %f] %s [%f, %f] == [%f, %f]\n", 
  148.         I1.lo(), I1.hi(), argv[3], I2.lo(), I2.hi(), 
  149.         Result.lo(), Result.hi());
  150.  
  151.     // Now check out += and integer-to-interval conversion
  152.     Result2 = Result;
  153.     Result2 += 4;
  154.     printf("[%f, %f] + 4 == [%f, %f]\n", 
  155.         Result.lo(), Result.hi(), Result2.lo(), Result2.hi());
  156.     }
  157. }
  158. SHAR_EOF
  159. if test 1020 -ne "`wc -c < 'Tinterval.c'`"
  160. then
  161.     echo shar: "error transmitting 'Tinterval.c'" '(should have been 1020 characters)'
  162. fi
  163. fi
  164. echo shar: "extracting 'interval.c'" '(4825 characters)'
  165. if test -f 'interval.c'
  166. then
  167.     echo shar: "will not over-write existing file 'interval.c'"
  168. else
  169. cat << \SHAR_EOF > 'interval.c'
  170. // This software is in the public domain.
  171. // Permission is granted to use this software for any
  172. // purpose commercial or non-commercial.
  173. // The implementer/distributer assume no liability for
  174. // any damages, incidental, consequential or otherwise,
  175. // arising from the use of this software.
  176.  
  177. // Implementation of class interval
  178. //
  179. // Most operators are trivial and could be inlined.
  180. //
  181. // Multiply is more interesting.  The complicated testing
  182. // may be more time consuming than the floating-point
  183. // operations it attempts to avoid.
  184. // How could this testing be streamlined?
  185. //
  186. // Note that divide blurts out error/warning messages using fprintf :-(
  187. // Strictly speaking, division is undefined if the divisor interval
  188. // contains zero.  This code is sloppy in allowing it, but it warns.
  189. //
  190. // No attempt is made to control rounding so accuracy is dubious.
  191. //
  192. // No attempt is made to detect/recover from overflow or underflow.
  193. // This is a serious deficiency.  Is there a 'portable' around it?
  194.  
  195. #include <CC/stdio.h>
  196. #include "interval.h"
  197.  
  198. interval interval::operator+(interval I)
  199. {
  200.     interval temp(lo_bound, hi_bound);
  201.     return temp += I;
  202. }
  203.  
  204. interval interval::operator-(interval I)
  205. {
  206.     interval temp(lo_bound, hi_bound);
  207.     return temp -= I;
  208. }
  209.  
  210. interval interval::operator*(interval I)
  211. {
  212.     interval temp(lo_bound, hi_bound);
  213.     return temp *= I;
  214. }
  215.  
  216. interval interval::operator/(interval I)
  217. {
  218.     interval temp(lo_bound, hi_bound);
  219.     return temp /= I;
  220. }
  221.  
  222. interval &interval::operator+=(interval I)
  223. {
  224.     lo_bound += I.lo_bound;
  225.     hi_bound += I.hi_bound;
  226.     return *this;
  227. }
  228.  
  229. interval &interval::operator-=(interval I)
  230. {
  231.     lo_bound -= I.lo_bound;
  232.     hi_bound -= I.hi_bound;
  233.     return *this;
  234. }
  235.  
  236. // Multiply is a bit complicated.
  237. // A naive version simply takes the minimum of all combinations
  238. // as the new lo_bound and the maximum as the new hi_bound.
  239. // But this involves four double precision floating multiples.
  240. // A quick examination of the signs of the intervals reveals
  241. // that there are nine possible cases of which only one
  242. // requires all four multiplications.  The other cases only
  243. // require two multiplications.  We assume that it is faster to
  244. // compute and dispatch to the right case than to perform
  245. // the extra multiplies.  If not, use the naive scheme.
  246. //
  247. // The cases can be described as follows:
  248. // Each interval is either non-negative (lo bound >= 0),
  249. // non-positive (hi bound <= 0), or crosses(includes) zero,
  250. // including some numbers on either side.
  251. // If we label the intervals, A and B, corresponding to
  252. // "this" and "I", we get the following matrix:
  253. //
  254. //        Bcross        Bnonneg        Bnonpos
  255. //
  256. //    Across
  257. //
  258. //    Anonneg
  259. //
  260. //    Anonpos
  261.  
  262. interval &interval::operator*=(interval I)
  263. {
  264.     enum possibility {
  265.     AcrossBcross, AcrossBnn, AcrossBnp,
  266.     AnnBcross, AnnBnn, AnnBnp,
  267.     AnpBcross, AnpBnn, AnpBnp
  268.     } choice;
  269.  
  270.     if (lo_bound >= 0.0) choice = AnnBcross;
  271.     else if (hi_bound <= 0.0) choice = AnpBcross;
  272.     else choice = AcrossBcross;
  273.  
  274.     if (I.lo_bound >= 0.0) choice += 1;
  275.     else if (I.hi_bound <= 0.0) choice += 2;
  276.  
  277.     switch (choice)
  278.     {
  279.     case AcrossBcross: {
  280.             double  HL = hi_bound*I.lo_bound,
  281.                 HH = hi_bound*I.hi_bound,
  282.                 LL = lo_bound*I.lo_bound,
  283.                 LH = lo_bound*I.hi_bound;
  284.             lo_bound = HL<LH?HL:LH;
  285.             hi_bound = LL>HH?LL:HH;
  286.         }
  287.         break;
  288.  
  289.     case AcrossBnn:
  290.         lo_bound *= I.hi_bound;
  291.         hi_bound *= I.hi_bound;
  292.         break;
  293.  
  294.     case AcrossBnp: {
  295.         double new_hi_bound = lo_bound * I.lo_bound;
  296.  
  297.         lo_bound = hi_bound*I.lo_bound;
  298.         hi_bound = new_hi_bound;
  299.         }
  300.         break;
  301.  
  302.     case AnnBcross:
  303.         lo_bound = hi_bound*I.lo_bound;
  304.         hi_bound *= I.hi_bound;
  305.         break;
  306.  
  307.     case AnnBnn:
  308.         lo_bound *= I.lo_bound;
  309.         hi_bound *= I.hi_bound;
  310.         break;
  311.  
  312.     case AnnBnp: {
  313.         double new_hi_bound = lo_bound * I.hi_bound;
  314.  
  315.         lo_bound = hi_bound*I.lo_bound;
  316.         hi_bound = new_hi_bound;
  317.         }
  318.         break;
  319.  
  320.     case AnpBcross:
  321.         hi_bound = lo_bound*I.lo_bound;
  322.         lo_bound *= I.hi_bound;
  323.         break;
  324.  
  325.     case AnpBnn:
  326.         lo_bound *= I.hi_bound;
  327.         hi_bound *= I.lo_bound;
  328.         break;
  329.  
  330.     case AnpBnp: {
  331.         double new_hi_bound = lo_bound * I.lo_bound;
  332.  
  333.         lo_bound = hi_bound*I.hi_bound;
  334.         hi_bound = new_hi_bound;
  335.         }
  336.         break;
  337.  
  338.     default:
  339.         fprintf(stderr, 
  340.             "Bad interval (low bound > hi bound) [%f, %f] * [%f, %f]\n",
  341.             lo_bound, hi_bound, I.lo_bound, I.hi_bound);
  342.         lo_bound = hi_bound = 0.0;
  343.         break;
  344.     }
  345.  
  346.     return *this;
  347. }
  348.  
  349. interval &interval::operator/=(interval I)
  350. {
  351.     if (I.lo_bound == 0.0 && I.hi_bound == 0)
  352.     fprintf(stderr, "Error: Division by zero attempted - [%f,%f]/[f,%f]\n",
  353.         lo_bound, hi_bound, I.lo_bound, I.hi_bound);
  354.     else {
  355.     if (I.lo_bound < 0 && I.hi_bound > 0)
  356.         fprintf(stderr, "Warning: Division by zero possible - [%f,%f]/[%f,%f]\n",
  357.         lo_bound, hi_bound, I.lo_bound, I.hi_bound);
  358.     interval inv(1.0/I.lo_bound, 1.0/I.hi_bound);
  359.     return inv *= *this;
  360.     }
  361. }
  362. SHAR_EOF
  363. if test 4825 -ne "`wc -c < 'interval.c'`"
  364. then
  365.     echo shar: "error transmitting 'interval.c'" '(should have been 4825 characters)'
  366. fi
  367. fi
  368. echo shar: "extracting 'interval.h'" '(1684 characters)'
  369. if test -f 'interval.h'
  370. then
  371.     echo shar: "will not over-write existing file 'interval.h'"
  372. else
  373. cat << \SHAR_EOF > 'interval.h'
  374. // Class interval
  375. //    A definition of intervals for interval arithmetic.
  376. //
  377. //    The representation of the interval is two doubles
  378. //    but conversion operators are defined for int and float.
  379. //
  380. //    Error handling is weak.  In an IEEE compliant implementation,
  381. //    quiet NANs could be used to indicate bounds that were
  382. //    invalid.  Otherwise, more state (and checking) would
  383. //    need to be added.
  384.  
  385. // To get the definition of MAXDOUBLE
  386. #include <values.h>
  387.  
  388. class interval {
  389.     double lo_bound, hi_bound;
  390. public:
  391.     double lo(void)        { return lo_bound; }
  392.     double hi(void)        { return hi_bound; }
  393.     double width(void)        { return hi_bound - lo_bound; }
  394.     interval()            { lo_bound = -MAXDOUBLE; hi_bound = MAXDOUBLE; }
  395.     interval(double l, double h){ lo_bound = l<h?l:h; hi_bound = l<h?h:l; }
  396.     interval(double f)        { lo_bound = f; hi_bound = f; }
  397.     interval(float l, float h)    { lo_bound = l<h?l:h; hi_bound = l<h?h:l; }
  398.     interval(float f)        { lo_bound = f; hi_bound = f; }
  399.     interval(int l, int h)    { lo_bound = l<h?l:h; hi_bound = l<h?h:l; }
  400.     interval(int i)        { lo_bound = i; hi_bound = i; }
  401.     interval operator+(interval I);
  402.     interval operator-(interval I);
  403.     interval operator*(interval I);
  404.     interval operator/(interval I);
  405.     interval &operator+=(interval I);
  406.     interval &operator-=(interval I);
  407.     interval &operator*=(interval I);
  408.     interval &operator/=(interval I);
  409.     int contains(interval I)    { return lo_bound <= I.lo_bound && 
  410.                      hi_bound >= I.hi_bound; }
  411.     int overlaps(interval I)    { return lo_bound <= I.hi_bound && 
  412.                      hi_bound >= I.lo_bound; }
  413.     int equal(interval I)    { return lo_bound == I.lo_bound && 
  414.                      hi_bound == I.hi_bound; }
  415. };
  416. SHAR_EOF
  417. if test 1684 -ne "`wc -c < 'interval.h'`"
  418. then
  419.     echo shar: "error transmitting 'interval.h'" '(should have been 1684 characters)'
  420. fi
  421. fi
  422. exit 0
  423. #    End of shell archive
  424. #
  425. #
  426. #
  427.  
  428.  
  429.