Click here to Skip to main content
15,867,921 members
Articles / Programming Languages / C#
Tip/Trick

Getting as Close to Rational as You Can

Rate me:
Please Sign up or sign in to vote.
5.00/5 (1 vote)
24 Jul 2011CPOL1 min read 10.8K   1   5
A C# method for finding ratio Q/R approximating real numbers.

When scaling to integer coordinates (especially pixel coordinates), using rounding to binary fractions, or rounding to powers of some approximate differential length, you can use GetRational() to return an array of candidate rational values Q and R such that Q/R offers the scaling you require. This function returns if Ndigits decimal digits to the right of the decimal point in the check-value calculated at each iteration are all zeros or all nines, which is normally the case when the ratio approaches close to integer Q and R. If that proves unsatisfactory, you can increase Ndigits.


The function arguments are:



  • X: Real value to be given a rational approximation.
  • ArrayItems: Length of array to be returned.
  • Ndigits: Number of digits expected.

The function return is: a D3tuple array, each element of which contains the fields Numer, Denom, and Ratio.


C#
public D3tuple[] GetRational(double X, int ArrayItems, int Ndigits)
{
  int d = ArrayItems+1;
  D3tuple[] res = new D3tuple[d]; res[0] = new D3tuple();
  int I;
  double A,E,F,G,W;
  double localEta = Math.Pow(10.0,-Math.Abs(Ndigits));
  I = 0; F = 0; G = 1;
  E = Math.Abs(X); A = E;
  if (X != E) { G = -1; }
  while ((I < d) & (F <= 23))
  {
    W = A - Math.Truncate(A);
    if ((W <= localEta) | ((1.0-W) < localEta))
    {
      return res;
    }
    else
    {
      A = 1 / W;
      if (A > Eta)
      {
        F = F + Math.Log(A);
        if (F > 23) { return res; }
        else if (I < d)
        {
          I++; res[0].Ratio = (double)I;
          res[I] = new D3tuple();
          res[I].Numer = G*Math.Exp(F + Math.Log(E)); res[I].Denom = Math.Exp(F);
          if (F > 0) res[I].Ratio = res[I].Numer / res[I].Denom;
        }
      }
      else { return res; }
    }
  }
  return res;
}

Only two things are missing above. The first is the constant Eta, which you may adjust as you see fit:


C#
const double Eta = 1.0E-15;

The second is the definition of the class which returns the elements of the function result array:


