Computer Science Canada :: View topic - Floating point operation accuracy

/*
* Vector.java
*
* Created on June 26, 2007, 3:01 PM
*
* To change this template, choose Tools | Template Manager
* and open the template in the editor.
*/
package org.aaziz.math.geom;

import java.util.Arrays;
import java.text.NumberFormat;

import java.util.Random;
import java.util.List;

/**
* Represents a geometric vector.
*
* 
* This class is abstract because certain methods (such as multiplication)
* are different for vectors of different dimension
* 
*
* 
* e.g. for 3-dimensionsal vectors there is a 'dot' product and a 'cross'
* product, while 2-dimensional vectors only have the 'dot' product operation.
* 
*
* 
* A vector {@code v} could be represented as {@code v = (a,b)}. In that case,
* the vector has 2 dimensions, such as {@code x} and {@code y} in a
* Cartesian plane.
* 
*
* 
* The component {@code a} would be the component of the first dimension,
* {@code x} for example. Likewise, {@code b} would be the component of the
* second dimension, {@code y} for example.
* 
*
* 
* So, for {@code v}, {@code |v|} represents the magnitude of the vector. That
* is, {@code |v|} is the length of the vector. The square of the magnitude
* is equal to the square of each component, summed.
* 
*
* 
* Subclasses that define a fixed dimension should override the constructor
* {@link #Vector(int)} to forbid the client to set a different number
* of dimensions other than intended.
* 
*
* @author ANTHONYA
*/;
public class Vector implements Dimensional
{
 /**
 * The components of this vector
 */
 private double[] components;

 /**
 * Creates a new, blank {@code Vector} of zero dimensions, thus has
 * no magnitude or direction.
 */
 public Vector()
 {
 this(0);
 }

 /**
 * Create a new {@code Vector} as a copy from another.
 *
 * The new {@code Vector} is given the same amount of dimensions as
 * {@code v}. The new {@code Vector} will be parellel to
 * {@code v} and have the same magnitude as well.
 *
 * @param v {@code Vector} to be copied
 */
 public Vector(Vector v)
 {
 this(v.getDimensions());

 for (int i = 0; i < this.getDimensions(); i++)
 this.setComponent(i, v.getComponent(i));
 }

 /**
 * Instantiate new {@code Vector} with {@code dimesions}
 * dimensions, with no magnitude (a zero Vector).
 *
 * @param dimensions the amount of dimensions this {@code Vector}
 * will have
 */
 public Vector(int dimensions)
 {
 this.components = new double[dimensions];

 for (int i = 0; i < this.getDimensions(); i++)
 this.setComponent(i, 0);

 }

 /**
 * Instantiate new {@code Vector} using the supplied
 * values in {@code components}.
 *
 * 
 * This object will have the same dimensions as there is components, or
 * otherwise: {@code this.getDimensions() == components.length}.
 * 
 *
 * @param components component values that this object will have
 */
 public Vector(double... components)
 {
 this(components.length);

 for (int i = 0; i < this.getDimensions(); i++)
 this.setComponent(i, components[i]);
 }

 /**
 * Returns the amount of dimensions this {@code Vector} has.
 *
 * @return the dimensions this {@code Vector} has.
 */
 public int getDimensions()
 {
 return components.length;
 }

 /**
 * Returns the component of the specified {@code dimension}.
 *
 * @param dimension dimension to get component of (such as x, y, etc)
 *
 * @return component (magnitude) of the specified {@code dimension}.
 * @throws InvalidDimensionException when
 * {@code dimension < 0 || dimension > getDimensions()}.
 */
 public double getComponent(int dimension) throws InvalidDimensionException
 {
 if (dimension < 0 || dimension > getDimensions())
 throw new InvalidDimensionException(dimension);

 return components[dimension];
 }

 /**
 * Gets all the components as an array of integers.
 *
 * 
 * The first index represents the first dimension, etc.
 * 
 *
 * @return list of components, in order
 */
 public double[] getComponents()
 {
 return Arrays.copyOf(components, components.length);
 }

