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-2025.
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    private ScalarMath() {
55      // Prevent instantiation.
56    }
57  
58    /**
59     * b2u
60     *
61     * @param b a double.
62     * @return a double.
63     */
64    public static double b2u(double b) {
65      return b / eightPi2;
66    }
67  
68    /**
69     * binomial
70     *
71     * @param n a long.
72     * @param k a long.
73     * @return Returns the binomial of (n,k).
74     */
75    public static long binomial(long n, long k) {
76      return factorial(n) / (factorial(n - k) * factorial(k));
77    }
78  
79    /**
80     * Returns n!! Precondition: n .GE. -1 Returning 1 for -1 input is analogous to Maple behavior.
81     *
82     * @param n long.
83     * @return Returns the n!!.
84     */
85    public static long doubleFactorial(long n) {
86      if (n < -1) {
87        throw new RuntimeException("Underflow error in doubleFactorial");
88      } else if (n == -1 || n == 0 || n == 1) {
89        return 1;
90      } else {
91        return n * doubleFactorial(n - 2);
92      }
93    }
94  
95    /**
96     * Returns n! <br> Precondition: n .GE. 0 and n .LE. 20 <br> Max long = 9223372036854775807 <br>
97     * 20! = 2432902008176640000 is ok. <br> 21! returns an overflow: -4249290049419214848
98     *
99     * @param n long.
100    * @return Returns n!.
101    */
102   public static long factorial(long n) {
103     if (n < 0) {
104       throw new RuntimeException("Underflow error in factorial");
105     } else if (n > 20) {
106       throw new RuntimeException("Overflow error in factorial");
107     } else if (n == 0) {
108       return 1;
109     } else {
110       return n * factorial(n - 1);
111     }
112   }
113 
114   /**
115    * Compute 1.0 / (1.0 + exp(x)).
116    *
117    * @param x Input.
118    * @return Returns 1.0 / (1.0 + exp(x)).
119    */
120   public static double fermiFunction(double x) {
121     return 1.0 / (1.0 + exp(x));
122   }
123 
124   /**
125    * This is an atypical mod function used by crystallography methods.
126    *
127    * <p>mod
128    *
129    * @param a Value to mod.
130    * @param b Value to mod by.
131    * @return Positive a % b.
132    */
133   public static double mod(double a, double b) {
134     var res = a % b;
135     if (res < 0.0) {
136       res += b;
137     }
138     return res;
139   }
140 
141   /**
142    * This is an atypical mod function used by crystallography methods.
143    *
144    * @param a an int.
145    * @param b an int.
146    * @return an int.
147    */
148   public static int mod(int a, int b) {
149     int res = a % b;
150     if (res < 0) {
151       res += b;
152     }
153     return res;
154   }
155 
156   /**
157    * Atypical mod function used to move a value into the range lb &lt;= value &lt; ub, assuming the
158    * domain is periodic with a period of (ub - lb).
159    *
160    * @param value Value to move between bounds.
161    * @param lb    Lower bound.
162    * @param ub    Upper bound.
163    * @return Returns periodic copy of value, in the range lb &lt;= value &lt; ub.
164    */
165   public static double modToRange(double value, double lb, double ub) {
166     value -= lb;
167     var range = ub - lb;
168     value = mod(value, range);
169     value += lb;
170     return value;
171   }
172 
173   /**
174    * u2b
175    *
176    * @param u a double.
177    * @return a double.
178    */
179   public static double u2b(double u) {
180     return u * eightPi2;
181   }
182 
183   /**
184    * quadForm
185    *
186    * @param v   an array of double.
187    * @param mat an array of double.
188    * @return a double.
189    */
190   public static double quadForm(double[] v, double[][] mat) {
191     return (v[0] * (v[0] * mat[0][0] + 2 * (v[1] * mat[0][1] + v[2] * mat[0][2]))
192         + v[1] * (v[1] * mat[1][1] + 2 * (v[2] * mat[1][2]))
193         + v[2] * v[2] * mat[2][2]);
194   }
195 
196   /**
197    * Reflect proposed angles to be within 0 and 180 degrees.
198    *
199    * @param angle Proposed angle in radians.
200    * @return value of mirrored angle in radians.
201    */
202   public static double mirrorRadians(double angle) {
203     double angleDegrees = angle * 180 / PI;
204     return mirrorDegrees(angleDegrees) * PI / 180;
205   }
206 
207   /**
208    * Reflect proposed angles to be within 0 and 180 degrees.
209    *
210    * @param angle Proposed angle in radians.
211    * @return value of mirrored angle in radians.
212    */
213   public static double mirrorDegrees(double angle) {
214     if (angle > 180.0) {
215       angle = 180.0 - (angle - 180.0);
216     } else if (angle < 0.0) {
217       angle = 0.0 - angle;
218     }
219     return angle;
220   }
221 }