Copied from rough/bigdecimal,documents & some sample programs added.

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@3625 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
This commit is contained in:
shigek 2003-03-28 05:00:21 +00:00
Родитель 8369db4b3b
Коммит 7744351708
17 изменённых файлов: 6325 добавлений и 0 удалений

16
ext/bigdecimal/MANIFEST Normal file
Просмотреть файл

@ -0,0 +1,16 @@
MANIFEST
README
bigdecimal.def
bigdecimal.c
bigdecimal.h
depend
extconf.rb
bigdecimal_en.html
bigdecimal_ja.html
lib/nlsolve.rb
lib/jacobian.rb
lib/newton.rb
lib/linear.rb
lib/ludcmp.rb
lib/pai.rb
lib/bigdecimal-rational.rb

60
ext/bigdecimal/README Normal file
Просмотреть файл

@ -0,0 +1,60 @@
Ruby BIGDECIMAL(Variable Precision) extension library.
Copyright (C) 1999 by Shigeo Kobayashi(shigeo@tinyforest.gr.jp)
BigDecimal is copyrighted free software by Shigeo Kobayashi <shigeo@tinyforest.gr.jp>.
You can redistribute it and/or modify it under either the terms of the GPL
(see COPYING file), 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) rename any non-standard executables so the names do not conflict
with standard executables, which must also be provided.
d) make other distribution arrangements with the author.
3. You may distribute the software in object code or executable
form, provided that you do at least ONE of the following:
a) distribute the executables 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 executables 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).
5. THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
PURPOSE.
* The Author
Feel free to send comments and bug reports to the author. Here is the
author's latest mail address:
shigeo@tinyforest.gr.jp
-------------------------------------------------------
created at: Thu Dec 22 1999

4198
ext/bigdecimal/bigdecimal.c Normal file

Разница между файлами не показана из-за своего большого размера Загрузить разницу

Просмотреть файл

@ -0,0 +1,2 @@
EXPORTS
Init_bigdecimal

190
ext/bigdecimal/bigdecimal.h Normal file
Просмотреть файл

@ -0,0 +1,190 @@
/*
*
* Ruby BigDecimal(Variable decimal precision) extension library.
*
* Copyright(C) 2002 by Shigeo Kobayashi(shigeo@tinyforest.gr.jp)
*
* You may distribute under the terms of either the GNU General Public
* License or the Artistic License, as specified in the README file
* of this BigDecimal distribution.
*
* 2002-12-06
* The sqrt() bug reported by Bret Jolly fixed.
*
*/
#ifndef ____BIG_DECIMAL__H____
#define ____BIG_DECIMAL__H____
#if defined(__cplusplus)
extern "C" {
#endif
/*
* #define VP_EXPORT other than static to let VP_ routines
* be called from outside of this module.
*/
#define VP_EXPORT static
#define U_LONG unsigned long
#define S_LONG long
#define U_INT unsigned int
#define S_INT int
/* Exception codes */
#define VP_EXCEPTION_ALL ((unsigned short)0xFFFF)
#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)0x0001) /* 0x0010) */
/* Following 2 exceptions cann't controlled by user */
#define VP_EXCEPTION_OP ((unsigned short)0x0020)
#define VP_EXCEPTION_MEMORY ((unsigned short)0x0040)
#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 */
/*
* VP representation
* r = 0.xxxxxxxxx *BASE**exponent
*/
typedef struct {
VALUE obj; /* Back pointer(VALUE) for Ruby object. */
U_LONG MaxPrec; /* Maximum precision size */
/* This is the actual size of pfrac[] */
/*(frac[0] to frac[MaxPrec] are available). */
U_LONG Prec; /* Current precision size. */
/* This indicates how much the. */
/* the array frac[] is actually used. */
S_INT 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. */
U_LONG frac[1]; /* Pointer to array of fraction part. */
} Real;
/*
* ------------------
* EXPORTables.
* ------------------
*/
VP_EXPORT Real *
VpNewRbClass(U_LONG mx,char *str,VALUE klass);
VP_EXPORT Real *VpCreateRbObject(U_LONG mx,char *str);
VP_EXPORT U_LONG VpBaseFig(void);
VP_EXPORT U_LONG VpDblFig(void);
VP_EXPORT U_LONG VpBaseVal(void);
/* 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 U_LONG VpGetPrecLimit(void);
VP_EXPORT U_LONG VpSetPrecLimit(U_LONG n);
VP_EXPORT int VpException(unsigned short f,char *str,int always);
VP_EXPORT int VpIsNegDoubleZero(double v);
VP_EXPORT U_LONG VpNumOfChars(Real *vp);
VP_EXPORT U_LONG VpInit(U_LONG BaseVal);
VP_EXPORT void *VpMemAlloc(U_LONG mb);
VP_EXPORT void VpFree(Real *pv);
VP_EXPORT Real *VpAlloc(U_LONG mx, char *szVal);
VP_EXPORT int VpAsgn(Real *c,Real *a,int isw);
VP_EXPORT int VpAddSub(Real *c,Real *a,Real *b,int operation);
VP_EXPORT int VpMult(Real *c,Real *a,Real *b);
VP_EXPORT int VpDivd(Real *c,Real *r,Real *a,Real *b);
VP_EXPORT int VpComp(Real *a,Real *b);
VP_EXPORT S_LONG VpExponent10(Real *a);
VP_EXPORT void VpSzMantissa(Real *a,char *psz);
VP_EXPORT void VpToString(Real *a,char *psz,int fFmt);
VP_EXPORT int VpCtoV(Real *a,char *int_chr,U_LONG ni,char *frac,U_LONG nf,char *exp_chr,U_LONG ne);
VP_EXPORT void VpVtoD(double *d,S_LONG *e,Real *m);
VP_EXPORT void VpDtoV(Real *m,double d);
VP_EXPORT void VpItoV(Real *m,S_INT ival);
VP_EXPORT int VpSqrt(Real *y,Real *x);
VP_EXPORT void VpRound(Real *y,Real *x,int sw,int f,int il);
VP_EXPORT void VpFrac(Real *y,Real *x);
VP_EXPORT int VpPower(Real *y,Real *x,S_INT n);
VP_EXPORT void VpPai(Real *y);
VP_EXPORT void VpExp1(Real *y);
VP_EXPORT void VpExp(Real *y,Real *x);
VP_EXPORT void VpSinCos(Real *psin,Real *pcos,Real *x);
VP_EXPORT int VPrint(FILE *fp,char *cntl_chr,Real *a);
/*
* ------------------
* 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 IsWhiteChar(ch) (((ch==' ')||(ch=='\n')||(ch=='\t')||(ch=='\b'))?1:0)
#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((S_LONG)(a)->sign);else (a)->sign=-(short)Abs((S_LONG)(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)->frac[0]=(a)->exponent=(a)->Prec=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) ( ((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) ( ((s)>0)?VpSetPosInf(a):VpSetNegInf(a) )
#ifdef _DEBUG
int VpVarCheck(Real * v);
#endif /* _DEBUG */
#if defined(__cplusplus)
} /* extern "C" { */
#endif
#endif //____BIG_DECIMAL__H____

Просмотреть файл