C#
public class D3tuple
{
  double xx; double yy; double zz;
  public D3tuple() { xx = 0; yy = 0; zz = 0; }
  public D3tuple(double x, double y, double z)
  { xx = x; yy = y; zz = z; }
  public void Set(double x, double y, double z)
  { xx = x; yy = y; zz = z; }
  public void Get(ref double x, ref double y, ref double z) 
  { x = xx; y = yy; z = zz; }
  //  X = a[0] = a4[1] = Numer
  public double X 
  { get { return xx; } set { xx = (double)value; } }
  // Y = a[1] = a4[2] = Denom
  public double Y 
  { get { return yy; } set { yy = (double)value; } }
  // Z = a[2] = a4[3] = Ratio
  public double Z 
  { get { return zz; } set { zz = (double)value; } }
  // Numer = X = a[0] = a4[1]
  public double Numer 
  { get { return xx; } set { xx = (double)value; } }
  // Denom = Y = a[1] = a4[2]
  public double Denom 
  { get { return yy; } set { yy = (double)value; } }
  // Ratio = Z = a[2] = a4[3]
  public double Ratio 
  { get { return zz; } set { zz = (double)value; } }
  // This accessor returns a component of the D3tuple.  
  // If D is a D3tuple, D[1] is the X-component, D[2] is the Y-component,
  // and D[3] is the Z-component.
  //  
  // Arg: i = Index of component (1..3)
  // Ret: X, Y, or Z when i = 1, 2, or 3
  public double this[int i]
  {
    get
    {
      double res = new double();
      if ((i >= 1) & (i <= 3)) res = a4[i];
      return res;
    }
  }
  // a = [X,Y,Z] = [Numer,Denom,Ratio] = [a4[1],a4[2],a4[3]]
  //
  // This accessor returns the three components as an array 
  // of three doubles, indexed 0..2, whereas a4 returns an 
  // array of length four indexed 1..3 with the zeroth element unused.
  public double[] a
  {
    get { double[] v = new double[3]; 
          v[0] = xx; v[1] = yy; v[2] = zz; return v; }
    set
    {
      double[] v = value;
      if (value.Length == 4) 
      { xx = v[1]; yy = v[2]; zz = v[3]; }
      else if (value.Length == 3) 
      { xx = v[0]; yy = v[1]; zz = v[2]; }
    }
  }
  // a4 = [0,X,Y,Z] = [0,Numer,Denom,Ratio] = [0,a[0],a[1],a[2]]
  //
  // The a4 accessor is provided to allow an array of 
  // length 4 with elements 1..3 containing the X,Y,Z values
  // to be used to get and set 3tuples.
  public double[] a4
  {
    get { double[] v = new double[4]; v[1] = xx; v[2] = yy; v[3] = zz; return v; }
    set
    {
      double[] v = value;
      if (value.Length == 4) { xx = v[1]; yy = v[2]; zz = v[3]; }
    }
  }
  public D3tuple Square()
  {
    D3tuple res = new D3tuple();
    res.X = X*X; res.Y = Y*Y; res.Z = Z*Z;
    return res;
  }
  public double Sum()
  {
    double res = new double();
    res = X+Y+Z;
    return res;
  }
  // ZeroRef() sets the components X,Y,Z of a D3tuple to zero.
  public void Zero() { xx = 0; yy = 0; zz = 0; }
  // IsZero() returns true if the current instance
  // is the zero vector.
  public bool IsZero() { return ((X*X+Y*Y+Z*Z) == 0); }
  // Eq() returns true if the components of a D3tuple are
  // equal to the components of the D3tuple argument.
  public bool Eq(D3tuple d)
  {
    return ((X==d.X)&&(Y==d.Y)&&(Z==d.Z));
  }
  // == Equality operator:  true if a.X=b.X, a.Y=b.Y, and a.Z=b.Z
  public static bool operator==(D3tuple d1, D3tuple d2)
  {
    return d1.Eq(d2);
  }
  // != Non-equality operator: true if a==b is false.
  public static bool operator!=(D3tuple d1, D3tuple d2)
  {
    return !d1.Eq(d2);
  }
  // This override provides object reference comparisons which
  // are normally provided by the == operator.  a.Equals(b) returns
  // true if the object reference a equals the object reference b.
  public override bool Equals(object obj)
  {
    return base.Equals(obj);
  }
  // This override must be present when the ==/!= operators are
  // subject to overrides in the class.
  public override int GetHashCode()
  {
    return base.GetHashCode();
  }
  // ~ Norm (Vector Length) Operator: ~(A) = |A|
  public static double operator~(D3tuple d1)
  {
    double d; d = Math.Sqrt(d1.X*d1.X+d1.Y*d1.Y+d1.Z*d1.Z);
    return d;
  }
  // ! Unit Vector Operator: !A = A/|A|
  public static D3tuple operator!(D3tuple d1)
  {
    D3tuple d = new D3tuple();
    double L = ~d1; 
    if (L != 0) { d.X = d1.X / L; d.Y = d1.Y / L; d.Z = d1.Z / L; }
    return d;
  }
  // - Pole Reversal Operator: A' = -A
  public static D3tuple operator-(D3tuple d1)
  {
    return new D3tuple(-d1.X, -d1.Y, -d1.Z);
  }
  // + Addition Operator C = A + B
  public static D3tuple operator+(D3tuple d1, D3tuple d2)
  {
    return new D3tuple(d1.X+d2.X, d1.Y+d2.Y, d1.Z+d2.Z);
  }
  // - Subtraction Operator C = B - A
  public static D3tuple operator-(D3tuple d1, D3tuple d2)
  {
    return new D3tuple(d1.X-d2.X, d1.Y-d2.Y, d1.Z-d2.Z);
  }
  // * Scalar (Dot) Product Operator: A*B = a1*b1+a2*b2+a3*b3
  public static double operator*(D3tuple d1, D3tuple d2)
  {
    double d = new double();
    d = d1.xx*d2.xx+d1.yy*d2.yy+d1.zz*d2.zz;
    return d;
  }
  // * Scalar Component Multiplication Operator: h*A = [h*a1,h*a2,h*a3]
  // h: Scalar multiplier h
  // d1: Vector A to be multiplied by h
  // Returns C = [h*d1.X,h*d1.Y,h*d1.Z]
  public static D3tuple operator*(double h, D3tuple d1)
  {
    D3tuple dr = new D3tuple(h*d1.X,h*d1.Y,h*d1.Z);
    return dr;
  }
  // / Length Division Operator: A/B = A / |B| or [0,0,0]
  // d1: Vector A whose components are to be divided
  // d2: Vector B whose length provides the divisor L
  // Returns: C = [d1.X/L,d1.Y/L,d1.Z/L]
  public static D3tuple operator/(D3tuple d1, D3tuple d2)
  {
    D3tuple d = new D3tuple(); 
    double L = Math.Sqrt(d2.X*d2.X+d2.Y*d2.Y+d2.Z*d2.Z);
    if (L != 0) { d.X = d1.X / L; d.Y = d1.Y / L; d.Z = d1.Z / L; }
    return d;
  }
  // / Scalar Divide Components Operator: A/h = [a1/h,a2/h,a3/h]
  // d1: Vector A whose components are to be divided
  // v: Scalar divisor v
  // Returns C = [a1/v,a2/v,a3/v]
  public static D3tuple operator/(D3tuple d1, double v)
  {
    D3tuple d = new D3tuple(d1.X,d1.Y,d1.Z);
    if (v != 0) { d.X = d1.X / v; d.Y = d1.Y / v; d.Z = d1.Z / v; }
    return d;
  }
  // ^ Vector Product Operator: 
  //     A^B = [a2b3-a3b2,a3b1-a1b3,a1b2-a2b1]
  // Operator ^ is the vector (cross) product operator
  // d1: Vector A
  // d2: Vector B
  // Returns: CrossProd(A,B)
  public static D3tuple operator^(D3tuple d1, D3tuple d2)
  {
    return new D3tuple(
      d1.yy*d2.zz-d2.yy*d1.zz, 
      d1.zz*d2.xx-d2.zz*d1.xx, 
      d1.xx*d2.yy-d2.xx*d1.yy);
  }
  // & Carrot Product Operator:
  //    A & B = [a1*b1,a2*b2,a3*b3]
  // Operator & is the Carrot Product operator,
  // which multiplies the components of vector A with
  // those of vector B.
  // d1: Vector A = [a1,a2,a3]
  // d2: Vector B = [b1,b2,b3]
  // Returns: [a1*b1,a2*b2,a3*b3]
  public static D3tuple operator&(D3tuple d1, D3tuple d2)
  {
    return new D3tuple(d1.X*d2.X, d1.Y*d2.Y, d1.Z*d2.Z);
  }
  // >> Right-shift Indexes Operator
  // The right shift operator shifts component indexing
  // to the right in the 3-component ring.  A shift count
  // of +1 will cause the Y component to be the first 
  // component of the result vector, with the Z and X 
  // components following it in cyclic order.  A shift count
  // of 2 would place the Z component in the first component
  // position in the result vector, with the X and Y components
  // following it.
  // Shift counts of 1 and -2 are equivalent, as are 2 and -1.
  // d1: The vector whose components are to be rotated.
  // k: The rotation count
  // Returns: The rotated vector
  public static D3tuple operator>>(D3tuple d1, int k)
  {
    int c = k % 3; D3tuple d = new D3tuple(d1.X,d1.Y,d1.Z);
    if ((c == 1)|(c == -2)) d.Set(d1.Y, d1.Z, d1.X);
    else if ((c == 2)|(c == -1)) d.Set(d1.Z, d1.X, d1.Y);
    return d;
  }
  // << Left-shift Indexes Operator
  // The left shift operator shifts component indexing to the
  // left in the 3-component ring.  A shift count of +1 would
  // cause the Z component to be the first component of the 
  // result vector, with the X and Y components following it
  // in cyclic order.
  public static D3tuple operator<<(D3tuple d1, int k)
  {
    int c = k % 3; D3tuple d = new D3tuple(d1.X,d1.Y,d1.Z);
    if ((c == -1)|(c == 2)) d.Set(d1.Y, d1.Z, d1.X);
    else if ((c == -2)|(c == 1)) d.Set(d1.Z, d1.X, d1.Y);
    return d;
  }
  bool WriteSquare = false;
  string sf(double d)
  {
    string res = ""; if (WriteSquare) d = d * d;
    res = string.Format("{0,18:F15}", d);
    int p = res.IndexOf(' '),pm=res.IndexOf('-'),L=res.Length;
    if (p >= 0)
    {
      while ((p < L) & (res[p] == ' ')) p++;
      if ((p > 0) & (p == pm)) { res.Remove(p-1, 1); }
      else if (pm < 0) { }
    }
    else if (pm < 0) { res = " " + res; }
    return res;
  }
  // The ToString() override produces the form 
  //   [xxxxxx,yyyyyy,zzzzzz] as a string containing the 
  // numeric values of the components of the D3tuple.
  public override string ToString()
  {
    return "["+sf(xx)+","+sf(yy)+","+sf(zz)+"]";
  }
  public string ToSqString()
  {
    WriteSquare = true;  string s=this.ToString(); WriteSquare = false;
    return s;
  }
}

