From 71df39e4900c80397611f03ba479fd5dbaccc9fd Mon Sep 17 00:00:00 2001 From: axexlck Date: Sat, 10 Aug 2024 12:06:45 +0200 Subject: [PATCH] WIP #1034 `NSum` for complex numbers - added *Complexx.java variants --- symja_android_library/doc/functions/NSum.md | 3 +- .../core/numerics/series/dp/Aitken.java | 174 +++--- .../core/numerics/series/dp/Cohen.java | 87 ++- .../core/numerics/series/dp/Direct.java | 28 +- .../core/numerics/series/dp/Ensemble.java | 470 ++++++++------- .../numerics/series/dp/IteratedTheta.java | 92 +-- .../core/numerics/series/dp/Richardson.java | 562 +++++++++--------- .../numerics/series/dp/SeriesAlgorithm.java | 88 +++ .../series/dp/complex/AitkenComplex.java | 117 ++++ .../dp/complex/BrezinskiThetaComplex.java | 95 +++ .../series/dp/complex/CohenComplex.java | 68 +++ .../series/dp/complex/DirectComplex.java | 26 + .../series/dp/complex/EnsembleComplex.java | 279 +++++++++ .../dp/complex/IteratedThetaComplex.java | 65 ++ .../series/dp/complex/LevinComplex.java | 164 +++++ .../series/dp/complex/RichardsonComplex.java | 315 ++++++++++ .../dp/complex/SeriesAlgorithmComplex.java | 256 ++++++++ .../series/dp/complex/WynnEpsilonComplex.java | 288 +++++++++ .../series/dp/complex/WynnRhoComplex.java | 109 ++++ .../core/reflection/system/NSum.java | 150 +++-- .../src/main/resources/doc/functions/NSum.md | 25 + .../org/matheclipse/core/system/SumTest.java | 12 +- 22 files changed, 2695 insertions(+), 778 deletions(-) create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/AitkenComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/BrezinskiThetaComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/CohenComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/DirectComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/EnsembleComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/IteratedThetaComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/LevinComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/RichardsonComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/SeriesAlgorithmComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/WynnEpsilonComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/WynnRhoComplex.java create mode 100644 symja_android_library/matheclipse-core/src/main/resources/doc/functions/NSum.md diff --git a/symja_android_library/doc/functions/NSum.md b/symja_android_library/doc/functions/NSum.md index b540a5a876..c0223c8d89 100644 --- a/symja_android_library/doc/functions/NSum.md +++ b/symja_android_library/doc/functions/NSum.md @@ -13,12 +13,13 @@ NSum(expr, {i, imin, imax, di}) > `i` ranges from `imin` to `imax` in steps `di` See +* [Wikipedia - Series (mathematics)](https://en.wikipedia.org/wiki/Series_(mathematics)) * [Wikipedia - Summation](https://en.wikipedia.org/wiki/Summation) * [Wikipedia - Convergence_tests](https://en.wikipedia.org/wiki/Convergence_tests) ### Examples ``` ->> Num(k, {k, 1, 10}) +>> NSum(k, {k, 1, 10}) 55 ``` diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Aitken.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Aitken.java index 3864399527..ab58d1c93b 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Aitken.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Aitken.java @@ -1,106 +1,106 @@ package org.matheclipse.core.numerics.series.dp; /** - * Implements a modified adaptive Aitken delta^2 process for estimating infinite - * series, based on the method in [1]. + * Implements a modified adaptive Aitken delta^2 process for estimating infinite series, based on + * the method in [1]. * *

* References: *

*