@ -0,0 +1,767 @@
<HTML>
<HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html">
<style type="text/css"><!--
body {
color: #3f0f0f;
background: #fefeff;
margin-left: 2em; margin-right: 2em;
}
h1 {
color: #ffffff;
background-color: #3939AD;
border-color: #FF00FF;
width: 100%;
border-style: solid;
border-top-width: 0.1em;
border-bottom-width: 0.1em;
border-right: none;
border-left: none;
padding: 0.1em;
font-weight: bold;
font-size: 160%;
text-align: center;
}
h2 {
color: #00007f;
background-color: #e7e7ff;
border-color: #000094;
width: 100%;
border-style: solid;
border-left: none;
border-right: none;
border-top-width: 0.1em;
border-bottom-width: 0.1em;
padding: 0.1em;
font-weight: bold;
font-size: 110%;
}
h3 {
color: #00007f;
padding: 0.2em;
font-size: 110%;
}
h4, h5 {
color: #000000;
padding: 0.2em;
font-size: 100%;
}
table {
margin-top: 0.2em; margin-bottom: 0.2em;
margin-left: 2em; margin-right: 2em;
}
caption {
color: #7f0000;
font-weight: bold;
}
th {
background: #e7e7ff;
padding-left: 0.2em; padding-right: 0.2em;
}
td {
background: #f3f7ff;
padding-left: 0.2em; padding-right: 0.2em;
}
code {
color: #0000df;
}
dt {
margin-top: 0.2em;
}
li {
margin-top: 0.2em;
}
pre
{
BACKGROUND-COLOR: #d0d0d0;
BORDER-BOTTOM: medium none;
BORDER-LEFT: medium none;
BORDER-RIGHT: medium none;
BORDER-TOP: medium none;
LINE-HEIGHT: 100%;
MARGIN: 12px 12px 12px 12px;
PADDING-BOTTOM: 12px;
PADDING-LEFT: 12px;
PADDING-RIGHT: 12px;
PADDING-TOP: 12px;
WHITE-SPACE: pre;
WIDTH: 100%
}
--></style>
<TITLE>BigDecimal:An extension library for Ruby</TITLE>
</HEAD>
<BODY BGCOLOR=#FFFFE0>
<H1>BigDecimal(Variable Precision Floating Library for Ruby)</H1>
<DIV align="right"><A HREF="./bigdecimal_ja.html">Japanese</A></DIV><BR>
BigDecimal is an extension library for the Ruby interpreter.
Using BigDecimal class, you can obtain any number of significant digits in computation.
For the details about Ruby see:<BR>
<UL>
<LI><A HREF="http://www.ruby-lang.org/en/">http://www.ruby-lang.org/en/</A>:Official Ruby page(English).</LI>
<LI><A HREF="http://ruby.freak.ne.jp/">http://ruby.freak.ne.jp/</A>:Ruby informations(Japanese).</LI>
<LI><A HREF="http://kahori.com/ruby/ring/">http://kahori.com/ruby/ring/</A>:Mutually linked pages relating to Ruby(Japanese).
</LI>
</UL>
NOTE:<BR>
This software is provided "AS IS" and without any express or
implied warranties,including,without limitation,the implied
warranties of merchantibility and fitness for a particular
purpose. For the details,see COPYING and README included in this
distribution.
<BR>
<hr>
<H2>Contents</H2>
<UL>
<LI><A HREF="#INTRO">Introduction</LI>
<LI><A HREF="#SPEC">Usage and methods</A></LI>
<LI><A HREF="#UNDEF">Infinity,NaN,Zero</A></LI>
<LI><A HREF="#STRUCT">Internal structure</A></LI>
<LI><A HREF="#BASE">Binary or decimal number representation</A></LI>
<LI><A HREF="#PREC">Resulting number of significant digits</A></LI>
</UL>
<HR>
<A NAME="#INTRO">
<H2>Introduction</H2>
Ruby already has builtin (variable length integer number) class Bignum. Using Bignum class,you can obtain
any integer value in magnitude. But, variable length floating number class is not yet built in.
This is why I made variable length floating class BigDecimal.
Feel free to send any comments or bug reports to me.
<A HREF="mailto:shigeo@tinyforest.gr.jp">shigeo@tinyforest.gr.jp</A>
I will try(but can't promise) to fix bugs reported.
<hr>
<H2>Installation</H2>
The Ruby latest version can be downloaded from <A HREF="http://www.ruby-lang.org/en/">Official Ruby page</A>.
Once decompress the downloaded Ruby archive,follow the normal installation procedures according to the
documents included.
<A NAME="#SPEC">
<H2>Usage and methods</H2>
Suppose you already know Ruby programming,
to create BigDecimal objects,the program would like:<BR>
<CODE><PRE>
require 'bigdecimal'
a=BigDecimal::new("0.123456789123456789")
b=BigDecimal::new("123456.78912345678",40)
c=a+b
</PRE></CODE>
<H3>List of methods</H3>
In the following explanations,n specifies the minimum number of resulting significant digits,
not exactly but slightly excess memories will be allocated to newly created object.
In 32 bits integer system,every 4 digits(in decimal) are computed simultaneously.
This means the number of significant digits in BigDecimal is always a multiple of 4.
<UL>
<LI><B>new</B></LI><BR>
"new" method creates a new BigDecimal object.<BR>
a=BigDecimal::new(s[,n])<BR>
where:<BR>
s: Initial value string.<BR>
n: Maximum number of significant digits of a. n must be a Fixnum object.
If n is omitted or is equal to 0,then the maximum number of significant digits of a is determined from the length of s.
<LI><B>double_fig</B></LI><BR>
double_fig is a class method which returns the number of digits
the Float class can have.
<CODE><PRE>
p BigDecimal::double_fig # ==> 20 (depends on the CPU etc.)
</PRE></CODE>
The equivalent C programs which calculates the value of
double_fig is:
<CODE><PRE>
double v = 1.0;
int double_fig = 0;
while(v + 1.0 > 1.0) {
++double_fig;
v /= 10;
}
</PRE></CODE>
<LI><B>prec</B></LI><BR>
r,m = a.prec<BR>
where r is the number of significant digits of a,
m is the maximum number of significant digits a can hold.<br>
<CODE><PRE>
require "bigdecimal"
a = BigDecimal.new("0.12345")
p a.prec # ==> [8, 12]
b = BigDecimal.new("0.1234500000000")
p b.prec # ==> [8, 20]
c = BigDecimal.new("0.12345",20)
p c.prec # ==> [8, 24]
</PRE></CODE>
r and m are always the multiple of log10(BigDecimal::BASE).
<LI><B>+</B></LI><BR>
addition(c = a + b)<BR>
For the resulting number of significant digits of c,see <A HREF="#PREC">Resulting number of significant digits</A>.
<LI><B>-</B></LI><BR>
subtraction (c = a - b) or negation (c = -a)<BR>
For the resulting number of significant digits of c,see <A HREF="#PREC">Resulting number of significant digits</A>.
<LI><B>*</B></LI><BR>
multiplication(c = a * b)<BR>
For the resulting number of significant digits of c,see <A HREF="#PREC">Resulting number of significant digits</A>.
<LI><B>/</B></LI><BR>
division(c = a / b)<BR>
For the resulting number of significant digits of c,see <A HREF="#PREC">Resulting number of significant digits</A>.
<LI><B>assign</B></LI><BR>
c = a.assign(n,f)<BR>
assigns the value of a to c.<BR>
n is the number of significant digits of resulting c.<BR>
If f > 0,then a is assigned to c.<BR>
If f < 0,then -a is assigned to c.<BR>
The absolute value of f (|f|) must be 1 or 2.
If |f|=2,then proper round operation over c is performed,when the maximum
number of significant digits of c is less than current
number of significant digits of a.
If |f|=1 then extra digits are discarded when the maximum
number of significant digits of c is less than current
number of significant digits of a.
<LI><B>add</B></LI><BR>
c = a.add(b,n)<BR>
c = a.add(b,n) performs c = a + b.
If n is less than the actual significant digits of a + b,
then c is rounded properly.
<LI><B>sub</B></LI><BR>
c = a.sub(b,n)<BR>
c = a.sub(b,n) performs c = a - b.
If n is less than the actual significant digits of a - b,
then c is rounded properly.
<LI><B>mult</B></LI><BR>
c = a.mult(b,n)<BR>
c = a.mult(b,n) performs c = a * b.
If n is less than the actual significant digits of a * b,
then c is rounded properly.
<LI><B>div</B></LI><BR>
c,r = a.div(b,n)<BR>
c,r = a.div(b,n) performs c = a / b, r is the residue of a / b.
If necessary,the divide operation continues to n digits which c
can hold.
Unlike the divmod method,c is not always an integer.
c is never rounded,and the equation a = c*b + r is always
valid unless c is NaN or Infinity.
<LI><B>%</B></LI><BR>
r = a%b <BR>
is the same as:<BR>
r = a-((a/b).floor)*b<BR>
<LI><B>fix</B></LI><BR>
c = a.fix<BR>
returns integer part of a.<BR>
<LI><B>frac</B></LI><BR>
c = a.frac<BR>
returns fraction part of a.<BR>
<LI><B>floor[(n)]</B></LI><BR>
c = a.floor<BR>
returns the maximum integer value (in BigDecimal) which is less than or equal to a.<BR>
As shown in the following example,an optional integer argument (n) specifying the position
of 'floor'ed digit can be given.
If n> 0,then the (n+1)th digit counted from the decimal point in fraction part is 'floor'ed.
If n<0,then the n-th digit counted from the decimal point in integer part is 'floor'ed.<BR>
c = BigDecimal::new("1.23456")<BR>
d = c.floor(4) # d = 1.2345<BR>
c = BigDecimal::new("15.23456")<BR>
d = c.floor(-1) # d = 10.0<BR>
<LI><B>ceil[(n)]</B></LI><BR>
c = a.ceil<BR>
returns the minimum integer value (in BigDecimal) which is greater than or equal to a.<BR>
As shown in the following example,an optional integer argument (n) specifying the position
of 'ceil'ed digit can be given.
If n>0,then the (n+1)th digit counted from the decimal point in fraction part is 'ceil'ed.
If n<0,then the n-th digit counted from the decimal point in integer part is 'ceil'ed.<BR>
c = BigDecimal::new("1.23456")<BR>
d = c.ceil(4) # d = 1.2346<BR>
c = BigDecimal::new("15.23456")<BR>
d = c.ceil(-1) # d = 20.0<BR>
<LI><B>round[(n)]</B></LI><BR>
c = a.round<BR>
round off a to the nearest 1<>D<BR>
As shown in the following example,an optional integer argument (n) specifying the position
of rounded digit can be given.
If n>0,then the (n+1)th digit counted from the decimal point in fraction part is rounded.
If n<0,then the n-th digit counted from the decimal point in integer part is rounded.<BR>
c = BigDecimal::new("1.23456")<BR>
d = c.round(4) # d = 1.235 <BR>
c = BigDecimal::new("15.23456")<BR>
d = c.round(-1) # d = 20.0<BR>
<LI><B>truncate[(n)]</B></LI><BR>
c = a.truncate<BR>
truncate a to the nearest 1<>D<BR>
As shown in the following example,an optional integer argument (n) specifying the position
of truncated digit can be given.
If n>0,then the (n+1)th digit counted from the decimal point in fraction part is truncated.
If n<0,then the n-th digit counted from the decimal point in integer part is truncated.<BR>
c = BigDecimal::new("1.23456")<BR>
d = c.truncate(4) # d = 1.2345<BR>
c = BigDecimal::new("15.23456")<BR>
d = c.truncate(-1) # d = 10.0<BR>
<LI><B>divmod</B></LI><BR>
c,r = a.divmod(b) # a = c*b + r<BR>
returns the quotient and remainder of a/b.<BR>
a = c * b + r is always satisfied.<BR>
where c is the integer sutisfying
c = (a/b).floor <BR>
and,therefore
r = a - c*b<BR>
<LI><B>remainder</B></LI><BR>
r=a.remainder(b)<BR>
returns the remainder of a/b.<BR>
where c is the integer sutisfying
c = (a/b).fix <BR>
and,therefore:
r = a - c*b<BR>
<LI><B>abs</B></LI><BR>
c = a.abs<BR>
returns an absolute value of a.<BR>
<LI><B>to_i</B></LI><BR>
changes a to an integer.<BR>
i = a.to_i<BR>
i becomes to Fixnum or Bignum.
IF a is Infinity or NaN,then i becomes to nil.
<LI><B>to_s[(n)]</B></LI><BR>
converts to string(results look like "0.xxxxxEn").<BR>
s = a.to_s<BR>
If n is given,then a space is inserted after every n digits for readability.<BR>
s = a.to_s(n)
<LI><B>exponent</B></LI><BR>
returns an integer holding exponent value of a.<BR>
n = a.exponent <BR>
means a = 0.xxxxxxx*10**n.
<LI><B>to_f</B></LI><BR>
same as dup method.
creates a new BigDecimal object having same value.
<LI><B>E</B></LI><BR>
e = BigDecimal::E(n)<BR>
where e(=2.718281828....) is the base value of natural logarithm.<BR>
n specifies the length of significant digits of e.
<LI><B>PI</B></LI><BR>
e = BigDecimal::PI(n)<BR>
returns at least n digits of the ratio of the circumference of a circle to its dirmeter
(pi=3.14159265358979....) using J.Machin's formula.<BR>
<LI><B>BASE</B></LI><BR>
Base value used in the BigDecimal calculation.
On 32 bit integer system,the value of BASE is 10000.<BR>
b = BigDecimal::BASE<BR>
<LI><B>mode</B></LI><BR>
mode method controls BigDecimal computation.
Following usage are defined.<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_NaN,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_INFINITY,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_UNDERFLOW,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_OVERFLOW,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_ZERODIVIDE,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_ALL,flag)<BR>
EXCEPTION_NaN controls the execution once computation results to NaN.
EXCEPTION_INFINITY controls the execution once computation results to Infinity(<28>}Infinity).
EXCEPTION_UNDERFLOW controls the execution once computation underflows.
EXCEPTION_OVERFLOW controls the execution once computation overflows.
EXCEPTION_ZERODIVIDE controls the execution once zero-division occures.
EXCEPTION_ALL controls the execution for any exception defined occures.
If the flag is true,then the relating exception is thrown.
No exception is thrown when the flag is false(default) and computation
continues with the result:<BR>
EXCEPTION_NaN results to NaN<BR>
EXCEPTION_INFINITY results to +Infinity or -Infinity<BR>
EXCEPTION_UNDERFLOW results to 0.<BR>
EXCEPTION_OVERFLOW results to +Infinity or -Infinity<BR>
EXCEPTION_ZERODIVIDE results to +Infinity or -Infinity<BR>
EXCEPTION_INFINITY,EXCEPTION_OVERFLOW, and EXCEPTION_ZERODIVIDE are
currently the same.<BR>
The return value of mode method is the value set.
Suppose the return value of the mode method is f,then
f & BigDecimal::EXCEPTION_NaN !=0 means EXCEPTION_NaN is set to on.
If the value of the argument flag is other than nil,true nor false then
current mode status is returned.
<LI><B>limit[(n)]</B></LI><BR>
Limits the maximum digits that the newly created BigDecimal objects can hold
never exceed n. Returns maximum value before set.
Zero,the default value,means no upper limit.<BR>
mf = BigDecimal::limit(n)<BR>
<LI><B>sign</B></LI><BR>
returns the 'attribute'.
n = a.sign <BR>
where the value of n means that a is:<BR>
n = BigDecimal::SIGN_NaN(0) : a is NaN<BR>
n = BigDecimal::SIGN_POSITIVE_ZERO(1) : a is +0<BR>
n = BigDecimal::SIGN_NEGATIVE_ZERO(-1) : a is -0<BR>
n = BigDecimal::SIGN_POSITIVE_FINITE(2) : a is positive<BR>
n = BigDecimal::SIGN_NEGATIVE_FINITE(-2) : a is negative<BR>
n = BigDecimal::SIGN_POSITIVE_INFINITE(3) : a is +Infinity<BR>
n = BigDecimal::SIGN_NEGATIVE_INFINITE(-3) : a is -Infinity<BR>
The value in () is the actual value,see (<A HREF="#STRUCT">Internal structure</A>.<BR>
<LI><B>nan?</B></LI><BR>
a.nan? returns True when a is NaN.
<LI><B>infinite?</B></LI><BR>
a.infinite? returns True when a is +<2B>‡ or -<2D>‡.
<LI><B>finite?</B></LI><BR>
a.finite? returns True when a is neither <20>‡ nor NaN.
<LI><B>to_parts</B></LI><BR>
decomposes a BigDecimal value to 4 parts.
All 4 parts are returned as an array.<BR>
Parts consist of a sign(0 when the value is NaN,+1 for positive and
-1 for negative value), a string representing fraction part,base value(always 10 currently),and an integer(Fixnum) for exponent respectively.
a=BigDecimal::new("3.14159265",10)<BR>
f,x,y,z = a.to_parts<BR>
where f=+1,x="314159265",y=10 and z=1<BR>
therefore,you can translate BigDecimal value to Float as:<BR>
s = "0."+x<BR>
b = f*(s.to_f)*(y**z)<BR>
<LI><B>inspect</B></LI><BR>
is used for debugging output.<BR>
p a=BigDecimal::new("3.14",10)<BR>
should produce output like "#&lt;0x112344:'0.314E1',4(12)%gt;".
where "0x112344" is the address,
'0.314E1' is the value,4 is the number of the significant digits,
and 12 is the maximum number of the significant digits
the object can hold.
<LI><B>dup</B></LI><BR>
creates a new BigDecimal object having same value.
<LI><B>sqrt</B></LI><BR>
c = a.sqrt(n)<BR>
computes square root value of a with significant digit number n at least.<BR>
<LI><B>sincos</B></LI><BR>
computes and returns sine and cosine value of a with significant digit number n at least.<BR>
sin,cos = a.sincos(n)<BR>
<LI><B>exp</B></LI><BR>
c = a.exp(n)<BR>
computes the base of natural logarithm value(e=2.718281828....) powered by a
with significant digit number n at least.<BR>
<LI><B>power</B></LI><BR>
c = a.power(n)<BR>
returns the value of a powered by n(c=a**n).
n must be an integer.<BR>
<LI><B>zero?</B></LI><BR>
c = a.zero?<BR>
returns true if a is equal to 0,otherwise returns false<BR>
<LI><B>nonzero?</B></LI><BR>
c = a.nonzero?<BR>
returns false if a is 0,otherwise returns a itself.<BR>
<LI><B>&lt;=&gt;</B></LI><BR>
c = a &lt;=&gt; b <BR>
returns 0 if a==b,1 if a &gt b,and returns -1 if a &lt b.<BR>
</UL>
Following methods need no explanation.<BR>
<UL>
<LI>==</LI>
<LI>===</LI>
same as ==,used in case statement.
<LI>!=</LI>
<LI>&lt;</LI>
<LI>&lt;=</LI>
<LI>&gt;</LI>
<LI>&gt;=</LI>
</UL>
<HR>
<H3>About 'coerce'</H3>
<B>For the binary operation like A op B:</B>
<DL>
<DT> 1.Both A and B are BigDecimal objects</DT>
<DD> A op B is normally performed.</DD>
<DT> 2.A is the BigDecimal object but B is other than BigDecimal object</DT>
<DD> Operation is performed,after B is translated to correcponding BigDecimal object(because BigDecimal supports coerce method).</DD>
<DT> 3.A is not the BigDecimal object but B is BigDecimal object</DT>
<DD>If A has coerce mthod,then B will translate A to corresponding
BigDecimal object and the operation is performed,otherwise an error occures.</DD>
</DL>
Attention must be paid when a String is to be translated to BigDecimal.
Translation stops without error at the character representing non digit.
For instance,"10XX" is translated to 10,"XXXX" is translated to 0.<BR>
String representing zero or infinity such as "Infinity","+Infinity","-Infinity",and "NaN" can also be translated to BigDecimal unless false is specified by mode method.<BR>
BigDecimal class supports coerce method(for the details about coerce method,see Ruby documentations). This means the most binary operation can be performed if the BigDecimal object is at the left hand side of the operation.<BR><BR>
For example:
<CODE><PRE>
a = BigDecimal.E(20)
c = a * "0.123456789123456789123456789" # A String is changed to BigDecimal object.
</PRE></CODE>
is performed normally.<BR>
But,because String does not have coerce method,the following example can not be performed.<BR>
<CODE><PRE>
a = BigDecimal.E(20)
c = "0.123456789123456789123456789" * a # ERROR
</PRE></CODE>
If you actually have any inconvenience about the error above.
You can define a new class derived from String class,
and define coerce method within the new class.<BR>
<hr>
<A NAME="#UNDEF">
<H2>Infinity,Not a Number(NaN),Zero</H2>
Infinite numbers and NaN can be represented by string writing "+Infinity"(or "Infinity"),"-Infinity",and "NaN" respectively in your program.
Infinite numbers can be obtained by 1.0/0.0(=Infinity) or -1.0/0.0(=-Infinity).
<BR><BR>
NaN(Not a number) can be obtained by undefined computation like 0.0/0.0
or Infinity-Infinity.
Any computation including NaN results to NaN.
Comparisons with NaN never become true,including comparison with NaN itself.
<BR><BR>
Zero has two different variations as +0.0 and -0.0.
But,still, +0.0==-0.0 is true.
<BR><BR>
Computation results including Infinity,NaN,+0.0 or -0.0 become complicated.
Run following program and comfirm the results.
Send me any incorrect result if you find.
<PRE><CODE>
require "bigdecimal"
aa = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
ba = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
opa = %w(+ - * / <=> > >= < == != <=)
for a in aa
for b in ba
for op in opa
x = BigDecimal::new(a)
y = BigDecimal::new(b)
eval("ans= x #{op} y;print a,' ',op,' ',b,' ==> ',ans.to_s,\"\n\"")
end
end
end
</CODE></PRE>
<hr>
<A NAME="#STRUCT">
<H2>Internal structure</H2>
BigDecimal number is defined by the structure Real in BigDecimal.h.
Digits representing a float number are kept in the array frac[] defined in the structure.
In the program,any floating number(BigDecimal number) is represented as:<BR>
<BigDecimal number> = 0.xxxxxxxxx*BASE**n<BR><BR>
where 'x' is any digit representing mantissa(kept in the array frac[]),
BASE is base value(=10000 in 32 bit integer system),
and n is the exponent value.<BR>
Larger BASE value enables smaller size of the array frac[],and increases computation speed.
The value of BASE is defined ind VpInit(). In 32 bit integer system,this value is
10000. In 64 bit integer system,the value becomes larger.
BigDecimal has not yet been compiled and tested on 64 bit integer system.
It will be very nice if anyone try to run BigDecimal on 64 bit system and
inform me the results.
When BASE is 10000,an element of the array frac[] can have vale of from 0 to 9999.
(up to 4 digits).<BR>
The structure Real is defined in bigdecimal.h as:<BR>
<CODE><PRE>
typedef struct {
VALUE obj; /* Back pointer(VALUE) for Ruby object. */
unsigned long MaxPrec; /* The size of the array frac[] */
unsigned long Prec; /* Current size of frac[] actually used. */
short sign; /* Attribute of the value. */
/* ==0 : NaN */
/* 1 : +0 */
/* -1 : -0 */
/* 2 : Positive number */
/* -2 : Negative number */
/* 3 : +Infinity */
/* -3 : -Infinity */
unsigned short flag; /* Control flag */
int exponent; /* Exponent value(0.xxxx*BASE**exponent) */
unsigned long frac[1]; /* An araay holding mantissa(Variable) */
} Real;
</CODE></PRE>
The decimal value 1234.56784321 is represented as(BASE=10000):<BR>
<PRE>
0.1234 5678 4321*(10000)**1
</PRE>
where frac[0]=1234,frac[1]=5678,frac[2]=4321,
Prec=3,sign=2,exponent=1. MaxPrec can be any value greater than or equal to
Prec.
<hr>
<A NAME="#BASE">
<H2>Binary or decimal number representation</H2>
I adopted decimal number representation for BigDecimal implementation.
Of cource,binary number representation is common on the most computers.
<H3>Advantages using decimal representation</H3>
The reason why I adopted decimal number representation for BigDecimal is:<BR>
<DL>
<DT>Easy for debugging
<DD>The floating number 1234.56784321 can be easily represented as:<BR>
frac[0]=1234,frac[1]=5678,frac[2]=4321,exponent=1,and sign=2.
<DT>Exact representation
<DD>Following program can add all numbers(in decimal) in a file
without any error(no round operation).<BR>
<PRE><CODE>
file = File::open(....,"r")
s = BigDecimal::new("0")
while line = file.gets
s = s + line
end
</CODE></PRE>
If the internal representation is binary,translation from decimal to
binary is required and the translation error is inevitable.
For example, 0.1 can not exactly be represented in binary.<BR>
0.1 => b1*2**(-1)+b1*2**(-2)+b3*2**(-3)+b4*2**(-4)....<BR>
where b1=0,b2=0,b3=0,b4=1...<BR>
bn(n=1,2,3,...) is infinite series of digit with value of 0 or 1,
and rounding operation is necessary but where we should round the series ?
Of cource,exact "0.1" is printed if the rouding operation is properly done,
<DT>Significant digit we can have is automatically determined
<DD>In binary representation,0.1 can not be represented in finite series of digit.
But we only need one element(frac[0]=1) in decimal representation.
This means that we can always determine the size of the array frac[] in Real
structure.
</DL>
<H3>Disadvantage of decimal representation</H3>
Advantages stated so far can also be disadvantages if the input from outside is
represented in binary.
Translation error from decimal to binary or vice versa is inevitable.
So,translation from Float(binary) to BigDecimal(decimal) is not alway done exactly.
<H4>Which is the first input?</H4>
Because most people uses decimal notatin for numeric data representation,
BigDecimal can handle numeric data without loss of translation error.
<hr>
<A NAME="#PREC">
<H2>Resulting number of significant digits</H2>
For the fundamental arithmetics such as addition,subtraction,
multiplication,and division,I prepared 2 group of methods<BR>
<H3>1. +,-,*,/</H3>
For the operation + - * /,you can not specify the resulting
number of significant digits.<BR>
Resulting number of significant digits are defined as:<BR>
1.1 For * and /,resulting number of significant digits is the sum of the
significant digits of both side of the operator.<BR>
1.2 For + and -,resulting number of significant digits is determined so that
no round operation is needed. <br>
For example, c has more than 100 siginificant digits if c is computed as:<BR>
c = 0.1+0.1*10**(-100)<br>
<BR>
As +,-,and * are always exact(no round operation is performed),
which means more momories are required to keep computation results.
As for the division as c = a/b,the significant digits of c is the same
as a*b. Division such as c=1.0/3.0 will be rounded.<BR>
<H3>2. assign,add,sub,mult,div</H3>
The length of the significant digits obtained from +,-,*,/
is always defined by that of right and left side of the operator.
To specify the length of the significant digits by your self,
use methos assign,add,sub,mult,div, or limit(class method).
Following example compute the ratio of the circumference of a circle to
its dirmeter(pi=3.14159265358979....) using J.Machin's formula.
<BR><BR>
<CODE><PRE>
#!/usr/local/bin/ruby
#
# pai.rb
# USAGE: ruby pai.rb n
# where n is the number of digits required.
# EX.: ruby pai.rb 1000
#
require "bigdecimal"
#
# Calculates 3.1415.... using J. Machin's formula.
#
def pai(sig) # sig: Number of significant figures
exp = -sig
pi = BigDecimal::new("0")
two = BigDecimal::new("2")
m25 = BigDecimal::new("-0.04")
m57121 = BigDecimal::new("-57121")
u = BigDecimal::new("1")
k = BigDecimal::new("1")
w = BigDecimal::new("1")
t = BigDecimal::new("-80")
while (u.exponent >= exp)
t = t*m25
u,r = t.div(k,sig)
pi = pi + u
k = k+two
end
u = BigDecimal::new("1")
k = BigDecimal::new("1")
w = BigDecimal::new("1")
t = BigDecimal::new("956")
while (u.exponent >= exp )
t,r = t.div(m57121,sig)
u,r = t.div(k,sig)
pi = pi + u
k = k+two
end
pi
end
if $0 == __FILE__
print "PAI("+ARGV[0]+"):\n"
p pai(ARGV[0].to_i)
end
</PRE></CODE>
<HR>
<FONT size=2>
<I>
<A HREF="http://www.tinyforest.gr.jp">
Shigeo Kobayashi
</A>
(E-Mail:<A HREF="mailto:shigeo@tinyforest.gr.jp">&lt;shigeo@tinyforest.gr.jp&gt;</U></A>)
</I>
</FONT>
</TD>
</TR>
</TABLE>
</BODY>
</HTML>

Просмотреть файл

@ -0,0 +1,706 @@
<HTML>
<HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=Shift_JIS">
<style type="text/css"><!--
body {
color: #3f0f0f;
background: #fefeff;
margin-left: 2em; margin-right: 2em;
}
h1 {
color: #ffffff;
background-color: #3939AD;
border-color: #FF00FF;
width: 100%;
border-style: solid;
border-top-width: 0.1em;
border-bottom-width: 0.1em;
border-right: none;
border-left: none;
padding: 0.1em;
font-weight: bold;
font-size: 160%;
text-align: center;
}
h2 {
color: #00007f;
background-color: #e7e7ff;
border-color: #000094;
width: 100%;
border-style: solid;
border-left: none;
border-right: none;
border-top-width: 0.1em;
border-bottom-width: 0.1em;
padding: 0.1em;
font-weight: bold;
font-size: 110%;
}
h3 {
color: #00007f;
padding: 0.2em;
font-size: 110%;
}
h4, h5 {
color: #000000;
padding: 0.2em;
font-size: 100%;
}
table {
margin-top: 0.2em; margin-bottom: 0.2em;
margin-left: 2em; margin-right: 2em;
}
caption {
color: #7f0000;
font-weight: bold;
}
th {
background: #e7e7ff;
padding-left: 0.2em; padding-right: 0.2em;
}
td {
background: #f3f7ff;
padding-left: 0.2em; padding-right: 0.2em;
}
code {
color: #0000df;
}
dt {
margin-top: 0.2em;
}
li {
margin-top: 0.2em;
}
pre
{
BACKGROUND-COLOR: #d0d0d0;
BORDER-BOTTOM: medium none;
BORDER-LEFT: medium none;
BORDER-RIGHT: medium none;
BORDER-TOP: medium none;
LINE-HEIGHT: 100%;
MARGIN: 12px 12px 12px 12px;
PADDING-BOTTOM: 12px;
PADDING-LEFT: 12px;
PADDING-RIGHT: 12px;
PADDING-TOP: 12px;
WHITE-SPACE: pre;
WIDTH: 100%
}
--></style>
<TITLE>BigDecimal:An extension library for Ruby</TITLE>
</HEAD>
<BODY BGCOLOR=#FFFFE0>
<H1>BigDecimal(可変長浮動少数点演算用拡張ライブラリ)</H1>
<DIV align="right"><A HREF="./bigdecimal_en.html">English</A></DIV><BR>
BigDecimal はオブジェクト指向の強力なスクリプト言語である Ruby に可変長浮動小数点
計算機能を追加するための拡張ライブラリです。
Ruby についての詳しい内容は以下のURLを参照してください。
<UL>
<LI><A HREF="http://www.ruby-lang.org/ja/">http://www.ruby-lang.org/ja/</A>Ruby公式ページ</LI>
<LI><A HREF="http://ruby.freak.ne.jp/">http://ruby.freak.ne.jp/</A>Rubyに関する情報ページ</LI>
<LI><A HREF="http://kahori.com/ruby/ring/">http://kahori.com/ruby/ring/</A>Rubyに関するページを辿れます</LI>
</UL>
<hr>
<H2>目次</H2>
<UL>
<LI><A HREF="#INTRO">はじめに</LI>
<LI><A HREF="#SPEC">使用方法とメソッドの一覧</A></LI>
<LI><A HREF="#UNDEF">無限、非数、ゼロの扱い</A></LI>
<LI><A HREF="#STRUCT">内部構造</A></LI>
<LI><A HREF="#BASE">2進と10進</A></LI>
<LI><A HREF="#PREC">計算精度について</A></LI>
</UL>
<HR>
<A NAME="#INTRO">
<H2>はじめに</H2>
Ruby には Bignum というクラスがあり、数百桁の整数でも計算することができます。
ただ、任意桁の浮動少数点演算用クラスが無いようです。そこで、
任意桁の浮動少数点演算用拡張ライブラリ BigDecimal を作成しました。
不具合や助言・提案がある場合どしどし、
<A HREF="mailto:shigeo@tinyforest.gr.jp">shigeo@tinyforest.gr.jp</A>
までお知らせください。不具合を直す気は大いにあります。ただ、時間などの関係で約束
はできません。また、結果についても保証できるものではありません。
予め、ご了承ください。
<BR><BR>
このプログラムは、自由に配布・改変して構いません。ただし、著作権は放棄していません。
配布・改変等の権利は Ruby のそれに準じます。詳しくは README を読んでください。
<hr>
<H2>インストールについて</H2>
BigDecimal を含む Ruby の最新版は<A HREF="http://www.ruby-lang.org/ja/">Ruby公式ページ</A>からダウンロードできます。
ダウンロードした最新版を解凍したら、通常のインストール手順を実行して下さい。
Ruby が正しくインストールされれば、同時に BigDecimal も利用できるようになるはずです。
ソースファイルは
bigdecimal.c,bigdecimal.h
の2個のみです。<BR>
<hr>
<A NAME="#SPEC">
<H2>使用方法とメソッドの一覧</H2>
「Rubyは既に書ける」という前提で、<br><br>
<CODE>
<PRE>
require 'bigdecimal'
a=BigDecimal::new("0.123456789123456789")
b=BigDecimal::new("123456.78912345678",40)
c=a+b
</PRE>
</CODE>
<br>
というような感じで使用します。
<H3>メソッド一覧</H3>
以下のようなメソッドが利用可能です。<BR>
記述上、BigDecimal オブジェクトを a,b,c,rで、String文字列オブジェクトを
s、整数を n で表記します。また、「有効桁数」とは BigDecimal が精度を保証する
桁数です。ぴったりではありません、若干の余裕を持って計算されます。また、
例えば32ビットのシステムでは10進で4桁毎に計算します。従って、現状では、
内部の「有効桁数」は4の倍数となっています。
<UL>
<LI>new</LI><BR>
新しい BigDecimal オブジェクトを生成します。<BR>
a=BigDecimal::new(s[,n])<BR>
s は初期値を文字列で指定します.
n は必要な有効桁数a の最大有効桁数)を整数で指定します。
n が 0 または省略されたときは、n の値は s の有効桁数とみなされます。
s の有効桁数より n が小さいときも n=0 のときと同じです。
a の最大有効桁数は n より若干大い値が採用されます。
<LI>+</LI><BR>
加算c = a + b<BR>
c の精度については「<A HREF="#PREC">計算精度について</A>」を参照してください。
<LI>-</LI><BR>
減算c = a - b、または符号反転c = -a<BR>
c の精度については「<A HREF="#PREC">計算精度について</A>」を参照してください。
<LI>*</LI><BR>
乗算(c = a * b)<BR>
cの精度は(aの精度)+(bの精度)程度です。<br>
詳しくは「<A HREF="#PREC">計算精度について</A>」を参照してください。
<LI>/</LI><BR>
除算(c = a / b)<BR>
c の精度については「<A HREF="#PREC">計算精度について</A>」を参照してください。
<LI>assign</LI><BR>
以下のように使用します。<BR>
c = a.assign(n,f)<BR>
f > 0 なら、a を c に、そのまま代入します。
f < 0 なら-a c に代入します
f の絶対値(|f|)は1か2を指定してください。
|f|=2 のときは、c の最大精度が a の実精度より小さいときには
丸められます。|f|=1 のときは切り捨てられます。
n は c の有効桁数ですn 桁以上の精度を持つ c が生成されます)。
<LI>add</LI><BR>
以下のように使用します。<BR>
c = a.add(b,n)<BR>
c = a + b を最大で n 桁まで計算します。
a + b の精度が n より大きいときは丸められます。
<LI>sub</LI><BR>
以下のように使用します。<BR>
c = a.sub(b,n)<BR>
c = a - b を最大で n 桁まで計算します。
a - b の精度が n より大きいときは丸められます。
<LI>mult</LI><BR>
以下のように使用します。<BR>
c = a.mult(b,n)<BR>
c = a * b を最大で n 桁まで計算します。
a * b の精度が n より大きいときは丸められます。
<LI>div</LI><BR>
以下のように使用します。<BR>
c,r = a.div(b,n)<BR>
c=a/b の計算をします。 r には剰余が代入されます。a/bは
必要ならn 桁まで計算されます。divmod メソッド
と異なり、c は整数とは限りません。
また、 c は丸められることはありません。
a = c*b + r の関係は成立します。
<LI>%</LI><BR>
r = a%b <BR>
a/b の余りを計算します。以下の計算と同じものです。<BR>
r = a-((a/b).floor)*b<BR>
<LI>fix</LI><BR>
a の小数点以下の切り捨て。<BR>
c = a.fix
<LI>frac</LI><BR>
a の整数部分の切り捨て。<BR>
c = a.frac
<LI>floor</LI><BR>
a 以下の最大整数を表す値BigDecimal 値)を返します。<BR>
c = a.floor<BR>
以下のように引数を与えて、小数点以下 n+1 位の数字を操作することもできます
(少数点以下を、最大 n 桁にします)。<BR>
c = BigDecimal("1.23456")<BR>
d = c.floor(4) # d = 1.2345 になります。<BR>
n が負のときは小数点以上 n 桁目を操作します。<BR>
c = BigDecimal("15.23456")<BR>
d = c.floor(-1) # d = 10.0 になります。<BR>
<LI>ceil</LI><BR>
a 以上の整数のうち、最も小さい整数を計算し、その値BigDecimal 値)を返します。<BR>
c = a.ceil<BR>
以下のように引数を与えて、小数点以下 n+1 位の数字を操作することもできます
(少数点以下を、最大 n 桁にします)。<BR>
c = BigDecimal::new("1.23456")<BR>
d = c.ceil(4) # d = 1.2346 になります。<BR>
n が負のときは小数点以上 n 桁目をを操作します。<BR>
c = BigDecimal::new("15.23456")<BR>
d = c.ceil(-1) # d = 20.0 になります。<BR>
<LI>round</LI><BR>
小数点以下第一位の数を四捨五入して整数BigDecimal 値)にします。<BR>
c = a.round<BR>
以下のように引数を与えて、小数点以下 n+1 位の数字を操作することもできます
(少数点以下を、最大 n 桁にします)。<BR>
n が正の時は、小数点以下 n+1 位の数字を四捨五入します。
c = BigDecimal::new("1.23456")<BR>
d = c.round(4) # d = 1.235 になります。<BR>
n が負のときは小数点以上 n 桁目をを操作します。<BR>
c = BigDecimal::new("15.23456")<BR>
d = c.round(-1) # d = 20.0 になります。<BR>
<LI>truncate</LI><BR>
小数点以下の数を切り捨てて整数BigDecimal 値)にします。<BR>
c = a.truncate<BR>
以下のように引数を与えて、小数点以下 n+1 位の数字を操作することもできます
(少数点以下を、最大 n 桁にします)。<BR>
n が正の時は、小数点以下 n+1 位の数字を切り捨てます。
c = BigDecimal::new("1.23456")<BR>
d = c.truncate(4) # d = 1.2345 になります。<BR>
n が負のときは小数点以上 n 桁目をを操作します。<BR>
c = BigDecimal::new("15.23456")<BR>
d = c.truncate(-1) # d = 10.0 になります。<BR>
<LI>divmod</LI><BR>
商と剰余の配列を返します。<BR>
c,r = a.divmod(b) # a = c*b + r<BR>
divmodメソッドは a = c * b + r となる a / b の浮動小数点型の商 c と剰余 r を
計算します。ここで c は整数(少数部分のない実数)になります。<BR>
c = (a/b).floor <BR>
r = a - c*b<BR>
で計算されます。
<LI>remainder</LI><BR>
r=a.remainder(b)<BR>
a/b の剰余 r を計算します。<BR>
c = (a/b).fix <BR>
r = a - c*b<BR>
で計算されます。
<LI>abs</LI><BR>
aの絶対値<BR>
c = a.abs<BR>
<LI>to_i</LI><BR>
少数点以下を切り捨てて整数に変換します。<BR>
i = a.to_i<BR>
i は値に応じて Fixnum か Bignum になります。
a が Infinity や NaN のとき、i は nil になります。
<LI>to_f</LI><BR>
dup と全く同じです。
同じ値の BigDecimal オブジェクトを生成します。
<LI>to_s</LI><BR>
文字列に変換します("0.xxxxxEn"の形になります)。<BR>
s = a.to_s
<LI>to_s2</LI><BR>
文字列に変換します。仮数部分を n 桁毎に空白で区切ります。<BR>
s = a.to_s2(n)
<LI>exponent</LI><BR>
指数部を整数値で返します。
n = a.exponent <BR>
は a の値が 0.xxxxxxx*10**n を意味します。
<LI>E</LI><BR>
自然対数の底e(=2.718281828....)を計算します(正直にテイラー展開で)。<BR>
e = BigDecimal::E(n)<BR>
nは必要な有効桁数を整数で指定します。
<LI>PI</LI><BR>
円周率(=3.14159265358979....)を計算しますJMachinの公式を用います<BR>
e = BigDecimal::PI(n)<BR>
n は必要な有効桁数を整数で指定します。
<LI>BASE</LI><BR>
内部で使用される基数の値です。整数が 32 ビットの処理系では10000です。<BR>
b = BigDecimal::BASE<BR>
<LI>mode</LI><BR>
BigDecimalの実行結果を制御します。以下の使用方法が定義されています。<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_NaN,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_INFINITY,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_UNDERFLOW,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_OVERFLOW,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_ZERODIVIDE,flag)<BR>
f = BigDecimal::mode(BigDecimal::EXCEPTION_ALL,flag)<BR>
EXCEPTION_NaN は結果が NaN になったときの指定です。
EXCEPTION_INFINITY は結果が無限大(±Infinity)
になったときの指定です。
EXCEPTION_UNDERFLOW は指数部がアンダーフローするときの指定です。
EXCEPTION_OVERFLOW は指数部がオーバーフローするときの指定です。
EXCEPTION_ZERODIVIDE はゼロによる割り算を実行したときの指定です。
EXCEPTION_ALL は、可能な全てに対して一括して設定するときに
使用します。
flag が true のときは、指定した状態になったときに例外を発行
するようになります。
flag が falseデフォルトなら、例外は発行されません。計算結果は
以下のようになります。<BR>
EXCEPTION_NaN のとき、非数(NaN)<BR>
EXCEPTION_INFINITY のとき、無限(+ or -Infinity)<BR>
EXCEPTION_UNDERFLOW のとき、ゼロ<BR>
EXCEPTION_OVERFLOW のとき、+Infinity か -Infinity<BR>
EXCEPTION_ZERODIVIDE のとき、+Infinity か -Infinity<BR>
EXCEPTION_INFINITY、EXCEPTION_OVERFLOW、EXCEPTION_ZERODIVIDE
は今のところ同じです。<BR>
戻り値は、設定後の値です。「値」の意味は、例えば
BigDecimal::EXCEPTION_NaNと「値」の & が ゼロ以外ならば
EXCEPTION_NaNが設定されているという意味です。
flag が nil、または、true と false 以外なら現在の設定値が返ります。
<LI>limit([n])</LI><BR>
生成されるBigDecimalオブジェクトの最大桁数をn桁に制限します。戻り値は
設定する前の値です。設定値のデフォルト値は0で、桁数無制限という意味です。
nを指定しない場合は、現状の最大桁数が返ります。<BR>
mf = BigDecimal::limit(n)<BR>
<LI>sign</LI><BR>
値の属性を返します。
n = a.sign <BR>
としたとき n の値は a が以下のときを意味します。<BR>
() の中の数字は、実際の値です(<A HREF="#STRUCT">「内部構造」</A>を参照)。<BR>
n = BigDecimal::SIGN_NaN(0) : a は NaN<BR>
n = BigDecimal::SIGN_POSITIVE_ZERO(1) : a は +0<BR>
n = BigDecimal::SIGN_NEGATIVE_ZERO(-1) : a は -0<BR>
n = BigDecimal::SIGN_POSITIVE_FINITE(2) : a は正の値<BR>
n = BigDecimal::SIGN_NEGATIVE_FINITE(-2) : a は負の値<BR>
n = BigDecimal::SIGN_POSITIVE_INFINITE(3) : a は+Infinity<BR>
n = BigDecimal::SIGN_NEGATIVE_INFINITE(-3) : a は-Infinity<BR>
<LI>nan?</LI><BR>
a.nan? は a がNaNのとき真を返します。
<LI>infinite?</LI><BR>
a.infinite? は a が+∞または-∞のとき真を返します。
<LI>finite?</LI><BR>
a.finite? は a が∞または NaN でないとき真を返します。
<LI>to_parts</LI><BR>
BigDecimal 値を 0.xxxxxxx*10**n と表現したときに、符号NaNのときは
0、それ以外は+1か-1になります
仮数部分の文字列("xxxxxxx"と、基数10、更に指数 n を配列で
返します。<BR>
a=BigDecimal::new("3.14159265",10)<BR>
f,x,y,z = a.to_parts<BR>
とすると、f=+1、x="314159265"、y=10、z=1になります。<BR>
従って、<BR>
s = "0."+x<BR>
b = f*(s.to_f)*(y**z)<BR>
で Float に変換することができます。
<LI>inspect</LI><BR>
デバッグ出力に使用されます。<BR>
p a=BigDecimal::new("3.14",10)<BR>
とすると、[0x112344:'0.314E1',4(12)]のように出力されます。
最初の16進数はオブジェクトのアドレス、次の '0.314E1' は値、
次の4は現在の有効桁数(表示より若干大きいことがあります)、
最後はオブジェクトが取り得る最大桁数になります。
<LI>dup</LI><BR>
同じ値の BigDecimal オブジェクトを生成します。
<LI>sqrt</LI><BR>
aの有効桁 n 桁の平方根n の平方根ではありません)。
これまた、正直にニュートン法で計算します。<BR>
c = a.sqrt(n)<BR>
<LI>sincos</LI><BR>
a の有効桁 n 桁の sin と cos を同時に(テイラー展開で)計算して、
sin と cos の配列を返します。
n は必要な有効桁数です( n の sin や cos を計算するわけではありません)。
<BR>
sin,cos = a.sincos(n)<BR>
<LI>exp</LI><BR>
自然対数の底e(=2.718281828....)の a 乗を計算します。<BR>
c = a.exp(n)<BR>
n は必要な有効桁数です。
<LI>power</LI><BR>
a の n 乗を計算します。nは整数。<BR>
c = a.power(n)<BR>
結果として c の有効桁は a の n 倍以上になるので注意。
<LI>zero?</LI><BR>
a が 0 なら true になります。<BR>
c = a.zero?<BR>
<LI>nonzero?</LI><BR>
a が 0 なら false、0 以外なら a そのものが返ります。<BR>
c = a.nonzero?<BR>
<LI>&lt=&gt</LI><BR>
a==b なら 0、a &gt b なら 1、a &lt b なら -1 になります。<BR>
c = a &lt=&gt b <BR>
</UL>
後は、読んで字の如くです。<BR>
<UL>
<LI>==</LI>
<LI>===</LI>
「==」と同じですが case 文で使用されます。
<LI>!=</LI>
<LI>&lt</LI>
<LI>&lt=</LI>
<LI>&gt</LI>
<LI>&gt=</LI>
</UL>
<H3>coerceについて</H3>
BigDecimal オブジェクトが算術演算子の左にあるときは、BigDecimal オブジェクトが
右にあるオブジェクトを(必要なら) BigDecimal に変換してから計算します。
従って、BigDecimal オブジェクト以外でも数値を意味するものなら右に置けば
演算は可能です。<BR><BR>
文字列で数値を与える場合は注意が必要です。数値に変換できない文字があると、
単に変換を止めるだけでエラーにはなりません。"10XX"なら10、"XXXX"は0
と扱われます。<BR>
<CODE><PRE>
a = BigDecimal.E(20)
c = a * "0.123456789123456789123456789" # 文字を BigDecimal に変換してから計算
</PRE></CODE>
無限大や非数を表す文字として、"Infinity"、"+Infinity"、"-Infinity"、"NaN"
も使用できます(大文字・小文字を区別します)。ただし、mode メソッドで false を
指定した場合は例外が発生します。
<BR>
また、BigDecimalクラスは coerceRuby本参照をサポートしています。
従って、BigDecimal オブジェクトが右にある場合も大抵は大丈夫です。
ただ、現在の Ruby インタプリタの仕様上、文字列が左にあると計算できません。<BR>
<CODE><PRE>
a = BigDecimal.E(20)
c = "0.123456789123456789123456789" * a # エラー
</PRE></CODE>
必要性があるとは思いませんが、どうしてもと言う人は
String オブジェクトを継承した新たなクラスを作成してから、
そのクラスで coerce をサポートしてください。
<hr>
<A NAME="#UNDEF">
<H2>無限、非数、ゼロの扱い</H2>
「無限」とは表現できないくらい大きな数です。特別に扱うために
+Infinity正の無限大や -Infinity負の無限大という
ように表記されます。
無限は 1.0/0.0 のようにゼロで割るような計算をしたときに生成されます。
<BR><BR>
「非数」は 0.0/0.0 や Infinity-Infinity 等の結果が定義できない
計算をしたときに生成されます。非数は NaNNot a Numberと表記されます。
NaN を含む計算は全て NaN になります。また NaN は自分も含めて、どんな数
とも一致しません。
<BR><BR>
ゼロは +0.0 と -0.0 が存在します。ただし、+0.0==-0.0 は true です。
<BR><BR>
Infinity、NaN、 +0.0 と -0.0 等を含んだ計算結果は組み合わせに
より複雑です。興味のある人は、以下のプログラムを実行して結果を
確認してください(結果について、疑問や間違いを発見された方は
お知らせ願います)。
<PRE>
<CODE>
require "bigdecimal"
aa = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
ba = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
opa = %w(+ - * / <=> > >= < == != <=)
for a in aa
for b in ba
for op in opa
x = BigDecimal::new(a)
y = BigDecimal::new(b)
eval("ans= x #{op} y;print a,' ',op,' ',b,' ==> ',ans.to_s,\"\n\"")
end
end
end
</CODE>
</PRE>
<hr>
<A NAME="#STRUCT">
<H2>内部構造</H2>
BigDecimal内部で浮動小数点は構造体(Real)で表現されます。
そのうち仮数部は unsigned long の配列(以下の構造体要素frac)で管理されます。
概念的には、以下のようになります。<BR><BR>
<浮動小数点数> = 0.xxxxxxxxx*BASE**n<BR><BR>
ここで、xは仮数部を表す数字、BASEは基数進なら、nは指数部を表す
整数値です。BASEが大きいほど、大きな数値が表現できます。つまり、配列のサイズを
少なくできます。BASEは大きいほど都合がよいわけですが、デバッグのやりやすさなどを
考慮して、10000になっていますBASEはVpInit()関数で自動的に計算します)。
これは、32ビット整数の場合です。64ビット整数の場合はもっと大きな値になります。
残念ながら、64ビット整数でのテストはまだやっていませんもし、やられた方がいれば
結果を教えていただければありがたいです)。
BASEが10000のときは、以下の仮数部の配列(frac)の各要素には最大で4桁の
数字が格納されます。<BR><BR>
浮動小数点構造体(Real)は以下のようになっています。
<BR>
<CODE><PRE>
typedef struct {
unsigned long MaxPrec; // 最大精度(frac[]の配列サイズ)
unsigned long Prec; // 精度(frac[]の使用サイズ)
short sign; // 以下のように符号等の状態を定義します。
// ==0 : NaN
// 1 : +0
// -1 : -0
// 2 : 正の値
// -2 : 負の値
// 3 : +Infinity
// -3 : -Infinity
unsigned short flag; // 各種の制御フラッグ
int exponent; // 指数部の値(仮数部*BASE**exponent)
unsigned long frac[1]; // 仮数部の配列(可変)
} Real;
</CODE></PRE>
例えば 1234.56784321 という数字は(BASE=10000なら)<BR>
<PRE>
0.1234 5678 4321*(10000)**1
</PRE>
ですから frac[0]=1234、frac[1]=5678、frac[2]=4321、
Prec=3、sign=2、exponent=1 となります。MaxPrecは
Prec より大きければいくつでもかまいません。flag の
使用方法は実装に依存して内部で使用されます。
<hr>
<A NAME="#BASE">
<H2>2進と10進</H2>
BigDecimal は <浮動小数点数> = 0.xxxxxxxxx*10**n という10進形式で数値を保持します。
しかし、計算機の浮動小数点数の内部表現は、言うまでもなく <浮動小数点数> = 0.bbbbbbbb*2**n という
2進形式が普通です(x は 0 から 9 まで、b は 0 か 1 の数字)。
BigDecimal がなぜ10進の内部表現形式を採用したのかを以下に説明します。
<H4>10進のメリット</H4>
<DL>
<DT>デバッグのしやすさ
<DD>まず、プログラム作成が楽です。frac[0]=1234、frac[1]=5678、frac[2]=4321、
exponent=1、sign=2 なら数値が 1234.56784321 であるのは見れば直ぐに分かります。
<DT>10進表記された数値なら確実に内部表現に変換できる
<DD>例えば、以下のようなプログラムは全く誤差無しで
計算することができます。以下の例は、一行に一つの数値
が書いてあるファイル file の合計数値を求めるものです。
<PRE><CODE>
file = File::open(....,"r")
s = BigDecimal::new("0")
while line = file.gets
s = s + line
end
</CODE></PRE>
この例を2進数でやると誤差が入り込む可能性があります。
例えば 0.1 を2進で表現すると 0.1 = b1*2**(-1)+b1*2**(-2)+b3*2**(-3)+b4*2**(-4)....
と無限に続いてしまいます(b1=0,b2=0,b3=0,b4=1...)。ここで bn(n=1,2,3,...) は
2進を表現する 0 か 1 の数字列です。従って、どこかで打ち切る必要があります。
ここで変換誤差が入ります。もちろん、これを再度10進表記にして印刷するような
場合は適切な丸め操作(四捨五入)によって再び "0.1" と表示されます。しかし、
内部では正確な 0.1 ではありません。
<DT>有効桁数は有限である(つまり自動決定できる)
<DD>0.1 を表現するための領域はたった一つの配列要素( frac[0]=1 )で済みます。
配列要素の数は10進数値から自動的に決定できます。これは、可変長浮動小数点演算では
大事なことです。逆に 0.1 を2進表現したときには2進の有効桁をいくつにするのか 0.1 を
見ただけでは決定できません。
</DL>
<H3>10進のデメリット</H3>
実は今までのメリットは、そのままデメリットにもなります。
そもそも、10進を2進、2進を10進に変換するような操作は変換誤差
を伴う場合を回避することはできません。
既に計算機内部に取り込まれた2進数値を BigDecimal の内部表現に
変換するときには誤差が避けられない場合があります。
<H3>最初は何か?</H3>
自分で計算するときにわざわざ2進数を使う人は極めてまれです。
計算機にデータを入力するときもほとんどの場合、
10進数で入力します。その結果、double 等の計算機内部
表現は最初から誤差が入っている場合があります。
BigDecimal はユーザ入力を誤差無しで取り込むことができます。
デバッグがしやすいのと、データ読みこみ時に誤差が入らない
というのが実際のメリットです。
<hr>
<A NAME="#PREC">
<H2>計算精度について</H2>
c = a op b という計算(op は + - * /)をしたときの動作は
以下のようになります。<BR><BR>
1.乗算と除算は(a の有効桁数)+(a の有効桁数)分の最大桁数(実際は、余裕を持って、
もう少し大きくなります)を持つ変数 c を新たに生成します。
加減算の場合は、誤差が出ないだけの精度を持つ c を生成します。例えば
c = 0.1+0.1*10**(-100) のような場合、c の精度は100桁以上の精度を
持つようになります。
<BR>
2.次に c = a op b の計算を実行します。<BR><BR>
このように、加減算と乗算での c は必ず「誤差が出ない」だけの精度を
持って生成されます。除算は(a の有効桁数)+(a の有効桁数)分の最大桁数
を持つ c が生成されますが、c = 1.0/3.0 のような計算で明らかなように、
c の最大精度を超えるところで計算が打ち切られる場合があります。<BR><BR>
いずれにせよ、c の最大精度は a や b より大きくなりますので c が必要とする
メモリー領域は大きくなることに注意して下さい。
<BR><BR>
注意:「+,-,*,/」では結果の精度(有効桁数)を自分で指定できません。
精度をコントロールしたい場合は、以下の add,sub 等のメソッド
を使用します。<BR>
<H3>自分で精度をコントロールしたい場合</H3>
自分で精度(有効桁数)をコントロールしたい場合は assign、add、sub、mult、div 等のメソッド
が使用できます。
以下の円周率を計算するプログラム例のように、
求める桁数は自分で指定することができます。
<BR><BR>
<CODE><PRE>
#!/usr/local/bin/ruby
#
# pai.rb
# USAGE: ruby pai.rb n
# where n is the number of digits required.
# EX.: ruby pai.rb 1000
#
require "bigdecimal"
#
# Calculates 3.1415.... using J. Machin's formula.
#
def pai(sig) # sig: Number of significant figures
exp = -sig
pi = BigDecimal::new("0")
two = BigDecimal::new("2")
m25 = BigDecimal::new("-0.04")
m57121 = BigDecimal::new("-57121")
u = BigDecimal::new("1")
k = BigDecimal::new("1")
w = BigDecimal::new("1")
t = BigDecimal::new("-80")
while (u.exponent >= exp)
t = t*m25
u,r = t.div(k,sig)
pi = pi + u
k = k+two
end
u = BigDecimal::new("1")
k = BigDecimal::new("1")
w = BigDecimal::new("1")
t = BigDecimal::new("956")
while (u.exponent >= exp )
t,r = t.div(m57121,sig)
u,r = t.div(k,sig)
pi = pi + u
k = k+two
end
pi
end
if $0 == __FILE__
print "PAI("+ARGV[0]+"):\n"
p pai(ARGV[0].to_i)
end
</PRE></CODE>
<HR>
<FONT size=2>
<I>
<A HREF="http://www.tinyforest.gr.jp">
小林 茂雄
</A>
(E-Mail:<A HREF="mailto:shigeo@tinyforest.gr.jp">&ltshigeo@tinyforest.gr.jp&gt</U></A>)
</I>
</FONT>
</TD>
</TR>
</TABLE>
</BODY>
</HTML>