 /**
 * Sets the component of the specified {@code dimension}.
 *
 * @param dimension dimension to set component for (such as x, y, etc)
 * @param value value to set to the component of the specified
 * {@code dimension}
 *
 * @throws InvalidDimensionException when
 * {@code dimension < 0 || dimension > getDimensions()}.
 */
 public void setComponent(int dimension, double value)
 throws InvalidDimensionException
 {
 if (dimension < 0 || dimension > getDimensions())
 throw new InvalidDimensionException(dimension);

 components[dimension] = value;
 }

 /**
 * Calculates and returns the magnitude of this {@code Vector}.
 *
 * 
 * Calculates the magnitute of this Vector, in unit, using an algorithm
 * based on Pythagorean's Theorem. The magnitude is always a positive
 * scalar value (i.e. no direction).
 * 
 *
 * 
 * The formula used is:
 * 
 *
 * <pre>
 * v = a vector
 * |v| = magnitude of v
 * v_1 = component of v in 1st dimension
 * v_2 = component of v in 2nd dimension
 * .
 * .
 * v_n = component of v in nth dimension
 *
 * (^ represents exponent operator)
 *
 * |v| = sqrt(v_1^2 + v_2^2 + .. + v_n^n)
 * </pre>
 *
 * @return magnitute of this {@code Vector} in units
 */
 public double getMagnitude()
 {
 // Need to get the sum of each component, squared
 // eg. x^2 + y^2 + z^2

 int sumOfSquares = 0;
 for (double component : components)
 sumOfSquares += Math.pow(component, 2);

 return Math.sqrt(sumOfSquares);
 }

 /**
 * Calculates the angle between this {@code Vector} and {@code other}.
 *
 * 
 * The angle between two vectors is the angle when they are joined at the
 * tails. The vectors can be of any dimensions, as long as they have the
 * same amount of dimensions. The analytical formula for finding the angle
 * between any two vectors is:
 * 
 *
 * <pre>
 * let v and u be vectors of the same dimension
 *
 * let |v| and |u| represent the magnitudes of v and u, respectively
 *
 * let A represent the angle between v and u
 *
 * v @ u where @ represents dot product
 * cos A = -------
 * |v|*|u|
 * </pre>
 *
 * @param other another {@code Vector}
 *
 * @return the angle, in degrees, between this {@code Vector} and
 * {@code other}.
 *
 * @throws InvalidDimensionException if
 * {@code this.getDimensions() != other.getDimensions()}
 */
 public double angleBetween(Vector other) throws InvalidDimensionException
 {
 if (this.getDimensions() != other.getDimensions())
 {
 throw new InvalidDimensionException(
 "Vectors must be of the same dimension to" +
 "find the angle between them.");
 }

 // Calculate right hand side of equation (in javadoc comments)
 double rhs = (this.dot(other))
 / (this.getMagnitude() * other.getMagnitude());

 // Use inverse of cos to find angle
 // Trigonometric functions of the Math class give angles in radians
 double angle = Math.toDegrees(Math.acos(rhs));

 //
 return angle ;
 }

 /**
 * Returns a {@code Vector} that is opposite to this one.
 *
 * 
 * If {@code v} is a vector, than the opposite is denoted as {@code -v}.
 * This opposite vector has the same magnitude as the original, and is in
 * an opposite direction. Each component in {@code -v} is the negative of
 * its respective component in {@code v}.
 * 
 *
 * 
 * Since this new vector has the same magnitude as the old one, it is true
 * that {@code v.getMagnitude() == v.opposite().getMagnitude()}.
 * 
 *
 * This is the same as calling {@code scale(-1)}
 *
 * @return a {@code Vector} opposite to this one
 *
 * @see #scale(double)
 */
 public Vector opposite()
 {
 return scale(-1);
 }

 public UnitVector unitVector()
 {
 return new UnitVector(this);
 }

 /**
 * Returns a {@code Vector} that is opposite to this one.
 *
 * 
 * If {@code v} is a vector, than the opposite is denoted as {@code -v}.
 * This opposite vector has the same magnitude as the original, and is in
 * an opposite direction. Each component in {@code -v} is the negative of
 * its respective component in {@code v}.
 * 
 *
 * 
 * Since this new vector has the same magnitude as the old one, it is true
 * that {@code v.getMagnitude() == v.opposite().getMagnitude()}.
 * 
 *
 * This is the same as calling {@code scale(-1)}
 *
 * @param factor factor by which to scale this {@code Vector}
 *
 * @return a {@code Vector} that is parallel to this one, whose magnitude is
 * a scale of
 *
 * @see #scale(double)
 */
 public Vector scale(double factor)
 {
 // Create a new set of components and scale each one
 double[] scaledComps = getComponents();
 for (double c : scaledComps)
 c *= factor;

 return new Vector(scaledComps);
 }

