.gitignore: add bigdecimal (native gem)

This commit is contained in:
Bo Anderson 2024-05-23 14:30:09 +01:00
parent d3e8e0b47a
commit d70e5d4071
No known key found for this signature in database
8 changed files with 1 additions and 737 deletions

1
.gitignore vendored
View File

@ -71,6 +71,7 @@
# Ignore dependencies we don't wish to vendor
**/vendor/bundle/ruby/*/gems/ast-*/
**/vendor/bundle/ruby/*/gems/bigdecimal-*/
**/vendor/bundle/ruby/*/gems/bootsnap-*/
**/vendor/bundle/ruby/*/gems/bundler-*/
**/vendor/bundle/ruby/*/gems/byebug-*/

View File

@ -1,56 +0,0 @@
Ruby is copyrighted free software by Yukihiro Matsumoto <matz@netlab.jp>.
You can redistribute it and/or modify it under either the terms of the
2-clause BSDL (see the file BSDL), or the conditions below:
1. You may make and give away verbatim copies of the source form of the
software without restriction, provided that you duplicate all of the
original copyright notices and associated disclaimers.
2. You may modify your copy of the software in any way, provided that
you do at least ONE of the following:
a) place your modifications in the Public Domain or otherwise
make them Freely Available, such as by posting said
modifications to Usenet or an equivalent medium, or by allowing
the author to include your modifications in the software.
b) use the modified software only within your corporation or
organization.
c) give non-standard binaries non-standard names, with
instructions on where to get the original software distribution.
d) make other distribution arrangements with the author.
3. You may distribute the software in object code or binary form,
provided that you do at least ONE of the following:
a) distribute the binaries and library files of the software,
together with instructions (in the manual page or equivalent)
on where to get the original distribution.
b) accompany the distribution with the machine-readable source of
the software.
c) give non-standard binaries non-standard names, with
instructions on where to get the original software distribution.
d) make other distribution arrangements with the author.
4. You may modify and include the part of the software into any other
software (possibly commercial). But some files in the distribution
are not written by the author, so that they are not under these terms.
For the list of those files and their copying conditions, see the
file LEGAL.
5. The scripts and library files supplied as input to or produced as
output from the software do not automatically fall under the
copyright of the software, but belong to whomever generated them,
and may be sold commercially, and may be aggregated with this
software.
6. THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE.

View File

@ -1,5 +0,0 @@
if RUBY_ENGINE == 'jruby'
JRuby::Util.load_ext("org.jruby.ext.bigdecimal.BigDecimalLibrary")
else
require 'bigdecimal.so'
end

View File

@ -1,90 +0,0 @@
# frozen_string_literal: false
require 'bigdecimal'
# require 'bigdecimal/jacobian'
#
# Provides methods to compute the Jacobian matrix of a set of equations at a
# point x. In the methods below:
#
# f is an Object which is used to compute the Jacobian matrix of the equations.
# It must provide the following methods:
#
# f.values(x):: returns the values of all functions at x
#
# f.zero:: returns 0.0
# f.one:: returns 1.0
# f.two:: returns 2.0
# f.ten:: returns 10.0
#
# f.eps:: returns the convergence criterion (epsilon value) used to determine whether two values are considered equal. If |a-b| < epsilon, the two values are considered equal.
#
# x is the point at which to compute the Jacobian.
#
# fx is f.values(x).
#
module Jacobian
module_function
# Determines the equality of two numbers by comparing to zero, or using the epsilon value
def isEqual(a,b,zero=0.0,e=1.0e-8)
aa = a.abs
bb = b.abs
if aa == zero && bb == zero then
true
else
if ((a-b)/(aa+bb)).abs < e then
true
else
false
end
end
end
# Computes the derivative of +f[i]+ at +x[i]+.
# +fx+ is the value of +f+ at +x+.
def dfdxi(f,fx,x,i)
nRetry = 0
n = x.size
xSave = x[i]
ok = 0
ratio = f.ten*f.ten*f.ten
dx = x[i].abs/ratio
dx = fx[i].abs/ratio if isEqual(dx,f.zero,f.zero,f.eps)
dx = f.one/f.ten if isEqual(dx,f.zero,f.zero,f.eps)
until ok>0 do
deriv = []
nRetry += 1
if nRetry > 100
raise "Singular Jacobian matrix. No change at x[" + i.to_s + "]"
end
dx = dx*f.two
x[i] += dx
fxNew = f.values(x)
for j in 0...n do
if !isEqual(fxNew[j],fx[j],f.zero,f.eps) then
ok += 1
deriv <<= (fxNew[j]-fx[j])/dx
else
deriv <<= f.zero
end
end
x[i] = xSave
end
deriv
end
# Computes the Jacobian of +f+ at +x+. +fx+ is the value of +f+ at +x+.
def jacobian(f,fx,x)
n = x.size
dfdx = Array.new(n*n)
for i in 0...n do
df = dfdxi(f,fx,x,i)
for j in 0...n do
dfdx[j*n+i] = df[j]
end
end
dfdx
end
end