1
ext/bigdecimal/depend Normal file
Просмотреть файл

@ -0,0 +1 @@
bigdecimal.o: bigdecimal.c bigdecimal.h $(hdrdir)/ruby.h

Просмотреть файл

@ -0,0 +1,2 @@
require 'mkmf'
create_makefile('bigdecimal')

Просмотреть файл

@ -0,0 +1,30 @@
#
# BigDecimal <-> Rational
#
class BigDecimal
# Convert BigDecimal to Rational
def to_r
sign,digits,base,power = self.to_parts
numerator = sign*digits.to_i
denomi_power = power - digits.size # base is always 10
if denomi_power < 0
denominator = base ** (-denomi_power)
else
denominator = base ** denomi_power
end
Rational.new(numerator,denominator)
end
end
class Rational
# Convert Rational to BigDecimal
# to_d returns an array [quotient,residue]
def to_d(nFig=0)
num = self.numerator.to_s
if nFig<=0
nFig = BigDecimal.double_fig*2+1
end
BigDecimal.new(num).div(self.denominator,nFig)
end
end

7
ext/bigdecimal/lib/delcr Normal file
Просмотреть файл

@ -0,0 +1,7 @@
tr -d '\r' < bigdecimal-rational.rb > unix/bigdecimal-rational.rb
tr -d '\r' <jacobian.rb > unix/jacobian.rb
tr -d '\r' <linear.rb > unix/linear.rb
tr -d '\r' <ludcmp.rb > unix/ludcmp.rb
tr -d '\r' <newton.rb > unix/newton.rb
tr -d '\r' <nlsolve.rb > unix/nlsolve.rb
tr -d '\r' <pai.rb > unix/pai.rb