 /**
 * Returns the sum of this {@code Vector} added to another.
 *
 * 
 * This method adds the vectors by adding each component to the
 * corresponding component in {@code other}. It is calculated using the
 * following equations:
 * 
 *
 * <pre>
 * v = a vector
 * u = a vector of same dimensions of v
 *
 * v_1 = component of v in 1st dimension
 * v_2 = component of v in 2nd dimension
 * .
 * .
 * v_n = component of v in nth dimension
 *
 * u_1 = component of u in 1st dimension
 * u_2 = component of u in 2nd dimension
 * .
 * .
 * u_n = component of u in nth dimension
 *
 * where v = (v_1, v_2, .. , v_n)
 * and u = (u_1, u_2, .. , u_n)
 *
 * v + u = (v_1, v_2, .. , v_n) + (u_1, u_2, .. , u_n)
 *
 * if t = v + u then t is a vector such that
 *
 * t = (v_1 + u_1, v_2 + u_2, .. , v_n + u_n)
 * </pre>
 *
 * 
 * Furthermore, say there is a vector {@code a} and a vector {@code b}, and
 * {@code c = a + b}. If {@code b} is placed so that it's tail is at the
 * head of {@code a} (or vice versa), then the tail of {@code c} would be
 * at the tail of {@code a} and the head of {@code c} would be at the head
 * of {@code b}. This is called the parallelogram law:
 * 
 *
 * 
 * <img src="doc-files/parallelogram_law.gif">
 * 
 *
 * 
 * Note: this operation is purely functional and does not modify
 * the original {@code Vector} in any way. Only vectors of the same
 * dimensions can be added.
 * 
 *
 * @param other {@code Vector} to add to this one
 *
 * @return a {@code Vector} that is the sum of this {@code Vector} and
 * {@code other}
 *
 * @throws InvalidDimensionException if
 * {@code this.getDimensions() != other.getDimensions()}
 */
 public Vector add(Vector other) throws InvalidDimensionException
 {
 if (this.getDimensions() != other.getDimensions())
 {
 throw new InvalidDimensionException(
 "Only Vectors of same dimension can be added");
 }

 double[] otherComponents = other.getComponents();

 // Add components from this object and the specified one into a
 // new array of components
 double[] newComponents = new double[components.length];
 for (int i = 0; i < newComponents.length; i++)
 {
 newComponents[i] = components[i] + otherComponents[i];
 }

 Vector v = new Vector(newComponents);

 return v;
 }

 /**
 * Returns the dot product of this {@code Vector} and another.
 *
 * 
 * The dot product (or scalar product) of two vectors. It is calculated
 * using the following equation:
 * 
 *
 *
 * <pre>
 * v = a vector
 * u = a vector of same dimensions of v
 *
 * v_1 = component of v in 1st dimension
 * v_2 = component of v in 2nd dimension
 * .
 * .
 * v_n = component of v in nth dimension
 *
 * u_1 = component of u in 1st dimension
 * u_2 = component of u in 2nd dimension
 * .
 * .
 * u_n = component of u in nth dimension
 *
 * where v = (v_1, v_2, .. , v_n)
 * and u = (u_1, u_2, .. , u_n)
 *
 * u • v = (v_1, v_2, .. , v_n) • (u_1, u_2, .. , u_n)
 *
 * u • v = (v_1)(u_1) + (v_2)(u_2) + .. + (v_n)(u_n)
 * </pre>
 *
 * 
 * Note: this operation is purely functional and does not modify
 * the original {@code Vector} in any way. The dot product can only be
 * calculated on vectors of the same dimension.
 * 
 *
 * @param other {@code Vector} to perform dot product operation with
 *
 * @return an scalar ({@code double}) that is dot product of this
 * {@code Vector} and {@code other}
 *
 * @throws InvalidDimensionException if
 * {@code this.getDimensions() != other.getDimensions()}
 */
 public double dot(Vector other) throws InvalidDimensionException
 {
 if (this.getDimensions() != other.getDimensions())
 {
 throw new InvalidDimensionException(
 "Only Vectors of same dimension can be added");
 }

 double[] otherComponents = other.getComponents();

 // Multiply each component in this vector by its respective component
 // in other, and sum the result (using a running total)

 double dotProduct = 0;
 for (int i = 0; i < components.length; i++)
 {
 dotProduct += components[i] * otherComponents[i];
 }

 return dotProduct;
 }

