home *** CD-ROM | disk | FTP | other *** search
/ Computer Shopper 275 / DPCS0111DVD.ISO / Toolkit / Audio-Visual / VirtualDub / Source / VirtualDub-1.9.10-src.7z / src / system / source / Fraction.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  2009-09-14  |  7.5 KB  |  328 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. VDFraction& VDFraction::operator*=(VDFraction b) {
  176.     return *this = reduce((uint64)hi * b.hi, (uint64)lo * b.lo);
  177. }
  178.  
  179. VDFraction& VDFraction::operator/=(VDFraction b) {
  180.     return *this = reduce((uint64)hi * b.lo, (uint64)lo * b.hi);
  181. }
  182.  
  183. VDFraction& VDFraction::operator*=(unsigned long b) {
  184.     return *this = reduce((uint64)hi * b, lo);
  185. }
  186.  
  187. VDFraction& VDFraction::operator/=(unsigned long b) {
  188.     return *this = reduce(hi, (uint64)lo * b);
  189. }
  190.  
  191. ///////////////////////////////////////////////////////////////////////////
  192.  
  193. sint64 VDFraction::scale64t(sint64 v) const {
  194.     uint32 r;
  195.     return v<0 ? -VDFractionScale64(-v, hi, lo, r) : VDFractionScale64(v, hi, lo, r);
  196. }
  197.  
  198. sint64 VDFraction::scale64u(sint64 v) const {
  199.     uint32 r;
  200.     if (v<0) {
  201.         v = -VDFractionScale64(-v, hi, lo, r);
  202.         return v;
  203.     } else {
  204.         v = +VDFractionScale64(+v, hi, lo, r);
  205.         return v + (r > 0);
  206.     }
  207. }
  208.  
  209. sint64 VDFraction::scale64r(sint64 v) const {
  210.     uint32 r;
  211.     if (v<0) {
  212.         v = -VDFractionScale64(-v, hi, lo, r);
  213.         return v - (r >= (lo>>1) + (lo&1));
  214.     } else {
  215.         v = +VDFractionScale64(+v, hi, lo, r);
  216.         return v + (r >= (lo>>1) + (lo&1));
  217.     }
  218. }
  219.  
  220. sint64 VDFraction::scale64it(sint64 v) const {
  221.     uint32 r;
  222.     return v<0 ? -VDFractionScale64(-v, lo, hi, r) : +VDFractionScale64(+v, lo, hi, r);
  223. }
  224.  
  225. sint64 VDFraction::scale64ir(sint64 v) const {
  226.     uint32 r;
  227.     if (v<0) {
  228.         v = -VDFractionScale64(-v, lo, hi, r);
  229.         return v - (r >= (hi>>1) + (hi&1));
  230.     } else {
  231.         v = +VDFractionScale64(+v, lo, hi, r);
  232.         return v + (r >= (hi>>1) + (hi&1));
  233.     }
  234. }
  235.  
  236. sint64 VDFraction::scale64iu(sint64 v) const {
  237.     uint32 r;
  238.     if (v<0) {
  239.         v = -VDFractionScale64(-v, lo, hi, r);
  240.         return v;
  241.     } else {
  242.         v = +VDFractionScale64(+v, lo, hi, r);
  243.         return v + (r > 0);
  244.     }
  245. }
  246.  
  247. ///////////////////////////////////////////////////////////////////////////
  248.  
  249. double VDFraction::asDouble() const {
  250.     return (double)hi / (double)lo;
  251. }
  252.  
  253. double VDFraction::AsInverseDouble() const {
  254.     return (double)lo / (double)hi;
  255. }
  256.  
  257. unsigned long VDFraction::roundup32ul() const {
  258.     return (hi + (lo-1)) / lo;
  259. }
  260.  
  261. ///////////////////////////////////////////////////////////////////////////
  262.  
  263. bool VDFraction::Parse(const char *s) {
  264.     char c;
  265.  
  266.     // skip whitespace
  267.     while((c = *s) && (c == ' ' || c == '\t'))
  268.         ++s;
  269.  
  270.     // accumulate integer digits
  271.     uint64 x = 0;
  272.     uint64 y = 1;
  273.  
  274.     while(c = *s) {
  275.         uint32 offset = (uint32)c - '0';
  276.  
  277.         if (offset >= 10)
  278.             break;
  279.  
  280.         x = (x * 10) + offset;
  281.  
  282.         // check for overflow
  283.         if (x >> 32)
  284.             return false;
  285.  
  286.         ++s;
  287.     }
  288.  
  289.     if (c == '.') {
  290.         ++s;
  291.  
  292.         while(c = *s) {
  293.             uint32 offset = (uint32)c - '0';
  294.  
  295.             if (offset >= 10)
  296.                 break;
  297.  
  298.             if (x >= 100000000000000000 ||
  299.                 y >= 100000000000000000) {
  300.                 if (offset >= 5)
  301.                     ++x;
  302.                 while((c = *s) && (unsigned)(c - '0') < 10)
  303.                     ++s;
  304.                 break;
  305.             }
  306.  
  307.             x = (x * 10) + offset;
  308.             y *= 10;
  309.             ++s;
  310.         }
  311.     }
  312.  
  313.     while(c == ' ' || c == '\t')
  314.         c = *++s;
  315.  
  316.     // check for trailing garbage
  317.     if (c)
  318.         return false;
  319.  
  320.     // check for overflow
  321.     if (!(y >> 32) && ((uint64)(uint32)y << 32) <= x)
  322.         return false;
  323.  
  324.     // reduce fraction and return success
  325.     *this = reduce(x, y);
  326.     return true;
  327. }
  328.