1 package MathUtil;
   2 #
   3 # File: MathUtil.pm
   4 # Author: Manish Sud <msud@san.rr.com>
   5 #
   6 # Copyright (C) 2024 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