Просмотреть файл

@ -0,0 +1,63 @@
#
# jacobian.rb
#
# Computes Jacobian matrix of f at x
#
module Jacobian
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
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
s = f.zero
deriv = []
if(nRetry>100) then
raize "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
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,46 @@
#!/usr/local/bin/ruby
#
# 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.
#
require "bigdecimal"
require "ludcmp"
include LUSolve
def rd_order
printf("Number of equations ?")
n = gets().chomp.to_i
end
zero = BigDecimal::new("0.0")
one = BigDecimal::new("1.0")
while (n=rd_order())>0
a = []
as= []
b = []
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 = gets
a <<=BigDecimal::new(s);
as<<=BigDecimal::new(s);
end
printf("Contatant vector element b[%d] ? ",i);b<<=BigDecimal::new(gets);
end
printf "ANS="
x = lusolve(a,b,ludecomp(a,n,zero,one),zero)
p x
printf "A*x-b\n"
for i in 0...n do
s = zero
for j in 0...n do
s = s + as[i*n+j]*x[j]
end
p s-b[i]
end
end

Просмотреть файл

@ -0,0 +1,75 @@
#
# ludcmp.rb
#
module LUSolve
def ludecomp(a,n,zero=0.0,one=1.0)
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/nrmrow
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]/pivot
if mult!=zero then
pskn = ps[k]*n
for j in (k+1)...n do
a[psin+j] -= mult*a[pskn+j]
end
end
end
end
raise "Singular matrix" if a[ps[n1]*n+n1] == zero
ps
end
def lusolve(a,b,ps,zero=0.0)
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]*x[j] + 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]*x[j] + dot
end
x[i] = (x[i]-dot)/a[psin+i]
end
x
end
end