View File

@ -1,89 +0,0 @@
# frozen_string_literal: false
require 'bigdecimal'
#
# Solves a*x = b for x, using LU decomposition.
#
module LUSolve
module_function
# Performs LU decomposition of the n by n matrix a.
def ludecomp(a,n,zero=0,one=1)
prec = BigDecimal.limit(nil)
ps = []
scales = []
for i in 0...n do # pick up largest(abs. val.) element in each row.
ps <<= i
nrmrow = zero
ixn = i*n
for j in 0...n do
biggst = a[ixn+j].abs
nrmrow = biggst if biggst>nrmrow
end
if nrmrow>zero then
scales <<= one.div(nrmrow,prec)
else
raise "Singular matrix"
end
end
n1 = n - 1
for k in 0...n1 do # Gaussian elimination with partial pivoting.
biggst = zero;
for i in k...n do
size = a[ps[i]*n+k].abs*scales[ps[i]]
if size>biggst then
biggst = size
pividx = i
end
end
raise "Singular matrix" if biggst<=zero
if pividx!=k then
j = ps[k]
ps[k] = ps[pividx]
ps[pividx] = j
end
pivot = a[ps[k]*n+k]
for i in (k+1)...n do
psin = ps[i]*n
a[psin+k] = mult = a[psin+k].div(pivot,prec)
if mult!=zero then
pskn = ps[k]*n
for j in (k+1)...n do
a[psin+j] -= mult.mult(a[pskn+j],prec)
end
end
end
end
raise "Singular matrix" if a[ps[n1]*n+n1] == zero
ps
end
# Solves a*x = b for x, using LU decomposition.
#
# a is a matrix, b is a constant vector, x is the solution vector.
#
# ps is the pivot, a vector which indicates the permutation of rows performed
# during LU decomposition.
def lusolve(a,b,ps,zero=0.0)
prec = BigDecimal.limit(nil)
n = ps.size
x = []
for i in 0...n do
dot = zero
psin = ps[i]*n
for j in 0...i do
dot = a[psin+j].mult(x[j],prec) + dot
end
x <<= b[ps[i]] - dot
end
(n-1).downto(0) do |i|
dot = zero
psin = ps[i]*n
for j in (i+1)...n do
dot = a[psin+j].mult(x[j],prec) + dot
end
x[i] = (x[i]-dot).div(a[psin+i],prec)
end
x
end
end

View File

