#!/usr/bin/env --split-string=${JDK_HOME}/bin/java --source ${JDK_VERSION}
package scratch;

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;

import static java.lang.Math.min;
import static java.lang.Integer.parseInt;
import static java.lang.System.err;
import static java.lang.System.exit;
import static java.lang.System.out;
import static java.math.BigDecimal.ONE;
import static java.math.BigDecimal.ZERO;
import static java.math.RoundingMode.HALF_EVEN;


/** A calculator of the likelihood (P) of an event of decaying probability.  The calculation is a sum
  * over a large number of trials (t = 1, 2, 3, …) using two parameters: the initial per-trial
  * probability 0 ≤ p ≤ 1 of the event; and the per-trial decay rate 0 ≤ δ ≤ 1 of the probability
  * (applied after each trial), where δ = 1 yields a constant probability.
  *
  * <p>Calculations are done using numbers either of a given precision (number of digits), or unlimited
  * precision (the default).  Results are rounded to 128 digits for output.</p>*//*
  *
  *   https://math.stackexchange.com/questions/4692693/bernoulli-series-with-decaying-probability-p
  */
class PPrecise {


    public static void main( final String[] arguments ) {
        exitOnDemand( arguments );
        final int n = arguments.length;
        final int precision;
        if( n == 3 ) precision = parseInt( arguments[2] );
        else {
            precision = 0/*unlimited*/;
            if( n != 2 ) {
                err.println( commandName + ": Expecting 2 or 3 arguments, found " + n );
                exitWithUsage( 1 ); }}
        new PPrecise().run( /*p*/new BigDecimal(arguments[0]),
                                     /*δ*/new BigDecimal(arguments[1]), precision ); }



////  P r i v a t e  ////////////////////////////////////////////////////////////////////////////////////


    static final String commandName = "P-precise";



    static void exitOnDemand( final String[] arguments ) {
        for( final String arg: arguments ) {
            if( arg.equals("-?") || arg.equals("-help") ) exitWithUsage( 0 ); }}



    static void exitWithUsage( final int status ) {
        out.println( commandName +
          ": Likelihood of an event of initial probability p, decaying at rate δ" );
        out.println( "Usage: " + commandName + " p δ [precision]" );
        out.println( "       " + commandName + " -help | -?" );
        exit( status ); }



    static final int maxPrintScale = 128;



    /** Cumulative likelihood of occurence.
      */
    BigDecimal pCumulative = ZERO;



    static final RoundingMode rounding = HALF_EVEN;



    /** Per-trial probability at last trial.
      */
    BigDecimal pLast = ZERO;



    /** Cumulative likelihood of failure.
      */
    BigDecimal qCumulative = ONE;



    /** @param p Per-trial probability, initial value.
      * @param δ Per-trial decay rate of `p`.
      */
    void run( BigDecimal p, final BigDecimal δ, final int precision ) {
        testParameter( p );
        testParameter( δ );
        final MathContext c = new MathContext( precision, rounding);
        final BigDecimal δInverse = ONE.subtract( δ );
        do {
            qCumulative = qCumulative.multiply( ONE.subtract( pLast ));
            final BigDecimal pJustNow/*earlier failure *and* occurence on the present trial*/ =
              qCumulative/*likelihood of earlier failure*/.multiply( p/*per-trial probability*/, c );
            pCumulative = pCumulative.add( pJustNow );
            if( t == tReport ) {
                final int scale = min( maxPrintScale, pCumulative.scale() ); // This does nothing.
                out.println( "t " + t + "     \t" + pCumulative.setScale( scale, rounding ));
                if( tReport < 20 ) tReport += 1;
                else if (tReport == 20 ) tReport = trialReportInterval;
                else tReport += trialReportInterval; }
            pLast = p;
            p = p.multiply( δInverse, c ); }
            while( ++t <= trialsMaximum ); }



    /** Trial number.
      */
    int t = 1;



    static void testParameter( final BigDecimal d ) {
        if( d.signum() < 0 || d.compareTo(ONE) > 0 ) {
            throw new IllegalArgumentException( "Parameter outside range 0 to 1" ); }}



    /** Trial number for the next report.
      */
    int tReport = t;



    static final int trialReportInterval = 1000;



    static final int trialsMaximum = trialReportInterval * 12; }