This commit is contained in:
@@ -0,0 +1,56 @@
|
||||
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.
|
||||
@@ -0,0 +1,57 @@
|
||||
# coding: utf-8
|
||||
|
||||
name = File.basename(__FILE__, '.*')
|
||||
source_version = ["", "ext/#{name}/"].find do |dir|
|
||||
begin
|
||||
break File.foreach(File.join(__dir__, "#{dir}#{name}.c")) {|line|
|
||||
break $1.sub("-", ".") if /^#define\s+#{name.upcase}_VERSION\s+"(.+)"/o =~ line
|
||||
}
|
||||
rescue Errno::ENOENT
|
||||
end
|
||||
end or raise "can't find #{name.upcase}_VERSION"
|
||||
|
||||
Gem::Specification.new do |s|
|
||||
s.name = name
|
||||
s.version = source_version
|
||||
s.authors = ["Kenta Murata", "Zachary Scott", "Shigeo Kobayashi"]
|
||||
s.email = ["mrkn@mrkn.jp"]
|
||||
|
||||
s.summary = "Arbitrary-precision decimal floating-point number library."
|
||||
s.description = "This library provides arbitrary-precision decimal floating-point number class."
|
||||
s.homepage = "https://github.com/ruby/bigdecimal"
|
||||
s.licenses = ["Ruby", "BSD-2-Clause"]
|
||||
|
||||
s.require_paths = %w[lib]
|
||||
s.files = %w[
|
||||
LICENSE
|
||||
bigdecimal.gemspec
|
||||
lib/bigdecimal.rb
|
||||
lib/bigdecimal/jacobian.rb
|
||||
lib/bigdecimal/ludcmp.rb
|
||||
lib/bigdecimal/math.rb
|
||||
lib/bigdecimal/newton.rb
|
||||
lib/bigdecimal/util.rb
|
||||
sample/linear.rb
|
||||
sample/nlsolve.rb
|
||||
sample/pi.rb
|
||||
]
|
||||
if Gem::Platform === s.platform and s.platform =~ 'java' or RUBY_ENGINE == 'jruby'
|
||||
s.platform = 'java'
|
||||
else
|
||||
s.extensions = %w[ext/bigdecimal/extconf.rb]
|
||||
s.files += %w[
|
||||
ext/bigdecimal/bigdecimal.c
|
||||
ext/bigdecimal/bigdecimal.h
|
||||
ext/bigdecimal/bits.h
|
||||
ext/bigdecimal/feature.h
|
||||
ext/bigdecimal/missing.c
|
||||
ext/bigdecimal/missing.h
|
||||
ext/bigdecimal/missing/dtoa.c
|
||||
ext/bigdecimal/static_assert.h
|
||||
]
|
||||
end
|
||||
|
||||
s.required_ruby_version = Gem::Requirement.new(">= 2.5.0")
|
||||
|
||||
s.metadata["changelog_uri"] = s.homepage + "/blob/master/CHANGES.md"
|
||||
end
|
||||
@@ -0,0 +1,270 @@
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# V=0 quiet, V=1 verbose. other values don't work.
|
||||
V = 0
|
||||
V0 = $(V:0=)
|
||||
Q1 = $(V:1=)
|
||||
Q = $(Q1:0=@)
|
||||
ECHO1 = $(V:1=@ :)
|
||||
ECHO = $(ECHO1:0=@ echo)
|
||||
NULLCMD = :
|
||||
|
||||
#### Start of system configuration section. ####
|
||||
|
||||
srcdir = .
|
||||
topdir = /usr/include/ruby-3.2.0
|
||||
hdrdir = $(topdir)
|
||||
arch_hdrdir = /usr/include/x86_64-linux-gnu/ruby-3.2.0
|
||||
PATH_SEPARATOR = :
|
||||
VPATH = $(srcdir):$(arch_hdrdir)/ruby:$(hdrdir)/ruby
|
||||
prefix = $(DESTDIR)/usr
|
||||
rubysitearchprefix = $(sitearchlibdir)/$(RUBY_BASE_NAME)
|
||||
rubyarchprefix = $(archlibdir)/$(RUBY_BASE_NAME)
|
||||
rubylibprefix = $(libdir)/$(RUBY_BASE_NAME)
|
||||
exec_prefix = $(prefix)
|
||||
vendorarchhdrdir = $(sitearchincludedir)/$(RUBY_VERSION_NAME)/vendor_ruby
|
||||
sitearchhdrdir = $(sitearchincludedir)/$(RUBY_VERSION_NAME)/site_ruby
|
||||
rubyarchhdrdir = $(archincludedir)/$(RUBY_VERSION_NAME)
|
||||
vendorhdrdir = $(rubyhdrdir)/vendor_ruby
|
||||
sitehdrdir = $(rubyhdrdir)/site_ruby
|
||||
rubyhdrdir = $(includedir)/$(RUBY_VERSION_NAME)
|
||||
vendorarchdir = $(rubysitearchprefix)/vendor_ruby/$(ruby_version)
|
||||
vendorlibdir = $(vendordir)/$(ruby_version)
|
||||
vendordir = $(rubylibprefix)/vendor_ruby
|
||||
sitearchdir = $(DESTDIR)/usr/local/lib/x86_64-linux-gnu/site_ruby
|
||||
sitelibdir = $(sitedir)/$(ruby_version)
|
||||
sitedir = $(DESTDIR)/usr/local/lib/site_ruby
|
||||
rubyarchdir = $(rubyarchprefix)/$(ruby_version)
|
||||
rubylibdir = $(rubylibprefix)/$(ruby_version)
|
||||
sitearchincludedir = $(includedir)/$(sitearch)
|
||||
archincludedir = $(includedir)/$(arch)
|
||||
sitearchlibdir = $(libdir)/$(sitearch)
|
||||
archlibdir = $(libdir)/$(arch)
|
||||
ridir = $(datarootdir)/$(RI_BASE_NAME)
|
||||
mandir = $(datarootdir)/man
|
||||
localedir = $(datarootdir)/locale
|
||||
libdir = $(exec_prefix)/lib
|
||||
psdir = $(docdir)
|
||||
pdfdir = $(docdir)
|
||||
dvidir = $(docdir)
|
||||
htmldir = $(docdir)
|
||||
infodir = $(datarootdir)/info
|
||||
docdir = $(datarootdir)/doc/$(PACKAGE)
|
||||
oldincludedir = $(DESTDIR)/usr/include
|
||||
includedir = $(prefix)/include
|
||||
runstatedir = $(DESTDIR)/var/run
|
||||
localstatedir = $(DESTDIR)/var
|
||||
sharedstatedir = $(prefix)/com
|
||||
sysconfdir = $(DESTDIR)/etc
|
||||
datadir = $(datarootdir)
|
||||
datarootdir = $(prefix)/share
|
||||
libexecdir = $(exec_prefix)/libexec
|
||||
sbindir = $(exec_prefix)/sbin
|
||||
bindir = $(exec_prefix)/bin
|
||||
archdir = $(rubyarchdir)
|
||||
|
||||
|
||||
CC_WRAPPER =
|
||||
CC = x86_64-linux-gnu-gcc
|
||||
CXX = x86_64-linux-gnu-g++
|
||||
LIBRUBY = $(LIBRUBY_SO)
|
||||
LIBRUBY_A = lib$(RUBY_SO_NAME)-static.a
|
||||
LIBRUBYARG_SHARED = -l$(RUBY_SO_NAME)
|
||||
LIBRUBYARG_STATIC = -l$(RUBY_SO_NAME)-static $(MAINLIBS)
|
||||
empty =
|
||||
OUTFLAG = -o $(empty)
|
||||
COUTFLAG = -o $(empty)
|
||||
CSRCFLAG = $(empty)
|
||||
|
||||
RUBY_EXTCONF_H =
|
||||
cflags = $(optflags) $(debugflags) $(warnflags)
|
||||
cxxflags =
|
||||
optflags = -O3 -fno-fast-math
|
||||
debugflags = -ggdb3
|
||||
warnflags = -Wall -Wextra -Wdeprecated-declarations -Wdiv-by-zero -Wduplicated-cond -Wimplicit-function-declaration -Wimplicit-int -Wmisleading-indentation -Wpointer-arith -Wwrite-strings -Wold-style-definition -Wimplicit-fallthrough=0 -Wmissing-noreturn -Wno-cast-function-type -Wno-constant-logical-operand -Wno-long-long -Wno-missing-field-initializers -Wno-overlength-strings -Wno-packed-bitfield-compat -Wno-parentheses-equality -Wno-self-assign -Wno-tautological-compare -Wno-unused-parameter -Wno-unused-value -Wsuggest-attribute=format -Wsuggest-attribute=noreturn -Wunused-variable -Wundef
|
||||
cppflags =
|
||||
CCDLFLAGS = -fPIC
|
||||
CFLAGS = $(CCDLFLAGS) -g -O2 -fno-omit-frame-pointer -mno-omit-leaf-frame-pointer -ffile-prefix-map=BUILDDIR=. -fstack-protector-strong -fstack-clash-protection -Wformat -Werror=format-security -fcf-protection -fdebug-prefix-map=BUILDDIR=/usr/src/ruby3.2-3.2.3-1ubuntu0.24.04.8 -fPIC $(ARCH_FLAG)
|
||||
INCFLAGS = -I. -I$(arch_hdrdir) -I$(hdrdir)/ruby/backward -I$(hdrdir) -I$(srcdir)
|
||||
DEFS =
|
||||
CPPFLAGS = -DHAVE_BUILTIN___BUILTIN_CLZ -DHAVE_BUILTIN___BUILTIN_CLZL -DHAVE_BUILTIN___BUILTIN_CLZLL -DHAVE_FLOAT_H -DHAVE_MATH_H -DHAVE_STDBOOL_H -DHAVE_STDLIB_H -DHAVE_X86INTRIN_H -DHAVE_LABS -DHAVE_LLABS -DHAVE_FINITE -DHAVE_RUBY_ATOMIC_H -DHAVE_RUBY_INTERNAL_HAS_BUILTIN_H -DHAVE_RUBY_INTERNAL_STATIC_ASSERT_H -DHAVE_RB_RATIONAL_NUM -DHAVE_RB_RATIONAL_DEN -DHAVE_RB_COMPLEX_REAL -DHAVE_RB_COMPLEX_IMAG -DHAVE_RB_OPTS_EXCEPTION_P -DHAVE_RB_CATEGORY_WARN -DHAVE_CONST_RB_WARN_CATEGORY_DEPRECATED -Wdate-time -D_FORTIFY_SOURCE=3 $(DEFS) $(cppflags)
|
||||
CXXFLAGS = $(CCDLFLAGS) -g -O2 -fno-omit-frame-pointer -mno-omit-leaf-frame-pointer -ffile-prefix-map=BUILDDIR=. -fstack-protector-strong -fstack-clash-protection -Wformat -Werror=format-security -fcf-protection -fdebug-prefix-map=BUILDDIR=/usr/src/ruby3.2-3.2.3-1ubuntu0.24.04.8 $(ARCH_FLAG)
|
||||
ldflags = -L. -Wl,-Bsymbolic-functions -Wl,-z,relro -Wl,-z,now -fstack-protector-strong -rdynamic -Wl,-export-dynamic -Wl,--no-as-needed
|
||||
dldflags = -Wl,-Bsymbolic-functions -Wl,-z,relro -Wl,-z,now
|
||||
ARCH_FLAG =
|
||||
DLDFLAGS = $(ldflags) $(dldflags) $(ARCH_FLAG)
|
||||
LDSHARED = $(CC) -shared
|
||||
LDSHAREDXX = $(CXX) -shared
|
||||
AR = x86_64-linux-gnu-gcc-ar
|
||||
EXEEXT =
|
||||
|
||||
RUBY_INSTALL_NAME = $(RUBY_BASE_NAME)3.2
|
||||
RUBY_SO_NAME = ruby-3.2
|
||||
RUBYW_INSTALL_NAME =
|
||||
RUBY_VERSION_NAME = $(RUBY_BASE_NAME)-$(ruby_version)
|
||||
RUBYW_BASE_NAME = rubyw
|
||||
RUBY_BASE_NAME = ruby
|
||||
|
||||
arch = x86_64-linux-gnu
|
||||
sitearch = $(arch)
|
||||
ruby_version = 3.2.0
|
||||
ruby = $(bindir)/$(RUBY_BASE_NAME)3.2
|
||||
RUBY = $(ruby)
|
||||
BUILTRUBY = $(bindir)/$(RUBY_BASE_NAME)3.2
|
||||
ruby_headers = $(hdrdir)/ruby.h $(hdrdir)/ruby/backward.h $(hdrdir)/ruby/ruby.h $(hdrdir)/ruby/defines.h $(hdrdir)/ruby/missing.h $(hdrdir)/ruby/intern.h $(hdrdir)/ruby/st.h $(hdrdir)/ruby/subst.h $(arch_hdrdir)/ruby/config.h
|
||||
|
||||
RM = rm -f
|
||||
RM_RF = rm -fr
|
||||
RMDIRS = rmdir --ignore-fail-on-non-empty -p
|
||||
MAKEDIRS = /bin/mkdir -p
|
||||
INSTALL = /usr/bin/install -c
|
||||
INSTALL_PROG = $(INSTALL) -m 0755
|
||||
INSTALL_DATA = $(INSTALL) -m 644
|
||||
COPY = cp
|
||||
TOUCH = exit >
|
||||
|
||||
#### End of system configuration section. ####
|
||||
|
||||
preload =
|
||||
libpath = . $(archlibdir)
|
||||
LIBPATH = -L. -L$(archlibdir)
|
||||
DEFFILE =
|
||||
|
||||
CLEANFILES = mkmf.log
|
||||
DISTCLEANFILES =
|
||||
DISTCLEANDIRS =
|
||||
|
||||
extout =
|
||||
extout_prefix =
|
||||
target_prefix =
|
||||
LOCAL_LIBS =
|
||||
LIBS = $(LIBRUBYARG_SHARED) -lm -lpthread -lc
|
||||
ORIG_SRCS = bigdecimal.c missing.c
|
||||
SRCS = $(ORIG_SRCS)
|
||||
OBJS = bigdecimal.o missing.o
|
||||
HDRS = $(srcdir)/bigdecimal.h $(srcdir)/bits.h $(srcdir)/feature.h $(srcdir)/missing.h $(srcdir)/static_assert.h
|
||||
LOCAL_HDRS =
|
||||
TARGET = bigdecimal
|
||||
TARGET_NAME = bigdecimal
|
||||
TARGET_ENTRY = Init_$(TARGET_NAME)
|
||||
DLLIB = $(TARGET).so
|
||||
EXTSTATIC =
|
||||
STATIC_LIB =
|
||||
|
||||
TIMESTAMP_DIR = .
|
||||
BINDIR = $(bindir)
|
||||
RUBYCOMMONDIR = $(sitedir)$(target_prefix)
|
||||
RUBYLIBDIR = $(sitelibdir)$(target_prefix)
|
||||
RUBYARCHDIR = $(sitearchdir)$(target_prefix)
|
||||
HDRDIR = $(sitehdrdir)$(target_prefix)
|
||||
ARCHHDRDIR = $(sitearchhdrdir)$(target_prefix)
|
||||
TARGET_SO_DIR =
|
||||
TARGET_SO = $(TARGET_SO_DIR)$(DLLIB)
|
||||
CLEANLIBS = $(TARGET_SO) false
|
||||
CLEANOBJS = $(OBJS) *.bak
|
||||
TARGET_SO_DIR_TIMESTAMP = $(TIMESTAMP_DIR)/.sitearchdir.time
|
||||
BIGDECIMAL_RB = $(srcdir)/../../lib/bigdecimal.rb
|
||||
|
||||
all: $(DLLIB)
|
||||
static: $(STATIC_LIB)
|
||||
.PHONY: all install static install-so install-rb
|
||||
.PHONY: clean clean-so clean-static clean-rb
|
||||
|
||||
clean-static::
|
||||
clean-rb-default::
|
||||
clean-rb::
|
||||
clean-so::
|
||||
clean: clean-so clean-static clean-rb-default clean-rb
|
||||
-$(Q)$(RM_RF) $(CLEANLIBS) $(CLEANOBJS) $(CLEANFILES) .*.time
|
||||
|
||||
distclean-rb-default::
|
||||
distclean-rb::
|
||||
distclean-so::
|
||||
distclean-static::
|
||||
distclean: clean distclean-so distclean-static distclean-rb-default distclean-rb
|
||||
-$(Q)$(RM) Makefile $(RUBY_EXTCONF_H) conftest.* mkmf.log
|
||||
-$(Q)$(RM) core ruby$(EXEEXT) *~ $(DISTCLEANFILES)
|
||||
-$(Q)$(RMDIRS) $(DISTCLEANDIRS) 2> /dev/null || true
|
||||
|
||||
realclean: distclean
|
||||
install: install-so install-rb
|
||||
|
||||
install-so: $(DLLIB) $(TARGET_SO_DIR_TIMESTAMP)
|
||||
$(INSTALL_PROG) $(DLLIB) $(RUBYARCHDIR)
|
||||
clean-static::
|
||||
-$(Q)$(RM) $(STATIC_LIB)
|
||||
install-rb: pre-install-rb do-install-rb install-rb-default
|
||||
install-rb-default: pre-install-rb-default do-install-rb-default
|
||||
pre-install-rb: Makefile
|
||||
pre-install-rb-default: Makefile
|
||||
do-install-rb:
|
||||
do-install-rb-default:
|
||||
pre-install-rb-default:
|
||||
@$(NULLCMD)
|
||||
$(TARGET_SO_DIR_TIMESTAMP):
|
||||
$(Q) $(MAKEDIRS) $(@D) $(RUBYARCHDIR)
|
||||
$(Q) $(TOUCH) $@
|
||||
|
||||
site-install: site-install-so site-install-rb
|
||||
site-install-so: install-so
|
||||
site-install-rb: install-rb
|
||||
|
||||
.SUFFIXES: .c .m .cc .mm .cxx .cpp .o .S
|
||||
|
||||
.cc.o:
|
||||
$(ECHO) compiling $(<)
|
||||
$(Q) $(CXX) $(INCFLAGS) $(CPPFLAGS) $(CXXFLAGS) $(COUTFLAG)$@ -c $(CSRCFLAG)$<
|
||||
|
||||
.cc.S:
|
||||
$(ECHO) translating $(<)
|
||||
$(Q) $(CXX) $(INCFLAGS) $(CPPFLAGS) $(CXXFLAGS) $(COUTFLAG)$@ -S $(CSRCFLAG)$<
|
||||
|
||||
.mm.o:
|
||||
$(ECHO) compiling $(<)
|
||||
$(Q) $(CXX) $(INCFLAGS) $(CPPFLAGS) $(CXXFLAGS) $(COUTFLAG)$@ -c $(CSRCFLAG)$<
|
||||
|
||||
.mm.S:
|
||||
$(ECHO) translating $(<)
|
||||
$(Q) $(CXX) $(INCFLAGS) $(CPPFLAGS) $(CXXFLAGS) $(COUTFLAG)$@ -S $(CSRCFLAG)$<
|
||||
|
||||
.cxx.o:
|
||||
$(ECHO) compiling $(<)
|
||||
$(Q) $(CXX) $(INCFLAGS) $(CPPFLAGS) $(CXXFLAGS) $(COUTFLAG)$@ -c $(CSRCFLAG)$<
|
||||
|
||||
.cxx.S:
|
||||
$(ECHO) translating $(<)
|
||||
$(Q) $(CXX) $(INCFLAGS) $(CPPFLAGS) $(CXXFLAGS) $(COUTFLAG)$@ -S $(CSRCFLAG)$<
|
||||
|
||||
.cpp.o:
|
||||
$(ECHO) compiling $(<)
|
||||
$(Q) $(CXX) $(INCFLAGS) $(CPPFLAGS) $(CXXFLAGS) $(COUTFLAG)$@ -c $(CSRCFLAG)$<
|
||||
|
||||
.cpp.S:
|
||||
$(ECHO) translating $(<)
|
||||
$(Q) $(CXX) $(INCFLAGS) $(CPPFLAGS) $(CXXFLAGS) $(COUTFLAG)$@ -S $(CSRCFLAG)$<
|
||||
|
||||
.c.o:
|
||||
$(ECHO) compiling $(<)
|
||||
$(Q) $(CC) $(INCFLAGS) $(CPPFLAGS) $(CFLAGS) $(COUTFLAG)$@ -c $(CSRCFLAG)$<
|
||||
|
||||
.c.S:
|
||||
$(ECHO) translating $(<)
|
||||
$(Q) $(CC) $(INCFLAGS) $(CPPFLAGS) $(CFLAGS) $(COUTFLAG)$@ -S $(CSRCFLAG)$<
|
||||
|
||||
.m.o:
|
||||
$(ECHO) compiling $(<)
|
||||
$(Q) $(CC) $(INCFLAGS) $(CPPFLAGS) $(CFLAGS) $(COUTFLAG)$@ -c $(CSRCFLAG)$<
|
||||
|
||||
.m.S:
|
||||
$(ECHO) translating $(<)
|
||||
$(Q) $(CC) $(INCFLAGS) $(CPPFLAGS) $(CFLAGS) $(COUTFLAG)$@ -S $(CSRCFLAG)$<
|
||||
|
||||
$(TARGET_SO): $(OBJS) Makefile
|
||||
$(ECHO) linking shared-object $(DLLIB)
|
||||
-$(Q)$(RM) $(@)
|
||||
$(Q) $(LDSHARED) -o $@ $(OBJS) $(LIBPATH) $(DLDFLAGS) $(LOCAL_LIBS) $(LIBS)
|
||||
|
||||
|
||||
|
||||
$(OBJS): $(HDRS) $(ruby_headers)
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,313 @@
|
||||
/*
|
||||
*
|
||||
* Ruby BigDecimal(Variable decimal precision) extension library.
|
||||
*
|
||||
* Copyright(C) 2002 by Shigeo Kobayashi(shigeo@tinyforest.gr.jp)
|
||||
*
|
||||
*/
|
||||
|
||||
#ifndef RUBY_BIG_DECIMAL_H
|
||||
#define RUBY_BIG_DECIMAL_H 1
|
||||
|
||||
#define RUBY_NO_OLD_COMPATIBILITY
|
||||
#include "ruby/ruby.h"
|
||||
#include "missing.h"
|
||||
|
||||
#ifdef HAVE_FLOAT_H
|
||||
# include <float.h>
|
||||
#endif
|
||||
|
||||
#ifdef HAVE_INT64_T
|
||||
# define DECDIG uint32_t
|
||||
# define DECDIG_DBL uint64_t
|
||||
# define DECDIG_DBL_SIGNED int64_t
|
||||
# define SIZEOF_DECDIG 4
|
||||
# define PRI_DECDIG_PREFIX ""
|
||||
# ifdef PRI_LL_PREFIX
|
||||
# define PRI_DECDIG_DBL_PREFIX PRI_LL_PREFIX
|
||||
# else
|
||||
# define PRI_DECDIG_DBL_PREFIX "l"
|
||||
# endif
|
||||
#else
|
||||
# define DECDIG uint16_t
|
||||
# define DECDIG_DBL uint32_t
|
||||
# define DECDIG_DBL_SIGNED int32_t
|
||||
# define SIZEOF_DECDIG 2
|
||||
# define PRI_DECDIG_PREFIX "h"
|
||||
# define PRI_DECDIG_DBL_PREFIX ""
|
||||
#endif
|
||||
|
||||
#define PRIdDECDIG PRI_DECDIG_PREFIX"d"
|
||||
#define PRIiDECDIG PRI_DECDIG_PREFIX"i"
|
||||
#define PRIoDECDIG PRI_DECDIG_PREFIX"o"
|
||||
#define PRIuDECDIG PRI_DECDIG_PREFIX"u"
|
||||
#define PRIxDECDIG PRI_DECDIG_PREFIX"x"
|
||||
#define PRIXDECDIG PRI_DECDIG_PREFIX"X"
|
||||
|
||||
#define PRIdDECDIG_DBL PRI_DECDIG_DBL_PREFIX"d"
|
||||
#define PRIiDECDIG_DBL PRI_DECDIG_DBL_PREFIX"i"
|
||||
#define PRIoDECDIG_DBL PRI_DECDIG_DBL_PREFIX"o"
|
||||
#define PRIuDECDIG_DBL PRI_DECDIG_DBL_PREFIX"u"
|
||||
#define PRIxDECDIG_DBL PRI_DECDIG_DBL_PREFIX"x"
|
||||
#define PRIXDECDIG_DBL PRI_DECDIG_DBL_PREFIX"X"
|
||||
|
||||
#if SIZEOF_DECDIG == 4
|
||||
# define BIGDECIMAL_BASE ((DECDIG)1000000000U)
|
||||
# define BIGDECIMAL_COMPONENT_FIGURES 9
|
||||
/*
|
||||
* The number of components required for a 64-bit integer.
|
||||
*
|
||||
* INT64_MAX: 9_223372036_854775807
|
||||
* UINT64_MAX: 18_446744073_709551615
|
||||
*/
|
||||
# define BIGDECIMAL_INT64_MAX_LENGTH 3
|
||||
|
||||
#elif SIZEOF_DECDIG == 2
|
||||
# define BIGDECIMAL_BASE ((DECDIG)10000U)
|
||||
# define BIGDECIMAL_COMPONENT_FIGURES 4
|
||||
/*
|
||||
* The number of components required for a 64-bit integer.
|
||||
*
|
||||
* INT64_MAX: 922_3372_0368_5477_5807
|
||||
* UINT64_MAX: 1844_6744_0737_0955_1615
|
||||
*/
|
||||
# define BIGDECIMAL_INT64_MAX_LENGTH 5
|
||||
|
||||
#else
|
||||
# error Unknown size of DECDIG
|
||||
#endif
|
||||
|
||||
#define BIGDECIMAL_DOUBLE_FIGURES (1+DBL_DIG)
|
||||
|
||||
#if defined(__cplusplus)
|
||||
extern "C" {
|
||||
#if 0
|
||||
} /* satisfy cc-mode */
|
||||
#endif
|
||||
#endif
|
||||
|
||||
extern VALUE rb_cBigDecimal;
|
||||
|
||||
/*
|
||||
* NaN & Infinity
|
||||
*/
|
||||
#define SZ_NaN "NaN"
|
||||
#define SZ_INF "Infinity"
|
||||
#define SZ_PINF "+Infinity"
|
||||
#define SZ_NINF "-Infinity"
|
||||
|
||||
/*
|
||||
* #define VP_EXPORT other than static to let VP_ routines
|
||||
* be called from outside of this module.
|
||||
*/
|
||||
#define VP_EXPORT static
|
||||
|
||||
/* Exception mode */
|
||||
#define VP_EXCEPTION_ALL ((unsigned short)0x00FF)
|
||||
#define VP_EXCEPTION_INFINITY ((unsigned short)0x0001)
|
||||
#define VP_EXCEPTION_NaN ((unsigned short)0x0002)
|
||||
#define VP_EXCEPTION_UNDERFLOW ((unsigned short)0x0004)
|
||||
#define VP_EXCEPTION_OVERFLOW ((unsigned short)0x0001) /* 0x0008) */
|
||||
#define VP_EXCEPTION_ZERODIVIDE ((unsigned short)0x0010)
|
||||
|
||||
/* Following 2 exceptions can't controlled by user */
|
||||
#define VP_EXCEPTION_OP ((unsigned short)0x0020)
|
||||
|
||||
#define BIGDECIMAL_EXCEPTION_MODE_DEFAULT 0U
|
||||
|
||||
/* This is used in BigDecimal#mode */
|
||||
#define VP_ROUND_MODE ((unsigned short)0x0100)
|
||||
|
||||
/* Rounding mode */
|
||||
#define VP_ROUND_UP RBD_ROUND_UP
|
||||
#define VP_ROUND_DOWN RBD_ROUND_DOWN
|
||||
#define VP_ROUND_HALF_UP RBD_ROUND_HALF_UP
|
||||
#define VP_ROUND_HALF_DOWN RBD_ROUND_HALF_DOWN
|
||||
#define VP_ROUND_CEIL RBD_ROUND_CEIL
|
||||
#define VP_ROUND_FLOOR RBD_ROUND_FLOOR
|
||||
#define VP_ROUND_HALF_EVEN RBD_ROUND_HALF_EVEN
|
||||
|
||||
enum rbd_rounding_mode {
|
||||
RBD_ROUND_UP = 1,
|
||||
RBD_ROUND_DOWN = 2,
|
||||
RBD_ROUND_HALF_UP = 3,
|
||||
RBD_ROUND_HALF_DOWN = 4,
|
||||
RBD_ROUND_CEIL = 5,
|
||||
RBD_ROUND_FLOOR = 6,
|
||||
RBD_ROUND_HALF_EVEN = 7,
|
||||
|
||||
RBD_ROUND_DEFAULT = RBD_ROUND_HALF_UP,
|
||||
RBD_ROUND_TRUNCATE = RBD_ROUND_DOWN,
|
||||
RBD_ROUND_BANKER = RBD_ROUND_HALF_EVEN,
|
||||
RBD_ROUND_CEILING = RBD_ROUND_CEIL
|
||||
};
|
||||
|
||||
#define BIGDECIMAL_ROUNDING_MODE_DEFAULT VP_ROUND_HALF_UP
|
||||
|
||||
/* Sign flag */
|
||||
#define VP_SIGN_NaN 0 /* NaN */
|
||||
#define VP_SIGN_POSITIVE_ZERO 1 /* Positive zero */
|
||||
#define VP_SIGN_NEGATIVE_ZERO -1 /* Negative zero */
|
||||
#define VP_SIGN_POSITIVE_FINITE 2 /* Positive finite number */
|
||||
#define VP_SIGN_NEGATIVE_FINITE -2 /* Negative finite number */
|
||||
#define VP_SIGN_POSITIVE_INFINITE 3 /* Positive infinite number */
|
||||
#define VP_SIGN_NEGATIVE_INFINITE -3 /* Negative infinite number */
|
||||
|
||||
/* The size of fraction part array */
|
||||
#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
|
||||
#define FLEXIBLE_ARRAY_SIZE /* */
|
||||
#elif defined(__GNUC__) && !defined(__STRICT_ANSI__)
|
||||
#define FLEXIBLE_ARRAY_SIZE 0
|
||||
#else
|
||||
#define FLEXIBLE_ARRAY_SIZE 1
|
||||
#endif
|
||||
|
||||
/*
|
||||
* VP representation
|
||||
* r = 0.xxxxxxxxx *BASE**exponent
|
||||
*/
|
||||
typedef struct {
|
||||
VALUE obj; /* Back pointer(VALUE) for Ruby object. */
|
||||
size_t MaxPrec; /* Maximum precision size */
|
||||
/* This is the actual size of frac[] */
|
||||
/*(frac[0] to frac[MaxPrec] are available). */
|
||||
size_t Prec; /* Current precision size. */
|
||||
/* This indicates how much the */
|
||||
/* array frac[] is actually used. */
|
||||
SIGNED_VALUE exponent; /* Exponent part. */
|
||||
short sign; /* Attributes of the value. */
|
||||
/*
|
||||
* ==0 : NaN
|
||||
* 1 : Positive zero
|
||||
* -1 : Negative zero
|
||||
* 2 : Positive number
|
||||
* -2 : Negative number
|
||||
* 3 : Positive infinite number
|
||||
* -3 : Negative infinite number
|
||||
*/
|
||||
short flag; /* Not used in vp_routines,space for user. */
|
||||
DECDIG frac[FLEXIBLE_ARRAY_SIZE]; /* Array of fraction part. */
|
||||
} Real;
|
||||
|
||||
/*
|
||||
* ------------------
|
||||
* EXPORTables.
|
||||
* ------------------
|
||||
*/
|
||||
|
||||
VP_EXPORT Real *VpNewRbClass(size_t mx, char const *str, VALUE klass, bool strict_p, bool raise_exception);
|
||||
|
||||
VP_EXPORT Real *VpCreateRbObject(size_t mx, const char *str, bool raise_exception);
|
||||
|
||||
#define VpBaseFig() BIGDECIMAL_COMPONENT_FIGURES
|
||||
#define VpDblFig() BIGDECIMAL_DOUBLE_FIGURES
|
||||
#define VpBaseVal() BIGDECIMAL_BASE
|
||||
|
||||
/* Zero,Inf,NaN (isinf(),isnan() used to check) */
|
||||
VP_EXPORT double VpGetDoubleNaN(void);
|
||||
VP_EXPORT double VpGetDoublePosInf(void);
|
||||
VP_EXPORT double VpGetDoubleNegInf(void);
|
||||
VP_EXPORT double VpGetDoubleNegZero(void);
|
||||
|
||||
/* These 2 functions added at v1.1.7 */
|
||||
VP_EXPORT size_t VpGetPrecLimit(void);
|
||||
VP_EXPORT size_t VpSetPrecLimit(size_t n);
|
||||
|
||||
/* Round mode */
|
||||
VP_EXPORT int VpIsRoundMode(unsigned short n);
|
||||
VP_EXPORT unsigned short VpGetRoundMode(void);
|
||||
VP_EXPORT unsigned short VpSetRoundMode(unsigned short n);
|
||||
|
||||
VP_EXPORT int VpException(unsigned short f,const char *str,int always);
|
||||
#if 0 /* unused */
|
||||
VP_EXPORT int VpIsNegDoubleZero(double v);
|
||||
#endif
|
||||
VP_EXPORT size_t VpNumOfChars(Real *vp,const char *pszFmt);
|
||||
VP_EXPORT size_t VpInit(DECDIG BaseVal);
|
||||
VP_EXPORT Real *VpAlloc(size_t mx, const char *szVal, int strict_p, int exc);
|
||||
VP_EXPORT size_t VpAsgn(Real *c, Real *a, int isw);
|
||||
VP_EXPORT size_t VpAddSub(Real *c,Real *a,Real *b,int operation);
|
||||
VP_EXPORT size_t VpMult(Real *c,Real *a,Real *b);
|
||||
VP_EXPORT size_t VpDivd(Real *c,Real *r,Real *a,Real *b);
|
||||
VP_EXPORT int VpComp(Real *a,Real *b);
|
||||
VP_EXPORT ssize_t VpExponent10(Real *a);
|
||||
VP_EXPORT void VpSzMantissa(Real *a, char *buf, size_t bufsize);
|
||||
VP_EXPORT int VpToSpecialString(Real *a, char *buf, size_t bufsize, int fPlus);
|
||||
VP_EXPORT void VpToString(Real *a, char *buf, size_t bufsize, size_t fFmt, int fPlus);
|
||||
VP_EXPORT void VpToFString(Real *a, char *buf, size_t bufsize, size_t fFmt, int fPlus);
|
||||
VP_EXPORT int VpCtoV(Real *a, const char *int_chr, size_t ni, const char *frac, size_t nf, const char *exp_chr, size_t ne);
|
||||
VP_EXPORT int VpVtoD(double *d, SIGNED_VALUE *e, Real *m);
|
||||
VP_EXPORT void VpDtoV(Real *m,double d);
|
||||
#if 0 /* unused */
|
||||
VP_EXPORT void VpItoV(Real *m,S_INT ival);
|
||||
#endif
|
||||
VP_EXPORT int VpSqrt(Real *y,Real *x);
|
||||
VP_EXPORT int VpActiveRound(Real *y, Real *x, unsigned short f, ssize_t il);
|
||||
VP_EXPORT int VpMidRound(Real *y, unsigned short f, ssize_t nf);
|
||||
VP_EXPORT int VpLeftRound(Real *y, unsigned short f, ssize_t nf);
|
||||
VP_EXPORT void VpFrac(Real *y, Real *x);
|
||||
VP_EXPORT int VpPowerByInt(Real *y, Real *x, SIGNED_VALUE n);
|
||||
#define VpPower VpPowerByInt
|
||||
|
||||
/* VP constants */
|
||||
VP_EXPORT Real *VpOne(void);
|
||||
|
||||
/*
|
||||
* ------------------
|
||||
* MACRO definitions.
|
||||
* ------------------
|
||||
*/
|
||||
#define Abs(a) (((a)>= 0)?(a):(-(a)))
|
||||
#define Max(a, b) (((a)>(b))?(a):(b))
|
||||
#define Min(a, b) (((a)>(b))?(b):(a))
|
||||
|
||||
#define VpMaxPrec(a) ((a)->MaxPrec)
|
||||
#define VpPrec(a) ((a)->Prec)
|
||||
#define VpGetFlag(a) ((a)->flag)
|
||||
|
||||
/* Sign */
|
||||
|
||||
/* VpGetSign(a) returns 1,-1 if a>0,a<0 respectively */
|
||||
#define VpGetSign(a) (((a)->sign>0)?1:(-1))
|
||||
/* Change sign of a to a>0,a<0 if s = 1,-1 respectively */
|
||||
#define VpChangeSign(a,s) {if((s)>0) (a)->sign=(short)Abs((ssize_t)(a)->sign);else (a)->sign=-(short)Abs((ssize_t)(a)->sign);}
|
||||
/* Sets sign of a to a>0,a<0 if s = 1,-1 respectively */
|
||||
#define VpSetSign(a,s) {if((s)>0) (a)->sign=(short)VP_SIGN_POSITIVE_FINITE;else (a)->sign=(short)VP_SIGN_NEGATIVE_FINITE;}
|
||||
|
||||
/* 1 */
|
||||
#define VpSetOne(a) {(a)->Prec=(a)->exponent=(a)->frac[0]=1;(a)->sign=VP_SIGN_POSITIVE_FINITE;}
|
||||
|
||||
/* ZEROs */
|
||||
#define VpIsPosZero(a) ((a)->sign==VP_SIGN_POSITIVE_ZERO)
|
||||
#define VpIsNegZero(a) ((a)->sign==VP_SIGN_NEGATIVE_ZERO)
|
||||
#define VpIsZero(a) (VpIsPosZero(a) || VpIsNegZero(a))
|
||||
#define VpSetPosZero(a) ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_POSITIVE_ZERO)
|
||||
#define VpSetNegZero(a) ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_NEGATIVE_ZERO)
|
||||
#define VpSetZero(a,s) (void)(((s)>0)?VpSetPosZero(a):VpSetNegZero(a))
|
||||
|
||||
/* NaN */
|
||||
#define VpIsNaN(a) ((a)->sign==VP_SIGN_NaN)
|
||||
#define VpSetNaN(a) ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_NaN)
|
||||
|
||||
/* Infinity */
|
||||
#define VpIsPosInf(a) ((a)->sign==VP_SIGN_POSITIVE_INFINITE)
|
||||
#define VpIsNegInf(a) ((a)->sign==VP_SIGN_NEGATIVE_INFINITE)
|
||||
#define VpIsInf(a) (VpIsPosInf(a) || VpIsNegInf(a))
|
||||
#define VpIsDef(a) ( !(VpIsNaN(a)||VpIsInf(a)) )
|
||||
#define VpSetPosInf(a) ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_POSITIVE_INFINITE)
|
||||
#define VpSetNegInf(a) ((a)->frac[0]=0,(a)->Prec=1,(a)->sign=VP_SIGN_NEGATIVE_INFINITE)
|
||||
#define VpSetInf(a,s) (void)(((s)>0)?VpSetPosInf(a):VpSetNegInf(a))
|
||||
#define VpHasVal(a) (a->frac[0])
|
||||
#define VpIsOne(a) ((a->Prec==1)&&(a->frac[0]==1)&&(a->exponent==1))
|
||||
#define VpExponent(a) (a->exponent)
|
||||
#ifdef BIGDECIMAL_DEBUG
|
||||
int VpVarCheck(Real * v);
|
||||
#endif /* BIGDECIMAL_DEBUG */
|
||||
|
||||
#if defined(__cplusplus)
|
||||
#if 0
|
||||
{ /* satisfy cc-mode */
|
||||
#endif
|
||||
} /* extern "C" { */
|
||||
#endif
|
||||
#endif /* RUBY_BIG_DECIMAL_H */
|
||||
@@ -0,0 +1,141 @@
|
||||
#ifndef BIGDECIMAL_BITS_H
|
||||
#define BIGDECIMAL_BITS_H
|
||||
|
||||
#include "feature.h"
|
||||
#include "static_assert.h"
|
||||
|
||||
#if defined(__x86_64__) && defined(HAVE_X86INTRIN_H)
|
||||
# include <x86intrin.h> /* for _lzcnt_u64, etc. */
|
||||
#elif defined(_MSC_VER) && defined(HAVE_INTRIN_H)
|
||||
# include <intrin.h> /* for the following intrinsics */
|
||||
#endif
|
||||
|
||||
#if defined(_MSC_VER) && defined(__AVX2__)
|
||||
# pragma intrinsic(__lzcnt)
|
||||
# pragma intrinsic(__lzcnt64)
|
||||
#endif
|
||||
|
||||
#define numberof(array) ((int)(sizeof(array) / sizeof((array)[0])))
|
||||
#define roomof(x, y) (((x) + (y) - 1) / (y))
|
||||
#define type_roomof(x, y) roomof(sizeof(x), sizeof(y))
|
||||
|
||||
#define MUL_OVERFLOW_SIGNED_INTEGER_P(a, b, min, max) ( \
|
||||
(a) == 0 ? 0 : \
|
||||
(a) == -1 ? (b) < -(max) : \
|
||||
(a) > 0 ? \
|
||||
((b) > 0 ? (max) / (a) < (b) : (min) / (a) > (b)) : \
|
||||
((b) > 0 ? (min) / (a) < (b) : (max) / (a) > (b)))
|
||||
|
||||
#ifdef HAVE_UINT128_T
|
||||
# define bit_length(x) \
|
||||
(unsigned int) \
|
||||
(sizeof(x) <= sizeof(int32_t) ? 32 - nlz_int32((uint32_t)(x)) : \
|
||||
sizeof(x) <= sizeof(int64_t) ? 64 - nlz_int64((uint64_t)(x)) : \
|
||||
128 - nlz_int128((uint128_t)(x)))
|
||||
#else
|
||||
# define bit_length(x) \
|
||||
(unsigned int) \
|
||||
(sizeof(x) <= sizeof(int32_t) ? 32 - nlz_int32((uint32_t)(x)) : \
|
||||
64 - nlz_int64((uint64_t)(x)))
|
||||
#endif
|
||||
|
||||
static inline unsigned nlz_int32(uint32_t x);
|
||||
static inline unsigned nlz_int64(uint64_t x);
|
||||
#ifdef HAVE_UINT128_T
|
||||
static inline unsigned nlz_int128(uint128_t x);
|
||||
#endif
|
||||
|
||||
static inline unsigned int
|
||||
nlz_int32(uint32_t x)
|
||||
{
|
||||
#if defined(_MSC_VER) && defined(__AVX2__) && defined(HAVE___LZCNT)
|
||||
/* Note: It seems there is no such thing like __LZCNT__ predefined in MSVC.
|
||||
* AMD CPUs have had this instruction for decades (since K10) but for
|
||||
* Intel, Haswell is the oldest one. We need to use __AVX2__ for maximum
|
||||
* safety. */
|
||||
return (unsigned int)__lzcnt(x);
|
||||
|
||||
#elif defined(__x86_64__) && defined(__LZCNT__) && defined(HAVE__LZCNT_U32)
|
||||
return (unsigned int)_lzcnt_u32(x);
|
||||
|
||||
#elif defined(_MSC_VER) && defined(HAVE__BITSCANREVERSE)
|
||||
unsigned long r;
|
||||
return _BitScanReverse(&r, x) ? (31 - (int)r) : 32;
|
||||
|
||||
#elif __has_builtin(__builtin_clz)
|
||||
STATIC_ASSERT(sizeof_int, sizeof(int) * CHAR_BIT == 32);
|
||||
return x ? (unsigned int)__builtin_clz(x) : 32;
|
||||
|
||||
#else
|
||||
uint32_t y;
|
||||
unsigned n = 32;
|
||||
y = x >> 16; if (y) {n -= 16; x = y;}
|
||||
y = x >> 8; if (y) {n -= 8; x = y;}
|
||||
y = x >> 4; if (y) {n -= 4; x = y;}
|
||||
y = x >> 2; if (y) {n -= 2; x = y;}
|
||||
y = x >> 1; if (y) {return n - 2;}
|
||||
return (unsigned int)(n - x);
|
||||
#endif
|
||||
}
|
||||
|
||||
static inline unsigned int
|
||||
nlz_int64(uint64_t x)
|
||||
{
|
||||
#if defined(_MSC_VER) && defined(__AVX2__) && defined(HAVE___LZCNT64)
|
||||
return (unsigned int)__lzcnt64(x);
|
||||
|
||||
#elif defined(__x86_64__) && defined(__LZCNT__) && defined(HAVE__LZCNT_U64)
|
||||
return (unsigned int)_lzcnt_u64(x);
|
||||
|
||||
#elif defined(_WIN64) && defined(_MSC_VER) && defined(HAVE__BITSCANREVERSE64)
|
||||
unsigned long r;
|
||||
return _BitScanReverse64(&r, x) ? (63u - (unsigned int)r) : 64;
|
||||
|
||||
#elif __has_builtin(__builtin_clzl) && __has_builtin(__builtin_clzll) && !(defined(__sun) && defined(__sparc))
|
||||
if (x == 0) {
|
||||
return 64;
|
||||
}
|
||||
else if (sizeof(long) * CHAR_BIT == 64) {
|
||||
return (unsigned int)__builtin_clzl((unsigned long)x);
|
||||
}
|
||||
else if (sizeof(long long) * CHAR_BIT == 64) {
|
||||
return (unsigned int)__builtin_clzll((unsigned long long)x);
|
||||
}
|
||||
else {
|
||||
/* :FIXME: Is there a way to make this branch a compile-time error? */
|
||||
__builtin_unreachable();
|
||||
}
|
||||
|
||||
#else
|
||||
uint64_t y;
|
||||
unsigned int n = 64;
|
||||
y = x >> 32; if (y) {n -= 32; x = y;}
|
||||
y = x >> 16; if (y) {n -= 16; x = y;}
|
||||
y = x >> 8; if (y) {n -= 8; x = y;}
|
||||
y = x >> 4; if (y) {n -= 4; x = y;}
|
||||
y = x >> 2; if (y) {n -= 2; x = y;}
|
||||
y = x >> 1; if (y) {return n - 2;}
|
||||
return (unsigned int)(n - x);
|
||||
|
||||
#endif
|
||||
}
|
||||
|
||||
#ifdef HAVE_UINT128_T
|
||||
static inline unsigned int
|
||||
nlz_int128(uint128_t x)
|
||||
{
|
||||
uint64_t y = (uint64_t)(x >> 64);
|
||||
|
||||
if (x == 0) {
|
||||
return 128;
|
||||
}
|
||||
else if (y == 0) {
|
||||
return (unsigned int)nlz_int64(x) + 64;
|
||||
}
|
||||
else {
|
||||
return (unsigned int)nlz_int64(y);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif /* BIGDECIMAL_BITS_H */
|
||||
@@ -0,0 +1,64 @@
|
||||
# frozen_string_literal: false
|
||||
require 'mkmf'
|
||||
|
||||
def have_builtin_func(name, check_expr, opt = "", &b)
|
||||
checking_for checking_message(name.funcall_style, nil, opt) do
|
||||
if try_compile(<<SRC, opt, &b)
|
||||
int foo;
|
||||
int main() { #{check_expr}; return 0; }
|
||||
SRC
|
||||
$defs.push(format("-DHAVE_BUILTIN_%s", name.tr_cpp))
|
||||
true
|
||||
else
|
||||
false
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
have_builtin_func("__builtin_clz", "__builtin_clz(0)")
|
||||
have_builtin_func("__builtin_clzl", "__builtin_clzl(0)")
|
||||
have_builtin_func("__builtin_clzll", "__builtin_clzll(0)")
|
||||
|
||||
have_header("float.h")
|
||||
have_header("math.h")
|
||||
have_header("stdbool.h")
|
||||
have_header("stdlib.h")
|
||||
|
||||
if have_header("x86intrin.h")
|
||||
have_func("_lzcnt_u32", "x86intrin.h")
|
||||
have_func("_lzcnt_u64", "x86intrin.h")
|
||||
end
|
||||
|
||||
if have_header("intrin.h")
|
||||
have_func("__lzcnt", "intrin.h")
|
||||
have_func("__lzcnt64", "intrin.h")
|
||||
have_func("_BitScanReverse", "intrin.h")
|
||||
have_func("_BitScanReverse64", "intrin.h")
|
||||
end
|
||||
|
||||
have_func("labs", "stdlib.h")
|
||||
have_func("llabs", "stdlib.h")
|
||||
have_func("finite", "math.h")
|
||||
have_func("isfinite", "math.h")
|
||||
|
||||
have_header("ruby/atomic.h")
|
||||
have_header("ruby/internal/has/builtin.h")
|
||||
have_header("ruby/internal/static_assert.h")
|
||||
|
||||
have_func("rb_rational_num", "ruby.h")
|
||||
have_func("rb_rational_den", "ruby.h")
|
||||
have_func("rb_complex_real", "ruby.h")
|
||||
have_func("rb_complex_imag", "ruby.h")
|
||||
have_func("rb_opts_exception_p", "ruby.h")
|
||||
have_func("rb_category_warn", "ruby.h")
|
||||
have_const("RB_WARN_CATEGORY_DEPRECATED", "ruby.h")
|
||||
|
||||
if File.file?(File.expand_path('../lib/bigdecimal.rb', __FILE__))
|
||||
bigdecimal_rb = "$(srcdir)/lib/bigdecimal.rb"
|
||||
else
|
||||
bigdecimal_rb = "$(srcdir)/../../lib/bigdecimal.rb"
|
||||
end
|
||||
|
||||
create_makefile('bigdecimal') {|mf|
|
||||
mf << "BIGDECIMAL_RB = #{bigdecimal_rb}\n"
|
||||
}
|
||||
@@ -0,0 +1,68 @@
|
||||
#ifndef BIGDECIMAL_HAS_FEATURE_H
|
||||
#define BIGDECIMAL_HAS_FEATURE_H
|
||||
|
||||
/* ======== __has_feature ======== */
|
||||
|
||||
#ifndef __has_feature
|
||||
# define __has_feature(_) 0
|
||||
#endif
|
||||
|
||||
/* ======== __has_extension ======== */
|
||||
|
||||
#ifndef __has_extension
|
||||
# define __has_extension __has_feature
|
||||
#endif
|
||||
|
||||
/* ======== __has_builtin ======== */
|
||||
|
||||
#ifdef HAVE_RUBY_INTERNAL_HAS_BUILTIN_H
|
||||
# include <ruby/internal/has/builtin.h>
|
||||
#endif
|
||||
|
||||
#ifdef RBIMPL_HAS_BUILTIN
|
||||
# define BIGDECIMAL_HAS_BUILTIN(...) RBIMPL_HAS_BUILTIN(__VA_ARGS__)
|
||||
|
||||
#else
|
||||
# /* The following section is copied from CRuby's builtin.h */
|
||||
#
|
||||
# ifdef __has_builtin
|
||||
# if defined(__INTEL_COMPILER)
|
||||
# /* :TODO: Intel C Compiler has __has_builtin (since 19.1 maybe?), and is
|
||||
# * reportedly broken. We have to skip them. However the situation can
|
||||
# * change. They might improve someday. We need to revisit here later. */
|
||||
# elif defined(__GNUC__) && ! __has_builtin(__builtin_alloca)
|
||||
# /* FreeBSD's <sys/cdefs.h> defines its own *broken* version of
|
||||
# * __has_builtin. Cygwin copied that content to be a victim of the
|
||||
# * broken-ness. We don't take them into account. */
|
||||
# else
|
||||
# define HAVE___HAS_BUILTIN 1
|
||||
# endif
|
||||
# endif
|
||||
#
|
||||
# if defined(HAVE___HAS_BUILTIN)
|
||||
# define BIGDECIMAL_HAS_BUILTIN(_) __has_builtin(_)
|
||||
#
|
||||
# elif defined(__GNUC__)
|
||||
# define BIGDECIMAL_HAS_BUILTIN(_) BIGDECIMAL_HAS_BUILTIN_ ## _
|
||||
# if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 6))
|
||||
# define BIGDECIMAL_HAS_BUILTIN___builtin_clz 1
|
||||
# define BIGDECIMAL_HAS_BUILTIN___builtin_clzl 1
|
||||
# else
|
||||
# define BIGDECIMAL_HAS_BUILTIN___builtin_clz 0
|
||||
# define BIGDECIMAL_HAS_BUILTIN___builtin_clzl 0
|
||||
# endif
|
||||
# elif defined(_MSC_VER)
|
||||
# define BIGDECIMAL_HAS_BUILTIN(_) 0
|
||||
#
|
||||
# else
|
||||
# define BIGDECIMAL_HAS_BUILTIN(_) BIGDECIMAL_HAS_BUILTIN_ ## _
|
||||
# define BIGDECIMAL_HAS_BUILTIN___builtin_clz HAVE_BUILTIN___BUILTIN_CLZ
|
||||
# define BIGDECIMAL_HAS_BUILTIN___builtin_clzl HAVE_BUILTIN___BUILTIN_CLZL
|
||||
# endif
|
||||
#endif /* RBIMPL_HAS_BUILTIN */
|
||||
|
||||
#ifndef __has_builtin
|
||||
# define __has_builtin(...) BIGDECIMAL_HAS_BUILTIN(__VA_ARGS__)
|
||||
#endif
|
||||
|
||||
#endif /* BIGDECIMAL_HAS_FEATURE_H */
|
||||
@@ -0,0 +1,28 @@
|
||||
#include <ruby/ruby.h>
|
||||
|
||||
#ifdef HAVE_RUBY_ATOMIC_H
|
||||
# include <ruby/atomic.h>
|
||||
#endif
|
||||
|
||||
#ifdef RUBY_ATOMIC_PTR_CAS
|
||||
# define ATOMIC_PTR_CAS(var, old, new) RUBY_ATOMIC_PTR_CAS(var, old, new)
|
||||
#endif
|
||||
|
||||
#if defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6))
|
||||
/* GCC warns about unknown sanitizer, which is annoying. */
|
||||
# undef NO_SANITIZE
|
||||
# define NO_SANITIZE(x, y) \
|
||||
_Pragma("GCC diagnostic push") \
|
||||
_Pragma("GCC diagnostic ignored \"-Wattributes\"") \
|
||||
__attribute__((__no_sanitize__(x))) y; \
|
||||
_Pragma("GCC diagnostic pop") \
|
||||
y
|
||||
#endif
|
||||
|
||||
#undef strtod
|
||||
#define strtod BigDecimal_strtod
|
||||
#undef dtoa
|
||||
#define dtoa BigDecimal_dtoa
|
||||
#undef hdtoa
|
||||
#define hdtoa BigDecimal_hdtoa
|
||||
#include "missing/dtoa.c"
|
||||
@@ -0,0 +1,196 @@
|
||||
#ifndef MISSING_H
|
||||
#define MISSING_H 1
|
||||
|
||||
#if defined(__cplusplus)
|
||||
extern "C" {
|
||||
#if 0
|
||||
} /* satisfy cc-mode */
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#ifdef HAVE_STDLIB_H
|
||||
# include <stdlib.h>
|
||||
#endif
|
||||
|
||||
#ifdef HAVE_MATH_H
|
||||
# include <math.h>
|
||||
#endif
|
||||
|
||||
#ifndef RB_UNUSED_VAR
|
||||
# if defined(_MSC_VER) && _MSC_VER >= 1911
|
||||
# define RB_UNUSED_VAR(x) x [[maybe_unused]]
|
||||
|
||||
# elif defined(__has_cpp_attribute) && __has_cpp_attribute(maybe_unused)
|
||||
# define RB_UNUSED_VAR(x) x [[maybe_unused]]
|
||||
|
||||
# elif defined(__has_c_attribute) && __has_c_attribute(maybe_unused)
|
||||
# define RB_UNUSED_VAR(x) x [[maybe_unused]]
|
||||
|
||||
# elif defined(__GNUC__)
|
||||
# define RB_UNUSED_VAR(x) x __attribute__ ((unused))
|
||||
|
||||
# else
|
||||
# define RB_UNUSED_VAR(x) x
|
||||
# endif
|
||||
#endif /* RB_UNUSED_VAR */
|
||||
|
||||
#if defined(_MSC_VER) && _MSC_VER >= 1310
|
||||
# define HAVE___ASSUME 1
|
||||
|
||||
#elif defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1300
|
||||
# define HAVE___ASSUME 1
|
||||
#endif
|
||||
|
||||
#ifndef UNREACHABLE
|
||||
# if __has_builtin(__builtin_unreachable)
|
||||
# define UNREACHABLE __builtin_unreachable()
|
||||
|
||||
# elif defined(HAVE___ASSUME)
|
||||
# define UNREACHABLE __assume(0)
|
||||
|
||||
# else
|
||||
# define UNREACHABLE /* unreachable */
|
||||
# endif
|
||||
#endif /* UNREACHABLE */
|
||||
|
||||
/* bool */
|
||||
|
||||
#if defined(__bool_true_false_are_defined)
|
||||
# /* Take that. */
|
||||
|
||||
#elif defined(HAVE_STDBOOL_H)
|
||||
# include <stdbool.h>
|
||||
|
||||
#else
|
||||
typedef unsigned char _Bool;
|
||||
# define bool _Bool
|
||||
# define true ((_Bool)+1)
|
||||
# define false ((_Bool)-1)
|
||||
# define __bool_true_false_are_defined
|
||||
#endif
|
||||
|
||||
/* abs */
|
||||
|
||||
#ifndef HAVE_LABS
|
||||
static inline long
|
||||
labs(long const x)
|
||||
{
|
||||
if (x < 0) return -x;
|
||||
return x;
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifndef HAVE_LLABS
|
||||
static inline LONG_LONG
|
||||
llabs(LONG_LONG const x)
|
||||
{
|
||||
if (x < 0) return -x;
|
||||
return x;
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifdef vabs
|
||||
# undef vabs
|
||||
#endif
|
||||
#if SIZEOF_VALUE <= SIZEOF_INT
|
||||
# define vabs abs
|
||||
#elif SIZEOF_VALUE <= SIZEOF_LONG
|
||||
# define vabs labs
|
||||
#elif SIZEOF_VALUE <= SIZEOF_LONG_LONG
|
||||
# define vabs llabs
|
||||
#endif
|
||||
|
||||
/* finite */
|
||||
|
||||
#ifndef HAVE_FINITE
|
||||
static int
|
||||
finite(double)
|
||||
{
|
||||
return !isnan(n) && !isinf(n);
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifndef isfinite
|
||||
# ifndef HAVE_ISFINITE
|
||||
# define HAVE_ISFINITE 1
|
||||
# define isfinite(x) finite(x)
|
||||
# endif
|
||||
#endif
|
||||
|
||||
/* dtoa */
|
||||
char *BigDecimal_dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve);
|
||||
|
||||
/* rational */
|
||||
|
||||
#ifndef HAVE_RB_RATIONAL_NUM
|
||||
static inline VALUE
|
||||
rb_rational_num(VALUE rat)
|
||||
{
|
||||
#ifdef RRATIONAL
|
||||
return RRATIONAL(rat)->num;
|
||||
#else
|
||||
return rb_funcall(rat, rb_intern("numerator"), 0);
|
||||
#endif
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifndef HAVE_RB_RATIONAL_DEN
|
||||
static inline VALUE
|
||||
rb_rational_den(VALUE rat)
|
||||
{
|
||||
#ifdef RRATIONAL
|
||||
return RRATIONAL(rat)->den;
|
||||
#else
|
||||
return rb_funcall(rat, rb_intern("denominator"), 0);
|
||||
#endif
|
||||
}
|
||||
#endif
|
||||
|
||||
/* complex */
|
||||
|
||||
#ifndef HAVE_RB_COMPLEX_REAL
|
||||
static inline VALUE
|
||||
rb_complex_real(VALUE cmp)
|
||||
{
|
||||
#ifdef RCOMPLEX
|
||||
return RCOMPLEX(cmp)->real;
|
||||
#else
|
||||
return rb_funcall(cmp, rb_intern("real"), 0);
|
||||
#endif
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifndef HAVE_RB_COMPLEX_IMAG
|
||||
static inline VALUE
|
||||
rb_complex_imag(VALUE cmp)
|
||||
{
|
||||
# ifdef RCOMPLEX
|
||||
return RCOMPLEX(cmp)->imag;
|
||||
# else
|
||||
return rb_funcall(cmp, rb_intern("imag"), 0);
|
||||
# endif
|
||||
}
|
||||
#endif
|
||||
|
||||
/* st */
|
||||
|
||||
#ifndef ST2FIX
|
||||
# undef RB_ST2FIX
|
||||
# define RB_ST2FIX(h) LONG2FIX((long)(h))
|
||||
# define ST2FIX(h) RB_ST2FIX(h)
|
||||
#endif
|
||||
|
||||
/* warning */
|
||||
|
||||
#if !defined(HAVE_RB_CATEGORY_WARN) || !defined(HAVE_CONST_RB_WARN_CATEGORY_DEPRECATED)
|
||||
# define rb_category_warn(category, ...) rb_warn(__VA_ARGS__)
|
||||
#endif
|
||||
|
||||
#if defined(__cplusplus)
|
||||
#if 0
|
||||
{ /* satisfy cc-mode */
|
||||
#endif
|
||||
} /* extern "C" { */
|
||||
#endif
|
||||
|
||||
#endif /* MISSING_H */
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,54 @@
|
||||
#ifndef BIGDECIMAL_STATIC_ASSERT_H
|
||||
#define BIGDECIMAL_STATIC_ASSERT_H
|
||||
|
||||
#include "feature.h"
|
||||
|
||||
#ifdef HAVE_RUBY_INTERNAL_STATIC_ASSERT_H
|
||||
# include <ruby/internal/static_assert.h>
|
||||
#endif
|
||||
|
||||
#ifdef RBIMPL_STATIC_ASSERT
|
||||
# define STATIC_ASSERT RBIMPL_STATIC_ASSERT
|
||||
#endif
|
||||
|
||||
#ifndef STATIC_ASSERT
|
||||
# /* The following section is copied from CRuby's static_assert.h */
|
||||
|
||||
# if defined(__cplusplus) && defined(__cpp_static_assert)
|
||||
# /* https://isocpp.org/std/standing-documents/sd-6-sg10-feature-test-recommendations */
|
||||
# define BIGDECIMAL_STATIC_ASSERT0 static_assert
|
||||
|
||||
# elif defined(__cplusplus) && defined(_MSC_VER) && _MSC_VER >= 1600
|
||||
# define BIGDECIMAL_STATIC_ASSERT0 static_assert
|
||||
|
||||
# elif defined(__INTEL_CXX11_MODE__)
|
||||
# define BIGDECIMAL_STATIC_ASSERT0 static_assert
|
||||
|
||||
# elif defined(__cplusplus) && __cplusplus >= 201103L
|
||||
# define BIGDECIMAL_STATIC_ASSERT0 static_assert
|
||||
|
||||
# elif defined(__cplusplus) && __has_extension(cxx_static_assert)
|
||||
# define BIGDECIMAL_STATIC_ASSERT0 __extension__ static_assert
|
||||
|
||||
# elif defined(__STDC_VERSION__) && __has_extension(c_static_assert)
|
||||
# define BIGDECIMAL_STATIC_ASSERT0 __extension__ _Static_assert
|
||||
|
||||
# elif defined(__STDC_VERSION__) && defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6))
|
||||
# define BIGDECIMAL_STATIC_ASSERT0 __extension__ _Static_assert
|
||||
#endif
|
||||
|
||||
# if defined(__DOXYGEN__)
|
||||
# define STATIC_ASSERT static_assert
|
||||
|
||||
# elif defined(BIGDECIMAL_STATIC_ASSERT0)
|
||||
# define STATIC_ASSERT(name, expr) \
|
||||
BIGDECIMAL_STATIC_ASSERT0(expr, #name ": " #expr)
|
||||
|
||||
# else
|
||||
# define STATIC_ASSERT(name, expr) \
|
||||
typedef int static_assert_ ## name ## _check[1 - 2 * !(expr)]
|
||||
# endif
|
||||
#endif /* STATIC_ASSERT */
|
||||
|
||||
|
||||
#endif /* BIGDECIMAL_STATIC_ASSERT_H */
|
||||
@@ -0,0 +1,5 @@
|
||||
if RUBY_ENGINE == 'jruby'
|
||||
JRuby::Util.load_ext("org.jruby.ext.bigdecimal.BigDecimalLibrary")
|
||||
else
|
||||
require 'bigdecimal.so'
|
||||
end
|
||||
Executable
BIN
Binary file not shown.
@@ -0,0 +1,90 @@
|
||||
# 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
|
||||
@@ -0,0 +1,89 @@
|
||||
# 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
|
||||
@@ -0,0 +1,232 @@
|
||||
# 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
|
||||
@@ -0,0 +1,80 @@
|
||||
# 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
|
||||
@@ -0,0 +1,185 @@
|
||||
# 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 error
|
||||
|
||||
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
|
||||
@@ -0,0 +1,74 @@
|
||||
#!/usr/local/bin/ruby
|
||||
# frozen_string_literal: false
|
||||
|
||||
#
|
||||
# linear.rb
|
||||
#
|
||||
# Solves linear equation system(A*x = b) by LU decomposition method.
|
||||
# where A is a coefficient matrix,x is an answer vector,b is a constant vector.
|
||||
#
|
||||
# USAGE:
|
||||
# ruby linear.rb [input file solved]
|
||||
#
|
||||
|
||||
# :stopdoc:
|
||||
require "bigdecimal"
|
||||
require "bigdecimal/ludcmp"
|
||||
|
||||
#
|
||||
# NOTE:
|
||||
# Change following BigDecimal.limit() if needed.
|
||||
BigDecimal.limit(100)
|
||||
#
|
||||
|
||||
include LUSolve
|
||||
def rd_order(na)
|
||||
printf("Number of equations ?") if(na <= 0)
|
||||
n = ARGF.gets().to_i
|
||||
end
|
||||
|
||||
na = ARGV.size
|
||||
zero = BigDecimal("0.0")
|
||||
one = BigDecimal("1.0")
|
||||
|
||||
while (n=rd_order(na))>0
|
||||
a = []
|
||||
as= []
|
||||
b = []
|
||||
if na <= 0
|
||||
# Read data from console.
|
||||
printf("\nEnter coefficient matrix element A[i,j]\n")
|
||||
for i in 0...n do
|
||||
for j in 0...n do
|
||||
printf("A[%d,%d]? ",i,j); s = ARGF.gets
|
||||
a << BigDecimal(s)
|
||||
as << BigDecimal(s)
|
||||
end
|
||||
printf("Contatant vector element b[%d] ? ",i)
|
||||
b << BigDecimal(ARGF.gets)
|
||||
end
|
||||
else
|
||||
# Read data from specified file.
|
||||
printf("Coefficient matrix and constant vector.\n")
|
||||
for i in 0...n do
|
||||
s = ARGF.gets
|
||||
printf("%d) %s",i,s)
|
||||
s = s.split
|
||||
for j in 0...n do
|
||||
a << BigDecimal(s[j])
|
||||
as << BigDecimal(s[j])
|
||||
end
|
||||
b << BigDecimal(s[n])
|
||||
end
|
||||
end
|
||||
x = lusolve(a,b,ludecomp(a,n,zero,one),zero)
|
||||
printf("Answer(x[i] & (A*x-b)[i]) follows\n")
|
||||
for i in 0...n do
|
||||
printf("x[%d]=%s ",i,x[i].to_s)
|
||||
s = zero
|
||||
for j in 0...n do
|
||||
s = s + as[i*n+j]*x[j]
|
||||
end
|
||||
printf(" & %s\n",(s-b[i]).to_s)
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,40 @@
|
||||
#!/usr/local/bin/ruby
|
||||
# frozen_string_literal: false
|
||||
|
||||
#
|
||||
# nlsolve.rb
|
||||
# An example for solving nonlinear algebraic equation system.
|
||||
#
|
||||
|
||||
require "bigdecimal"
|
||||
require "bigdecimal/newton"
|
||||
include Newton
|
||||
|
||||
class Function # :nodoc: all
|
||||
def initialize()
|
||||
@zero = BigDecimal("0.0")
|
||||
@one = BigDecimal("1.0")
|
||||
@two = BigDecimal("2.0")
|
||||
@ten = BigDecimal("10.0")
|
||||
@eps = BigDecimal("1.0e-16")
|
||||
end
|
||||
def zero;@zero;end
|
||||
def one ;@one ;end
|
||||
def two ;@two ;end
|
||||
def ten ;@ten ;end
|
||||
def eps ;@eps ;end
|
||||
def values(x) # <= defines functions solved
|
||||
f = []
|
||||
f1 = x[0]*x[0] + x[1]*x[1] - @two # f1 = x**2 + y**2 - 2 => 0
|
||||
f2 = x[0] - x[1] # f2 = x - y => 0
|
||||
f <<= f1
|
||||
f <<= f2
|
||||
f
|
||||
end
|
||||
end
|
||||
|
||||
f = BigDecimal.limit(100)
|
||||
f = Function.new
|
||||
x = [f.zero,f.zero] # Initial values
|
||||
n = nlsolve(f,x)
|
||||
p x
|
||||
@@ -0,0 +1,21 @@
|
||||
#!/usr/local/bin/ruby
|
||||
# frozen_string_literal: false
|
||||
|
||||
#
|
||||
# pi.rb
|
||||
#
|
||||
# Calculates 3.1415.... (the number of times that a circle's diameter
|
||||
# will fit around the circle) using J. Machin's formula.
|
||||
#
|
||||
|
||||
require "bigdecimal"
|
||||
require "bigdecimal/math.rb"
|
||||
|
||||
include BigMath
|
||||
|
||||
if ARGV.size == 1
|
||||
print "PI("+ARGV[0]+"):\n"
|
||||
p PI(ARGV[0].to_i)
|
||||
else
|
||||
print "TRY: ruby pi.rb 1000 \n"
|
||||
end
|
||||
Reference in New Issue
Block a user