001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.rng.examples.quadrature;
018
019import org.apache.commons.rng.UniformRandomProvider;
020import org.apache.commons.rng.simple.RandomSource;
021
022/**
023 * Computation of \( \pi \) using Monte-Carlo integration.
024 *
025 * The computation estimates the value by computing the probability that
026 * a point \( p = (x, y) \) will lie in the circle of radius \( r = 1 \)
027 * inscribed in the square of side \( r = 1 \).
028 * The probability could be computed by \[ area_{circle} / area_{square} \],
029 * where \( area_{circle} = \pi * r^2 \) and \( area_{square} = 4 r^2 \).
030 * Hence, the probability is \( \frac{\pi}{4} \).
031 *
032 * The Monte Carlo simulation will produce \( N \) points.
033 * Defining \( N_c \) as the number of point that satisfy \( x^2 + y^2 \le 1 \),
034 * we will have \( \frac{N_c}{N} \approx \frac{\pi}{4} \).
035 */
036public class ComputePi extends MonteCarloIntegration {
037    /** Expected number of arguments. */
038    private static final int EXPECTED_ARGUMENTS = 2;
039    /** Domain dimension. */
040    private static final int DIMENSION = 2;
041
042    /**
043     * Create an instance.
044     *
045     * @param rng RNG.
046     */
047    public ComputePi(UniformRandomProvider rng) {
048        super(rng, DIMENSION);
049    }
050
051    /**
052     * Program entry point.
053     *
054     * @param args Arguments.
055     * The order is as follows:
056     * <ol>
057     *  <li>
058     *   Number of random 2-dimensional points to generate.
059     *  </li>
060     *  <li>
061     *   {@link RandomSource Random source identifier}.
062     *  </li>
063     * </ol>
064     */
065    public static void main(String[] args) {
066        if (args.length != EXPECTED_ARGUMENTS) {
067            throw new IllegalStateException("Require arguments: [points] [RNG name]");
068        }
069
070        final long numPoints = Long.parseLong(args[0]);
071        final RandomSource randomSource = RandomSource.valueOf(args[1]);
072
073        final ComputePi piApp = new ComputePi(randomSource.create());
074        final double piMC = piApp.compute(numPoints);
075
076        //CHECKSTYLE: stop all
077        System.out.printf("After generating %d random numbers, the error on |𝛑 - %s| is %s%n",
078                          DIMENSION * numPoints, piMC, Math.abs(piMC - Math.PI));
079        //CHECKSTYLE: resume all
080    }
081
082    /**
083     * Compute the value of pi.
084     *
085     * @param numPoints Number of random points to generate.
086     * @return the approximate value of pi.
087     */
088    public double compute(long numPoints) {
089        return 4 * integrate(numPoints);
090    }
091
092    /** {@inheritDoc} */
093    @Override
094    protected boolean isInside(double... rand) {
095        final double r2 = rand[0] * rand[0] + rand[1] * rand[1];
096        return r2 <= 1;
097    }
098}