 /**
 * Equates if two vectors are equal (parallel and same magnitude).
 *
 * 
 * Two {@code Vector}s are considered equal if they are parallel
 * (have the same direction) as well as have the same magnitude.
 * 
 *
 * 
 * This can also be equated by comparing that the components in each
 * dimension are the same.
 * 
 *
 * @return
 * - {@code true} if: 
 * {@code obj} is of a non-null object of type {@code Vector}
 * ; and 
 * {@code obj.getDimensions() == this.getDimensions()}
 * ; and 
 * {@code obj} has the same magnitude/direction as this
 * {@code Vector}
 * 
 * 
 * - {@code false} otherwise
 */
 @Override
 public boolean equals(Object obj)
 {
 if (obj == this)
 return true;

 if (obj == null)
 return false;

 if (getClass() != obj.getClass())
 return false;

 // Cast obj to Vector to compare
 final Vector other = (Vector) obj;

 // Compare amount of dimensions
 if (this.getDimensions() != other.getDimensions())
 return false;

 // Compare each component
 for (int i = 0; i < this.getDimensions(); i++)
 {
 if (this.getComponent(i) != other.getComponent(i))
 return false;
 }

 return true;
 }

 /**
 * Returns a hashcode value of this {@code Vector}.
 *
 * @return hashcode of this object.
 */
 @Override
 public int hashCode()
 {
 // Multiple each component by the magnitude and a mix-up number,
 // then add it to the hash value

 int hash = 7;

 double magnitude = getMagnitude();
 for (double i : components)
 {
 hash = 31 * hash + (int)Math.round(i * getMagnitude());
 }

 return hash;
 }

 /**
 * Returns a {@code String} representation of this {@code Vector}.
 *
 * 
 * The string representation of a vector is its components ordered in
 * brackets. For example, {@code new Vector(2, 4).toString()} will return
 * "{@code (2,4)}".
 * 
 *
 * @return string representation of this {@code Vector}.
 */
 @Override
 public String toString()
 {
 NumberFormat nf = NumberFormat.getInstance();
 nf.setMinimumFractionDigits(0);

 StringBuilder vComps = new StringBuilder();

 // Don't add a comma after last component - (1,2,3,4)
 for (int i = 0; i < components.length - 1; i ++)
 {
 vComps.append(
 nf.format( components[i] )
 ).append(",");
 }
 vComps.append(
 nf.format( components[components.length - 1] )
 );

 StringBuilder builder = new StringBuilder();
 builder.append("(").append(vComps).append(")");

 return builder.toString();
 }

 /**
 * Test
 */
 public static void test()
 {
 Vector i = new Vector(2, 2);
 Vector j = new Vector(2, 2);

 double angle = i.angleBetween(j);

 System.out.format("Angle is %s%n", angle);

 Random r = new Random();

 double x = r.nextDouble() * r.nextInt(Integer.MAX_VALUE);
 double y = r.nextDouble() * r.nextInt(Integer.MAX_VALUE);
 double z = r.nextDouble() * r.nextInt(Integer.MAX_VALUE);

 System.out.format("x is %s; y is %s; z is %s%n", x, y, z);
 System.out.format(
 "Math.ulp(x) is %s; Math.ulp(y) is %s; Math.ulp(z) is %s%n",
 Math.ulp(x), Math.ulp(y), Math.ulp(z));
 }

