       1 package MathUtil;
       2 #
       3 # File: MathUtil.pm
       4 # Author: Manish Sud <msud@san.rr.com>
       5 #
       6 # Copyright (C) 2018 Manish Sud. All rights reserved.
       7 #
       8 # This file is part of MayaChemTools.
       9 #
      10 # MayaChemTools is free software; you can redistribute it and/or modify it under
      11 # the terms of the GNU Lesser General Public License as published by the Free
      12 # Software Foundation; either version 3 of the License, or (at your option) any
      13 # later version.
      14 #
      15 # MayaChemTools is distributed in the hope that it will be useful, but without
      16 # any warranty; without even the implied warranty of merchantability of fitness
      17 # for a particular purpose.  See the GNU Lesser General Public License for more
      18 # details.
      19 #
      20 # You should have received a copy of the GNU Lesser General Public License
      21 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
      22 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
      23 # Boston, MA, 02111-1307, USA.
      24 #
      25
      26 use strict;
      27 use Exporter;
      28 use Constants;
      29 use Math::Trig ();
      30 use POSIX ();
      31
      32 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
      33
      34 @ISA = qw(Exporter);
      35 @EXPORT = qw(acos asin atan tan ceil floor log10 min max srandom random round GeneratePrimeNumbersUpToLimit GeneratePrimeNumbersUpToCount);
      36 @EXPORT_OK = qw();
      37
      38 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]
      39                );
      40
      41
      42 # Return next largest integer...
      43 sub ceil ($) {
      44   my($Value) = @_;
      45
      46   return POSIX::ceil($Value);
      47 }
      48
      49 # Return previous smallest integer...
      50 sub floor ($) {
      51   my($Value) = @_;
      52
      53   return POSIX::floor($Value);
      54 }
      55
      56 # Calculate log value using base 10...
      57 sub log10 ($) {
      58   my($Value) = @_;
      59
      60   return CORE::log($Value)/CORE::log(10);
      61 }
      62
      63 # Return the smaller of two numbers...
      64 sub min ($$) {
      65   my($Value1, $Value2) = @_;
      66
      67   return ($Value1 <= $Value2) ? $Value1 : $Value2;
      68 }
      69
      70 # Return the larger of two numbers...
      71 sub max ($$) {
      72   my($Value1, $Value2) = @_;
      73
      74   return ($Value1 >= $Value2) ? $Value1 : $Value2;
      75 }
      76
      77 # The random number generator implemented in MayaChemTools is a variant of linear
      78 # congruential generator (LCG) as described by Miller et al. [ Ref 120 ]. It is
      79 # also referred to as Lehmer random number generator or Park-Miller random number
      80 # generator.
      81 #
      82 # Unlike Perl's core random number generator function rand, the random number
      83 # generator implemented in MayaChemTools generates consistent random values
      84 # across different platforms - Windows, CygWin, Linux, Unix - for a specific random
      85 # seed.
      86 #
      87
      88 # $RandomModulus = 2**31 - 1;
      89 # $RandomMultiplier = 16807;
      90 # $RandomQuotient = $RandomModulus / $RandomMultiplier;
      91 # $RandomRemainder = $RandomModulus % $RandomMultiplier
      92 #
      93 # $MaxRandomSeed = 2*31 -2
      94 #
      95 my($MaxRandomSeed, $RandomSeed, $RandomModulus, $RandomMultiplier, $RandomQuotient, $RandomRemainder);
      96
      97 $MaxRandomSeed = 2147483646;
      98 $RandomSeed = 123456789;
      99
     100 $RandomModulus = 2147483647;
     101 $RandomMultiplier = 16807;
     102 $RandomQuotient = 127773;
     103 $RandomRemainder = 2836;
     104
     105 # Set random number seed...
     106 #
     107 # The intial value of random number seed is recommeded to be an integer between 1
     108 # and 2**31 - 2 [Ref 120] which translates to be 1 and 2147483646
     109 #
     110 sub srandom ($) {
     111   my($Seed) = @_;
     112
     113   if ($Seed <= 0 ) {
     114     die "Error: srandom: Specified seed value must be greater than 0...";
     115   }
     116
     117   $RandomSeed = ($Seed > $MaxRandomSeed) ? ($Seed % $MaxRandomSeed) : $Seed;
     118
     119   return $RandomSeed;
     120 }
     121
     122 # Retrun a random number between 0 and less than 1 or specified size...
     123 #
     124 sub random (;$) {
     125   my($Size) = @_;
     126   my($Value, $LowValue, $HighValue);
     127
     128   $Size = defined $Size ? $Size : 1.0;
     129
     130   $HighValue = $RandomSeed / $RandomQuotient;
     131   $LowValue = $RandomSeed % $RandomQuotient;
     132
     133   $Value = $RandomMultiplier * $LowValue - $RandomRemainder * $HighValue;
     134
     135   $RandomSeed = ($Value > 0) ? $Value : ($Value + $RandomModulus);
     136
     137   return ($RandomSeed / $RandomModulus) * $Size;
     138 }
     139
     140 # Round a integer/real number to:
     141 # . A nearest integer
     142 # . Specified number of decimal places
     143 #
     144 sub round ($;$) {
     145   my($Value, $DecimalPlaces) = @_;
     146   my($RoundedValue);
     147
     148   if (defined($DecimalPlaces) && $DecimalPlaces > 0) {
     149     $RoundedValue = sprintf "%.${DecimalPlaces}f", $Value;
     150   }
     151   else {
     152     if ($Value < 0) {
     153       $RoundedValue = int($Value - 0.5);
     154     }
     155     else {
     156       $RoundedValue = int($Value + 0.5);
     157     }
     158   }
     159   return $RoundedValue;
     160 }
     161
     162 # Return tangent of an angle expressed in radians.
     163 sub tan {
     164   my($Value) = @_;
     165
     166   return (CORE::sin($Value)/CORE::cos($Value));
     167 }
     168
     169 # Return inverse sine of an angle expressed in radians.
     170 #
     171 # For a right angle triangle defined by sides X and Y in a unit circle, Pythagorean theorem implies
     172 # X**2 + Y**2 = 1 and sin value corresponds to Y. So asin is equivalent to atan2(Y, sqrt(1-Y**2)).
     173 # However, taking sqrt of negative numbers is problematic; Math::Trig::asin handles it using complex
     174 # numbers.
     175 #
     176 sub asin ($) {
     177   my($Value) = @_;
     178
     179   return Math::Trig::asin($Value);
     180 }
     181
     182 # Return inverse cosine of an angle expressed in radians.
     183 #
     184 # For a right angle triangle defined by sides X and Y in a unit circle, Pythagorean theorem implies
     185 # X**2 + Y**2 = 1 and cos value corresponds to X. So asin is equivalent to atan2(sqrt(1-X**2), X)
     186 # However, taking sqrt of negative numbers is problematic; Math::Trig::acos handles it using complex
     187 # numbers.
     188 #
     189 sub acos ($) {
     190   my($Value) = @_;
     191
     192   return Math::Trig::acos($Value);
     193 }
     194
     195 # Generate prime numbers up to a specified limit and return a reference to an
     196 # array containing the prime numbers.
     197 #
     198 # By default, the first 1000 prime numbers are generated. The 1000th prime
     199 # number is 7919 and that's why default limit is set to 7920.
     200 #
     201 sub GeneratePrimeNumbersUpToLimit (;$) {
     202   my($Limit) = @_;
     203
     204   $Limit = defined $Limit ? $Limit : 7920;
     205
     206   return _GeneratePrimeNumbers('ByLimit', $Limit)
     207 }
     208
     209 # Generate prime numbers up to specified count of prime numbers and return a
     210 # reference to an array containing the prime numbers.
     211 #
     212 # By default, the first 1000 prime numbers are generated. The 1000th prime
     213 # number is 7919.
     214 #
     215 sub GeneratePrimeNumbersUpToCount (;$) {
     216   my($Count) = @_;
     217
     218   $Count = defined $Count ? $Count : 1000;
     219
     220   return _GeneratePrimeNumbers('ByCount', $Count)
     221 }
     222
     223 # Generate prime numbers up to specified limit or count and return a reference
     224 # to an array containing the prime numbers.
     225 #
     226 # The algorithm to generate prime numbers is a modification of  Sieve of Erastothenes
     227 # prime number generator.
     228 #
     229 sub _GeneratePrimeNumbers {
     230   my($Mode, $Value) = @_;
     231   my($ByLimit, $PrimeNumber, $Number, $SqrtOfNumber, $NumberIsPrime, @PrimeNumbers);
     232
     233   $ByLimit = ($Mode =~ /^ByLimit$/i) ? 1 : 0;
     234
     235   @PrimeNumbers = (2, 3);
     236   $Number = 3;
     237
     238   # while ($Number <= $Limit) {
     239   while ($ByLimit ? ($Number < $Value) : (@PrimeNumbers < $Value)) {
     240     $Number += 2;
     241     $SqrtOfNumber = sqrt $Number;
     242
     243     $NumberIsPrime = 1;
     244     PRIMENUMBER: for $PrimeNumber (@PrimeNumbers) {
     245       if ($PrimeNumber > $SqrtOfNumber) {
     246         last PRIMENUMBER;
     247       }
     248       if (!($Number % $PrimeNumber)) {
     249         $NumberIsPrime = 0;
     250         last PRIMENUMBER;
     251       }
     252     }
     253     if ($NumberIsPrime) {
     254       push @PrimeNumbers, $Number;
     255     }
     256   }
     257   return \@PrimeNumbers;
     258 }
     259