@ -1,232 +0,0 @@
# frozen_string_literal: false
require 'bigdecimal'
#
#--
# Contents:
# sqrt(x, prec)
# sin (x, prec)
# cos (x, prec)
# atan(x, prec) Note: |x|<1, x=0.9999 may not converge.
# PI (prec)
# E (prec) == exp(1.0,prec)
#
# where:
# x ... BigDecimal number to be computed.
# |x| must be small enough to get convergence.
# prec ... Number of digits to be obtained.
#++
#
# Provides mathematical functions.
#
# Example:
#
# require "bigdecimal/math"
#
# include BigMath
#
# a = BigDecimal((PI(100)/2).to_s)
# puts sin(a,100) # => 0.99999999999999999999......e0
#
module BigMath
module_function
# call-seq:
# sqrt(decimal, numeric) -> BigDecimal
#
# Computes the square root of +decimal+ to the specified number of digits of
# precision, +numeric+.
#
# BigMath.sqrt(BigDecimal('2'), 16).to_s
# #=> "0.1414213562373095048801688724e1"
#
def sqrt(x, prec)
x.sqrt(prec)
end
# call-seq:
# sin(decimal, numeric) -> BigDecimal
#
# Computes the sine of +decimal+ to the specified number of digits of
# precision, +numeric+.
#
# If +decimal+ is Infinity or NaN, returns NaN.
#
# BigMath.sin(BigMath.PI(5)/4, 5).to_s
# #=> "0.70710678118654752440082036563292800375e0"
#
def sin(x, prec)
raise ArgumentError, "Zero or negative precision for sin" if prec <= 0
return BigDecimal("NaN") if x.infinite? || x.nan?
n = prec + BigDecimal.double_fig
one = BigDecimal("1")
two = BigDecimal("2")
x = -x if neg = x < 0
if x > (twopi = two * BigMath.PI(prec))
if x > 30
x %= twopi
else
x -= twopi while x > twopi
end
end
x1 = x
x2 = x.mult(x,n)
sign = 1
y = x
d = y
i = one
z = one
while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
m = BigDecimal.double_fig if m < BigDecimal.double_fig
sign = -sign
x1 = x2.mult(x1,n)
i += two
z *= (i-one) * i
d = sign * x1.div(z,m)
y += d
end
neg ? -y : y
end
# call-seq:
# cos(decimal, numeric) -> BigDecimal
#
# Computes the cosine of +decimal+ to the specified number of digits of
# precision, +numeric+.
#
# If +decimal+ is Infinity or NaN, returns NaN.
#
# BigMath.cos(BigMath.PI(4), 16).to_s
# #=> "-0.999999999999999999999999999999856613163740061349e0"
#
def cos(x, prec)
raise ArgumentError, "Zero or negative precision for cos" if prec <= 0
return BigDecimal("NaN") if x.infinite? || x.nan?
n = prec + BigDecimal.double_fig
one = BigDecimal("1")
two = BigDecimal("2")
x = -x if x < 0
if x > (twopi = two * BigMath.PI(prec))
if x > 30
x %= twopi
else
x -= twopi while x > twopi
end
end
x1 = one
x2 = x.mult(x,n)
sign = 1
y = one
d = y
i = BigDecimal("0")
z = one
while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
m = BigDecimal.double_fig if m < BigDecimal.double_fig
sign = -sign
x1 = x2.mult(x1,n)
i += two
z *= (i-one) * i
d = sign * x1.div(z,m)
y += d
end
y
end
# call-seq:
# atan(decimal, numeric) -> BigDecimal
#
# Computes the arctangent of +decimal+ to the specified number of digits of
# precision, +numeric+.
#
# If +decimal+ is NaN, returns NaN.
#
# BigMath.atan(BigDecimal('-1'), 16).to_s
# #=> "-0.785398163397448309615660845819878471907514682065e0"
#
def atan(x, prec)
raise ArgumentError, "Zero or negative precision for atan" if prec <= 0
return BigDecimal("NaN") if x.nan?
pi = PI(prec)
x = -x if neg = x < 0
return pi.div(neg ? -2 : 2, prec) if x.infinite?
return pi / (neg ? -4 : 4) if x.round(prec) == 1
x = BigDecimal("1").div(x, prec) if inv = x > 1
x = (-1 + sqrt(1 + x**2, prec))/x if dbl = x > 0.5
n = prec + BigDecimal.double_fig
y = x
d = y
t = x
r = BigDecimal("3")
x2 = x.mult(x,n)
while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
m = BigDecimal.double_fig if m < BigDecimal.double_fig
t = -t.mult(x2,n)
d = t.div(r,m)
y += d
r += 2
end
y *= 2 if dbl
y = pi / 2 - y if inv
y = -y if neg
y
end
# call-seq:
# PI(numeric) -> BigDecimal
#
# Computes the value of pi to the specified number of digits of precision,
# +numeric+.
#
# BigMath.PI(10).to_s
# #=> "0.3141592653589793238462643388813853786957412e1"
#
def PI(prec)
raise ArgumentError, "Zero or negative precision for PI" if prec <= 0
n = prec + BigDecimal.double_fig
zero = BigDecimal("0")
one = BigDecimal("1")
two = BigDecimal("2")
m25 = BigDecimal("-0.04")
m57121 = BigDecimal("-57121")
pi = zero
d = one
k = one
t = BigDecimal("-80")
while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0)
m = BigDecimal.double_fig if m < BigDecimal.double_fig
t = t*m25
d = t.div(k,m)
k = k+two
pi = pi + d
end
d = one
k = one
t = BigDecimal("956")
while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0)
m = BigDecimal.double_fig if m < BigDecimal.double_fig
t = t.div(m57121,n)
d = t.div(k,m)
pi = pi + d
k = k+two
end
pi
end
# call-seq:
# E(numeric) -> BigDecimal
#
# Computes e (the base of natural logarithms) to the specified number of
# digits of precision, +numeric+.
#
# BigMath.E(10).to_s
# #=> "0.271828182845904523536028752390026306410273e1"
#
def E(prec)
raise ArgumentError, "Zero or negative precision for E" if prec <= 0
BigMath.exp(1, prec)
end
end

