InverseGaussian.java

// http://ccs.adnu.edu.ph/acm-icpc2010/prob09_inverse_gaussian.pdf

import java.io.*;
import java.util.*;

public class InverseGaussian {
    
    private static double accuracyNormal = 0.000000000001;
    private static boolean diagnostic = false;
    
    public static void main(String[] args) throws Exception {
        (new InverseGaussian()).run();
        //System.out.println((new InverseGaussian()).G(0.121393675927));
        //System.out.println((new InverseGaussian()).InvG(0.1208));
        //System.out.println((new InverseGaussian()).G(1.224137298045));
    }
    
    private void run() throws Exception {
        BufferedReader br = new BufferedReader(new FileReader("prob09.in"));
        
        String input;
        int caseCount = Integer.parseInt(br.readLine());
        for (int caseNo = 1; caseNo <= caseCount; caseNo++) {
            input = br.readLine();
            double x = Double.parseDouble(input);
            double ans = InvG(x);
            System.out.println("Problem Set " + caseNo + ": " + String.format("%.12f", ans));
        }
    }
    
    private double InvG(double A) throws Exception {
        if (diagnostic) System.out.println("InvG(" + A + ")");
        
        double lower = 0.0001;
        double lowerG = G(lower);
        double upper = 0.8862;
        double upperG = G(upper);
        double x = 0;
        double prevX = 0;
        
        double bisectionPoint = 0.5;
        double intervalLength;
        double xG;
        double delta = 1;
        int execLimit = 0;
        do {
            if (execLimit++ > 500) {
                if (diagnostic) System.out.println("EXECUTION LIMIT REACHED");
                break;
            }
            prevX = x;
            
            intervalLength = upper - lower;
            x = bisectionPoint * intervalLength + lower;
            
            xG = G(x);
            
            if (diagnostic) System.out.println("   G(" + x + ") = " + xG + " >> lower=" + lower + ", upper=" + upper + ", bisectionPoint=" + bisectionPoint);
            
            delta = A - xG;
            if (diagnostic) System.out.println("delta = " + delta);
            //bisectionPoint = (x - lower) / intervalLength;
            if (delta < 0) {
                if (diagnostic) System.out.println("LOWER!");
                upper = x;
                continue;
            }
            else if (delta > 0) {
                if (diagnostic) System.out.println("HIGHER!");
                lower = x;
                continue;                   
            }
            else {
                break;
            }
        } while (Math.abs(delta) > accuracyNormal || Math.abs(x - prevX) > accuracyNormal);
        
        if (diagnostic) System.out.println("Loops: " + execLimit);
        
        return x;
    }
    
    // Converge the bisection point parabolically downward.
                //bisectionPoint = Math.sqrt(x - lower) / intervalLength;
    // Converge the bisection point exponentially upward.
                /*double amount = Math.pow(upper - x, 2.0) / intervalLength;
                if (amount > 1) {
                    bisectionPoint = 0.95;
                }
                else {
                    bisectionPoint = amount;
                }
                */
    
    private double G(double x) throws Exception {
        //System.out.println("G(" + x + ")");
        
        // Get integration of integrand e^(-u^2) from 0 to x.
        // (accurate to 12 decimal places)
        int integrandTermIndex = 0;
        double ret = 0;
        double termValue;
        do {
            // Compute for F(x).
            Pair<Double, Double> term = GIntegrandTermIntegrated(integrandTermIndex);
            termValue = term.Item1 * Math.pow(x, term.Item2);
            ret += termValue;
            //System.out.println("ret=" + ret + ", index=" + integrandTermIndex + ", termValue=" + termValue);
            
            integrandTermIndex++;
        } while (Math.abs(termValue) > accuracyNormal);
        if (ret == Double.NaN) {
            throw new Exception();
        }
        return ret;
    }
    
    private Pair<Double, Double> GIntegrandTermIntegrated(int index) throws IndexOutOfBoundsException {
        if (index < 0) {
            throw new IndexOutOfBoundsException();
        }
        double power = index * 2;
        double constant = 0.0;
        try {
             constant = 1 / FactorialSolver.factorial(index);    
        }
        catch (Exception ex) {
        }
        
        constant *= (index % 2 == 0) ? 1 : -1;
        
        double integratedPower = power + 1;
        double integratedMultiplier = 1 / (power + 1);
        
        double integratedConstant = constant * integratedMultiplier;
        
        // item 1 - constant
        // item 2 - power
        Pair<Double, Double> ret = new Pair<Double, Double>(integratedConstant, integratedPower);
        return ret;
    }
    
    private class Pair<T1, T2> {
        public T1 Item1;
        public T2 Item2;
        
        public Pair() {
        }
        
        public Pair(T1 item1, T2 item2) {
            Item1 = item1;
            Item2 = item2;
        }
    }
}