There's no sense in including the D3tuple class without showing how its major features can be used. So here's a demonstration:


C#
// This demonstration shows the use of the operator overrides
// defined in the D3tuple class involving the operators ^ * + -
// ! / & and == as well as the bit-shift operators and 
// explains the meaning of each.
public void AssignmentsExample()
{
  // Values can be assigned to the X,Y,Z fields 
  // when the D3tuple is instantiated.
  D3tuple Diag = new D3tuple(1,2,3);
  double[] v = new double[3];
  // The design of the D3tuple class allows you to set and 
  // get values from an array, either by passing a reference
  // to an array (here, v) or to its individual elements; i.e., v[1].
  v = Diag.a;
  v[0] = 0.1; v[1] = 0.2; v[2] = 0.3;  
  // Note above that v is indexed from 0, as usual
  Diag.a = v;
  v[1] = Diag.a[1]; 
  // Note above that the D3tuple is indexed from 0 also
  // using the .a accessor
  Diag.a[2] = v[2];
  // Now some vectors are defined for use in what follows.
  D3tuple 
    x1=new D3tuple(1,-2,3), 
    x2=new D3tuple(2,3,-4), 
    x3=new D3tuple(1,0,0);
  // Calculate the vector product of x1 and x2. 
  // Note that it's screwy to have to use indexes 0..2
  // rather than the mathematical indexes 1..3 which sane 
  // and sober mathematicians have used for centuries.
  Diag.Set(
    x1.a[1]*x2.a[2]-x1.a[2]*x2.a[1], 
    x1.a[2]*x2.a[0]-x1.a[0]*x2.a[2], 
    x1.a[0]*x2.a[1]-x1.a[1]*x2.a[0]);
  // This uses an operator overload for the ^ operator,
  // which calculates the vector product, and is non-
  // commutative; u.e., x1^x2 = -(x2^x1).  Much simpler 
  // than the previous statement.
  Diag = x1^x2;
  // This uses an operator overload for the * operator, 
  // which calculates the scalar product, and is commutative.
  v[0] = x1*x2;
  // This uses two operator overloads to calculate the 
  // scalar triple product of three vectors.
  // The parentheses are required.
  v[0] = Diag*(x1^x2);
  // This uses an operator overload for the & operator, 
  // which calculates the "carrot" product, and is commutative.
  Diag = x1 & x2;
  // The * operator overload can be used to pre-multiply a 
  // vector by a scalar, and this usage is not commutative.
  Diag = 5.0 * x1;
  // The / operator overload can be used to divide the components
  // of a vector by a scalar, and is not commutative.
  Diag = x1 / 2;
  // The / operator overload can be used to divide the components
  // of a vector by the length of a second vector following the /
  // operator, and is not commutative.
  Diag = x1 / x2;
  // The ~ operator overload returns the length of a vector, 
  // which is a scalar.
  v[0] = ~x1;
  // The ! operator overload returns a unit vector by dividing
  // the operand by its length. Note that !x1 equals x1 / x1; 
  // If v is any vector, u / !v is equal to u.
  // What is done here is to multiply the unit vector !x1 by
  // the length of x1, which is x1.
  Diag = !x1;
  if (!x1 == (x1 / x1)) { x1 = (~x1) * !x1; }
  // Here, the reciprocal vectors are calculated from x1,x2,x3
  // by dividing the cross product by the scalar triple product.
  // Note that the operand following the / is equal in all 3
  // expressions.
  Diag = (x1^x2) / (x3 * (x1^x2));
  Diag = (x2^x3) / (x1 * (x2^x3));
  Diag = (x3^x1) / (x2 * (x3^x1));
  // The following shows three methods of calculating the 
  // vector cross product. The second uses the >> operator overload,
  // which shifts the components such that a component with a higher
  // index is shifted into each position, with the indexes being 
  // evaluated modulo 3 before the assignment is made, and the third uses
  // both component-shift operators.
  Diag = x1^x2;
  Diag = ((x1>>1)&(x2>>2)) - ((x1>>2)&(x2>>1));
  Diag = ((x1>>1)&(x2>>2)) - ((x1<<1)&(x2<<2));
}

I hope you have fun with it!

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
United States United States
Writer, designer, inventor, musician, observer, and critic with many years of experience in each of these areas. I can do without the approval of "experts" because I believe candid statements and penetrating analysis bring greater rewards than a "pat on the back". And if I have something to say when you're not listening, I tell someone else about it.

Comments and Discussions

 
GeneralReason for my vote of 5 A quite useful tool when porting to ... Pin
YvesDaoust31-Jul-11 4:41
YvesDaoust31-Jul-11 4:41 
GeneralWhy do you need this vector class ? Pin
YvesDaoust25-Jul-11 20:56
YvesDaoust25-Jul-11 20:56 
GeneralRe: When working with un-normalized oblique coordinate systems, ... Pin
GAMerritt31-Jul-11 3:17
GAMerritt31-Jul-11 3:17 
GeneralWould have been instructive to know what is the method. Not ... Pin
YvesDaoust25-Jul-11 20:56
YvesDaoust25-Jul-11 20:56 
GeneralRe: That's a sharp observation; but continued fractions are even... Pin
GAMerritt31-Jul-11 3:31
GAMerritt31-Jul-11 3:31 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.