View File

@ -1,80 +0,0 @@
# frozen_string_literal: false
require "bigdecimal/ludcmp"
require "bigdecimal/jacobian"
#
# newton.rb
#
# Solves the nonlinear algebraic equation system f = 0 by Newton's method.
# This program is not dependent on BigDecimal.
#
# To call:
# n = nlsolve(f,x)
# where n is the number of iterations required,
# x is the initial value vector
# f is an Object which is used to compute the values of the equations to be solved.
# It must provide the following methods:
#
# f.values(x):: returns the values of all functions at x
#
# f.zero:: returns 0.0
# f.one:: returns 1.0
# f.two:: returns 2.0
# f.ten:: returns 10.0
#
# f.eps:: returns the convergence criterion (epsilon value) used to determine whether two values are considered equal. If |a-b| < epsilon, the two values are considered equal.
#
# On exit, x is the solution vector.
#
module Newton
include LUSolve
include Jacobian
module_function
def norm(fv,zero=0.0) # :nodoc:
s = zero
n = fv.size
for i in 0...n do
s += fv[i]*fv[i]
end
s
end
# See also Newton
def nlsolve(f,x)
nRetry = 0
n = x.size
f0 = f.values(x)
zero = f.zero
one = f.one
two = f.two
p5 = one/two
d = norm(f0,zero)
minfact = f.ten*f.ten*f.ten
minfact = one/minfact
e = f.eps
while d >= e do
nRetry += 1
# Not yet converged. => Compute Jacobian matrix
dfdx = jacobian(f,f0,x)
# Solve dfdx*dx = -f0 to estimate dx
dx = lusolve(dfdx,f0,ludecomp(dfdx,n,zero,one),zero)
fact = two
xs = x.dup
begin
fact *= p5
if fact < minfact then
raise "Failed to reduce function values."
end
for i in 0...n do
x[i] = xs[i] - dx[i]*fact
end
f0 = f.values(x)
dn = norm(f0,zero)
end while(dn>=d)
d = dn
end
nRetry
end
end

View File

