home *** CD-ROM | disk | FTP | other *** search
/ Computer Shopper 242 / Issue 242 - April 2008 - DPCS0408DVD.ISO / Software Money Savers / VirtualDub / Source / VirtualDub-1.7.7-src.7z / src / system / source / Fraction.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  2007-10-23  |  7.2 KB  |  316 lines
  1. //    VirtualDub - Video processing and capture application
  2. //    System library component
  3. //    Copyright (C) 1998-2006 Avery Lee, All Rights Reserved.
  4. //
  5. //    Beginning with 1.6.0, the VirtualDub system library is licensed
  6. //    differently than the remainder of VirtualDub.  This particular file is
  7. //    thus licensed as follows (the "zlib" license):
  8. //
  9. //    This software is provided 'as-is', without any express or implied
  10. //    warranty.  In no event will the authors be held liable for any
  11. //    damages arising from the use of this software.
  12. //
  13. //    Permission is granted to anyone to use this software for any purpose,
  14. //    including commercial applications, and to alter it and redistribute it
  15. //    freely, subject to the following restrictions:
  16. //
  17. //    1.    The origin of this software must not be misrepresented; you must
  18. //        not claim that you wrote the original software. If you use this
  19. //        software in a product, an acknowledgment in the product
  20. //        documentation would be appreciated but is not required.
  21. //    2.    Altered source versions must be plainly marked as such, and must
  22. //        not be misrepresented as being the original software.
  23. //    3.    This notice may not be removed or altered from any source
  24. //        distribution.
  25.  
  26. #include "stdafx.h"
  27. #include <math.h>
  28.  
  29. #include <vd2/system/fraction.h>
  30. #include <vd2/system/vdtypes.h>
  31. #include <vd2/system/math.h>
  32.  
  33. VDFraction::VDFraction(double d) {
  34.     int xp;
  35.     double mant = frexp(d, &xp);
  36.  
  37.     if (xp >= 33) {
  38.         hi = 0xFFFFFFFF;
  39.         lo = 1;
  40.     } else if (xp < -31) {
  41.         hi = 0;
  42.         lo = 1;
  43.     } else if (xp >= 0) {
  44.         *this = reduce((uint64)(0.5 + ldexp(mant, 62)), 1i64<<(62-xp));
  45.     } else {
  46.         // This is not quite accurate for very tiny numbers.
  47.         VDFraction t(1.0 / d);
  48.         lo = t.hi;
  49.         hi = t.lo;
  50.     }
  51. }
  52.  
  53. VDFraction VDFraction::reduce(uint64 hi, uint64 lo) {
  54.  
  55.     // Check for undefined.
  56.  
  57.     if (!lo)
  58.         return VDFraction(0,0);
  59.  
  60.     // Check for zero.
  61.  
  62.     if (!hi) {
  63.         return VDFraction(0,1);
  64.     }
  65.  
  66.     // Check for infinity.
  67.  
  68.     if (!((uint64)lo>>32) && (uint64)hi > ((uint64)lo<<32)-lo)
  69.         return VDFraction(0xFFFFFFFFUL, 1);
  70.  
  71.     // Algorithm from Wikipedia, Continued Fractions:
  72.     uint64 n0 = 0;
  73.     uint64 d0 = 1;
  74.     uint32 n1 = 1;
  75.     uint32 d1 = 0;
  76.     uint64 fp = 0;
  77.  
  78.     uint32 n_best;
  79.     uint32 d_best;
  80.  
  81.     for(;;) {
  82.         uint64 a = hi/lo;            // next continued fraction term
  83.         uint64 f = hi%lo;            // remainder
  84.  
  85.         uint64 n2 = n0 + n1*a;        // next convergent numerator
  86.         uint64 d2 = d0 + d1*a;        // next convergent denominator
  87.  
  88.         uint32 n_overflow = (uint32)(n2 >> 32);
  89.         uint32 d_overflow = (uint32)(d2 >> 32);
  90.  
  91.         if (n_overflow | d_overflow) {
  92.             uint64 a2 = a;
  93.  
  94.             // reduce last component until numerator and denominator are within range
  95.             if (n_overflow)
  96.                 a2 = (0xFFFFFFFF - n0) / n1;
  97.  
  98.             if (d_overflow) {
  99.                 uint64 a3 = (0xFFFFFFFF - d0) / d1;
  100.                 if (a2 > a3)
  101.                     a2 = a3;
  102.             }
  103.  
  104.             // check if new term is better
  105.             // 1/2a_k admissibility test
  106.             if (a2*2 < a || (a2*2 == a && d0*fp <= f*d1))
  107.                 return VDFraction((uint32)n_best, (uint32)d_best);
  108.  
  109.             return VDFraction((uint32)(n0 + n1*a2), (uint32)(d0 + d1*a2));
  110.         }
  111.  
  112.         n_best = (uint32)n2;
  113.         d_best = (uint32)d2;
  114.  
  115.         // if fraction is exact, we're done.
  116.         if (!f)
  117.             return VDFraction((uint32)n_best, (uint32)d_best);
  118.  
  119.         n0 = n1;
  120.         n1 = (uint32)n2;
  121.         d0 = d1;
  122.         d1 = (uint32)d2;
  123.         fp = f;
  124.  
  125.         hi = lo;
  126.         lo = f;
  127.     }
  128. }
  129.  
  130. // a (cond) b
  131. // a-b (cond) 0
  132. // aH*bL - aL*bh (cond) 0
  133. // aH*bL (cond) aL*bH
  134.  
  135. bool VDFraction::operator==(VDFraction b) const {
  136.     return (uint64)hi * b.lo == (uint64)lo * b.hi;
  137. }
  138.  
  139. bool VDFraction::operator!=(VDFraction b) const {
  140.     return (uint64)hi * b.lo != (uint64)lo * b.hi;
  141. }
  142.  
  143. bool VDFraction::operator< (VDFraction b) const {
  144.     return (uint64)hi * b.lo < (uint64)lo * b.hi;
  145. }
  146.  
  147. bool VDFraction::operator<=(VDFraction b) const {
  148.     return (uint64)hi * b.lo <= (uint64)lo * b.hi;
  149. }
  150.  
  151. bool VDFraction::operator> (VDFraction b) const {
  152.     return (uint64)hi * b.lo > (uint64)lo * b.hi;
  153. }
  154.  
  155. bool VDFraction::operator>=(VDFraction b) const {
  156.     return (uint64)hi * b.lo >= (uint64)lo * b.hi;
  157. }
  158.  
  159. VDFraction VDFraction::operator*(VDFraction b) const {
  160.     return reduce((uint64)hi * b.hi, (uint64)lo * b.lo);
  161. }
  162.  
  163. VDFraction VDFraction::operator/(VDFraction b) const {
  164.     return reduce((uint64)hi * b.lo, (uint64)lo * b.hi);
  165. }
  166.  
  167. VDFraction VDFraction::operator*(unsigned long b) const {
  168.     return reduce((uint64)hi * b, lo);
  169. }
  170.  
  171. VDFraction VDFraction::operator/(unsigned long b) const {
  172.     return reduce(hi, (uint64)lo * b);
  173. }
  174.  
  175. ///////////////////////////////////////////////////////////////////////////
  176.  
  177. sint64 VDFraction::scale64t(sint64 v) const {
  178.     uint32 r;
  179.     return v<0 ? -VDFractionScale64(-v, hi, lo, r) : VDFractionScale64(v, hi, lo, r);
  180. }
  181.  
  182. sint64 VDFraction::scale64u(sint64 v) const {
  183.     uint32 r;
  184.     if (v<0) {
  185.         v = -VDFractionScale64(-v, hi, lo, r);
  186.         return v;
  187.     } else {
  188.         v = +VDFractionScale64(+v, hi, lo, r);
  189.         return v + (r > 0);
  190.     }
  191. }
  192.  
  193. sint64 VDFraction::scale64r(sint64 v) const {
  194.     uint32 r;
  195.     if (v<0) {
  196.         v = -VDFractionScale64(-v, hi, lo, r);
  197.         return v - (r >= (lo>>1) + (lo&1));
  198.     } else {
  199.         v = +VDFractionScale64(+v, hi, lo, r);
  200.         return v + (r >= (lo>>1) + (lo&1));
  201.     }
  202. }
  203.  
  204. sint64 VDFraction::scale64it(sint64 v) const {
  205.     uint32 r;
  206.     return v<0 ? -VDFractionScale64(-v, lo, hi, r) : +VDFractionScale64(+v, lo, hi, r);
  207. }
  208.  
  209. sint64 VDFraction::scale64ir(sint64 v) const {
  210.     uint32 r;
  211.     if (v<0) {
  212.         v = -VDFractionScale64(-v, lo, hi, r);
  213.         return v - (r >= (hi>>1) + (hi&1));
  214.     } else {
  215.         v = +VDFractionScale64(+v, lo, hi, r);
  216.         return v + (r >= (hi>>1) + (hi&1));
  217.     }
  218. }
  219.  
  220. sint64 VDFraction::scale64iu(sint64 v) const {
  221.     uint32 r;
  222.     if (v<0) {
  223.         v = -VDFractionScale64(-v, lo, hi, r);
  224.         return v;
  225.     } else {
  226.         v = +VDFractionScale64(+v, lo, hi, r);
  227.         return v + (r > 0);
  228.     }
  229. }
  230.  
  231. ///////////////////////////////////////////////////////////////////////////
  232.  
  233. VDFraction::operator long() const {
  234.     return (long)((hi + lo/2) / lo);
  235. }
  236.  
  237. VDFraction::operator unsigned long() const {
  238.     return (hi + lo/2) / lo;
  239. }
  240.  
  241. VDFraction::operator double() const {
  242.     return (double)hi / (double)lo;
  243. }
  244.  
  245. unsigned long VDFraction::roundup32ul() const {
  246.     return (hi + (lo-1)) / lo;
  247. }
  248.  
  249. ///////////////////////////////////////////////////////////////////////////
  250.  
  251. bool VDFraction::Parse(const char *s) {
  252.     char c;
  253.  
  254.     // skip whitespace
  255.     while((c = *s) && (c == ' ' || c == '\t'))
  256.         ++s;
  257.  
  258.     // accumulate integer digits
  259.     uint64 x = 0;
  260.     uint64 y = 1;
  261.  
  262.     while(c = *s) {
  263.         uint32 offset = (uint32)c - '0';
  264.  
  265.         if (offset >= 10)
  266.             break;
  267.  
  268.         x = (x * 10) + offset;
  269.  
  270.         // check for overflow
  271.         if (x >> 32)
  272.             return false;
  273.  
  274.         ++s;
  275.     }
  276.  
  277.     if (c == '.') {
  278.         ++s;
  279.  
  280.         while(c = *s) {
  281.             uint32 offset = (uint32)c - '0';
  282.  
  283.             if (offset >= 10)
  284.                 break;
  285.  
  286.             if (x >= 100000000000000000 ||
  287.                 y >= 100000000000000000) {
  288.                 if (offset >= 5)
  289.                     ++x;
  290.                 while((c = *s) && (unsigned)(c - '0') < 10)
  291.                     ++s;
  292.                 break;
  293.             }
  294.  
  295.             x = (x * 10) + offset;
  296.             y *= 10;
  297.             ++s;
  298.         }
  299.     }
  300.  
  301.     while(c == ' ' || c == '\t')
  302.         c = *++s;
  303.  
  304.     // check for trailing garbage
  305.     if (c)
  306.         return false;
  307.  
  308.     // check for overflow
  309.     if (!(y >> 32) && ((uint64)(uint32)y << 32) <= x)
  310.         return false;
  311.  
  312.     // reduce fraction and return success
  313.     *this = reduce(x, y);
  314.     return true;
  315. }
  316.