Package gov.nih.mipav.model.algorithms

Class ConfluentHypergeometric

java.lang.Object
gov.nih.mipav.model.algorithms.ConfluentHypergeometric
public class ConfluentHypergeometric extends Object

This code calculates the confluent hypergeometric function of the first and second kinds For the confluent hypergeometric function of the first kind a typical usage for the routine requiring real parameters and a real argument would be: 

 double result[] = new double[1];
 ConfluentHypergeometric ch = new ConfluentHypergeometric( -0.5, 1, 1.0, result);
 ch.run();
 Preferences.debug("Confluent hypergeomtric result = " + result[0] + "\n");
 UI.setDataText("Confluent hypergeometric result = " + result[0] + "\n");

For the confluent hypergeometric function of the first kind a typcial usage for the routine requiring real parameters and a complex argument would be: 

 ConfluentHypergeometric cf = new ConfluentHypergeometric( -0.5, 1.0, realZ, imagZ, realResult, imagResult);
 ch.run();
 System.out.println("realResult[0] = " + realResult[0]);
 System.out.println("imagResult[0] = " + imagResult[0]);

For the confluent hypergeometric function of the first kind a typical usage for the routine allowing complex parameters and a complex argument would be: 

 ConfluentHypergeometric cf;
 int Lnchf = 0;
 int ip = 700;
 double realResult[] = new double[1];
 double imagResult[] = new double[1];
 cf = new ConfluentHypergeometric( -0.5, 0.0, 1.0, 0.0, 10.0, 0.0, Lnchf, ip, realResult, imagResult);
 cf.run();
 System.out.println("realResult[0] = " + realResult[0]);
 System.out.println("imagResult[0] = " + imagResult[0]);

Since this complex routine allows large magnitude x, it should generally be used instead of the real version routine.