@ -1,185 +0,0 @@
# frozen_string_literal: false
#
#--
# bigdecimal/util extends various native classes to provide the #to_d method,
# and provides BigDecimal#to_d and BigDecimal#to_digits.
#++
require 'bigdecimal'
class Integer < Numeric
# call-seq:
# int.to_d -> bigdecimal
#
# Returns the value of +int+ as a BigDecimal.
#
# require 'bigdecimal'
# require 'bigdecimal/util'
#
# 42.to_d # => 0.42e2
#
# See also Kernel.BigDecimal.
#
def to_d
BigDecimal(self)
end
end
class Float < Numeric
# call-seq:
# float.to_d -> bigdecimal
# float.to_d(precision) -> bigdecimal
#
# Returns the value of +float+ as a BigDecimal.
# The +precision+ parameter is used to determine the number of
# significant digits for the result. When +precision+ is set to +0+,
# the number of digits to represent the float being converted is determined
# automatically.
# The default +precision+ is +0+.
#
# require 'bigdecimal'
# require 'bigdecimal/util'
#
# 0.5.to_d # => 0.5e0
# 1.234.to_d # => 0.1234e1
# 1.234.to_d(2) # => 0.12e1
#
# See also Kernel.BigDecimal.
#
def to_d(precision=0)
BigDecimal(self, precision)
end
end
class String
# call-seq:
# str.to_d -> bigdecimal
#
# Returns the result of interpreting leading characters in +str+
# as a BigDecimal.
#
# require 'bigdecimal'
# require 'bigdecimal/util'
#
# "0.5".to_d # => 0.5e0
# "123.45e1".to_d # => 0.12345e4
# "45.67 degrees".to_d # => 0.4567e2
#
# See also Kernel.BigDecimal.
#
def to_d
BigDecimal.interpret_loosely(self)
end
end
class BigDecimal < Numeric
# call-seq:
# a.to_digits -> string
#
# Converts a BigDecimal to a String of the form "nnnnnn.mmm".
# This method is deprecated; use BigDecimal#to_s("F") instead.
#
# require 'bigdecimal/util'
#
# d = BigDecimal("3.14")
# d.to_digits # => "3.14"
#
def to_digits
if self.nan? || self.infinite? || self.zero?
self.to_s
else
i = self.to_i.to_s
_,f,_,z = self.frac.split
i + "." + ("0"*(-z)) + f
end
end
# call-seq:
# a.to_d -> bigdecimal
#
# Returns self.
#
# require 'bigdecimal/util'
#
# d = BigDecimal("3.14")
# d.to_d # => 0.314e1
#
def to_d
self
end
end
class Rational < Numeric
# call-seq:
# rat.to_d(precision) -> bigdecimal
#
# Returns the value as a BigDecimal.
#
# The required +precision+ parameter is used to determine the number of
# significant digits for the result.
#
# require 'bigdecimal'
# require 'bigdecimal/util'
#
# Rational(22, 7).to_d(3) # => 0.314e1
#
# See also Kernel.BigDecimal.
#
def to_d(precision)
BigDecimal(self, precision)
end
end
class Complex < Numeric
# call-seq:
# cmp.to_d -> bigdecimal
# cmp.to_d(precision) -> bigdecimal
#
# Returns the value as a BigDecimal.
#
# The +precision+ parameter is required for a rational complex number.
# This parameter is used to determine the number of significant digits
# for the result.
#
# require 'bigdecimal'
# require 'bigdecimal/util'
#
# Complex(0.1234567, 0).to_d(4) # => 0.1235e0
# Complex(Rational(22, 7), 0).to_d(3) # => 0.314e1
#
# See also Kernel.BigDecimal.
#
def to_d(*args)
BigDecimal(self) unless self.imag.zero? # to raise eerror
if args.length == 0
case self.real
when Rational
BigDecimal(self.real) # to raise error
end
end
self.real.to_d(*args)
end
end
class NilClass
# call-seq:
# nil.to_d -> bigdecimal
#
# Returns nil represented as a BigDecimal.
#
# require 'bigdecimal'
# require 'bigdecimal/util'
#
# nil.to_d # => 0.0
#
def to_d
BigDecimal(0)
end
end