*/ public final class Aitken extends SeriesAlgorithm { - private int status; - private double xx, dx, dd, xp, alpha; - private final double[] x; - private final double[][] s, ds, ts, dts; + private int status; + private double xx, dx, dd, xp, alpha; + private final double[] x; + private final double[][] s, ds, ts, dts; - public Aitken(final double tolerance, final int maxIters, final int patience) { - super(tolerance, maxIters, patience); - x = new double[maxIters]; - s = new double[2][maxIters + 1]; - ds = new double[2][maxIters + 1]; - ts = new double[2][maxIters + 1]; - dts = new double[2][maxIters + 1]; - } + public Aitken(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + x = new double[maxIters]; + s = new double[2][maxIters + 1]; + ds = new double[2][maxIters + 1]; + ts = new double[2][maxIters + 1]; + dts = new double[2][maxIters + 1]; + } - @Override - public final double next(final double e, final double term) { - if (myIndex == 0) { - xx = dx = dd = 0.0; - } - final double xxp = xx; - x[myIndex] = xx = term; - if (myIndex == 0) { - dx = xx; - } else if (myIndex == 1) { - final double dxp = dx; - dx = xx - xxp; - dd = dx / (dx - dxp); - } else { - final double dxp = dx; - final double ddp = dd; - dx = xx - xxp; - dd = dx / (dx - dxp); - alpha = 1.0 / (dd - ddp) + 1.0; - alpha = aitken(alpha, ts, dts, myIndex - 1, -2.0); - } - xp = xx; - if (myIndex >= 2) { - for (int m = 1; m <= myIndex + 1; ++m) { - xp = aitken(x[m - 1], s, ds, m, alpha); - } - } - ++myIndex; - if (status == 1) { - return Double.NaN; - } else { - return xp; - } + @Override + public final double next(final double e, final double term) { + if (myIndex == 0) { + xx = dx = dd = 0.0; } - - private final double aitken(final double xx, final double[][] s, final double[][] ds, final int n, - final double theta) { - status = 0; - if (n != 1) { - final int kend = (n - 1) >> 1; - System.arraycopy(s[1], 0, s[0], 0, kend + 1); - if ((n & 1) == 0) { - System.arraycopy(ds[1], 0, ds[0], 0, kend); - } else { - System.arraycopy(ds[1], 0, ds[0], 0, kend + 1); - } - } - s[1][0] = xx; - if (n == 1) { - ds[1][0] = xx; - return xx; - } - ds[1][0] = xx - s[0][0]; - for (int k = 1; k <= (n >> 1); ++k) { - final double w1 = ds[0][k - 1] * ds[1][k - 1]; - final double w2 = ds[1][k - 1] - ds[0][k - 1]; - if (Math.abs(w2) < TINY) { - status = 1; - return xx; - } - final int twok = k << 1; - final double coef = ((twok - 1) - theta) / (twok - 2 - theta); - s[1][k] = s[0][k - 1] - coef * w1 / w2; - if (n != twok - 1) { - ds[1][k] = s[1][k] - s[0][k]; - } - } - final int kopt = n >> 1; - return s[1][kopt]; + final double xxp = xx; + x[myIndex] = xx = term; + if (myIndex == 0) { + dx = xx; + } else if (myIndex == 1) { + final double dxp = dx; + dx = xx - xxp; + dd = dx / (dx - dxp); + } else { + final double dxp = dx; + final double ddp = dd; + dx = xx - xxp; + dd = dx / (dx - dxp); + alpha = 1.0 / (dd - ddp) + 1.0; + alpha = aitken(alpha, ts, dts, myIndex - 1, -2.0); + } + xp = xx; + if (myIndex >= 2) { + for (int m = 1; m <= myIndex + 1; ++m) { + xp = aitken(x[m - 1], s, ds, m, alpha); + } + } + ++myIndex; + if (status == 1) { + return Double.NaN; + } else { + return xp; } + } - @Override - public final String getName() { - return "Modified Aitken"; + private final double aitken(final double xx, final double[][] s, final double[][] ds, final int n, + final double theta) { + status = 0; + if (n != 1) { + final int kend = (n - 1) >> 1; + System.arraycopy(s[1], 0, s[0], 0, kend + 1); + if ((n & 1) == 0) { + System.arraycopy(ds[1], 0, ds[0], 0, kend); + } else { + System.arraycopy(ds[1], 0, ds[0], 0, kend + 1); + } } + s[1][0] = xx; + if (n == 1) { + ds[1][0] = xx; + return xx; + } + ds[1][0] = xx - s[0][0]; + for (int k = 1; k <= (n >> 1); ++k) { + final double w1 = ds[0][k - 1] * ds[1][k - 1]; + final double w2 = ds[1][k - 1] - ds[0][k - 1]; + if (Math.abs(w2) < TINY) { + status = 1; + return xx; + } + final int twok = k << 1; + final double coef = ((twok - 1) - theta) / (twok - 2 - theta); + s[1][k] = s[0][k - 1] - coef * w1 / w2; + if (n != twok - 1) { + ds[1][k] = s[1][k] - s[0][k]; + } + } + final int kopt = n >> 1; + return s[1][kopt]; + } + + @Override + public final String getName() { + return "Modified Aitken"; + } } diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Cohen.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Cohen.java index 0ba1c1aca7..de2573ce30 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Cohen.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Cohen.java @@ -4,65 +4,64 @@ import org.matheclipse.core.numerics.utils.SimpleMath; /** - * This is an implementation of Algorithm 1 in [1] for evaluating infinite - * series whose terms have alternating signs. + * This is an implementation of Algorithm 1 in [1] for evaluating infinite series whose terms have + * alternating signs. * *

* References: *

*

*/ public final class Cohen extends SeriesAlgorithm { - private final double[] myTab; - private double mySign0; + private final double[] myTab; + private double mySign0; - public Cohen(final double tolerance, final int maxIters, final int patience) { - super(tolerance, maxIters, patience); - myTab = new double[maxIters]; - } + public Cohen(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + myTab = new double[maxIters]; + } - @Override - public final double next(final double e, final double term) { + @Override + public final double next(final double e, final double term) { - // add next element - final int n = myIndex + 1; - myTab[myIndex] = Math.abs(e); + // add next element + final int n = myIndex + 1; + myTab[myIndex] = Math.abs(e); - // record the sign of the first term, since this method - // requires the first term of the sequence to be positive - if (myIndex == 0) { - mySign0 = Math.signum(e); - if (mySign0 == 0.0) { - mySign0 = 1.0; - } - } + // record the sign of the first term, since this method + // requires the first term of the sequence to be positive + if (myIndex == 0) { + mySign0 = Math.signum(e); + if (mySign0 == 0.0) { + mySign0 = 1.0; + } + } - // initialize d value - double d = SimpleMath.pow(3.0 + 2.0 * Constants.SQRT2, n); - d = 0.5 * (d + 1.0 / d); + // initialize d value + double d = SimpleMath.pow(3.0 + 2.0 * Constants.SQRT2, n); + d = 0.5 * (d + 1.0 / d); - // apply the Chebychef polynomial recursively (Algorithm 1) - double b = -1.0; - double c = -d; - double s = 0.0; - for (int k = 0; k < n; ++k) { - c = b - c; - s += c * myTab[k]; - final double numer = (k + n) * (k - n); - final double denom = (k + 0.5) * (k + 1); - b *= numer / denom; - } - ++myIndex; - return mySign0 * s / d; + // apply the Chebychef polynomial recursively (Algorithm 1) + double b = -1.0; + double c = -d; + double s = 0.0; + for (int k = 0; k < n; ++k) { + c = b - c; + s += c * myTab[k]; + final double numer = (k + n) * (k - n); + final double denom = (k + 0.5) * (k + 1); + b *= numer / denom; } + ++myIndex; + return mySign0 * s / d; + } - @Override - public final String getName() { - return "Cohen"; - } + @Override + public final String getName() { + return "Cohen"; + } } diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Direct.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Direct.java index 41543bd671..0415851cfc 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Direct.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Direct.java @@ -1,23 +1,23 @@ package org.matheclipse.core.numerics.series.dp; /** - * Naively computes the limit of a sequence or series by inspecting elements one - * a time. No transform of the original sequence is performed. + * Naively computes the limit of a sequence or series by inspecting elements one a time. No + * transform of the original sequence is performed. */ public final class Direct extends SeriesAlgorithm { - public Direct(final double tolerance, final int maxIters, final int patience) { - super(tolerance, maxIters, patience); - } + public Direct(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + } - @Override - public final double next(final double e, final double term) { - ++myIndex; - return term; - } + @Override + public final double next(final double e, final double term) { + ++myIndex; + return term; + } - @Override - public final String getName() { - return "Direct Sum"; - } + @Override + public final String getName() { + return "Direct Sum"; + } } diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Ensemble.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Ensemble.java index ba0faf47aa..e546c18408 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Ensemble.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Ensemble.java @@ -3,278 +3,276 @@ import java.util.ArrayList; import java.util.List; import org.matheclipse.core.numerics.series.dp.Levin.RemainderSequence; +import org.matheclipse.core.numerics.series.dp.SeriesAlgorithm.SeriesSolution; import org.matheclipse.core.numerics.series.dp.WynnEpsilon.ShanksMethod; /** - * An ensemble algorithm for evaluating the limits of sequences and series of - * real values. This algorithms works by initializing a number of convergence - * acceleration algorithms with different properties, and running them in - * parallel. The estimates of the sequence's or corresponding series limit are - * extracted from the instance that has converged first to a stable limit to the - * desired tolerance. + * An ensemble algorithm for evaluating the limits of sequences and series of real values. This + * algorithms works by initializing a number of convergence acceleration algorithms with different + * properties, and running them in parallel. The estimates of the sequence's or corresponding series + * limit are extracted from the instance that has converged first to a stable limit to the desired + * tolerance. */ public final class Ensemble extends SeriesAlgorithm { - private static final int PRINT_DIGITS = 18; + private static final int PRINT_DIGITS = 18; - private final int myPrint; - private final List myMethods; - private final List myForAlternating; - private final List myForSequence; + private final int myPrint; + private final List myMethods; + private final List myForAlternating; + private final List myForSequence; - private int myMethodCount, myBestMethod; - private boolean myIsAlternating, myCheckForSequenceOnly; - private double myOldSignum; - private boolean[] myExcept; - private double[] myPrevEst, myEst; + private int myMethodCount, myBestMethod; + private boolean myIsAlternating, myCheckForSequenceOnly; + private double myOldSignum; + private boolean[] myExcept; + private double[] myPrevEst, myEst; - /** - * Creates a new instance of an ensemble convergence acceleration algorithm. - * - * @param tolerance the smallest acceptable change in series evaluation in - * consecutive iterations that indicates the algorithm has - * converged - * @param maxIters the maximum number of sequence terms to evaluate before - * giving up - * @param patience the number of iterations used to assess whether or not a - * particular instance of convergence acceleration - * algorithm used by this algorithm has converged - * @param printProgress how often to print the progress of the current - * algorithm, e.g. the estimate of a sequence's or series's - * limit at each iteration according to each algorithm - * instance - */ - public Ensemble(final double tolerance, final int maxIters, final int patience, final int printProgress) { - super(tolerance, maxIters, patience); - myPrint = printProgress; - myMethods = new ArrayList<>(); - myForAlternating = new ArrayList<>(); - myForSequence = new ArrayList<>(); + /** + * Creates a new instance of an ensemble convergence acceleration algorithm. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of iterations used to assess whether or not a particular instance of + * convergence acceleration algorithm used by this algorithm has converged + * @param printProgress how often to print the progress of the current algorithm, e.g. the + * estimate of a sequence's or series's limit at each iteration according to each algorithm + * instance + */ + public Ensemble(final double tolerance, final int maxIters, final int patience, + final int printProgress) { + super(tolerance, maxIters, patience); + myPrint = printProgress; + myMethods = new ArrayList<>(); + myForAlternating = new ArrayList<>(); + myForSequence = new ArrayList<>(); - // only for alternating series - myMethods.add(new Cohen(myTol, myMaxIters, 1)); - myForAlternating.add(true); - myForSequence.add(false); - myMethods.add(new Levin(myTol, myMaxIters, 1, RemainderSequence.D)); - myForAlternating.add(true); - myForSequence.add(false); + // only for alternating series + myMethods.add(new Cohen(myTol, myMaxIters, 1)); + myForAlternating.add(true); + myForSequence.add(false); + myMethods.add(new Levin(myTol, myMaxIters, 1, RemainderSequence.D)); + myForAlternating.add(true); + myForSequence.add(false); - // for alternating and linear series - myMethods.add(new Aitken(myTol, myMaxIters, 1)); - myForAlternating.add(false); - myForSequence.add(true); - myMethods.add(new WynnEpsilon(myTol, myMaxIters, 1, ShanksMethod.WYNN)); - myForAlternating.add(false); - myForSequence.add(true); - myMethods.add(new Levin(myTol, myMaxIters, 1, RemainderSequence.T)); - myForAlternating.add(false); - myForSequence.add(false); + // for alternating and linear series + myMethods.add(new Aitken(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new WynnEpsilon(myTol, myMaxIters, 1, ShanksMethod.WYNN)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new Levin(myTol, myMaxIters, 1, RemainderSequence.T)); + myForAlternating.add(false); + myForSequence.add(false); - // for logarithmic series - myMethods.add(new Levin(myTol, myMaxIters, 1, RemainderSequence.U)); - myForAlternating.add(false); - myForSequence.add(false); - myMethods.add(new Levin(myTol, myMaxIters, 1, RemainderSequence.V)); - myForAlternating.add(false); - myForSequence.add(false); - myMethods.add(new Richardson(myTol, myMaxIters, 1)); - myForAlternating.add(false); - myForSequence.add(false); - myMethods.add(new BrezinskiTheta(myTol, myMaxIters, 1)); - myForAlternating.add(false); - myForSequence.add(true); - myMethods.add(new IteratedTheta(myTol, myMaxIters, 1)); - myForAlternating.add(false); - myForSequence.add(true); - myMethods.add(new WynnRho(myTol, myMaxIters, 1)); - myForAlternating.add(false); - myForSequence.add(true); - myMethods.add(new WynnEpsilon(myTol, myMaxIters, 1, ShanksMethod.ALTERNATING)); - myForAlternating.add(false); - myForSequence.add(true); - } - - public Ensemble(final double tolerance, final int maxIters, final int patience) { - this(tolerance, maxIters, patience, 0); - } + // for logarithmic series + myMethods.add(new Levin(myTol, myMaxIters, 1, RemainderSequence.U)); + myForAlternating.add(false); + myForSequence.add(false); + myMethods.add(new Levin(myTol, myMaxIters, 1, RemainderSequence.V)); + myForAlternating.add(false); + myForSequence.add(false); + myMethods.add(new Richardson(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(false); + myMethods.add(new BrezinskiTheta(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new IteratedTheta(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new WynnRho(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new WynnEpsilon(myTol, myMaxIters, 1, ShanksMethod.ALTERNATING)); + myForAlternating.add(false); + myForSequence.add(true); + } - @Override - public final double next(final double e, final double term) { + public Ensemble(final double tolerance, final int maxIters, final int patience) { + this(tolerance, maxIters, patience, 0); + } - // ******************************************************************** - // INITIALIZATION - // ******************************************************************** - if (myIndex == 0) { + @Override + public final double next(final double e, final double term) { - // initialize variables - myIsAlternating = true; - myMethodCount = myMethods.size(); - myOldSignum = Math.signum(e); - myPrevEst = new double[myMethodCount]; - myEst = new double[myMethodCount]; - myExcept = new boolean[myMethodCount]; + // ******************************************************************** + // INITIALIZATION + // ******************************************************************** + if (myIndex == 0) { - // disable methods for series only - if (myCheckForSequenceOnly) { - for (int m = 0; m < myMethodCount; ++m) { - if (!myForSequence.get(m)) { - myExcept[m] = true; - if (myPrint > 0) { - System.out.println("Disabling method " + myMethods.get(m).getName() - + " that is invalid for sequences"); - } - } - } - myCheckForSequenceOnly = false; - } + // initialize variables + myIsAlternating = true; + myMethodCount = myMethods.size(); + myOldSignum = Math.signum(e); + myPrevEst = new double[myMethodCount]; + myEst = new double[myMethodCount]; + myExcept = new boolean[myMethodCount]; - // initialize each method - for (int m = 0; m < myMethodCount; ++m) { - if (!myExcept[m]) { - final SeriesAlgorithm method = myMethods.get(m); - method.myIndex = 0; - myEst[m] = method.next(e, term); - } - } + // disable methods for series only + if (myCheckForSequenceOnly) { + for (int m = 0; m < myMethodCount; ++m) { + if (!myForSequence.get(m)) { + myExcept[m] = true; + if (myPrint > 0) { + System.out.println("Disabling method " + myMethods.get(m).getName() + + " that is invalid for sequences"); + } + } + } + myCheckForSequenceOnly = false; + } - // print header - if (myPrint > 0) { - String line = ""; - line += pad("Index", PRINT_DIGITS + 5); - for (final SeriesAlgorithm method : myMethods) { - line += "\t" + pad(method.getName(), PRINT_DIGITS + 5); - } - System.out.println(line); - } + // initialize each method + for (int m = 0; m < myMethodCount; ++m) { + if (!myExcept[m]) { + final SeriesAlgorithm method = myMethods.get(m); + method.myIndex = 0; + myEst[m] = method.next(e, term); + } + } - // print values - if (myPrint > 0) { - printRow(); - } + // print header + if (myPrint > 0) { + String line = ""; + line += pad("Index", PRINT_DIGITS + 5); + for (final SeriesAlgorithm method : myMethods) { + line += "\t" + pad(method.getName(), PRINT_DIGITS + 5); + } + System.out.println(line); + } - ++myIndex; - return term; - } + // print values + if (myPrint > 0) { + printRow(); + } - // ******************************************************************** - // SEQUENCE INSPECTION - // ******************************************************************** - // test for invalid value of series or sequence - if (!Double.isFinite(term) || !Double.isFinite(e)) { - if (myPrint > 0) { - System.out.println("Aborting due to NaN at iteration " + myIndex); - } - ++myIndex; - return Double.NaN; - } + ++myIndex; + return term; + } - // track the sign of the sequence - final double signum = Math.signum(e); - if (myIndex > 0 && myIsAlternating && signum * myOldSignum >= 0) { - myIsAlternating = false; - } - myOldSignum = signum; + // ******************************************************************** + // SEQUENCE INSPECTION + // ******************************************************************** + // test for invalid value of series or sequence + if (!Double.isFinite(term) || !Double.isFinite(e)) { + if (myPrint > 0) { + System.out.println("Aborting due to NaN at iteration " + myIndex); + } + ++myIndex; + return Double.NaN; + } - // exclude alternating series methods if the series does not alternate - if (!myIsAlternating) { - for (int m = 0; m < myMethodCount; ++m) { - if (!myExcept[m] && myForAlternating.get(m)) { - myExcept[m] = true; - if (myPrint > 0) { - System.out.println("Disabling method " + myMethods.get(m).getName() - + " because series is not alternating at iteration " + myIndex); - } - } - } - } + // track the sign of the sequence + final double signum = Math.signum(e); + if (myIndex > 0 && myIsAlternating && signum * myOldSignum >= 0) { + myIsAlternating = false; + } + myOldSignum = signum; - // ******************************************************************** - // MAIN UPDATE - // ******************************************************************** - myBestMethod = -1; - double bestError = Double.POSITIVE_INFINITY; - for (int m = 0; m < myMethodCount; ++m) { - if (myExcept[m]) { - continue; - } + // exclude alternating series methods if the series does not alternate + if (!myIsAlternating) { + for (int m = 0; m < myMethodCount; ++m) { + if (!myExcept[m] && myForAlternating.get(m)) { + myExcept[m] = true; + if (myPrint > 0) { + System.out.println("Disabling method " + myMethods.get(m).getName() + + " because series is not alternating at iteration " + myIndex); + } + } + } + } - // estimate the next terms - final SeriesAlgorithm method = myMethods.get(m); - myPrevEst[m] = myEst[m]; - myEst[m] = method.next(e, term); + // ******************************************************************** + // MAIN UPDATE + // ******************************************************************** + myBestMethod = -1; + double bestError = Double.POSITIVE_INFINITY; + for (int m = 0; m < myMethodCount; ++m) { + if (myExcept[m]) { + continue; + } - // if a method produces an invalid estimate it is excluded - if (!Double.isFinite(myEst[m])) { - if (myPrint > 0) { - System.out.println( - "Disabling method " + method.getName() + " due to instability at iteration " + myIndex); - } - myExcept[m] = true; - continue; - } + // estimate the next terms + final SeriesAlgorithm method = myMethods.get(m); + myPrevEst[m] = myEst[m]; + myEst[m] = method.next(e, term); - // estimate error and best method - final double error = Math.abs(myPrevEst[m] - myEst[m]); - if (error < bestError && myEst[m] != HUGE) { - myBestMethod = m; - bestError = error; - } - } + // if a method produces an invalid estimate it is excluded + if (!Double.isFinite(myEst[m])) { + if (myPrint > 0) { + System.out.println("Disabling method " + method.getName() + + " due to instability at iteration " + myIndex); + } + myExcept[m] = true; + continue; + } - // ******************************************************************** - // FINALIZE THIS ITERATION - // ******************************************************************** - if (myPrint > 0 && myIndex % myPrint == 0) { - printRow(); - } + // estimate error and best method + final double error = Math.abs(myPrevEst[m] - myEst[m]); + if (error < bestError && myEst[m] != HUGE) { + myBestMethod = m; + bestError = error; + } + } - // return the method with the best convergence so far - ++myIndex; - if (myBestMethod >= 0) { - return myEst[myBestMethod]; - } else { - return Double.NaN; - } + // ******************************************************************** + // FINALIZE THIS ITERATION + // ******************************************************************** + if (myPrint > 0 && myIndex % myPrint == 0) { + printRow(); } - @Override - public final SeriesSolution limit(final Iterable seq, final boolean series, final int extrapolateStart) { - myCheckForSequenceOnly = !series; - final SeriesSolution result = super.limit(seq, series, extrapolateStart); - if (myPrint > 0 && myBestMethod >= 0 && result.converged) { - System.out.println( - "Converged at iteration " + myIndex + " with method " + myMethods.get(myBestMethod).getName()); - } - return result; + // return the method with the best convergence so far + ++myIndex; + if (myBestMethod >= 0) { + return myEst[myBestMethod]; + } else { + return Double.NaN; } + } - @Override - public final String getName() { - return "Ensemble"; + @Override + public final SeriesSolution limit(final Iterable seq, final boolean series, + final int extrapolateStart) { + myCheckForSequenceOnly = !series; + final SeriesSolution result = super.limit(seq, series, extrapolateStart); + if (myPrint > 0 && myBestMethod >= 0 && result.converged) { + System.out.println("Converged at iteration " + myIndex + " with method " + + myMethods.get(myBestMethod).getName()); } + return result; + } - private final void printRow() { - String line = pad(myIndex + "", PRINT_DIGITS); - for (int m = 0; m < myMethodCount; ++m) { - final String str; - if (myExcept[m]) { - str = "-"; - } else { - str = Double.toString(myEst[m]); - } - line += "\t" + pad(str, PRINT_DIGITS); - } - System.out.println(line); + @Override + public final String getName() { + return "Ensemble"; + } + + private final void printRow() { + String line = pad(myIndex + "", PRINT_DIGITS); + for (int m = 0; m < myMethodCount; ++m) { + final String str; + if (myExcept[m]) { + str = "-"; + } else { + str = Double.toString(myEst[m]); + } + line += "\t" + pad(str, PRINT_DIGITS); } + System.out.println(line); + } - private static final String pad(final String str, final int len) { - if (str.length() >= len) { - return str; - } - String result = str; - for (int i = str.length() + 1; i <= len; ++i) { - result += " "; - } - return result; + private static final String pad(final String str, final int len) { + if (str.length() >= len) { + return str; + } + String result = str; + for (int i = str.length() + 1; i <= len; ++i) { + result += " "; } + return result; + } } diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/IteratedTheta.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/IteratedTheta.java index b1ca4f4f60..fb31c77670 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/IteratedTheta.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/IteratedTheta.java @@ -1,65 +1,67 @@ package org.matheclipse.core.numerics.series.dp; +import org.hipparchus.complex.Complex; + /** - * Implements the Iterated Theta algorithm for convergence of infinite series as - * described in [1]. + * Implements the Iterated Theta algorithm for convergence of infinite series as described in [1]. * *

* References: *

    - *
  • [1] Weniger, Ernst Joachim. "Nonlinear sequence transformations for the - * acceleration of convergence and the summation of divergent series." arXiv - * preprint math/0306302 (2003).
    • + *
  • [1] Weniger, Ernst Joachim. "Nonlinear sequence transformations for the acceleration of + * convergence and the summation of divergent series." arXiv preprint math/0306302 (2003).
    • *
*

*/ public final class IteratedTheta extends SeriesAlgorithm { - private final double[] myTab; - - public IteratedTheta(final double tolerance, final int maxIters, final int patience) { - super(tolerance, maxIters, patience); - myTab = new double[maxIters]; - } + private final double[] myTab; + private final Complex[] myComplexTab; - @Override - public final double next(final double e, final double term) { + public IteratedTheta(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + myTab = new double[maxIters]; + myComplexTab = new Complex[maxIters]; + } - // initial value for J - myTab[myIndex] = term; - if (myIndex < 3) { - ++myIndex; - return term; - } + @Override + public final double next(final double e, final double term) { - // higher-order J array - final int lmax = myIndex / 3; - int m = myIndex; - for (int l = 1; l <= lmax; ++l) { - m -= 3; - final double diff0 = myTab[m + 1] - myTab[m]; - final double diff1 = myTab[m + 2] - myTab[m + 1]; - final double diff2 = myTab[m + 3] - myTab[m + 2]; - final double dfsq1 = diff1 - diff0; - final double dfsq2 = diff2 - diff1; - final double denom = diff2 * dfsq1 - diff0 * dfsq2; - final double ratio = dfsq2 / denom; + // initial value for J + myTab[myIndex] = term; + if (myIndex < 3) { + ++myIndex; + return term; + } - // safeguarded against underflow - if (Math.abs(denom) < TINY) { - myTab[m] = HUGE; - } else { - myTab[m] = myTab[m + 1] - diff0 * diff1 * ratio; - } - } - ++myIndex; + // higher-order J array + final int lmax = myIndex / 3; + int m = myIndex; + for (int l = 1; l <= lmax; ++l) { + m -= 3; + final double diff0 = myTab[m + 1] - myTab[m]; + final double diff1 = myTab[m + 2] - myTab[m + 1]; + final double diff2 = myTab[m + 3] - myTab[m + 2]; + final double dfsq1 = diff1 - diff0; + final double dfsq2 = diff2 - diff1; + final double denom = diff2 * dfsq1 - diff0 * dfsq2; + final double ratio = dfsq2 / denom; - // return result - return myTab[myIndex % 3]; + // safeguarded against underflow + if (Math.abs(denom) < TINY) { + myTab[m] = HUGE; + } else { + myTab[m] = myTab[m + 1] - diff0 * diff1 * ratio; + } } + ++myIndex; - @Override - public final String getName() { - return "Iterated Theta"; - } + // return result + return myTab[myIndex % 3]; + } + + @Override + public final String getName() { + return "Iterated Theta"; + } } diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Richardson.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Richardson.java index e76ce8c3b5..d77b9de512 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Richardson.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/Richardson.java @@ -3,323 +3,309 @@ import java.util.Iterator; /** - * Implements a generalized form of the Richardson extrapolation for evaluating - * infinite series. This is implemented using a d-transform of the original - * series and is called the W(m) algorithm [1]. + * Implements a generalized form of the Richardson extrapolation for evaluating infinite series. + * This is implemented using a d-transform of the original series and is called the W(m) algorithm + * [1]. * *

- * The method in [1] defines an auxiliary increasing sequence of integers, R(l), - * defined by the recursion R(l + 1) = (int)(sigma * R(l)) + s, where sigma and - * s are two tuning parameters. The series are evaluated at each R(l) of the - * auxiliary sequence and forms a subsequence of the original partial sums that - * converges to the same limit. The d-transform is then applied to accelerate - * this sequence to infinity. + * The method in [1] defines an auxiliary increasing sequence of integers, R(l), defined by the + * recursion R(l + 1) = (int)(sigma * R(l)) + s, where sigma and s are two tuning parameters. The + * series are evaluated at each R(l) of the auxiliary sequence and forms a subsequence of the + * original partial sums that converges to the same limit. The d-transform is then applied to + * accelerate this sequence to infinity. *

* *

* References: *

    - *
  • [1] Ford, William & Sidi, Avram. (1987). An Algorithm for a - * Generalization of the Richardson Extrapolation Process. SIAM Journal on - * Numerical Analysis. 24. 10.1137/0724080.
    • + *
  • [1] Ford, William & Sidi, Avram. (1987). An Algorithm for a Generalization of the Richardson + * Extrapolation Process. SIAM Journal on Numerical Analysis. 24. 10.1137/0724080.
    • *
*

*/ public final class Richardson extends SeriesAlgorithm { - private final double mySigma; - private final int myIncr; - - private int myPrevIndex; - private final double[] myElements; - - private final int myM; - private final double[] myG; - private final double[][][] myPsiAI, myPsiQ, myPsiG; - private int evals; - - /** - * Creates an instance of the d-transform for a given auxiliary sequence. - * - * @param tolerance the smallest acceptable change in series evaluation in - * consecutive iterations that indicates the algorithm has - * converged - * @param maxIters the maximum number of sequence terms to evaluate before - * giving up - * @param patience the number of consecutive iterations required for the - * tolerance criterion to be satisfied to stop the algorithm - * @param m the parameter in the paper; must be >= 1 - * @param sigma the parameter sigma defining the auxiliary sequence R(l); - * must be > 1 - * @param s the parameter s defining the auxiliary sequence R(l); must - * be >= 1 - */ - public Richardson(final double tolerance, final int maxIters, final int patience, final int m, final double sigma, - final int s) { - super(tolerance, maxIters, patience); - myM = m; - mySigma = sigma; - myIncr = s; - - myElements = new double[maxIters + 1]; - myG = new double[myM]; - myPsiAI = new double[maxIters + 1][2][2]; - myPsiQ = new double[maxIters + 1][myM][2]; - myPsiG = new double[maxIters + 1][myM][2]; + private final double mySigma; + private final int myIncr; + + private int myPrevIndex; + private final double[] myElements; + + private final int myM; + private final double[] myG; + private final double[][][] myPsiAI, myPsiQ, myPsiG; + private int evals; + + /** + * Creates an instance of the d-transform for a given auxiliary sequence. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param m the parameter in the paper; must be >= 1 + * @param sigma the parameter sigma defining the auxiliary sequence R(l); must be > 1 + * @param s the parameter s defining the auxiliary sequence R(l); must be >= 1 + */ + public Richardson(final double tolerance, final int maxIters, final int patience, final int m, + final double sigma, final int s) { + super(tolerance, maxIters, patience); + myM = m; + mySigma = sigma; + myIncr = s; + + myElements = new double[maxIters + 1]; + myG = new double[myM]; + myPsiAI = new double[maxIters + 1][2][2]; + myPsiQ = new double[maxIters + 1][myM][2]; + myPsiG = new double[maxIters + 1][myM][2]; + } + + /** + * Creates an instance of the d-transform for a given auxiliary sequence with sigma = 1. The case + * s = 1 is most suitable for alternating series or quickly-converging series, while s > 1 is + * suitable for power series, trigonometric series, and other series representing analytic + * functions, for arguments close to the functions' singular values. [1] notes that s = 1 is OK + * for arguments far away from the functions' singular values. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param m the parameter in the paper; must be >= 1 + * @param s the parameter s defining the auxiliary sequence R(l); must be >= 1 + */ + public Richardson(final double tolerance, final int maxIters, final int patience, final int m, + final int s) { + this(tolerance, maxIters, patience, m, 1.0, s); + } + + /** + * Creates an instance of the d-transform for a given auxiliary sequence with s = 1. Most suitable + * for monotonic series that converge slowly. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param m the parameter in the paper; must be >= 1 + * @param sigma the parameter sigma defining the auxiliary sequence R(l); must be > 1 + */ + public Richardson(final double tolerance, final int maxIters, final int patience, final int m, + final double sigma) { + this(tolerance, maxIters, patience, m, sigma, 1); + } + + /** + * Creates an instance of the d-transform for the simplest auxiliary sequence R(l) = l. Most + * suitable for alternating series. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param m the parameter in the paper; must be >= 1 + */ + public Richardson(final double tolerance, final int maxIters, final int patience, final int m) { + this(tolerance, maxIters, patience, m, 1.0, 1); + } + + /** + * Creates an instance of the d-transform for the simplest auxiliary sequence R(l) = l. Most + * suitable for alternating series. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + */ + public Richardson(final double tolerance, final int maxIters, final int patience) { + this(tolerance, maxIters, patience, 2); + } + + @Override + public final double next(final double e, final double term) { + + // first m elements do not induce an update of g or psi + myElements[myIndex] = e; + ++myIndex; + if (myIndex < myM) { + return term; } - /** - * Creates an instance of the d-transform for a given auxiliary sequence with - * sigma = 1. The case s = 1 is most suitable for alternating series or - * quickly-converging series, while s > 1 is suitable for power series, - * trigonometric series, and other series representing analytic functions, for - * arguments close to the functions' singular values. [1] notes that s = 1 is OK - * for arguments far away from the functions' singular values. - * - * @param tolerance the smallest acceptable change in series evaluation in - * consecutive iterations that indicates the algorithm has - * converged - * @param maxIters the maximum number of sequence terms to evaluate before - * giving up - * @param patience the number of consecutive iterations required for the - * tolerance criterion to be satisfied to stop the algorithm - * @param m the parameter in the paper; must be >= 1 - * @param s the parameter s defining the auxiliary sequence R(l); must - * be >= 1 - */ - public Richardson(final double tolerance, final int maxIters, final int patience, final int m, final int s) { - this(tolerance, maxIters, patience, m, 1.0, s); + // main updates for psi: in the iterative mode we can only use R(l) = l + final int l = myIndex - myM; + final int Rl = l; + double al = 0.0; + for (int i = 0; i <= Rl; ++i) { + al += myElements[i]; } - - /** - * Creates an instance of the d-transform for a given auxiliary sequence with s - * = 1. Most suitable for monotonic series that converge slowly. - * - * @param tolerance the smallest acceptable change in series evaluation in - * consecutive iterations that indicates the algorithm has - * converged - * @param maxIters the maximum number of sequence terms to evaluate before - * giving up - * @param patience the number of consecutive iterations required for the - * tolerance criterion to be satisfied to stop the algorithm - * @param m the parameter in the paper; must be >= 1 - * @param sigma the parameter sigma defining the auxiliary sequence R(l); - * must be > 1 - */ - public Richardson(final double tolerance, final int maxIters, final int patience, final int m, final double sigma) { - this(tolerance, maxIters, patience, m, sigma, 1); + updateG(Rl); + return updatePsi(al, Rl, l); + } + + @Override + public final SeriesSolution limit(final Iterable seq, final boolean series, + final int extrapolateStart) { + + // the sequence will be consumed one element at time + final Iterator it = seq.iterator(); + myPrevIndex = -1; + myIndex = 0; + evals = 0; + + // read the first elements in the partial sum + double partialSum = 0.0; + for (int i = 0; i < extrapolateStart; ++i) { + partialSum += it.next(); + ++myIndex; + ++evals; } - /** - * Creates an instance of the d-transform for the simplest auxiliary sequence - * R(l) = l. Most suitable for alternating series. - * - * @param tolerance the smallest acceptable change in series evaluation in - * consecutive iterations that indicates the algorithm has - * converged - * @param maxIters the maximum number of sequence terms to evaluate before - * giving up - * @param patience the number of consecutive iterations required for the - * tolerance criterion to be satisfied to stop the algorithm - * @param m the parameter in the paper; must be >= 1 - */ - public Richardson(final double tolerance, final int maxIters, final int patience, final int m) { - this(tolerance, maxIters, patience, m, 1.0, 1); + // extract the first element of the transformed series + int Rl = computeR(0); + double al = computePartialSum(it, Rl); + if (Double.isNaN(al)) { + return new SeriesSolution(Double.NaN, Double.NaN, evals, false); } - - /** - * Creates an instance of the d-transform for the simplest auxiliary sequence - * R(l) = l. Most suitable for alternating series. - * - * @param tolerance the smallest acceptable change in series evaluation in - * consecutive iterations that indicates the algorithm has - * converged - * @param maxIters the maximum number of sequence terms to evaluate before - * giving up - * @param patience the number of consecutive iterations required for the - * tolerance criterion to be satisfied to stop the algorithm - */ - public Richardson(final double tolerance, final int maxIters, final int patience) { - this(tolerance, maxIters, patience, 2); + updateG(Rl); + updatePsi(al, Rl, 0); + + // main loop + double result = al; + double oldResult = result; + int convergeSteps = 0; + for (int l = 1; l <= myMaxIters; ++l) { + + // extract the next element of the transformed series + Rl = computeR(l); + al = computePartialSum(it, Rl); + if (Double.isNaN(al)) { + break; + } + updateG(Rl); + result = updatePsi(al, Rl, l); + + // check for convergence + if (l >= 2) { + final double error = Math.abs(result - oldResult); + if (error <= myTol) { + ++convergeSteps; + } else { + convergeSteps = 0; + } + if (convergeSteps >= myPatience) { + return new SeriesSolution(result + partialSum, error, evals, true); + } + } + oldResult = result; } - @Override - public final double next(final double e, final double term) { - - // first m elements do not induce an update of g or psi - myElements[myIndex] = e; - ++myIndex; - if (myIndex < myM) { - return term; - } - - // main updates for psi: in the iterative mode we can only use R(l) = l - final int l = myIndex - myM; - final int Rl = l; - double al = 0.0; - for (int i = 0; i <= Rl; ++i) { - al += myElements[i]; - } - updateG(Rl); - return updatePsi(al, Rl, l); - } + // could not converge + return new SeriesSolution(Double.NaN, Double.NaN, evals, false); + } - @Override - public final SeriesSolution limit(final Iterable seq, final boolean series, final int extrapolateStart) { - - // the sequence will be consumed one element at time - final Iterator it = seq.iterator(); - myPrevIndex = -1; - myIndex = 0; - evals = 0; - - // read the first elements in the partial sum - double partialSum = 0.0; - for (int i = 0; i < extrapolateStart; ++i) { - partialSum += it.next(); - ++myIndex; - ++evals; - } - - // extract the first element of the transformed series - int Rl = computeR(0); - double al = computePartialSum(it, Rl); - if (Double.isNaN(al)) { - return new SeriesSolution(Double.NaN, Double.NaN, evals, false); - } - updateG(Rl); - updatePsi(al, Rl, 0); - - // main loop - double result = al; - double oldResult = result; - int convergeSteps = 0; - for (int l = 1; l <= myMaxIters; ++l) { - - // extract the next element of the transformed series - Rl = computeR(l); - al = computePartialSum(it, Rl); - if (Double.isNaN(al)) { - break; - } - updateG(Rl); - result = updatePsi(al, Rl, l); - - // check for convergence - if (l >= 2) { - final double error = Math.abs(result - oldResult); - if (error <= myTol) { - ++convergeSteps; - } else { - convergeSteps = 0; - } - if (convergeSteps >= myPatience) { - return new SeriesSolution(result + partialSum, error, evals, true); - } - } - oldResult = result; - } - - // could not converge - return new SeriesSolution(Double.NaN, Double.NaN, evals, false); - } + @Override + public final String getName() { + return "W(" + myM + ") Algorithm"; + } + + private final double computePartialSum(final Iterator it, final int R) { - @Override - public final String getName() { - return "W(" + myM + ") Algorithm"; + // cache the sequence elements so they only need to be consumed once + final int newIndex = R + myM - 1; + if (newIndex >= myMaxIters) { + return Double.NaN; + } + for (int i = 1 + myPrevIndex; i <= newIndex; ++i) { + myElements[i] = it.next(); + ++myIndex; + ++evals; } + myPrevIndex = newIndex; - private final double computePartialSum(final Iterator it, final int R) { - - // cache the sequence elements so they only need to be consumed once - final int newIndex = R + myM - 1; - if (newIndex >= myMaxIters) { - return Double.NaN; - } - for (int i = 1 + myPrevIndex; i <= newIndex; ++i) { - myElements[i] = it.next(); - ++myIndex; - ++evals; - } - myPrevIndex = newIndex; - - // compute the transformed series element a(l) (B.3) - double al = 0.0; - for (int i = 0; i <= R; ++i) { - al += myElements[i]; - } - return al; + // compute the transformed series element a(l) (B.3) + double al = 0.0; + for (int i = 0; i <= R; ++i) { + al += myElements[i]; + } + return al; + } + + private final double updatePsi(final double term, final int R, final int l) { + final int ic = (l + 1) & 1; + final int ip = 1 - ic; + myPsiAI[0][0][ip] = term / myG[0]; + myPsiAI[0][1][ip] = 1.0 / myG[0]; + myPsiQ[0][0][ip] = R + 1.0; + for (int k = 1; k < myM; ++k) { + myPsiG[0][k - 1][ip] = myG[k] / myG[0]; + } + myPsiG[0][myM - 1][ip] = 1.0 / (R + 1.0); + + double sign = -1.0; + for (int p = 1; p <= l; ++p) { + double d = 0.0; + if (p <= myM) { + d = myPsiG[p - 1][p - 1][ip] - myPsiG[p - 1][p - 1][ic]; + for (int i = p + 2; i <= myM + 1; ++i) { + myPsiG[p][i - 2][ip] = (myPsiG[p - 1][i - 2][ip] - myPsiG[p - 1][i - 2][ic]) / d; + } + } + if (p < myM) { + myPsiQ[p][p][ip] = sign / myPsiG[p][myM - 1][ip]; + sign = -sign; + } + for (int q = 1; q <= p - 1 && q <= myM - 1; ++q) { + final double ps = myPsiQ[p - 2][q - 1][ic]; + final double dq = ps / myPsiQ[p - 1][q - 1][ic] - ps / myPsiQ[p - 1][q - 1][ip]; + myPsiQ[p][q][ip] = (myPsiQ[p - 1][q][ip] - myPsiQ[p - 1][q][ic]) / dq; + } + if (p > myM) { + final double ps = myPsiQ[p - 2][myM - 1][ic]; + d = ps / myPsiQ[p - 1][myM - 1][ic] - ps / myPsiQ[p - 1][myM - 1][ip]; + } + myPsiQ[p][0][ip] = (myPsiQ[p - 1][0][ip] - myPsiQ[p - 1][0][ic]) / d; + myPsiAI[p][0][ip] = (myPsiAI[p - 1][0][ip] - myPsiAI[p - 1][0][ic]) / d; + myPsiAI[p][1][ip] = (myPsiAI[p - 1][1][ip] - myPsiAI[p - 1][1][ic]) / d; } - private final double updatePsi(final double term, final int R, final int l) { - final int ic = (l + 1) & 1; - final int ip = 1 - ic; - myPsiAI[0][0][ip] = term / myG[0]; - myPsiAI[0][1][ip] = 1.0 / myG[0]; - myPsiQ[0][0][ip] = R + 1.0; - for (int k = 1; k < myM; ++k) { - myPsiG[0][k - 1][ip] = myG[k] / myG[0]; - } - myPsiG[0][myM - 1][ip] = 1.0 / (R + 1.0); - - double sign = -1.0; - for (int p = 1; p <= l; ++p) { - double d = 0.0; - if (p <= myM) { - d = myPsiG[p - 1][p - 1][ip] - myPsiG[p - 1][p - 1][ic]; - for (int i = p + 2; i <= myM + 1; ++i) { - myPsiG[p][i - 2][ip] = (myPsiG[p - 1][i - 2][ip] - myPsiG[p - 1][i - 2][ic]) / d; - } - } - if (p < myM) { - myPsiQ[p][p][ip] = sign / myPsiG[p][myM - 1][ip]; - sign = -sign; - } - for (int q = 1; q <= p - 1 && q <= myM - 1; ++q) { - final double ps = myPsiQ[p - 2][q - 1][ic]; - final double dq = ps / myPsiQ[p - 1][q - 1][ic] - ps / myPsiQ[p - 1][q - 1][ip]; - myPsiQ[p][q][ip] = (myPsiQ[p - 1][q][ip] - myPsiQ[p - 1][q][ic]) / dq; - } - if (p > myM) { - final double ps = myPsiQ[p - 2][myM - 1][ic]; - d = ps / myPsiQ[p - 1][myM - 1][ic] - ps / myPsiQ[p - 1][myM - 1][ip]; - } - myPsiQ[p][0][ip] = (myPsiQ[p - 1][0][ip] - myPsiQ[p - 1][0][ic]) / d; - myPsiAI[p][0][ip] = (myPsiAI[p - 1][0][ip] - myPsiAI[p - 1][0][ic]) / d; - myPsiAI[p][1][ip] = (myPsiAI[p - 1][1][ip] - myPsiAI[p - 1][1][ic]) / d; - } - - final double denom = myPsiAI[l][1][ip]; - if (Math.abs(denom) < TINY) { - return HUGE; - } else { - return myPsiAI[l][0][ip] / denom; - } + final double denom = myPsiAI[l][1][ip]; + if (Math.abs(denom) < TINY) { + return HUGE; + } else { + return myPsiAI[l][0][ip] / denom; } + } + + private final void updateG(final int R) { - private final void updateG(final int R) { - - // compute the elements of the modified differences g(l) (B.4) - for (int k = 0; k < myM; ++k) { - myG[k] = myElements[R + k]; - } - for (int i = 2; i <= myM; ++i) { - for (int j = myM; j >= i; --j) { - myG[j - 1] -= myG[j - 2]; - } - } - double p = R + 1.0; - for (int k = 0; k < myM; ++k) { - myG[k] *= p; - p *= (R + 1.0); - } + // compute the elements of the modified differences g(l) (B.4) + for (int k = 0; k < myM; ++k) { + myG[k] = myElements[R + k]; + } + for (int i = 2; i <= myM; ++i) { + for (int j = myM; j >= i; --j) { + myG[j - 1] -= myG[j - 2]; + } + } + double p = R + 1.0; + for (int k = 0; k < myM; ++k) { + myG[k] *= p; + p *= (R + 1.0); } + } - private final int computeR(final int l) { - int r = 0; - for (int i = 0; i < l; ++i) { - r = (int) (mySigma * r) + myIncr; - } - return r; + private final int computeR(final int l) { + int r = 0; + for (int i = 0; i < l; ++i) { + r = (int) (mySigma * r) + myIncr; } + return r; + } } diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/SeriesAlgorithm.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/SeriesAlgorithm.java index 7ba862d959..df61c97ab5 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/SeriesAlgorithm.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/SeriesAlgorithm.java @@ -2,6 +2,7 @@ import java.util.Iterator; import java.util.function.Function; +import org.hipparchus.complex.Complex; /** * The base class of all convergence acceleration algorithms. Provides functionality for evaluating @@ -31,6 +32,28 @@ public String toString() { } } + public static final class SeriesSolutionComplex { + + public final Complex limit; + public final double error; + public final int evaluations; + public final boolean converged; + + public SeriesSolutionComplex(final Complex est, final double err, final int evals, + final boolean success) { + limit = est; + error = err; + evaluations = evals; + converged = success; + } + + @Override + public String toString() { + return String.format("%.08f", limit) + " +- " + String.format("%.08f", error) + "\n" + + "evaluations: " + evaluations + "\n" + "converged: " + converged; + } + } + protected static final double HUGE = 1e60; protected static final double TINY = 1e-60; @@ -64,6 +87,10 @@ public SeriesAlgorithm(final double tolerance, final int maxIters, final int pat */ public abstract double next(double e, double term); + public Complex nextComplex(Complex e, Complex term) { + return null; + } + /** * Returns a string representation of the current algorithm. * @@ -143,6 +170,67 @@ public SeriesSolution limit(final Iterable seq, final boolean series, return new SeriesSolution(est + partial, Double.NaN, evals, false); } + // public SeriesSolutionComplex limitComplex(final Iterable seq, final boolean series) { + // return limitComplex(seq, series, 0); + // } + // + // public SeriesSolutionComplex limitComplex(final Iterable seq, final boolean series, + // final int extrapolateStart) { + // myIndex = 0; + // int convergeSteps = 0; + // int indexBeforeExtrap = 0; + // int evals = 0; + // Complex partial = Complex.ZERO; + // Complex term = Complex.ZERO; + // Complex est = Complex.NaN; + // + // for (final Complex e : seq) { + // ++evals; + // + // // check if extrapolation starts + // if (indexBeforeExtrap < extrapolateStart) { + // if (series) { + // partial = partial.add(e); + // } + // ++indexBeforeExtrap; + // continue; + // } + // + // // get the next term of the sequence + // if (series) { + // term = term.add(e); + // } else { + // term = e; + // } + // + // // estimate the next term + // final Complex oldest = est; + // est = nextComplex(e, term); + // if (est.isNaN()) { + // break; + // } + // + // // monitor convergence + // if (myIndex >= 2) { + // final double error = oldest.subtract(est).norm(); + // if (error <= myTol && est.getReal() != HUGE && est.getImaginary() != HUGE) { + // ++convergeSteps; + // } else { + // convergeSteps = 0; + // } + // if (convergeSteps >= myPatience) { + // return new SeriesSolutionComplex(est.add(partial), error, evals, true); + // } + // if (myIndex >= myMaxIters || !Double.isFinite(error)) { + // break; + // } + // } + // } + // + // // did not achieve the desired error + // return new SeriesSolutionComplex(est.add(partial), Double.NaN, evals, false); + // } + /** * Given a sequence represented as an Iterable object, numerically evaluates the limit of the * sequence or the corresponding series. diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/AitkenComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/AitkenComplex.java new file mode 100644 index 0000000000..ad78ac3878 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/AitkenComplex.java @@ -0,0 +1,117 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import org.hipparchus.complex.Complex; + +/** + * Implements a modified adaptive Aitken delta^2 process for estimating infinite series, based on + * the method in [1]. + * + *

+ * References: + *

    + *
  • [1] Osada, Naoki. Acceleration methods for slowly convergent sequences and their + * applications. Diss. PhD thesis, Nagoya University, 1993.
    • + *
+ *

+ */ +public final class AitkenComplex extends SeriesAlgorithmComplex { + + private int status; + private Complex xx, dx, dd, xp, alpha; + private final Complex[] x; + private final Complex[][] s, ds, ts, dts; + + public AitkenComplex(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + x = new Complex[maxIters]; + s = new Complex[2][maxIters + 1]; + ds = new Complex[2][maxIters + 1]; + ts = new Complex[2][maxIters + 1]; + dts = new Complex[2][maxIters + 1]; + for (int i = 0; i < 2; i++) { + for (int j = 0; j < maxIters + 1; j++) { + s[i][j] = Complex.ZERO; + ds[i][j] = Complex.ZERO; + ts[i][j] = Complex.ZERO; + dts[i][j] = Complex.ZERO; + } + } + } + + @Override + public final Complex next(final Complex e, final Complex term) { + if (myIndex == 0) { + xx = dx = dd = Complex.ZERO; + } + final Complex xxp = xx; + x[myIndex] = xx = term; + if (myIndex == 0) { + dx = xx; + } else if (myIndex == 1) { + final Complex dxp = dx; + dx = xx.subtract(xxp); + dd = dx.divide(dx.subtract(dxp)); + } else { + final Complex dxp = dx; + final Complex ddp = dd; + dx = xx.subtract(xxp); + dd = dx.divide(dx.subtract(dxp)); + alpha = dd.subtract(ddp).reciprocal().add(1.0); + alpha = aitken(alpha, ts, dts, myIndex - 1, new Complex(-2.0)); + } + xp = xx; + if (myIndex >= 2) { + for (int m = 1; m <= myIndex + 1; ++m) { + xp = aitken(x[m - 1], s, ds, m, alpha); + } + } + ++myIndex; + if (status == 1) { + return Complex.NaN; + } else { + return xp; + } + } + + private final Complex aitken(final Complex xx, final Complex[][] s, final Complex[][] ds, + final int n, final Complex theta) { + status = 0; + if (n != 1) { + final int kend = (n - 1) >> 1; + System.arraycopy(s[1], 0, s[0], 0, kend + 1); + if ((n & 1) == 0) { + System.arraycopy(ds[1], 0, ds[0], 0, kend); + } else { + System.arraycopy(ds[1], 0, ds[0], 0, kend + 1); + } + } + s[1][0] = xx; + if (n == 1) { + ds[1][0] = xx; + return xx; + } + ds[1][0] = xx.subtract(s[0][0]); + for (int k = 1; k <= (n >> 1); ++k) { + final Complex w1 = ds[0][k - 1].multiply(ds[1][k - 1]); + final Complex w2 = ds[1][k - 1].subtract(ds[0][k - 1]); + if (w2.norm() < TINY) { + status = 1; + return xx; + } + final int twok = k << 1; + final Complex coef = + (new Complex(twok - 1).subtract(theta)).divide(new Complex(twok - 2).subtract(theta)); + s[1][k] = s[0][k - 1].subtract(coef.multiply(w1).divide(w2)); + if (n != twok - 1) { + ds[1][k] = s[1][k].subtract(s[0][k]); + } + } + final int kopt = n >> 1; + return s[1][kopt]; + } + + @Override + public final String getName() { + return "Modified Aitken Complex"; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/BrezinskiThetaComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/BrezinskiThetaComplex.java new file mode 100644 index 0000000000..61c8dbddf6 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/BrezinskiThetaComplex.java @@ -0,0 +1,95 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import org.hipparchus.complex.Complex; + +/** + * Implements Brezinski's Theta algorithm for convergence of infinite series based on [1] and as + * described in [2]. + * + *

+ * References: + *

    + *
  • [1] Brezinski, Claude. Acc�l�ration de la convergence en analyse num�rique. Vol. 584. + * Springer, 2006.
    • + *
  • [2] Weniger, Ernst Joachim. "Nonlinear sequence transformations for the acceleration of + * convergence and the summation of divergent series." arXiv preprint math/0306302 (2003).
    • + *
+ *

+ */ +public final class BrezinskiThetaComplex extends SeriesAlgorithmComplex { + + private final Complex[] myTabA, myTabB; + + public BrezinskiThetaComplex(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + myTabA = new Complex[maxIters]; + myTabB = new Complex[maxIters]; + } + + @Override + public final Complex next(final Complex e, final Complex term) { + + // first entry of table + if (myIndex == 0) { + myTabA[0] = term; + ++myIndex; + return term; + } + + // swapping the A and B arrays + final Complex[] a, b; + if ((myIndex & 1) == 0) { + a = myTabA; + b = myTabB; + } else { + a = myTabB; + b = myTabA; + } + + // the case n >= 1 + final int jmax = ((myIndex << 1) + 1) / 3; + Complex aux1 = a[0]; + Complex aux2 = Complex.ZERO; + a[0] = term; + for (int j = 1; j <= jmax; ++j) { + final Complex aux3 = aux2; + aux2 = aux1; + if (j < jmax) { + aux1 = a[j]; + } + if ((j & 1) == 0) { + final Complex denom = a[j - 1].subtract(b[j - 1].add(b[j - 1])).add(aux2); + + // correct for underflow + if (denom.norm() <= TINY) { + a[j] = new Complex(HUGE, HUGE); + } else { + a[j] = aux3 + .add((b[j - 2].subtract(aux3)).multiply((a[j - 1].subtract(b[j - 1]))).divide(denom)); + } + } else { + final Complex diff = a[j - 1].subtract(b[j - 1]); + + // correct for underflow + if (diff.norm() <= TINY) { + a[j] = new Complex(HUGE, HUGE); + } else { + a[j] = aux3.add(diff.reciprocal()); + } + } + } + ++myIndex; + + // return result + if ((jmax & 1) == 0) { + return a[jmax]; + } else { + return a[jmax - 1]; + } + } + + @Override + public final String getName() { + return "Brezinski Theta"; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/CohenComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/CohenComplex.java new file mode 100644 index 0000000000..85fab20d63 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/CohenComplex.java @@ -0,0 +1,68 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import org.hipparchus.complex.Complex; +import org.matheclipse.core.numerics.utils.Constants; +import org.matheclipse.core.numerics.utils.SimpleMath; + +/** + * This is an implementation of Algorithm 1 in [1] for evaluating infinite series whose terms have + * alternating signs. + * + *

+ * References: + *

    + *
  • [1] Cohen, Henri, Fernando Rodriguez Villegas, and Don Zagier. "Convergence acceleration of + * alternating series." Experimental mathematics 9.1 (2000): 3-12.
    • + *
+ *

+ */ +public final class CohenComplex extends SeriesAlgorithmComplex { + + private final Complex[] myTab; + private Complex mySign0; + + public CohenComplex(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + myTab = new Complex[maxIters]; + } + + @Override + public final Complex next(final Complex e, final Complex term) { + + // add next element + final int n = myIndex + 1; + myTab[myIndex] = e.abs(); + + // record the sign of the first term, since this method + // requires the first term of the sequence to be positive + if (myIndex == 0) { + mySign0 = e.sign(); + if (mySign0.isZero()) { + mySign0 = Complex.ONE; + } + } + + // initialize d value + double d = SimpleMath.pow(3.0 + 2.0 * Constants.SQRT2, n); + d = 0.5 * (d + 1.0 / d); + + // apply the Chebychef polynomial recursively (Algorithm 1) + double b = -1.0; + double c = -d; + Complex s = Complex.ZERO; + for (int k = 0; k < n; ++k) { + c = b - c; + s = s.add(myTab[k].multiply(c)); + final double numer = (k + n) * (k - n); + final double denom = (k + 0.5) * (k + 1); + b *= numer / denom; + } + ++myIndex; + return mySign0.multiply(s).divide(d); + } + + @Override + public final String getName() { + return "Cohen"; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/DirectComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/DirectComplex.java new file mode 100644 index 0000000000..3d7cf012dd --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/DirectComplex.java @@ -0,0 +1,26 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import org.hipparchus.complex.Complex; + +/** + * Naively computes the limit of a sequence or series by inspecting elements one a time. No + * transform of the original sequence is performed. + */ +public final class DirectComplex extends SeriesAlgorithmComplex { + + public DirectComplex(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + } + + + @Override + public final Complex next(final Complex e, final Complex term) { + ++myIndex; + return term; + } + + @Override + public final String getName() { + return "Direct Complex Sum"; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/EnsembleComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/EnsembleComplex.java new file mode 100644 index 0000000000..705ac83e50 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/EnsembleComplex.java @@ -0,0 +1,279 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import java.util.ArrayList; +import java.util.List; +import org.hipparchus.complex.Complex; +import org.matheclipse.core.numerics.series.dp.complex.LevinComplex.RemainderSequence; +import org.matheclipse.core.numerics.series.dp.complex.WynnEpsilonComplex.ShanksMethod; + +/** + * An ensemble algorithm for evaluating the limits of sequences and series of complex values. This + * algorithms works by initializing a number of convergence acceleration algorithms with different + * properties, and running them in parallel. The estimates of the sequence's or corresponding series + * limit are extracted from the instance that has converged first to a stable limit to the desired + * tolerance. + */ +public final class EnsembleComplex extends SeriesAlgorithmComplex { + + private static final int PRINT_DIGITS = 18; + + private final int myPrint; + private final List myMethods; + private final List myForAlternating; + private final List myForSequence; + + private int myMethodCount, myBestMethod; + private boolean myIsAlternating, myCheckForSequenceOnly; + private double myOldSignum; + private boolean[] myExcept; + private Complex[] myPrevEst, myEst; + + /** + * Creates a new instance of an ensemble convergence acceleration algorithm. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of iterations used to assess whether or not a particular instance of + * convergence acceleration algorithm used by this algorithm has converged + * @param printProgress how often to print the progress of the current algorithm, e.g. the + * estimate of a sequence's or series's limit at each iteration according to each algorithm + * instance + */ + public EnsembleComplex(final double tolerance, final int maxIters, final int patience, + final int printProgress) { + super(tolerance, maxIters, patience); + myPrint = printProgress; + myMethods = new ArrayList<>(); + myForAlternating = new ArrayList<>(); + myForSequence = new ArrayList<>(); + + // only for alternating series + myMethods.add(new CohenComplex(myTol, myMaxIters, 1)); + myForAlternating.add(true); + myForSequence.add(false); + myMethods.add(new LevinComplex(myTol, myMaxIters, 1, RemainderSequence.D)); + myForAlternating.add(true); + myForSequence.add(false); + + // for alternating and linear series + myMethods.add(new AitkenComplex(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new WynnEpsilonComplex(myTol, myMaxIters, 1, ShanksMethod.WYNN)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new LevinComplex(myTol, myMaxIters, 1, RemainderSequence.T)); + myForAlternating.add(false); + myForSequence.add(false); + + // for logarithmic series + myMethods.add(new LevinComplex(myTol, myMaxIters, 1, RemainderSequence.U)); + myForAlternating.add(false); + myForSequence.add(false); + myMethods.add(new LevinComplex(myTol, myMaxIters, 1, RemainderSequence.V)); + myForAlternating.add(false); + myForSequence.add(false); + myMethods.add(new RichardsonComplex(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(false); + myMethods.add(new BrezinskiThetaComplex(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new IteratedThetaComplex(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new WynnRhoComplex(myTol, myMaxIters, 1)); + myForAlternating.add(false); + myForSequence.add(true); + myMethods.add(new WynnEpsilonComplex(myTol, myMaxIters, 1, ShanksMethod.ALTERNATING)); + myForAlternating.add(false); + myForSequence.add(true); + } + + public EnsembleComplex(final double tolerance, final int maxIters, final int patience) { + this(tolerance, maxIters, patience, 0); + } + + @Override + public final Complex next(final Complex e, final Complex term) { + + // ******************************************************************** + // INITIALIZATION + // ******************************************************************** + if (myIndex == 0) { + + // initialize variables + myIsAlternating = true; + myMethodCount = myMethods.size(); + myOldSignum = Math.signum(e.getReal()); + myPrevEst = new Complex[myMethodCount]; + myEst = new Complex[myMethodCount]; + myExcept = new boolean[myMethodCount]; + + // disable methods for series only + if (myCheckForSequenceOnly) { + for (int m = 0; m < myMethodCount; ++m) { + if (!myForSequence.get(m)) { + myExcept[m] = true; + if (myPrint > 0) { + System.out.println("Disabling method " + myMethods.get(m).getName() + + " that is invalid for sequences"); + } + } + } + myCheckForSequenceOnly = false; + } + + // initialize each method + for (int m = 0; m < myMethodCount; ++m) { + if (!myExcept[m]) { + final SeriesAlgorithmComplex method = myMethods.get(m); + method.myIndex = 0; + myEst[m] = method.next(e, term); + } + } + + // print header + if (myPrint > 0) { + String line = ""; + line += pad("Index", PRINT_DIGITS + 5); + for (final SeriesAlgorithmComplex method : myMethods) { + line += "\t" + pad(method.getName(), PRINT_DIGITS + 5); + } + System.out.println(line); + } + + // print values + if (myPrint > 0) { + printRow(); + } + + ++myIndex; + return term; + } + + // ******************************************************************** + // SEQUENCE INSPECTION + // ******************************************************************** + // test for invalid value of series or sequence + if (!term.isFinite() || !e.isFinite()) { + if (myPrint > 0) { + System.out.println("Aborting due to NaN at iteration " + myIndex); + } + ++myIndex; + return Complex.NaN; + } + + // track the sign of the sequence + final double signum = Math.signum(e.getReal()); + if (myIndex > 0 && myIsAlternating && signum * myOldSignum >= 0) { + myIsAlternating = false; + } + myOldSignum = signum; + + // exclude alternating series methods if the series does not alternate + if (!myIsAlternating) { + for (int m = 0; m < myMethodCount; ++m) { + if (!myExcept[m] && myForAlternating.get(m)) { + myExcept[m] = true; + if (myPrint > 0) { + System.out.println("Disabling method " + myMethods.get(m).getName() + + " because series is not alternating at iteration " + myIndex); + } + } + } + } + + // ******************************************************************** + // MAIN UPDATE + // ******************************************************************** + myBestMethod = -1; + double bestError = Double.POSITIVE_INFINITY; + for (int m = 0; m < myMethodCount; ++m) { + if (myExcept[m]) { + continue; + } + + // estimate the next terms + final SeriesAlgorithmComplex method = myMethods.get(m); + myPrevEst[m] = myEst[m]; + myEst[m] = method.next(e, term); + + // if a method produces an invalid estimate it is excluded + if (!myEst[m].isFinite()) { + if (myPrint > 0) { + System.out.println("Disabling method " + method.getName() + + " due to instability at iteration " + myIndex); + } + myExcept[m] = true; + continue; + } + + // estimate error and best method + final double error = myPrevEst[m].subtract(myEst[m]).norm(); // Math.abs(myPrevEst[m] - + // myEst[m]); + if (error < bestError && myEst[m].getReal() != HUGE && myEst[m].getImaginary() != HUGE) { + myBestMethod = m; + bestError = error; + } + } + + // ******************************************************************** + // FINALIZE THIS ITERATION + // ******************************************************************** + if (myPrint > 0 && myIndex % myPrint == 0) { + printRow(); + } + + // return the method with the best convergence so far + ++myIndex; + if (myBestMethod >= 0) { + return myEst[myBestMethod]; + } else { + return Complex.NaN; + } + } + + @Override + public final SeriesSolutionComplex limit(final Iterable seq, final boolean series, + final int extrapolateStart) { + myCheckForSequenceOnly = !series; + final SeriesSolutionComplex result = super.limit(seq, series, extrapolateStart); + if (myPrint > 0 && myBestMethod >= 0 && result.converged) { + System.out.println("Converged at iteration " + myIndex + " with method " + + myMethods.get(myBestMethod).getName()); + } + return result; + } + + @Override + public final String getName() { + return "Ensemble"; + } + + private final void printRow() { + String line = pad(myIndex + "", PRINT_DIGITS); + for (int m = 0; m < myMethodCount; ++m) { + final String str; + if (myExcept[m]) { + str = "-"; + } else { + str = myEst[m].toString(); + } + line += "\t" + pad(str, PRINT_DIGITS); + } + System.out.println(line); + } + + private static final String pad(final String str, final int len) { + if (str.length() >= len) { + return str; + } + String result = str; + for (int i = str.length() + 1; i <= len; ++i) { + result += " "; + } + return result; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/IteratedThetaComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/IteratedThetaComplex.java new file mode 100644 index 0000000000..7f2ae8ac59 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/IteratedThetaComplex.java @@ -0,0 +1,65 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import org.hipparchus.complex.Complex; + +/** + * Implements the Iterated Theta algorithm for convergence of infinite series as described in [1]. + * + *

+ * References: + *

    + *
  • [1] Weniger, Ernst Joachim. "Nonlinear sequence transformations for the acceleration of + * convergence and the summation of divergent series." arXiv preprint math/0306302 (2003).
    • + *
+ *

+ */ +public final class IteratedThetaComplex extends SeriesAlgorithmComplex { + + private final Complex[] myComplexTab; + + public IteratedThetaComplex(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + myComplexTab = new Complex[maxIters]; + } + + @Override + public final Complex next(final Complex e, final Complex term) { + + // initial value for J + myComplexTab[myIndex] = term; + if (myIndex < 3) { + ++myIndex; + return term; + } + + // higher-order J array + final int lmax = myIndex / 3; + int m = myIndex; + for (int l = 1; l <= lmax; ++l) { + m -= 3; + final Complex diff0 = myComplexTab[m + 1].subtract(myComplexTab[m]); + final Complex diff1 = myComplexTab[m + 2].subtract(myComplexTab[m + 1]); + final Complex diff2 = myComplexTab[m + 3].subtract(myComplexTab[m + 2]); + final Complex dfsq1 = diff1.subtract(diff0); + final Complex dfsq2 = diff2.subtract(diff1); + final Complex denom = diff2.multiply(dfsq1).subtract(diff0.multiply(dfsq2)); + final Complex ratio = dfsq2.divide(denom); + + // safeguarded against underflow + if (denom.norm() < TINY) { + myComplexTab[m] = new Complex(HUGE, HUGE); + } else { + myComplexTab[m] = myComplexTab[m + 1].subtract(diff0.multiply(diff1).multiply(ratio)); + } + } + ++myIndex; + + // return result + return myComplexTab[myIndex % 3]; + } + + @Override + public final String getName() { + return "Iterated Theta Complex"; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/LevinComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/LevinComplex.java new file mode 100644 index 0000000000..6071ceb074 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/LevinComplex.java @@ -0,0 +1,164 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import org.hipparchus.complex.Complex; +import org.matheclipse.core.numerics.utils.SimpleMath; + +/** + * Implements the family of Generalized Levin algorithms for convergence of infinite series based on + * remainder sequences described in [1, 2] and as outlined in [3]. + * + *

+ * References: + *

    + *
  • [1] Levin, David. "Development of non-linear transformations for improving convergence of + * sequences." International Journal of Computer Mathematics 3.1-4 (1972): 371-388.
    • + *
  • [2] Smith, David A., and William F. Ford. "Acceleration of linear and logarithmic + * convergence." SIAM Journal on Numerical Analysis 16.2 (1979): 223-240.
    • + *
  • [3] Weniger, Ernst Joachim. "Nonlinear sequence transformations for the acceleration of + * convergence and the summation of divergent series." arXiv preprint math/0306302 (2003).
    • + *
+ *

+ */ +public final class LevinComplex extends SeriesAlgorithmComplex { + + /** + * The type of remainder sequence to assume when accelerating the convergence of a series. + */ + public static enum RemainderSequence { + T, D, U, V + } + + private final String prefix; + private final double myBeta; + private final RemainderSequence myRemainder; + + private Complex myPrevE, myPrevTerm; + private final Complex[] myTabHi, myTabLo; + + /** + * Creates a new instance of the Generalized Levin algorithm with customizable remainder + * sequences. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param beta the constant described in [1, 3] for the Levin u-transform remainder sequence: if + * this sequence is not selected this argument is ignored + * @param remainder the type of assumed remainder sequence as described in [3] + */ + public LevinComplex(final double tolerance, final int maxIters, final int patience, + final double beta, final RemainderSequence remainder) { + super(tolerance, maxIters, patience); + myBeta = beta; + myRemainder = remainder; + prefix = remainder.name(); + + myTabHi = new Complex[maxIters]; + myTabLo = new Complex[maxIters]; + } + + /** + * Creates a new instance of the Generalized Levin algorithm with customizable remainder + * sequences. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param remainder the type of assumed remainder sequence as described in [3] + */ + public LevinComplex(final double tolerance, final int maxIters, final int patience, + final RemainderSequence remainder) { + this(tolerance, maxIters, patience, 1.0, remainder); + } + + public LevinComplex(final double tolerance, final int maxIters, final int patience) { + this(tolerance, maxIters, patience, RemainderSequence.U); + } + + @Override + public final Complex next(final Complex e, final Complex term) { + + // shift over elements by one by skipping first iteration so that + // the iteration is actually computing the term one over from the + // one used in the current computation + if (myIndex == 0) { + myPrevE = e; + myPrevTerm = term; + if (myRemainder == RemainderSequence.V || myRemainder == RemainderSequence.D) { + ++myIndex; + return myPrevTerm; + } + } + + // compute the remainder term + final int k; + final Complex rem; + final Complex currentTerm; + switch (myRemainder) { + case T: + k = myIndex; + rem = e; + currentTerm = term; + break; + case U: + k = myIndex; + rem = e.multiply(myBeta + k); + currentTerm = term; + break; + case V: + k = myIndex - 1; + rem = myPrevE.multiply(e).divide(myPrevE.subtract(e)); + currentTerm = myPrevTerm; + break; + case D: + k = myIndex - 1; + rem = e; + currentTerm = myPrevTerm; + break; + default: + return Complex.NaN; + } + if (rem.norm() < TINY) { + return Complex.NaN; + } + + // update table + myTabHi[k] = currentTerm.divide(rem); + myTabLo[k] = rem.reciprocal(); + if (k > 0) { + myTabHi[k - 1] = myTabHi[k].subtract(myTabHi[k - 1]); + myTabLo[k - 1] = myTabLo[k].subtract(myTabLo[k - 1]); + if (k > 1) { + final double base = (myBeta + k - 1) / (myBeta + k); + for (int j = 2; j <= k; ++j) { + final int kmj = k - j; + final double cjm2 = SimpleMath.pow(base, j - 2); + final double fac = (myBeta + kmj) * cjm2 / (myBeta + k); + myTabHi[kmj] = myTabHi[kmj + 1].subtract(myTabHi[kmj].multiply(fac)); + myTabLo[kmj] = myTabLo[kmj + 1].subtract(myTabLo[kmj].multiply(fac)); + } + } + } + + // increment counters and temporaries + myPrevE = e; + myPrevTerm = term; + ++myIndex; + + // correct for underflow + if (myTabLo[0].norm() < TINY) { + return new Complex(HUGE, HUGE); + } else { + return myTabHi[0].divide(myTabLo[0]); + } + } + + @Override + public final String getName() { + return "Levin" + " " + prefix; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/RichardsonComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/RichardsonComplex.java new file mode 100644 index 0000000000..e773f4f25f --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/RichardsonComplex.java @@ -0,0 +1,315 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import java.util.Iterator; +import org.hipparchus.complex.Complex; + +/** + * Implements a generalized form of the Richardson extrapolation for evaluating infinite series. + * This is implemented using a d-transform of the original series and is called the W(m) algorithm + * [1]. + * + *

+ * The method in [1] defines an auxiliary increasing sequence of integers, R(l), defined by the + * recursion R(l + 1) = (int)(sigma * R(l)) + s, where sigma and s are two tuning parameters. The + * series are evaluated at each R(l) of the auxiliary sequence and forms a subsequence of the + * original partial sums that converges to the same limit. The d-transform is then applied to + * accelerate this sequence to infinity. + *

+ * + *

+ * References: + *

    + *
  • [1] Ford, William & Sidi, Avram. (1987). An Algorithm for a Generalization of the Richardson + * Extrapolation Process. SIAM Journal on Numerical Analysis. 24. 10.1137/0724080.
    • + *
+ *

+ */ +public final class RichardsonComplex extends SeriesAlgorithmComplex { + + private final double mySigma; + private final int myIncr; + + private int myPrevIndex; + private final Complex[] myElements; + + private final int myM; + private final Complex[] myG; + private final Complex[][][] myPsiAI, myPsiQ, myPsiG; + private int evals; + + /** + * Creates an instance of the d-transform for a given auxiliary sequence. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param m the parameter in the paper; must be >= 1 + * @param sigma the parameter sigma defining the auxiliary sequence R(l); must be > 1 + * @param s the parameter s defining the auxiliary sequence R(l); must be >= 1 + */ + public RichardsonComplex(final double tolerance, final int maxIters, final int patience, + final int m, final double sigma, final int s) { + super(tolerance, maxIters, patience); + myM = m; + mySigma = sigma; + myIncr = s; + + myElements = new Complex[maxIters + 1]; + myG = new Complex[myM]; + myPsiAI = new Complex[maxIters + 1][2][2]; + myPsiQ = new Complex[maxIters + 1][myM][2]; + myPsiG = new Complex[maxIters + 1][myM][2]; + } + + /** + * Creates an instance of the d-transform for a given auxiliary sequence with sigma = 1. The case + * s = 1 is most suitable for alternating series or quickly-converging series, while s > 1 is + * suitable for power series, trigonometric series, and other series representing analytic + * functions, for arguments close to the functions' singular values. [1] notes that s = 1 is OK + * for arguments far away from the functions' singular values. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param m the parameter in the paper; must be >= 1 + * @param s the parameter s defining the auxiliary sequence R(l); must be >= 1 + */ + public RichardsonComplex(final double tolerance, final int maxIters, final int patience, + final int m, final int s) { + this(tolerance, maxIters, patience, m, 1.0, s); + } + + /** + * Creates an instance of the d-transform for a given auxiliary sequence with s = 1. Most suitable + * for monotonic series that converge slowly. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param m the parameter in the paper; must be >= 1 + * @param sigma the parameter sigma defining the auxiliary sequence R(l); must be > 1 + */ + public RichardsonComplex(final double tolerance, final int maxIters, final int patience, + final int m, final double sigma) { + this(tolerance, maxIters, patience, m, sigma, 1); + } + + /** + * Creates an instance of the d-transform for the simplest auxiliary sequence R(l) = l. Most + * suitable for alternating series. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param m the parameter in the paper; must be >= 1 + */ + public RichardsonComplex(final double tolerance, final int maxIters, final int patience, + final int m) { + this(tolerance, maxIters, patience, m, 1.0, 1); + } + + /** + * Creates an instance of the d-transform for the simplest auxiliary sequence R(l) = l. Most + * suitable for alternating series. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + */ + public RichardsonComplex(final double tolerance, final int maxIters, final int patience) { + this(tolerance, maxIters, patience, 2); + } + + @Override + public final Complex next(final Complex e, final Complex term) { + + // first m elements do not induce an update of g or psi + myElements[myIndex] = e; + ++myIndex; + if (myIndex < myM) { + return term; + } + + // main updates for psi: in the iterative mode we can only use R(l) = l + final int l = myIndex - myM; + final int Rl = l; + Complex al = Complex.ZERO; + for (int i = 0; i <= Rl; ++i) { + al = al.add(myElements[i]); + } + updateG(Rl); + return updatePsi(al, Rl, l); + } + + @Override + public final SeriesSolutionComplex limit(final Iterable seq, final boolean series, + final int extrapolateStart) { + + // the sequence will be consumed one element at time + final Iterator it = seq.iterator(); + myPrevIndex = -1; + myIndex = 0; + evals = 0; + + // read the first elements in the partial sum + Complex partialSum = Complex.ZERO; + for (int i = 0; i < extrapolateStart; ++i) { + partialSum = partialSum.add(it.next()); + ++myIndex; + ++evals; + } + + // extract the first element of the transformed series + int Rl = computeR(0); + Complex al = computePartialSum(it, Rl); + if (al.isNaN()) { + return new SeriesSolutionComplex(Complex.NaN, Double.NaN, evals, false); + } + updateG(Rl); + updatePsi(al, Rl, 0); + + // main loop + Complex result = al; + Complex oldResult = result; + int convergeSteps = 0; + for (int l = 1; l <= myMaxIters; ++l) { + + // extract the next element of the transformed series + Rl = computeR(l); + al = computePartialSum(it, Rl); + if (al.isNaN()) { + break; + } + updateG(Rl); + result = updatePsi(al, Rl, l); + + // check for convergence + if (l >= 2) { + final double error = result.subtract(oldResult).norm();// Math.abs(result - oldResult); + if (error <= myTol) { + ++convergeSteps; + } else { + convergeSteps = 0; + } + if (convergeSteps >= myPatience) { + return new SeriesSolutionComplex(result.add(partialSum), error, evals, true); + } + } + oldResult = result; + } + + // could not converge + return new SeriesSolutionComplex(Complex.NaN, Double.NaN, evals, false); + } + + @Override + public final String getName() { + return "W(" + myM + ") Algorithm"; + } + + private final Complex computePartialSum(final Iterator it, final int R) { + + // cache the sequence elements so they only need to be consumed once + final int newIndex = R + myM - 1; + if (newIndex >= myMaxIters) { + return Complex.NaN; + } + for (int i = 1 + myPrevIndex; i <= newIndex; ++i) { + myElements[i] = it.next(); + ++myIndex; + ++evals; + } + myPrevIndex = newIndex; + + // compute the transformed series element a(l) (B.3) + Complex al = Complex.ZERO; + for (int i = 0; i <= R; ++i) { + al = al.add(myElements[i]); + } + return al; + } + + private final Complex updatePsi(final Complex term, final int R, final int l) { + final int ic = (l + 1) & 1; + final int ip = 1 - ic; + myPsiAI[0][0][ip] = term.divide(myG[0]); + myPsiAI[0][1][ip] = myG[0].reciprocal(); + myPsiQ[0][0][ip] = new Complex(R + 1.0); + for (int k = 1; k < myM; ++k) { + myPsiG[0][k - 1][ip] = myG[k].divide(myG[0]); + } + myPsiG[0][myM - 1][ip] = new Complex(R + 1.0).reciprocal(); + + double sign = -1.0; + for (int p = 1; p <= l; ++p) { + Complex d = Complex.ZERO; + if (p <= myM) { + d = myPsiG[p - 1][p - 1][ip].subtract(myPsiG[p - 1][p - 1][ic]); + for (int i = p + 2; i <= myM + 1; ++i) { + myPsiG[p][i - 2][ip] = + (myPsiG[p - 1][i - 2][ip].subtract(myPsiG[p - 1][i - 2][ic]).divide(d)); + } + } + if (p < myM) { + myPsiQ[p][p][ip] = new Complex(sign).divide(myPsiG[p][myM - 1][ip]); + sign = -sign; + } + for (int q = 1; q <= p - 1 && q <= myM - 1; ++q) { + final Complex ps = myPsiQ[p - 2][q - 1][ic]; + final Complex dq = + ps.divide(myPsiQ[p - 1][q - 1][ic]).subtract(ps.divide(myPsiQ[p - 1][q - 1][ip])); + myPsiQ[p][q][ip] = myPsiQ[p - 1][q][ip].subtract(myPsiQ[p - 1][q][ic]).divide(dq); + } + if (p > myM) { + final Complex ps = myPsiQ[p - 2][myM - 1][ic]; + d = ps.divide(myPsiQ[p - 1][myM - 1][ic]).subtract(ps.divide(myPsiQ[p - 1][myM - 1][ip])); + } + myPsiQ[p][0][ip] = myPsiQ[p - 1][0][ip].subtract(myPsiQ[p - 1][0][ic]).divide(d); + myPsiAI[p][0][ip] = myPsiAI[p - 1][0][ip].subtract(myPsiAI[p - 1][0][ic]).divide(d); + myPsiAI[p][1][ip] = myPsiAI[p - 1][1][ip].subtract(myPsiAI[p - 1][1][ic]).divide(d); + } + + final Complex denom = myPsiAI[l][1][ip]; + if (denom.norm() < TINY) { + return new Complex(HUGE, HUGE); + } else { + return myPsiAI[l][0][ip].divide(denom); + } + } + + private final void updateG(final int R) { + + // compute the elements of the modified differences g(l) (B.4) + for (int k = 0; k < myM; ++k) { + myG[k] = myElements[R + k]; + } + for (int i = 2; i <= myM; ++i) { + for (int j = myM; j >= i; --j) { + myG[j - 1] = myG[j - 1].subtract(myG[j - 2]); + } + } + double p = R + 1.0; + for (int k = 0; k < myM; ++k) { + myG[k] = myG[k].multiply(p); + p *= (R + 1.0); + } + } + + private final int computeR(final int l) { + int r = 0; + for (int i = 0; i < l; ++i) { + r = (int) (mySigma * r) + myIncr; + } + return r; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/SeriesAlgorithmComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/SeriesAlgorithmComplex.java new file mode 100644 index 0000000000..7c69bb4010 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/SeriesAlgorithmComplex.java @@ -0,0 +1,256 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import java.util.Iterator; +import java.util.function.Function; +import org.hipparchus.complex.Complex; + +/** + * The base class of all convergence acceleration algorithms. Provides functionality for evaluating + * infinite series and limits of sequences of complex values. + */ +public abstract class SeriesAlgorithmComplex { + + + public static final class SeriesSolutionComplex { + + public final Complex limit; + public final double error; + public final int evaluations; + public final boolean converged; + + public SeriesSolutionComplex(final Complex est, final double err, final int evals, + final boolean success) { + limit = est; + error = err; + evaluations = evals; + converged = success; + } + + @Override + public String toString() { + return String.format("%.08f", limit) + " +- " + String.format("%.08f", error) + "\n" + + "evaluations: " + evaluations + "\n" + "converged: " + converged; + } + } + + protected static final double HUGE = 1e60; + protected static final double TINY = 1e-60; + + protected final double myTol; + protected final int myMaxIters; + protected final int myPatience; + + protected int myIndex; + + /** + * Creates a new instance of the current convergence acceleration algorithm. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param maxIters the maximum number of sequence terms to evaluate before giving up + */ + public SeriesAlgorithmComplex(final double tolerance, final int maxIters, final int patience) { + myTol = tolerance; + myMaxIters = maxIters; + myPatience = patience; + } + + /** + * Update the current algorithm after observing the specified element and partial sum. + * + * @param e the next element of the sequence + * @param term the partial sum of all observed elements of the sequence thus far + * @return the next estimate of the limit of the sequence + */ + public abstract Complex next(Complex e, Complex term); + + /** + * Returns a string representation of the current algorithm. + * + * @return a string representation of the current algorithm + */ + public abstract String getName(); + + // ========================================================================== + // LIMITS + // ========================================================================== + /** + * Given a sequence represented as an Iterable object, numerically evaluates the limit of the + * sequence or the corresponding series. + * + * @param seq an Iterable object of type Complex whose limit to evaluate + * @param series a boolean variable indicating whether to evaluate the limit of the series or the + * sequence + * @param extrapolateStart the number of terms to observe (and sum trivially) before starting + * extrapolation + * @return a Double that approximates the limit of the sequence or corresponding series. If the + * limit cannot be determined, returns NaN + */ + public SeriesSolutionComplex limit(final Iterable seq, final boolean series) { + return limit(seq, series, 0); + } + + public SeriesSolutionComplex limit(final Iterable seq, final boolean series, + final int extrapolateStart) { + myIndex = 0; + int convergeSteps = 0; + int indexBeforeExtrap = 0; + int evals = 0; + Complex partial = Complex.ZERO; + Complex term = Complex.ZERO; + Complex est = Complex.NaN; + + for (final Complex e : seq) { + ++evals; + + // check if extrapolation starts + if (indexBeforeExtrap < extrapolateStart) { + if (series) { + partial = partial.add(e); + } + ++indexBeforeExtrap; + continue; + } + + // get the next term of the sequence + if (series) { + term = term.add(e); + } else { + term = e; + } + + // estimate the next term + final Complex oldest = est; + est = next(e, term); + if (est.isNaN()) { + break; + } + + // monitor convergence + if (myIndex >= 2) { + final double error = oldest.subtract(est).norm(); + if (error <= myTol && est.getReal() != HUGE && est.getImaginary() != HUGE) { + ++convergeSteps; + } else { + convergeSteps = 0; + } + if (convergeSteps >= myPatience) { + return new SeriesSolutionComplex(est.add(partial), error, evals, true); + } + if (myIndex >= myMaxIters || !Double.isFinite(error)) { + break; + } + } + } + + // did not achieve the desired error + return new SeriesSolutionComplex(est.add(partial), Double.NaN, evals, false); + } + + /** + * Accelerates a sequence or series according to the transformation of the current instance, and + * returns the transformed series. + * + * @param seq an Iterable object of type Double to transform + * @param series a boolean variable indicating whether to evaluate the limit of the series or the + * sequence + * @return an Iterable of type Double representing the transformed series or sequence + */ + public Iterable transform(final Iterable seq, final boolean series) { + return () -> new Iterator<>() { + + private final Iterator it = seq.iterator(); + + private Complex term = Complex.ZERO; + private int index = 0; + + @Override + public final boolean hasNext() { + return index < SeriesAlgorithmComplex.this.myMaxIters && it.hasNext(); + } + + @Override + public final Complex next() { + final Complex e = it.next(); + if (series) { + term = term.add(e); + } else { + term = e; + } + ++index; + return SeriesAlgorithmComplex.this.next(e, term); + } + }; + } + + // ========================================================================== + // TRANSLATION OR REPRESENTATION OF SEQUENCES/SERIES + // ========================================================================== + /** + * The Van Wijngaarden transformation converts a convergent sequence of positive terms into a + * sequence of alternating terms that has the same limit. This trick is incredibly useful for + * accelerating convergence of slowly-increasing sequences of positive terms, because alternating + * series are generally easier to accelerate. + * + *

+ * Formally, given a sequence of terms a(1), a(2)... converts this into an alternating sequence + * where the kth term is + *

+ * a(k) + 2a(2k) + 4a(4k) + ... + *

+ * This method uses the current algorithm's implementation to evaluate the limit of the above + * expression. The resulting series are evalated lazily and represented internally as a function. + *

+ * + *

+ * References: + *

    + *
  • [1] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. + * 2007. Numerical Recipes 3rd Edition: The Art of Scientific Computing (3rd. ed.). Cambridge + * University Press, USA.
    • + *
+ *

+ * + * @param seq a sequence of non-negative terms represented as a Function + * @param start a long representing the starting index at which the sequence begins + * @return a Function representing the resulting alternating sequence obtained with Van + * Wijngaarden's transformation + */ + public Function toAlternatingSeries(final Function func, + final long start) { + return (k) -> { + + // create the condensation sequence with kth term 2^k a(2^k) + final Iterable condensed = () -> new Iterator<>() { + + private long i = k; + private double coeff = 1.0; + + @Override + public final boolean hasNext() { + return i <= (Long.MAX_VALUE >> 1L) && coeff <= (Double.MAX_VALUE / 2.0); + } + + @Override + public final Complex next() { + final Complex term = func.apply(i + start - 1); + final Complex result = term.multiply(coeff); + i <<= 1L; + coeff *= 2.0; + return result; + } + }; + + // determine the next term as the infinite series of this sequence + final SeriesSolutionComplex sol = limit(condensed, true); + final Complex term = sol.limit; + if (((k - 1) & 1L) == 0) { + return term; + } else { + return term.negate(); + } + }; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/WynnEpsilonComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/WynnEpsilonComplex.java new file mode 100644 index 0000000000..1ef25eb04d --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/WynnEpsilonComplex.java @@ -0,0 +1,288 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import org.hipparchus.complex.Complex; + +/** + * Implements the Wynn Epsilon algorithm and its variants for convergence of infinite series. + * Fundamentally, this implementation contains several methods for implementing the Shanks transform + * efficiently. There are several implementations: + *
    + *
  1. the classic moving lozenge technique due to Wynn [1, 2]
    2. + *
  2. the Shanks transform, implemented using the G-transform with a certain remainder sequence + * [3]
    3. + *
  3. the one-parameter generalized epsilon transform of [4]
    4. + *
  4. the alternating epsilon transform of [5]
    5. + *
+ * + *

+ * References: + *

    + *
  • [1] Wynn, Peter. "On a device for computing the e m (S n) transformation." Mathematical + * Tables and Other Aids to Computation (1956): 91-96.
    • + *
  • [2] Weniger, Ernst Joachim. "Nonlinear sequence transformations for the acceleration of + * convergence and the summation of divergent series." arXiv preprint math/0306302 (2003).
    • + *
  • [3] Varga, R.S. Extrapolation methods: theory and practice. Numer Algor 4, 305 (1993). + * https://doi.org/10.1007/BF02144109
    • + *
  • [4] Broeck, Jean-Marc Vanden, and Leonard W. Schwartz. "A one-parameter family of sequence + * transformations." SIAM Journal on Mathematical Analysis 10.3 (1979): 658-666.
    • + *
  • [5] Barber, Michael N., and C. J. Hamer. "Extrapolation of sequences using a generalized + * epsilon-algorithm." The ANZIAM Journal 23.3 (1982): 229-240.
    • + *
+ *

+ */ +public final class WynnEpsilonComplex extends SeriesAlgorithmComplex { + + public static enum ShanksMethod { + WYNN, SHANKS, GENERALIZED, ALTERNATING + } + + private final ShanksMethod myMethod; + private final double myAlpha; + + private Complex myTerm; + private final Complex[] myTab; + private final Complex[] myTab2; + private final Complex[][] myEps; + private final Complex[][] myF; + + /** + * Creates a new instance of the Shanks transform with customizable method. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param method the method to use in computing the Shanks transform + * @param alpha the parameter alpha to use for the one-parameter epsilon transform; ignored if + * method is not set appropriately + */ + public WynnEpsilonComplex(final double tolerance, final int maxIters, final int patience, + final ShanksMethod method, final double alpha) { + super(tolerance, maxIters, patience); + myMethod = method; + myAlpha = alpha; + + switch (method) { + case SHANKS: + myTab = new Complex[maxIters]; + myTab2 = new Complex[maxIters]; + myEps = myF = null; + break; + case GENERALIZED: + case ALTERNATING: + myTab = myTab2 = null; + myEps = new Complex[2][maxIters]; + myF = new Complex[2][maxIters]; + break; + case WYNN: + default: + myTab = new Complex[maxIters]; + myTab2 = null; + myEps = myF = null; + break; + } + } + + /** + * Creates a new instance of the Shanks transform with customizable method. + * + * @param tolerance the smallest acceptable change in series evaluation in consecutive iterations + * that indicates the algorithm has converged + * @param maxIters the maximum number of sequence terms to evaluate before giving up + * @param patience the number of consecutive iterations required for the tolerance criterion to be + * satisfied to stop the algorithm + * @param method the method to use in computing the Shanks transform + */ + public WynnEpsilonComplex(final double tolerance, final int maxIters, final int patience, + final ShanksMethod method) { + this(tolerance, maxIters, patience, method, -1.0); + } + + public WynnEpsilonComplex(final double tolerance, final int maxIters, final int patience) { + this(tolerance, maxIters, patience, ShanksMethod.WYNN); + } + + @Override + public final Complex next(final Complex e, final Complex term) { + switch (myMethod) { + case SHANKS: + return nextGTransform(e, term); + case GENERALIZED: + case ALTERNATING: + return nextAlternatingEpsilon(e, term); + case WYNN: + default: + return nextWynnEpsilon(e, term); + } + } + + private final Complex nextAlternatingEpsilon(final Complex e, final Complex term) { + + // i is the index of old diagonal, j is new + final int i = myIndex & 1; + final int j = (i + 1) & 1; + final Complex[] eps0 = myEps[i], eps1 = myEps[j]; + final Complex[] f0 = myF[i], f1 = myF[j]; + + // compute the new diagonals of the epsilon and f tables + Complex result = f1[0] = term; + for (int k = 0; k < myIndex; ++k) { + if ((k & 1) == 0) { + + // finite difference for epsilon values + final int idx = k >> 1; + final Complex diff = f1[idx].subtract(f0[idx]); + if (diff.norm() < TINY) { + return Complex.NaN; + } + + // determine the parameter alpha + final double alpha; + if (myMethod == ShanksMethod.ALTERNATING) { + alpha = -1.0 * (idx & 1); + } else { + alpha = myAlpha; + } + + // update epsilon + if (idx == 0) { + eps1[idx] = diff.reciprocal(); + } else { + eps1[idx] = eps0[idx - 1].multiply(alpha).add(diff.reciprocal()); + } + } else { + + // finite difference for f values + final int idx = (k >> 1) + 1; + final Complex diff = eps1[idx - 1].subtract(eps0[idx - 1]); + if (diff.norm() < TINY) { + return Complex.NaN; + } + + // update f + result = f1[idx] = f0[idx - 1].add(diff.reciprocal()); + } + } + ++myIndex; + return result; + } + + private final Complex nextWynnEpsilon(final Complex e, final Complex term) { + + // initialization + myTab[myIndex] = term; + if (myIndex == 0) { + ++myIndex; + return term; + } + + // main loop of Wynn Epsilon algorithm + Complex aux = Complex.ZERO; + for (int j = myIndex; j >= 1; --j) { + + // compute the next epsilon + final Complex temp = aux; + aux = myTab[j - 1]; + Complex diff = myTab[j].subtract(aux); + + // correct denominators for underflow + if (diff.norm() <= TINY) { + diff = new Complex(HUGE, HUGE); + } else { + + diff = temp.add(diff.reciprocal()); + } + myTab[j - 1] = diff; + } + + // prepare result + final Complex result = myTab[myIndex & 1]; + ++myIndex; + return result; + } + + private final Complex nextGTransform(final Complex e, final Complex term) { + + // wait to compute series difference + if (myIndex == 0) { + ++myIndex; + myTerm = term; + return term; + } + + // initialize table on the first element and difference + if (myIndex == 1) { + myTab2[0] = term.subtract(myTerm); + myTab[0] = myTerm; + ++myIndex; + myTerm = term; + return term; + } + + // process the next element + ++myIndex; + final int idx = Math.min(myIndex - 1, myMaxIters) - 1; + Complex tmp1 = term.subtract(myTerm); + Complex tmp2 = Complex.ZERO, tmp3 = Complex.ZERO; + Complex extrap = myTerm; + for (int i = 0; i <= idx - 1; ++i) { + if (myTab2[i].norm() < TINY) { + myTerm = term; + return Complex.NaN; + } + final Complex rs = tmp1.divide(myTab2[i]).subtract(1.0); + final Complex tmp0; + if (i == 0) { + tmp0 = rs; + } else { + tmp0 = rs.multiply(myTab2[i - 1]); + } + if ((i & 1) == 0) { + if (rs.norm() < TINY) { + myTerm = term; + return Complex.NaN; + } + if (i == 0) { + tmp3 = myTab[0].subtract(extrap.subtract(myTab[0])).divide(rs); + myTab[0] = extrap; + if (myIndex <= 4) { + extrap = tmp3; + } + } else { + extrap = myTab[i - 1].subtract(myTab[i].subtract(myTab[i - 1])).divide(rs); + myTab[i - 1] = myTab[i]; + myTab[i] = tmp3; + tmp3 = extrap; + } + } + if (i > 0) { + myTab2[i - 1] = tmp2; + } + tmp2 = tmp1; + tmp1 = tmp0; + } + myTab2[idx - 1] = tmp2; + if (idx < myMaxIters) { + myTab2[idx] = tmp1; + } + myTab[idx] = extrap; + myTerm = term; + return extrap; + } + + @Override + public final String getName() { + switch (myMethod) { + case SHANKS: + return "Shanks"; + case GENERALIZED: + return "Gen. Epsilon(" + myAlpha + ")"; + case ALTERNATING: + return "Alt. Epsilon"; + case WYNN: + default: + return "Wynn Epsilon"; + } + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/WynnRhoComplex.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/WynnRhoComplex.java new file mode 100644 index 0000000000..4703b1a896 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/numerics/series/dp/complex/WynnRhoComplex.java @@ -0,0 +1,109 @@ +package org.matheclipse.core.numerics.series.dp.complex; + +import org.hipparchus.complex.Complex; + +/** + * Implements the Wynn Rho algorithm for convergence of infinite series introduced in [1] and as + * described in [2]. + * + *

+ * References: + *

    + *
  • [1] Wynn, Peter. "On a procrustean technique for the numerical transformation of slowly + * convergent sequences and series." Mathematical Proceedings of the Cambridge Philosophical + * Society. Vol. 52. No. 4. Cambridge University Press, 1956.
    • + *
  • [2] Osada, Naoki. Acceleration methods for slowly convergent sequences and their + * applications. Diss. PhD thesis, Nagoya University, 1993.
    • + *
+ *

+ */ +public final class WynnRhoComplex extends SeriesAlgorithmComplex { + + private int status; + private Complex xx, dx, dd, xp, alpha; + private final Complex[] x; + private final Complex[][] myTab, myTab2; + + public WynnRhoComplex(final double tolerance, final int maxIters, final int patience) { + super(tolerance, maxIters, patience); + x = new Complex[maxIters + 1]; + myTab = new Complex[2][maxIters + 1]; + myTab2 = new Complex[2][maxIters + 1]; + } + + @Override + public final Complex next(final Complex e, final Complex term) { + if (myIndex == 0) { + xx = dx = dd = Complex.ZERO; + } + final Complex xxp = xx; + x[myIndex] = xx = term; + if (myIndex == 0) { + dx = xx; + } else if (myIndex == 1) { + final Complex dxp = dx; + dx = xx.subtract(xxp); + dd = dx.divide(dx.subtract(dxp)); + } else { + final Complex dxp = dx; + final Complex ddp = dd; + dx = xx.subtract(xxp); + dd = dx.divide(dx.subtract(dxp)); + alpha = dd.subtract(ddp).reciprocal().add(1.0); + alpha = wynnrho(alpha, myTab2, myIndex - 1, new Complex(-2.0)); + } + xp = xx; + if (myIndex >= 2) { + for (int m = 1; m <= myIndex + 1; ++m) { + xp = wynnrho(x[m - 1], myTab, m, alpha); + } + } + ++myIndex; + if (status == 1) { + return Complex.NaN; + } else { + return xp; + } + } + + private final Complex wynnrho(final Complex xx, final Complex[][] tab, final int n, + final Complex theta) { + status = 0; + if (n != 1) { + System.arraycopy(tab[1], 0, tab[0], 0, n); + } + tab[1][0] = xx; + if (n == 1) { + tab[1][1] = theta.negate().divide(xx); + return xx; + } + Complex drho = tab[1][0].subtract(tab[0][0]); + if (drho.norm() < TINY) { + status = 1; + return xx; + } + tab[1][1] = theta.negate().divide(drho); + for (int k = 2; k <= n; ++k) { + drho = tab[1][k - 1].subtract(tab[0][k - 1]); + if (drho.norm() < TINY) { + if ((k & 1) == 0) { + status = 1; + return xx; + } else { + return tab[1][k - 1]; + } + } + tab[1][k] = tab[0][k - 2].add((new Complex(k - 1.0).subtract(theta)).divide(drho)); + } + if ((n & 1) == 0) { + return tab[1][n]; + } else { + return tab[1][n - 1]; + } + } + + @Override + public final String getName() { + return "Wynn Rho"; + } +} diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/NSum.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/NSum.java index 1d2d452bd0..6015912d9b 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/NSum.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/NSum.java @@ -5,13 +5,15 @@ import org.matheclipse.core.eval.EvalEngine; import org.matheclipse.core.expression.F; import org.matheclipse.core.expression.ImplementationStatus; +import org.matheclipse.core.expression.S; import org.matheclipse.core.interfaces.IAST; import org.matheclipse.core.interfaces.IASTMutable; import org.matheclipse.core.interfaces.IExpr; import org.matheclipse.core.interfaces.ISymbol; -import org.matheclipse.core.numerics.series.dp.Ensemble; -import org.matheclipse.core.numerics.series.dp.SeriesAlgorithm; -import org.matheclipse.core.numerics.series.dp.SeriesAlgorithm.SeriesSolution; +import org.matheclipse.core.numerics.series.dp.complex.DirectComplex; +import org.matheclipse.core.numerics.series.dp.complex.EnsembleComplex; +import org.matheclipse.core.numerics.series.dp.complex.SeriesAlgorithmComplex; +import org.matheclipse.core.numerics.series.dp.complex.SeriesAlgorithmComplex.SeriesSolutionComplex; import org.matheclipse.core.numerics.utils.Sequences; public class NSum extends Sum { @@ -67,36 +69,46 @@ public IExpr evaluate(IAST ast, EvalEngine engine) { return arg1.mapThread(ast, 1); } - if (ast.isAST2() && ast.arg2().isList3() || ast.arg2().isList4()) { - // if (Summations.isConvergent(ast)) { - // ast = ast.apply(S.Sum); - // IAST preevaledSum = engine.preevalForwardBackwardAST(ast, 1); - // - // IExpr temp = evaluateSum(preevaledSum, engine); - // if (temp.isPresent()) { - // if (temp.isNumericFunction(true)) { - // return engine.evalN(temp); - // } - // if (!temp.isFree(S.Sum, true)) { - // temp = F.subst(temp, x -> x == S.Sum ? S.NSum : F.NIL); - // } - // return temp; - // } - // } - - if (ast.isAST2() && ast.arg2().isList3()) { - IExpr function = arg1; - IAST limits = (IAST) engine.evaluate(ast.arg2()); - IExpr variable = limits.arg1(); - IExpr lowerLimit = limits.arg2(); - IExpr upperLimit = limits.arg3(); - IASTMutable copy = ast.copy(); - copy.set(1, function); - copy.set(2, limits); - IExpr temp = nsum(function, variable, lowerLimit, upperLimit, copy); + if (ast.isAST2() && ast.arg2().isList3()) { + IAST limits = (IAST) engine.evaluate(ast.arg2()); + IExpr variable = limits.arg1(); + IExpr lowerLimit = limits.arg2(); + IExpr upperLimit = limits.arg3(); + long start = lowerLimit.toLongDefault(); + long end = upperLimit.toLongDefault() + 1; + boolean preevaluateSymbolic = true; + if (start != Long.MIN_VALUE || end != Long.MIN_VALUE) { + long range = end - start; + if (range > 10000 || range < 10000) { + preevaluateSymbolic = false; + } + } + if (preevaluateSymbolic) { + // if (Summations.isConvergent(ast)) { + ast = ast.apply(S.Sum); + IAST preevaledSum = engine.preevalForwardBackwardAST(ast, 1); + + IExpr temp = evaluateSum(preevaledSum, engine); if (temp.isPresent()) { + if (temp.isNumericFunction(true)) { + return engine.evalN(temp); + } + if (!temp.isFree(S.Sum, true)) { + temp = F.subst(temp, x -> x == S.Sum ? S.NSum : F.NIL); + } return temp; } + // } + } + + IExpr function = arg1; + + IASTMutable copy = ast.copy(); + copy.set(1, function); + copy.set(2, limits); + IExpr temp = nsum(function, variable, lowerLimit, upperLimit, copy); + if (temp.isPresent()) { + return temp; } } @@ -150,41 +162,65 @@ private static IExpr sumStartToInfinity(IExpr function, IExpr variable, IExpr lo long start = lowerLimit.toLongDefault(); if (start != Long.MIN_VALUE) { // if (Summations.isConvergent(ast)) { - LongDoubleFunction longDoubleFunction = new LongDoubleFunction(function, variable); - Iterable iter = Sequences.toIterable(longDoubleFunction, start); - SeriesAlgorithm alg = new Ensemble(1e-8, 1000, 5); - SeriesSolution limit = alg.limit(iter, true); - return F.num(limit.limit); + LongComplexFunction longComplexFunction = new LongComplexFunction(function, variable); + Iterable seq = Sequences.toIterable(longComplexFunction, start); + SeriesAlgorithmComplex alg = new EnsembleComplex(1e-8, 1000, 5); + SeriesSolutionComplex limit = alg.limit(seq, true); + return F.inexactNum(limit.limit); - // LongComplexFunction longComplexFunction = new LongComplexFunction(function, variable); + // LongDoubleFunction longDoubleFunction = new LongDoubleFunction(function, variable); // - // Iterable iter = Sequences.toIterable(longComplexFunction, start); + // Iterable seq = Sequences.toIterable(longDoubleFunction, start); // SeriesAlgorithm alg = new Ensemble(1e-8, 1000, 5); - // SeriesSolution limit = alg.limit(iter, true); + // SeriesSolution limit = alg.limit(seq, true); // return F.num(limit.limit); - // } } } else { long start = lowerLimit.toLongDefault(); - long end = upperLimit.toLongDefault(); - if (start != Long.MIN_VALUE && end != Long.MIN_VALUE) { - LongDoubleFunction longDoubleFunction = new LongDoubleFunction(function, variable); - - Iterable iter = Sequences.toIterable(longDoubleFunction, start, end); - SeriesAlgorithm alg = new Ensemble(1e-8, 1000, 5); - SeriesSolution limit = alg.limit(iter, true); - return F.num(limit.limit); - - // LongComplexFunction longComplexFunction = new LongComplexFunction(function, variable); - // Iterator iter = Sequences.toIterable(longComplexFunction, start, - // end).iterator(); - // Complex result = Complex.ZERO; - // for (long i = start; i < end; i++) { - // Complex c = iter.next(); - // result = result.add(c); - // } - // return F.complexNum(result); + long end = upperLimit.toLongDefault() + 1; + if (start != Long.MIN_VALUE && end != Long.MIN_VALUE && start < end) { + + + long range = end + 1 - start; + if (range < Integer.MAX_VALUE) { + LongComplexFunction longComplexFunction = new LongComplexFunction(function, variable); + Iterable seq = Sequences.toIterable(longComplexFunction, start, end); + SeriesAlgorithmComplex alg = new DirectComplex(1e-8, (int) range, (int) range); + SeriesSolutionComplex limit = alg.limit(seq, true); + Complex result = limit.limit; + if (result.isNaN()) { + alg = new EnsembleComplex(1e-8, 1000, 5); + limit = alg.limit(seq, true); + return F.inexactNum(limit.limit); + } + return F.inexactNum(result); + // Complex result = Complex.ZERO; + // for (final Complex e : seq) { + // result = result.add(e); + // } + // return F.complexNum(result); + + // LongDoubleFunction longDoubleFunction = new LongDoubleFunction(function, variable); + // Iterable seq = Sequences.toIterable(longDoubleFunction, start, end); + // double result = 0.0; + // for (final double e : seq) { + // result += e; + // } + // return F.num(result); + } else { + LongComplexFunction longComplexFunction = new LongComplexFunction(function, variable); + Iterable seq = Sequences.toIterable(longComplexFunction, start, end); + SeriesAlgorithmComplex alg = new EnsembleComplex(1e-8, 1000, 5); + SeriesSolutionComplex limit = alg.limit(seq, true); + return F.inexactNum(limit.limit); + + // LongDoubleFunction longDoubleFunction = new LongDoubleFunction(function, variable); + // Iterable seq = Sequences.toIterable(longDoubleFunction, start, end); + // SeriesAlgorithm alg = new Ensemble(1e-8, 1000, 5); + // SeriesSolution limit = alg.limit(seq, true); + // return F.num(limit.limit); + } } } return F.NIL; diff --git a/symja_android_library/matheclipse-core/src/main/resources/doc/functions/NSum.md b/symja_android_library/matheclipse-core/src/main/resources/doc/functions/NSum.md new file mode 100644 index 0000000000..c0223c8d89 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/main/resources/doc/functions/NSum.md @@ -0,0 +1,25 @@ +## NSum + +``` +NSum(expr, {i, imin, imax}) +``` + +> evaluates the numerical approximated sum of `expr` with `i` ranging from `imin` to `imax`. + +``` +NSum(expr, {i, imin, imax, di}) +``` + +> `i` ranges from `imin` to `imax` in steps `di` + +See +* [Wikipedia - Series (mathematics)](https://en.wikipedia.org/wiki/Series_(mathematics)) +* [Wikipedia - Summation](https://en.wikipedia.org/wiki/Summation) +* [Wikipedia - Convergence_tests](https://en.wikipedia.org/wiki/Convergence_tests) + +### Examples + +``` +>> NSum(k, {k, 1, 10}) +55 +``` diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SumTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SumTest.java index a4557449b4..7f64808ce5 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SumTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SumTest.java @@ -452,11 +452,12 @@ public void testSumPolyLogInifinity() { @Test public void testNSum001() { - // TODO multi dimensional case // check("approxSum = NSum((-1)^n*(2/n)^k/k^2, {n,2,Infinity}, {k,1,Infinity})", // // ""); + check("NSum(k, {k, 1, 10})", // + "55.0"); check("approxSum = NSum(1/i^2, {i,100,Infinity})", // "0.0100502"); check("restSum = NSum(1/i^2, {i, 10^6,Infinity})", // @@ -464,12 +465,12 @@ public void testNSum001() { check("approxSum - restSum", // "0.0100497"); check("NSum(1/i^2, {i, 100, 10^6})", // - "0.0100502"); + "0.0100492"); check("NSum((-5)^i/i!, {i, 0, Infinity})", // "0.00673795"); // - // check("NSum(Log(x)/x^(2+2*I), {x, 1,Infinity})", // - // "-0.181854+I*(-0.13788)"); + check("NSum(Log(x)/x^(2+2*I), {x, 1,Infinity})", // + "-0.182175+I*(-0.136618)"); } @Test @@ -480,7 +481,6 @@ public void testNSum002() { check("NSum(1/(n^2*Log(n)), {n,2,Infinity})", // "0.605523"); - // reciprocal power checkNumeric("NSum((-1)^(k-1)/k^1.5, {k,1,Infinity})", // "0.7651470246254165"); @@ -533,6 +533,6 @@ public void testNSum003() { // check("NSum(E^x*Piecewise({{1/(E*x!), x >= 0}}, 0),{x,-Infinity,Infinity})", // // ""); checkNumeric("NSum(1/E^(k^2), {k,-Infinity,Infinity})", // - "1.7726372048266523"); + "1.7726372048266525"); } }