View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.numerics.math;
39  
40  import static org.apache.commons.math3.util.FastMath.PI;
41  import static org.apache.commons.math3.util.FastMath.exp;
42  
43  /**
44   * The ScalarMath class is a simple math library that operates on single variables
45   *
46   * @author Jacob M. Litman
47   * @author Michael J. Schnieders
48   * @since 1.0
49   */
50  public class ScalarMath {
51  
52    private static final double eightPi2 = 8.0 * PI * PI;
53  
54    /**
55     * b2u
56     *
57     * @param b a double.
58     * @return a double.
59     */
60    public static double b2u(double b) {
61      return b / eightPi2;
62    }
63  
64    /**
65     * binomial
66     *
67     * @param n a long.
68     * @param k a long.
69     * @return Returns the binomial of (n,k).
70     */
71    public static long binomial(long n, long k) {
72      return factorial(n) / (factorial(n - k) * factorial(k));
73    }
74  
75    /**
76     * Returns n!! Precondition: n .GE. -1 Returning 1 for -1 input is analogous to Maple behavior.
77     *
78     * @param n long.
79     * @return Returns the n!!.
80     */
81    public static long doubleFactorial(long n) {
82      if (n < -1) {
83        throw new RuntimeException("Underflow error in doubleFactorial");
84      } else if (n == -1 || n == 0 || n == 1) {
85        return 1;
86      } else {
87        return n * doubleFactorial(n - 2);
88      }
89    }
90  
91    /**
92     * Returns n! <br> Precondition: n .GE. 0 and n .LE. 20 <br> Max long = 9223372036854775807 <br>
93     * 20! = 2432902008176640000 is ok. <br> 21! returns an overflow: -4249290049419214848
94     *
95     * @param n long.
96     * @return Returns n!.
97     */
98    public static long factorial(long n) {
99      if (n < 0) {
100       throw new RuntimeException("Underflow error in factorial");
101     } else if (n > 20) {
102       throw new RuntimeException("Overflow error in factorial");
103     } else if (n == 0) {
104       return 1;
105     } else {
106       return n * factorial(n - 1);
107     }
108   }
109 
110   /**
111    * Compute 1.0 / (1.0 + exp(x)).
112    *
113    * @param x Input.
114    * @return Returns 1.0 / (1.0 + exp(x)).
115    */
116   public static double fermiFunction(double x) {
117     return 1.0 / (1.0 + exp(x));
118   }
119 
120   /**
121    * This is an atypical mod function used by crystallography methods.
122    *
123    * <p>mod
124    *
125    * @param a Value to mod.
126    * @param b Value to mod by.
127    * @return Positive a % b.
128    */
129   public static double mod(double a, double b) {
130     var res = a % b;
131     if (res < 0.0) {
132       res += b;
133     }
134     return res;
135   }
136 
137   /**
138    * This is an atypical mod function used by crystallography methods.
139    *
140    * @param a an int.
141    * @param b an int.
142    * @return an int.
143    */
144   public static int mod(int a, int b) {
145     int res = a % b;
146     if (res < 0) {
147       res += b;
148     }
149     return res;
150   }
151 
152   /**
153    * Atypical mod function used to move a value into the range lb &lt;= value &lt; ub, assuming the
154    * domain is periodic with a period of (ub - lb).
155    *
156    * @param value Value to move between bounds.
157    * @param lb Lower bound.
158    * @param ub Upper bound.
159    * @return Returns periodic copy of value, in the range lb &lt;= value &lt; ub.
160    */
161   public static double modToRange(double value, double lb, double ub) {
162     value -= lb;
163     var range = ub - lb;
164     value = mod(value, range);
165     value += lb;
166     return value;
167   }
168 
169   /**
170    * u2b
171    *
172    * @param u a double.
173    * @return a double.
174    */
175   public static double u2b(double u) {
176     return u * eightPi2;
177   }
178 
179   /**
180    * quadForm
181    *
182    * @param v an array of double.
183    * @param mat an array of double.
184    * @return a double.
185    */
186   public static double quadForm(double[] v, double[][] mat) {
187     return (v[0] * (v[0] * mat[0][0] + 2 * (v[1] * mat[0][1] + v[2] * mat[0][2]))
188         + v[1] * (v[1] * mat[1][1] + 2 * (v[2] * mat[1][2]))
189         + v[2] * v[2] * mat[2][2]);
190   }
191 
192   /**
193    * Reflect proposed angles to be within 0 and 180 degrees.
194    *
195    * @param angle Proposed angle in radians.
196    * @return value of mirrored angle in radians.
197    */
198   public static double mirrorRadians(double angle) {
199     double angleDegrees = angle * 180 / PI;
200     return mirrorDegrees(angleDegrees) * PI / 180;
201   }
202 
203   /**
204    * Reflect proposed angles to be within 0 and 180 degrees.
205    *
206    * @param angle Proposed angle in radians.
207    * @return value of mirrored angle in radians.
208    */
209   public static double mirrorDegrees(double angle) {
210     if (angle > 180.0) {
211       angle = 180.0 - (angle - 180.0);
212     } else if (angle < 0.0) {
213       angle = 0.0 - angle;
214     }
215     return angle;
216   }
217 }