Просмотреть файл

@ -0,0 +1,75 @@
#
# newton.rb
#
# Solves 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 solution vector.
# f is the object to be solved which must have following methods.
#
# f ... Object to compute Jacobian matrix of the equation systems.
# [Methods required for f]
# f.values(x) returns values of all functions at x.
# f.zero returns 0.0
# f.one returns 1.0
# f.two returns 1.0
# f.ten returns 10.0
# f.eps convergence criterion
# x ... initial values
#
require "ludcmp"
require "jacobian"
module Newton
include LUSolve
include Jacobian
def norm(fv,zero=0.0)
s = zero
n = fv.size
for i in 0...n do
s += fv[i]*fv[i]
end
s
end
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
raize "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,38 @@
#!/usr/local/bin/ruby
#
# nlsolve.rb
# An example for solving nonlinear algebraic equation system.
#
require "bigdecimal"
require "newton"
include Newton
class Function
def initialize()
@zero = BigDecimal::new("0.0")
@one = BigDecimal::new("1.0")
@two = BigDecimal::new("2.0")
@ten = BigDecimal::new("10.0")
@eps = BigDecimal::new("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

49
ext/bigdecimal/lib/pai.rb Normal file
Просмотреть файл

@ -0,0 +1,49 @@
#!/usr/local/bin/ruby
#
# pai.rb
#
require "bigdecimal"
#
# Calculates 3.1415.... using J. Machin's formula.
#
def pai(sig) # sig: Number of significant figures
exp = -sig
pi = BigDecimal::new("0")
two = BigDecimal::new("2")
m25 = BigDecimal::new("-0.04")
m57121 = BigDecimal::new("-57121")
u = BigDecimal::new("1")
k = BigDecimal::new("1")
w = BigDecimal::new("1")
t = BigDecimal::new("-80")
while (u.exponent >= exp)
t = t*m25
u,r = t.div(k,sig)
pi = pi + u
k = k+two
end
u = BigDecimal::new("1")
k = BigDecimal::new("1")
w = BigDecimal::new("1")
t = BigDecimal::new("956")
while (u.exponent >= exp )
t,r = t.div(m57121,sig)
u,r = t.div(k,sig)
pi = pi + u
k = k+two
end
pi
end
if $0 == __FILE__
if ARGV.size == 1
print "PAI("+ARGV[0]+"):\n"
p pai(ARGV[0].to_i)
else
print "TRY: ruby pai.rb 1000 \n"
end
end