 /**
 * Driver method
 *
 * @param args takes several args: 
 * (1) operation. Is one of {@code add}, {@code dot}, {@code angle},
 * @code none} 
 * (2) number of dimensions each vector will have 
 * (3) first vector, each component separated by a space, so for a
 * vector (2,3,4) you would type
 * {@code 2 3 4} 
 * (4) second vector, same as above 
 */
 public static void main(String[] args)
 {
 List<String> argList = Arrays.asList(args);

 // Perform no-arg test if run w/o args
 if (argList.isEmpty())
 {
 test();
 }

 // Needs at least 4 arguments
 if (argList.size() < 4)
 {
 System.out.format("Not enough arguments. Terminated.%n");
 System.exit(1);
 }

 final int NONE = 0;
 final int ADD = 1;
 final int DOT = 2;
 final int ANGLE = 3;

 int op = NONE;

 // Parse operation argument
 if (argList.get(0).equalsIgnoreCase("none"))
 {
 op = NONE;
 }
 else if (argList.get(0).equalsIgnoreCase("add"))
 {
 op = ADD;
 }
 else if (argList.get(0).equalsIgnoreCase("dot"))
 {
 op = DOT;
 }
 else if (argList.get(0).equalsIgnoreCase("angle"))
 {
 op = ANGLE;
 }
 else
 {
 System.out.format("Operation must be one of none, add, dot, or angle");
 System.exit(1);
 }

 // Determine amount of dimensions
 int dims = 0;
 try
 {
 dims = Integer.parseInt(argList.get(1));
 }
 catch (NumberFormatException e)
 {
 System.out.format("Second parameter must be an integer." +
 "You supplied %s.%n", argList.get(1));
 System.exit(1);
 }

 // Make sure there are d * 2 + 2 arguments (that is, 2 vectors
 // (each has 'dims' dimensions) and first 2 args)
 if (argList.size() < dims * 2 + 2)
 {
 System.out.format(
 "Needs %s arguments. Each vector's component is a separate argument." +
 "You supplied %s arguments.", dims * 2 + 2, args.length);
 System.exit(1);
 }

 // Read components into two lists
 double[] compsA = new double[dims];
 double[] compsB = new double[dims];

 for (int i = 2; i < argList.size(); i++)
 {
 //Try to parse component
 double compArg = 0;
 try
 {
 compArg = Double.parseDouble(argList.get(i));
 }
 catch (NumberFormatException e)
 {
 System.out.format(
 "Component %s is not valid. It must be a decimal.",
 argList.get(i));
 System.exit(1);
 }

 // Test if component is for first or second vector
 if (i - 2 < dims)
 {
 compsA[i - 2] = compArg;
 }
 else
 {
 compsB[i - dims - 2] = compArg;
 }
 }

 Vector a = new Vector(compsA);
 Vector b = new Vector(compsB);

 Object result = new String("Error!");

 switch (op)
 {
 case NONE:
 result = new String("No operation");
 break;
 case ADD:
 result = a.add(b);
 break;
 case DOT:
 result = a.dot(b);
 break;
 case ANGLE:
 result = a.angleBetween(b);
 break;
 default:
 break;
 }

 System.out.format(
 "Vector a is %s. It's magnitude is %s and it's UnitVector is %s.%n",
 a, a.getMagnitude(), a.unitVector());
 System.out.format(
 "Vector b is %s. It's magnitude is %s and it's UnitVector is %s.%n",
 b, b.getMagnitude(), b.unitVector());
 System.out.format("Operation result: %s.%n", result);
 }
}

Author:	DIIST [ Fri Jun 29, 2007 4:54 pm ]
Post subject:	Re: Floating point operation accuracy
What language? I know Java has a issues with floating point numbers. Post the code it might help.

Author:	md [ Fri Jun 29, 2007 7:21 pm ]
Post subject:	RE:Floating point operation accuracy
It doesn't matter the language; it's actually quite clear what the problem is. Read http://en.wikipedia.org/wiki/Floating_point (see the accuracy section). The solution would be to incorporate some error checking code. If you know that the answer is really small, then you can safely assume it's 0. Likewise 45.00000000000001 is really close to 45.00 so you know it should be 45.00. Unfortunately that's only really possible for specific cases, fortunately it's only really needed in specific cases. 0, 90, 180, 270. The rest of the angles would have errors so small as to be un-noticeable. The other solution which would make the errors smaller is to use doubles... but they are twice as large and perhaps not as fast.

Author:	wtd [ Fri Jun 29, 2007 7:38 pm ]
Post subject:	RE:Floating point operation accuracy
Use true fractional types to get accuracy.