For the confluent hypergeometric function of the second kind a typical usage for the routine requiring real parameters and a real argument would be: 

   double result[] = new double[1];
   int method[] = new int[1];
   ConfluentHypergeometric ch = new ConfluentHypergeometric(-0.5, 1, 1.0, result, method);
   ch.run();
   Preferences.debug("Confluent hypergeomtric result = " + result[0] +  " method = " + 
                       method[0] + \n");

For 1F1(-1/2, 1, x) tested with Shanjie Zhang and Jianming Jin Computation of Special Functions CHGM routine and ACM Algorithm 707 conhyp routine by Mark Nardin, W. F. Perger, and Atul Bhalla the results were:

  • From x = -706 to x = +184 the results of the 2 routines matched to at least 1 part in 1E5.
  • For x = -3055 to x = +184 the ACM routine seemed to give valid results.
  • For x invalid input: '<'= -3056 the ACM routine results oscillated wildly and at x = -3200 the result was frozen at 1.0E75.
  • For x >= +185 the ACM result was frozen at 1.0E75.
  • The log result was also seen to oscillate for x invalid input: '<'= -3056.
  • For x = -706 to x = +781 the CHGM routine seemed to give valid results.
  • For x = -707 to x = -745 the CHGM routine gave infinity and for x invalid input: '<'= -746 the CHGM routine gave NaN.
  • For x >= +783 the CHGM routine gave a result of -infinity.
  • For 1F1(-1/2, 2, x) the results were very similar:
  • From x = -706 to x = +189 the results of the 2 routines matched to at least 1 part in 1E5.
  • From x = -3058 to x = +189 the ACM routine seemed to give valid results.
  • At x invalid input: '<'= -3059 the ACM results oscillated wildly.
  • For x > +189 the ACM result was frozen at 1.0E75.
  • For x = -706 to x = +791 the CHGM routine seemed to give valid results.
  • For x = -707 to x = -745 the CHGM routine gave infinity and for x invalid input: '<'= -746 the CHGM routine gave NaN.
  • For x >= +792 the CHGM routine gave a result of -infinity.

From the ACM website: 

 ACM Software License Agreement
 
  All software, both binary and source published by the Association for Computing Machinery (hereafter, Software) is copyrighted by the Association (hereafter, ACM) and ownership of all right, title and interest in and to the Software remains with ACM. By using or copying the Software, User agrees to abide by the terms of this Agreement.
  Noncommercial Use
 
  The ACM grants to you (hereafter, User) a royalty-free, nonexclusive right to execute, copy, modify and distribute both the binary and source code solely for academic, research and other similar noncommercial uses, subject to the following conditions:
 
  1. User acknowledges that the Software is still in the development stage and that it is being supplied "as is," without any support services from ACM. Neither ACM nor the author makes any representations or warranties, express or implied, including, without limitation, any representations or warranties of the merchantability or fitness for any particular purpose, or that the application of the software, will not infringe on any patents or other proprietary rights of others.
  2. ACM shall not be held liable for direct, indirect, incidental or consequential damages arising from any claim by User or any third party with respect to uses allowed under this Agreement, or from any use of the Software.
  3. User agrees to fully indemnify and hold harmless ACM and/or the author(s) of the original work from and against any and all claims, demands, suits, losses, damages, costs and expenses arising out of the User's use of the Software, including, without limitation, arising out of the User's modification of the Software.
  4. User may modify the Software and distribute that modified work to third parties provided that: (a) if posted separately, it clearly acknowledges that it contains material copyrighted by ACM (b) no charge is associated with such copies, (c) User agrees to notify ACM and the Author(s) of the distribution, and (d) User clearly notifies secondary users that such modified work is not the original Software.
  5. User agrees that ACM, the authors of the original work and others may enjoy a royalty-free, non-exclusive license to use, copy, modify and redistribute these modifications to the Software made by the User and distributed to third parties as a derivative work under this agreement.
  6. This agreement will terminate immediately upon User's breach of, or non-compliance with, any of its terms. User may be held liable for any copyright infringement or the infringement of any other proprietary rights in the Software that is caused or facilitated by the User's failure to abide by the terms of this agreement.
  7. This agreement will be construed and enforced in accordance with the law of the state of New York applicable to contracts performed entirely within the State. The parties irrevocably consent to the exclusive jurisdiction of the state or federal courts located in the City of New York for all disputes concerning this agreement. 
 
  Commercial Use
 
  Any User wishing to make a commercial use of the Software must contact ACM at permissions@acm.org to arrange an appropriate license. Commercial use includes (1) integrating or incorporating all or part of the source code into a product for sale or license by, or on behalf of, User to third parties, or (2) distribution of the binary or source code to third parties for use with a commercial product sold or licensed by, or on behalf of, User.
 
  Revised 6/98

Excerpt from Appendix C of Computation of Special Functions by Shanjie Zhang and Jianming Jin:

All the programs and subroutines contained in this diskette are copyrighted. However, we give permission to the reader who purchases this book to incorporate any of these programs into his or her programs provided that the copyright is acknowledged. 
 DISCLAIMER OF WARRANTY
 
 Although we have made a great effort to test and validate the 
 computer programs, we make no warranties, express or implied, that 
 these  programs  are  free  of  error,  or are consistent with any 
 particular  standard  of  merchantability,  or that they will meet 
 your requirements for any particular application.  They should not 
 be relied on for  solving problems  whose incorrect solution could 
 result in  injury to  a person or loss of property.  If you do use 
 the programs in such a manner, it is at your own risk. The authors 
 and publisher  disclaim all liability for  direct or consequential 
 damages resulting from your use of the programs.

  • Field Summary

    Fields
    Modifier and Type
    Field
    Description
    private double
    a
    Input parameter
    private double
    b
    Input parameter
    private int
    bit
     
    static final int
    COMPLEX_VERSION
     
    static final int
    CONFLUENT_HYPERGEOMETRIC_FIRST_KIND
    Confluent Hypergeometric Function of the First Kind.
    static final int
    CONFLUENT_HYPERGEOMETRIC_SECOND_KIND
    Confluent Hypergeometric Function of the Second Kind
    private double
    imagA
     
    private double
    imagB
     
    private double[]
    imagResult
     
    private double
    imagZ
     
    private int
    ip
    Specifies how many array positions are desired (usually 10 is sufficient).
    private int
    kind
    Tells whether confluent hypergeometric function is of first or second kind
    private int
    L
    Size of the arrays
    private final int
    length
    Create arrays of size length+2
    private int
    Lnchf
    Tells how the result should be represented A '0' will return the result in standard exponential form.
    private int[]
    method
    Method used in calculating the confluent hypergeometric function of the second kind.
    static final int
    REAL_VERSION
     
    private double
    realA
    Input parameter
    private double
    realB
    Input parameter
    static final int
    REALPARAM_COMPLEXARG_VERSION
     
    private double[]
    realResult
    outputted result
    private double
    realZ
    Input argument
    private double[]
    result
    Outputted result
    private double
    rmax
    Number of digits required to represent the numbers with the required accuracy
    private int
    version
    Tells whether real number of complex number version
    private double
    x
    Input argument

  • Constructor Summary

    Constructors
    Constructor
    Description
    ConfluentHypergeometric()
     
    ConfluentHypergeometric(double a, double b, double x, double[] result, int[] method)
     
    ConfluentHypergeometric(double a, double b, double realZ, double imagZ, double[] realResult, double[] imagResult)
     
    ConfluentHypergeometric(double realA, double imagA, double realB, double imagB, double realZ, double imagZ, int Lnchf, int ip, double[] realResult, double[] imagResult)
     
    ConfluentHypergeometric(int kind, double a, double b, double x, double[] result)
     

  • Method Summary

    Modifier and Type
    Method
    Description
    private void
    aradd(double[] a, double[] b, double[] c)
    Accepts two arrays of numbers, a and b, and returns the sum of the array c.
    private void
    armult(double[] a, double b, double[] c)
    Accepts array a and scalar b, and returns the product array c.
    private void
    arsub(double[] a, double[] b, double[] c)
    Accepts two arrays and subtracts each element in the second array b from the element in the first array a and returns the solution c
    private void
    arydiv(double[] ar, double[] ai, double[] br, double[] bi)
    Returns the complex number resulting from the division of four arrays, representing two complex numbers.
    private int
    bits()
    Determines the number of significant figures of machine precision to arrive at the size of the array the numbers must be stored in to get the accuracy of the solution.
    private void
    chgubi(int[] id)
     
    private void
    chguit(int[] id)
     
    private void
    chgul(int[] id)
     
    private void
    chgus(int[] id)
     
    private void
    cmpadd(double[] ar, double[] ai, double[] br, double[] bi, double[] cr, double[] ci)
    Takes two arrays representing one real and one imaginary part, and adds two arrays representing another complex number and returns two arrays holding the complex sum.
    private void
    cmpmul(double[] ar, double[] ai, double br, double bi, double[] cr, double[] ci)
    Takes 2 arrays, ar and ai representing one real and one imaginary part, and multiplies it with br and bi, representing a complex number and returns the complex product
    void
    complexTestFirstKind()
     
    private void
    conv12(double realCN, double imagCN, double[][] cae)
    Converts a real and imaginary number pair to a form of a 2x2 real array.
    private void
    conv21(double[][] cae)
    Converts a number represented in a 2x2 real array to the form of a real and imaginary pair.
    private void
    eadd(double n1, double e1, double n2, double e2, double[] nf, double[] ef)
    Returns the sum of two numbers in the form of a base and an exponent.
    private void
    ecpdiv(double[][] a, double[][] b, double[][] c)
    Divides two numbers and returns the solution.
    private void
    ecpmul(double[][] a, double[][] b, double[][] c)
    Multiplies two numbers which are each represented in the form of a two by two array and returns the solution in the same form
    private void
    ediv(double n1, double e1, double n2, double e2, double[] nf, double[] ef)
    Returns the solution in the form of the base and exponent of the division of two exponential numbers.
    private void
    emult(double n1, double e1, double n2, double e2, double[] nf, double[] ef)
    Takes one base and exponent and multplies it by another number's base and exponent to give the product in the form of base and exponent
    private void
    esub(double n1, double e1, double n2, double e2, double[] nf, double[] ef)
    Returns the solution to the subtraction of two numbers in the form of base and exponent.
    void
    finalize()
    Cleanup memory.
    private void
    firstKindComplex()
    Port of Algorithm 707, Collected Algorithms from ACM.
    private void
    firstKindComplexArgument()
    This code is a port of the FORTRAN routine CCHG from the book Computation of Special Functions by Shanjie Zhang and Jianming Jin, John Wiley invalid input: '&' Sons, Inc., 1996, pp. 400-402.
    private void
    firstKindRealArgument()
    This code is a port of the FORTRAN routine CHGM from the book Computation of Special Functions by Shanjie Zhang and Jianming Jin, John Wiley invalid input: '&' Sons, Inc., 1996, pp. 398-400.
    void
    realTestFirstKind()
     
    void
    realTestSecondKind()
     
    void
    run()
     
    private void
    secondKindRealArgument()
    This code is a port of the FORTRAN routine CHGU from the book Computation of Special Functions by Shanjie Zhang and Jianming Jin, John Wiley invalid input: '&' Sons, Inc., 1996, pp. 403-405.
    private double
    store(double x)
    This function forces its argument x to be stored in a memory location, thus providing a means of determining floating point number characteristics (such as the machine precision) when it is necessary to avoid computation in high precision registers.

    Methods inherited from class java.lang.Object

    clone, equals, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

  • Field Details

    • CONFLUENT_HYPERGEOMETRIC_FIRST_KIND

      public static final int CONFLUENT_HYPERGEOMETRIC_FIRST_KIND
      Confluent Hypergeometric Function of the First Kind.
      See Also:

    • CONFLUENT_HYPERGEOMETRIC_SECOND_KIND

      public static final int CONFLUENT_HYPERGEOMETRIC_SECOND_KIND
      Confluent Hypergeometric Function of the Second Kind
      See Also:

    • REAL_VERSION

      public static final int REAL_VERSION
      See Also:

    • COMPLEX_VERSION

      public static final int COMPLEX_VERSION
      See Also:

    • REALPARAM_COMPLEXARG_VERSION

      public static final int REALPARAM_COMPLEXARG_VERSION
      See Also:

    • kind

      private int kind
      Tells whether confluent hypergeometric function is of first or second kind

    • version

      private int version
      Tells whether real number of complex number version

    • a

      private double a
      Input parameter

    • b

      private double b
      Input parameter

    • x

      private double x
      Input argument

    • result

      private double[] result
      Outputted result

    • method

      private int[] method
      Method used in calculating the confluent hypergeometric function of the second kind.

    • realA

      private double realA
      Input parameter

    • imagA

      private double imagA

    • realB

      private double realB
      Input parameter

    • imagB

      private double imagB

    • realZ

      private double realZ
      Input argument

    • imagZ

      private double imagZ

    • Lnchf

      private int Lnchf
      Tells how the result should be represented A '0' will return the result in standard exponential form. A '1' will return the log of the result.

    • ip

      private int ip
      Specifies how many array positions are desired (usually 10 is sufficient). More difficult cases may require this parameter to be 100 or more, and it is not always trivial to predict which cases these might be. One pragmatic approach is to simply run the evaluator using a value of 10, then increase ip, run the evaluator again and compare the returned results. Overwriting memory may occur if ip exceeds 776. Setting IP = 0 causes the program to estimate the number of array positions

    • realResult

      private double[] realResult
      outputted result

    • imagResult

      private double[] imagResult

    • length

      private final int length
      Create arrays of size length+2
      See Also:

    • L

      private int L
      Size of the arrays

    • rmax

      private double rmax
      Number of digits required to represent the numbers with the required accuracy

    • bit

      private int bit

  • Constructor Details

    • ConfluentHypergeometric

      public ConfluentHypergeometric()

    • ConfluentHypergeometric

      public ConfluentHypergeometric(int kind, double a, double b, double x, double[] result)
      Parameters:
      a - input parameter
      b - input parameter
      x - input argument
      result - outputted result

    • ConfluentHypergeometric

      public ConfluentHypergeometric(double a, double b, double x, double[] result, int[] method)
      Parameters:
      a - input parameter
      b - input parameter
      x - input argument, x > 0
      result - outputted result
      method - method code used

    • ConfluentHypergeometric

      public ConfluentHypergeometric(double a, double b, double realZ, double imagZ, double[] realResult, double[] imagResult)
      Parameters:
      a - input real parameter
      b - input real parameter
      realZ - input real part of argument
      imagZ - input imaginary part of argument
      realResult - real part of outputted result
      imagResult - imaginary part of outputted result

    • ConfluentHypergeometric

      public ConfluentHypergeometric(double realA, double imagA, double realB, double imagB, double realZ, double imagZ, int Lnchf, int ip, double[] realResult, double[] imagResult)
      Parameters:
      realA - Real part of first input parameter
      imagA - Imaginary part of first input parameter
      realB - Real part of second input parameter
      imagB - Imaginary part of second input parameter
      realZ - Real part of input argument
      imagZ - Imaginary part of input argument
      Lnchf - = 0 for standard result, = 1 for log of result
      ip - Number of array positions desired
      realResult - Real part of outputted result
      imagResult - Imaginary part of outputted result

  • Method Details

    • finalize

      public void finalize() throws Throwable
      Cleanup memory.
      Overrides:
      finalize in class Object
      Throws:
      Throwable - DOCUMENT ME!

    • run

      public void run()

    • firstKindComplex

      private void firstKindComplex()
      Port of Algorithm 707, Collected Algorithms from ACM. Original Work published in Transactions on Mathematical Software, Vol. 18. No. 3, September, 1992, pp. 345-349. This is a port of Solution to the Confluent Hypergeometric Function by Mark Nardin, W. F. Perger, and Atul Bhalla, Michigan Technological University, Copyright 1989. Description: A numerical evaluator for the confluent hypergeometric function for complex arguments with large magnitudes using a direct summation of the Kummer series. The method used allows an accuracy of up to thirteen decimal places through the use of large arrays and a single final division. The confluent hypergeometric function of the first kind is the solution to the equation: zf"(z) + (a-z)f'(z) - bf(z) = 0 Reference: Numerical evaluation of the confluent hypergeometric function for complex arguments of large magnitudes by Mark Nardin, W. F. Perger, and Atul Bhalla, Journal of Computational and Applied Mathematics, Vol. 39, 1992, pp. 193-200.

    • arydiv

      private void arydiv(double[] ar, double[] ai, double[] br, double[] bi)
      Returns the complex number resulting from the division of four arrays, representing two complex numbers. The number returned will be in one of two different forms. Either standard scientific or as the log of the number.
      Parameters:
      ar -
      ai -
      br -
      bi -

    • conv21

      private void conv21(double[][] cae)
      Converts a number represented in a 2x2 real array to the form of a real and imaginary pair.
      Parameters:
      cae -

    • emult

      private void emult(double n1, double e1, double n2, double e2, double[] nf, double[] ef)
      Takes one base and exponent and multplies it by another number's base and exponent to give the product in the form of base and exponent
      Parameters:
      n1 -
      e1 -
      n2 -
      e2 -
      nf -
      ef -

    • ediv

      private void ediv(double n1, double e1, double n2, double e2, double[] nf, double[] ef)
      Returns the solution in the form of the base and exponent of the division of two exponential numbers.
      Parameters:
      n1 -
      e1 -
      n2 -
      e2 -
      nf -
      ef -

    • eadd

      private void eadd(double n1, double e1, double n2, double e2, double[] nf, double[] ef)
      Returns the sum of two numbers in the form of a base and an exponent.
      Parameters:
      n1 -
      e1 -
      n2 -
      e2 -
      nf -
      ef -

    • esub

      private void esub(double n1, double e1, double n2, double e2, double[] nf, double[] ef)
      Returns the solution to the subtraction of two numbers in the form of base and exponent.
      Parameters:
      n1 -
      e1 -
      n2 -
      e2 -
      nf -
      ef -

    • ecpmul

      private void ecpmul(double[][] a, double[][] b, double[][] c)
      Multiplies two numbers which are each represented in the form of a two by two array and returns the solution in the same form
      Parameters:
      a -
      b -
      c -

    • ecpdiv

      private void ecpdiv(double[][] a, double[][] b, double[][] c)
      Divides two numbers and returns the solution. All numbers are represented by a 2x2 array.
      Parameters:
      a -
      b -
      c -

    • conv12

      private void conv12(double realCN, double imagCN, double[][] cae)
      Converts a real and imaginary number pair to a form of a 2x2 real array.
      Parameters:
      realCN -
      imagCN -
      cae -

    • cmpadd

      private void cmpadd(double[] ar, double[] ai, double[] br, double[] bi, double[] cr, double[] ci)
      Takes two arrays representing one real and one imaginary part, and adds two arrays representing another complex number and returns two arrays holding the complex sum. (CR, CI) = (AR+BR, AI+BI)
      Parameters:
      ar -
      ai -
      br -
      bi -
      cr -
      ci -

    • bits

      private int bits()
      Determines the number of significant figures of machine precision to arrive at the size of the array the numbers must be stored in to get the accuracy of the solution.

    • store

      private double store(double x)
      This function forces its argument x to be stored in a memory location, thus providing a means of determining floating point number characteristics (such as the machine precision) when it is necessary to avoid computation in high precision registers.
      Parameters:
      x - The value to be stored
      Returns:
      The value of x after it has been stored and possibly truncated or rounded to the double precision word length.

    • cmpmul

      private void cmpmul(double[] ar, double[] ai, double br, double bi, double[] cr, double[] ci)
      Takes 2 arrays, ar and ai representing one real and one imaginary part, and multiplies it with br and bi, representing a complex number and returns the complex product
      Parameters:
      ar -
      ai -
      br -
      bi -
      cr -
      ci -

    • armult

      private void armult(double[] a, double b, double[] c)
      Accepts array a and scalar b, and returns the product array c.
      Parameters:
      a -
      b -
      c -

    • aradd

      private void aradd(double[] a, double[] b, double[] c)
      Accepts two arrays of numbers, a and b, and returns the sum of the array c.
      Parameters:
      a -
      b -
      c -

    • arsub

      private void arsub(double[] a, double[] b, double[] c)
      Accepts two arrays and subtracts each element in the second array b from the element in the first array a and returns the solution c
      Parameters:
      a -
      b -
      c -

    • complexTestFirstKind

      public void complexTestFirstKind()

    • realTestFirstKind

      public void realTestFirstKind()

    • firstKindRealArgument

      private void firstKindRealArgument()
      This code is a port of the FORTRAN routine CHGM from the book Computation of Special Functions by Shanjie Zhang and Jianming Jin, John Wiley invalid input: '&' Sons, Inc., 1996, pp. 398-400. It computes confluent hypergeometric functions of the first kind for real parameters and argument. It works well except for the case of a invalid input: '<'invalid input: '<' -1 when x > 0 and a >> 1 when x invalid input: '<' 0. In this exceptional case, the evaluation involves the summation of a partially alternating series which results in a loss of significant digits.

    • firstKindComplexArgument

      private void firstKindComplexArgument()
      This code is a port of the FORTRAN routine CCHG from the book Computation of Special Functions by Shanjie Zhang and Jianming Jin, John Wiley invalid input: '&' Sons, Inc., 1996, pp. 400-402. It computes confluent hypergeometric functions of the first kind for real parameters and a complex argument. It works well except for the case of a invalid input: '<'invalid input: '<' -1 when x > 0 and a >> 1 when x invalid input: '<' 0. In this exceptional case, the evaluation involves the summation of a partially alternating series which results in a loss of significant digits.

    • realTestSecondKind

      public void realTestSecondKind()

    • secondKindRealArgument

      private void secondKindRealArgument()
      This code is a port of the FORTRAN routine CHGU from the book Computation of Special Functions by Shanjie Zhang and Jianming Jin, John Wiley invalid input: '&' Sons, Inc., 1996, pp. 403-405. It computes confluent hypergeometric functions of the second kind for real parameters and argument. It works well except for the case of a invalid input: '<'invalid input: '<' -1 because in this case the evaluation involves a summation of a partially alternating series which results in a loss of significant digits. Routines called: chgus for small x method = 1 chgul for large x method = 2 chgubi for integer b method = 3 chguit for numerical integration method = 4

    • chgus

      private void chgus(int[] id)

    • chgul

      private void chgul(int[] id)

    • chgubi

      private void chgubi(int[] id)

    • chguit

      private void chguit(int[] id)