Class: AppMath::R

Inherits:
Numeric
  • Object
show all
Defined in:
lib/rnum.rb

Overview

Class of real numbers which implements arbitrary precision arithmetic together with the standard elementary transcendental functions.

Implementation is based on classes BigDecimal and module BigMath created by Shigeo Kobayashi.

The term ‘precision’ is here used in one of its technical meanings: here it is the number of decimal digits used for representing the ‘significand’ of a number. Recall that each real number can be represented uniquely as

sign * significand * 10^exponent

where 

sign = +/- 
significand = 0.d_1 d_2 d_3 ...., where the d_1, d_2, are
  decimal digits in 0,1,2,3,4,5,6,7,8,9, but d_1 != 0
exponent is an integer number

For any natural number n, a real number ‘of precision n’ is one for which the digits d_i vanish for all i > n.

Notice that specifying the precision defines an infinite subset of the rational numbers. The infinity is due to the infinity of choices for the exponent. In the class BigDecimal, there is the restriction that the exponent be a Fixnum (and not a Bignum) so that there is a well-defined large but finite number of possible exponents in BigDecimal and so there is a well-defined finite set of BigDecimals for any specified value of the precision.

Unlike class BigDecimal, the present class R (recall the mathematical standard symbols N, Z, Q, R, C, H for natural, integer, rational, real, complex, quaternionic numbers) treats precision not as a parameter of the arithmetic operations but as a class variable: The computational representation of the real number world as a whole depends on the precision parameter. Class R makes it easy to write programs in a ‘precision independent style’: With a single statement, e. g.

R.prec = 100

one selects the precision of all numbers (here as 100 decimal places), and all arithmetic operations and transcendental functions take the required precision automaticaly int account. There is an important point here: Obviously the exact product of two numbers of precsion p is normally a number of precision 2p. Class BigDecimal allows to obtain the exact value of a product if the precision parameter of multiplication remains unspecified (or is specified as 0). This exact arithmetic, however, gives longer and longer numbers and slower and slower execution in virtually all algorithms. Thus cutting all arithmetic results back to precission p is essential for algorithms of real analysis to remain executable in a practical sense.

Since arbitrary precision arithmetic, not only for very large precision, is considerably slower than Float arithmetic (I found a factor 3.6 for a singular value decomposition of a 20 times 20 matrix between precision 17 computation time and Float computation time. Precision 100 led to a computation time of 5.39 s, which was 2.7 times longer than for precision 17.) it is desirable to allow switching to Float arithmetic with the same syntax:

R.prec = 0

and 

R.prec = "float"

both have this effect. In order for this to work, one has to introduce real numbers into the program by a simple specific syntax relying on R: Instead of 

x = 1.456e-6;  y = 7.7;  z = 1000.0;  n = 128

a R-based program may use the conventional mehod ‘new’ and write

x = R.new("1.456E-6");  y = R.new(7.7);  z = R.new(1000);  n = 128

Instead, we also may use the class method c (‘c’ for ‘converter’) and write (with or without argument brackets)

x = R.c 1.456e-6;  y = R.c(7.7);  z = R.c 1000;  n = 128

Written in this form, the variables x,y,z refer to R-objects for precision set as a positive integer, and to Float-objects as a consequence of the statement 

R.prec = "float"

Assume that the program continues as 

u = (x.sin + y.exp) * z

This works for R-objects since R implements as methods all the functions which module Math povides as module functions. When, however, we have the statement R.prec = “float” in action, the variables x and y refer to Float objects. Then, the terms x.sin and y.exp are not defined, unless we prevented the problem via

require 'float_ext'

which adds to the interface of Float all the methods of R so that any program hat makes sense for R, also makes sense for Float. (This is a harmless change of Float that breaks no code which relied on the unmodified class Float. If Float has been modified already for some reason, one has to make sure that this extension does not conflict with the previous one.) Instead of setting precision for the whole program, we may vary it for direct comparison by code like this:

for i in 0..8
  R.prec = 20*i
  # some computation, e. g.
  t1 = Time.now
  x = [R.c2, R.i2, R.i(0.33333),R.ran(137), R.tob(i), R.pi] 
  x.each{ |x|
    y = x.sin**2 + x.cos**2 - 1
    puts "x = #{x.round(6)}, y = #{y}"
  }
  t2 = Time.now
  puts "computation time was #{t2 - t1} for precision #{R.prec}"
end

Here all precision independent creation modes for real numbers are exemplified:

R.c2 equivalent to the yet known R.c(2), hence 2
R.i2 inverse of 2, i.e. 0.5
R.i(0.33333) inverse of 0.33333 i.e. approximately 3
R.ran(137) sine-floor random generator invoked for argument 137
R.tob(i) test object generator invoked for agument i
R.pi number pi = 3.14159... (also R.e = 2.71828...)

Since the method ‘coerce’, which determines the meaning of binary operator expressions, is defined appropriately in R, one may for a R-object r write 1/r instead of R.c1/r and r + 1 instead of r + R.c1. This allows us to avoid the appearance of R in all formulas except those where real variables are initialized, as in the part

x = [R.c2, R.i2, R.i(0.33333),R.ran(137), R.tob(3), R.pi]

of our example computation.

Class R is truly consitent with Ruby’s numerical concepts: R is derived from Ruby class Numeric so that R, Float, Fixnum, and Bignum have a common superclass. Notice that class BigDecimal cannot have this property since most of its functions need the precision parameter as input. Further, R implements as member funcions all functions which are defined in Ruby module Math and thus form the core of Ruby’s mathematics. It would therefore be natural to integrate class R into the BigMath module. It would also be natural to include float_ext.rb since this would enable the ‘precision-independent coding style’ scetched above within all number-oriented Ruby projects. One would then have a framework in which valuable algorithms such as singular value decomposition of matrices could be coded with permanent validity. I achieved to transform real matrix SVD code from ‘Numerical Recipes in C’ by Press et al. to precision-independent Ruby code and did not succeed so far with the same task for complex matrices. Inspecting existing sources for “Algorithm 358” shows how far we are away fom presenting algorithms of typical complexity in a form that is ready for usage at arbitrary precision. In order for this task to make sense, one has to have complex numbers, vectors and matrices in precision-independent Ruby style. As far as I can see it, the many mathematical tool boxes provided by the Ruby community don’t contain the tools which are needed here. The Ruby language makes it pure fun, however, to build these tools according to ones needs, and I did it with the results available from www.ulrichmutze.de. I think that it is extremely important to have an agreed framework for coding mathematical and physical algorithms. Only then one can hope that people will start to add their work and use the work of others. So far, I have collected all the algoritms that I ever had the opportunity to use in my CPM class system, which is written in C++, using a programming style that I call C+-. This is presented on my website.

Having now a tiny part of the CPM system carried over (not at all verbatim) to Ruby, I expect that a Ruby version of the whole system would be much smaller and more compact. I’m quite sure that I’ll lack enthusiasm (and even power) to create such an extended algorithm collection in Ruby myself. Would it not be great, if any algorithm (e.g. eigen-vectors/values of matrices, FFT, min and root finding) that finds its way into some Ruby library, would be written in a way that allows it to be executed in arbitrary precision, even if Float precision is the preferred option? The present class R and also the other classes to be collected in the present module AppMath, should add to the experience which is needed to come up with a working method to achieve this goal.

We know from the beginning that in BigDecimal and hence in R the exponents of numbers are restricted to the size which can be represented by a Fixnum. Although computations in physics and engineering are far from coming close to this limit, schematically applied mathematical functions easily transcend it:

R.prec = 60
a = R.c 2.2
for i in 0...9
  a = a.exp
  puts "a = " + a.to_s
end

yields 

a = 0.902501349943412092647177716688866402972021659669817926079804E1
a = 0.830832663077249493655084378868900432568369546441921929731276E4
a = 0.182141787499134800567191386180195980368655533517234407096707E3609
a = 0.0
a = 0.1E1
a = 0.271828182845904523536028747135266249775724709369995957496696E1
a = 0.151542622414792641897604302726299119055285485368561397691404E2
a = 0.381427910476022059220921959409820357102394053622666607552533E7
a = 0.233150439900719546228968991101213766633201742896351361434871E1656521

where only the first three results are correct and the following ones are determined by the overfow behavior of Fixnum.

Motivation and rational

Experimentation with class R provided many experiences which sharpened my long-standing diffuse ideas concerning a coding framework for mathematics.

In mathematics we have ‘two worlds’: the discrete one and the ‘continuous’ one. For the first world, Ruby is fit without modifications. Its automatic switch from Fixnum representation of integer numbers to a Bignum representation eliminates all limitations to the coding of discrete mathematics problems for which a coding approach is reasonable from a logical point of view.

The ‘continuous world’ is the one for which the real numbers lay the ground. The prime example of this world is ‘real analysis’ for which the concept of convergence is considered central. When a computationally oriented scientist describes to a pure mathematician his experience that all mathematically and physically relevant structures seem to have natural codable counterparts, this pure mathematician will probably admit that this holds for the trivial part of the story, but he will probably insist that all deeper questions, those concerning convergence, closedness, and completeness are a priori outside the scope of numerical methods.

According to my understanding, this mathematician’s point of view is misleading. It is, however, suggested by what mathematicians actually do. For them, it is natural, when having to work with a number the square of which is known to be 2, to ‘construct’ this number as a limit of objects (e.g. finite decimal fractions) with the intention to use this ‘exact solution of the equation x^2 = 2’ in further constructions and deductions within the framework of real analysis.

For a computationally oriented scientist an alterative view is more natural and more promising: Do everything (i.e solving x^2 = 2, and the further constructions and deductions mentioned above) with finite precision and consider this ‘whole thing’ (and not only x) as a function of this precision. When we consider the behavior for growing precision of the ‘whole thing’ we have a single limit (if we need to consider a limit at all), and the question never arises whether two limits are allowed to be interchanged. Such questions are typically not easy and a major part of the technical scills of mathematicians is devoted to them. I’m quite sure that I spent more than a year of my live struggling with such questions.

To have a realistic example, let us consider that the ‘whole thing’, mentioned above, is the task to produce all tables, figures, diagrams, and animations that will go into a presentation or publication. Then it is natural to consider a large structured Ruby program as the means to perform this task. (My experience is restricted to C++ programs instead of Ruby programs for this situation.) Let us assume that all data that normally would be represented as Float objects, now are represented as R objects. As discussed already, this means that e.g. instead of

x = 2.0

we hat to write 

x = R.c(2.0)

or 

x = R.c2

or we had to add to the statement

x = 2.0

the conversion statement

x = R.c(x)

No modification or conversion would be needed for the integer numbers! Now consider having an initial statement

R.prec = "float"

in the main program flow. By executing the program we get curves in diagrams, moving particles in animations, etc. If some of these curves are more jagged than expected, or some table values turn out to be NaN or Infinite we may try

R.prec = 40

and 

R.prec = 80
  ...

If the results stabilize somewhere, practitioners are sure that ‘the result is now independent of numerical errors’. Of course, the mathematician’s objection that behavior for a few finite values says nothing about the limit, also applies here. But we are in a very comfortable position to cope with this objection: Since we deal with the solution of a task, we are allowed to assume that we know the experience-based conventions and ‘best practices’ concerning solutions of tasks in the problem field under consideration. So, the judgement whether the behavior of the solution as a function of precision supports a specific conclusion or not has a firm basis when the context is taken into account. It is an illusionary hope that some abstract framework such as ‘Mathematics’ could replace the guidance provided by problem-field related knowledge and experience.

Let me finally describe an experience which suggested to me the concept of precision as a parameter of whole task (project) instead of a parameter of individual arithmetic operations:

In an industrial project I had to assess the influence of lens aberrations to the efficiency of coupling laser light into the core of an optical fiber. Although the situation is clearly a wave-optical one, the idea was to test whether a ray-tracing simulation would reproduce the few available measured data. In case of success one would rely on ray-tracing for unexpensive optimization and tolerance analysis. At those times, in my company all computation was done in FORTRAN and floating point number representation was by 4 bytes (just as float in C). The simulation reproduced the measured curve not really but it wiggled arround it in a remarkably symmetrical manner. Just at this time our FORTRAN compiler was upgraded to provide 8-byte numbers (corresonding to C’s double). What I had suspected turned out to be true: the ray-trace simulation now reproduced the measurements with magical precision. Congratulations to the optical lab for having got such highly consistent data! Soon I turned to programming in C and enjoyed many advantages over FORTRAN, but one paradise was lost: There was no longer a meaningful comparison between 4 byte and 8 byte precision since the available C-compiler worked with 8 bytes internally in both cases. When later we got the type long double in additon, this was an disappointment since it was identical to double for the MS-compiler and only 10 bytes (?) for GNU. So my desire was to have a C++ compiler which implemented, and consistenly used

float: 4 bytes
double: 8 bytes 
long double: 16 bytes

I then would write all my programs with a floating point type named R which gets its real meaning in a flexible manner by a statement like

typedef long double R;

I was rather sure for a long time that I would never meet a practical situation in which a simulation would show objectionable numerical errors when done with 16 byte numbers. I had to learn, however, that this was a silly idea: Let us consider a system of polyspherical elastic particles (the shape of each particle is a union of overlapping spheres) which are placed in a mirror-symmetrical container and let initial positions (placements) and velocities be arranged symmetrically (with respect to the same mirror). Then the exact motion can be seen to preserve this symmetry. However, each simulation will loose this symmetry after a few particle to particle collisions. (The particles start to rotate after collisions, which influences further collisions much more than for mono-spherical particles.) Increasing the number of bytes per number will only allow a few more symmetrical collisions, and the computation time needed to see these will increase. Of coarse, this deviation from symmetry is not objectionable from a ‘real world’ point of view. It teaches the positive fact that tolerances for making the particles, placing and boosting them, are unrealisticly tight for realizing a symmetrical motion of the system. This shows that there are computational tasks in which not all computed system properties will stabilize with increasing computational precision. With class R such computational phenomena are within the scope of Ruby programming!

Preferred command line for documentation:

rdoc -a -N -S --op doc_r r.rb

Constant Summary collapse

@@prec =

number of stored decimals

40
@@pi =

value of number pi = 3.14159…

BigDecimal("0.3141592653589793238462643383279502884197E1")
@@e =

value of Euler’s number e = 2.718281828…

BigDecimal("0.2718281828459045235360287471352662497757E1")
@@lge =

value of log10(e)

BigDecimal("0.4342944819032518276511289189166050822944E0")
@@ln10 =

value of log(10)

BigDecimal("0.2302585092994045684017991454684364207601E1")

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(aStringOrNumber = "0.0") ⇒ R

Only for @@prec > 0, the function works, otherwise it fails. The argument is assumed to be either a BigDecimal, or any object x for which x.to_s determines a BigDecimal bd via bd = BigDecimal(x.to_s). This allows to input number-related character strings and numerical literals, and numberical variables in a variety of styles:

i = 10000000
j = 123456789
a = R.new(i * j)
b = R.new(1e-23)
c = R.new("-117.07e99")
d = R.new(1e66)
e = R.new("1E666")

The most universal mehod for generating large numbers is the one from String objects as in the definition of the variables c and e above. According to the design of BigDecimal, the exponent (given by the digits after the ‘E’) has to be representable by a Fixnum.



464
465
466
467
468
469
470
471
472
473
474
475
 
# File 'lib/rnum.rb', line 464

def initialize(aStringOrNumber="0.0")
  fail "R.new makes sense only for @@prec > 0" unless @@prec > 0
  if aStringOrNumber.kind_of?(BigDecimal)
    @g = aStringOrNumber # in order to take the reuts from BigDecimal
    # arithmetic directly, speeds programs up efficiently
  else
    x = BigDecimal(aStringOrNumber.to_s) 
    fail "Unable to define a BigDecimal from the argument" unless 
      x.kind_of?(BigDecimal)
    @g = x
  end
end

Instance Attribute Details

#g ⇒ Object (readonly)

Returns the value of attribute g.



382
383
384
 
# File 'lib/rnum.rb', line 382

def g
  @g
end

Class Method Details

.aa(a) ⇒ Object

Adjust argument a so that its data type fits @@prec



772
773
774
775
776
777
778
 
# File 'lib/rnum.rb', line 772

def R.aa(a)
  if a.class == R
    a.g
  else
    BigDecimal(a.to_s)
  end
end

.add_prec(anInteger) ⇒ Object

Changes the number of digits by adding the argument to the present value. The argument may be negative too, so that one easily sets the value back by means of the same function.



442
443
444
445
 
# File 'lib/rnum.rb', line 442

def R.add_prec(anInteger)
  ai = anInteger.to_i
  @@prec += ai
end

.c(aNumber) ⇒ Object

Number generator which yiels a Float for R.prec ==0 and a R else. This is the recommended method for introducing real numbers in programs since, when followed consequently, the whole program can be switched from Float precission to arbitrary precision by a single R.prec= statement. Usage:

x = R.c("3.45E-45");  y = R.c 2;  z = R.c(1.25e13);  u = R.c "40.1e-1"



508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
 
# File 'lib/rnum.rb', line 508

def R.c(aNumber)
  if @@prec < 1 # yield a Float
    nf = aNumber.to_f
    if nf.class == Float
      nf
    else
      nff = aNumber.to_s.to_f
      if nff.class == Float
        nff
      else
        fail "Unable to define a Float from the argument"
      end
    end
  else # yield a R with the correct number of digits
    R.new(aNumber)
  end
end

.c0 ⇒ Object

The constant 0.



550
551
552
553
 
# File 'lib/rnum.rb', line 550

def R.c0
  return 0.0 if @@prec < 1
  R.new
end

.c1 ⇒ Object

The constant 1.



555
556
557
558
 
# File 'lib/rnum.rb', line 555

def R.c1
  return 1.0 if @@prec < 1
  R.new("1.0")
end

.c10 ⇒ Object

The constant 10.



600
601
602
603
 
# File 'lib/rnum.rb', line 600

def R.c10
  return 10.0 if @@prec < 1
  R.new("10.0")
end

.c2 ⇒ Object

The constant 2.



560
561
562
563
 
# File 'lib/rnum.rb', line 560

def R.c2
  return 2.0 if @@prec < 1
  R.new("2.0")
end

.c3 ⇒ Object

The constant 3.



565
566
567
568
 
# File 'lib/rnum.rb', line 565

def R.c3
  return 3.0 if @@prec < 1
  R.new("3.0")
end

.c4 ⇒ Object

The constant 4.



570
571
572
573
 
# File 'lib/rnum.rb', line 570

def R.c4
  return 4.0 if @@prec < 1
  R.new("4.0")
end

.c5 ⇒ Object

The constant 5.



575
576
577
578
 
# File 'lib/rnum.rb', line 575

def R.c5
  return 5.0 if @@prec < 1
  R.new("5.0")
end

.c6 ⇒ Object

The constant 6.



580
581
582
583
 
# File 'lib/rnum.rb', line 580

def R.c6
  return 6.0 if @@prec < 1
  R.new("6.0")
end

.c7 ⇒ Object

The constant 7.



585
586
587
588
 
# File 'lib/rnum.rb', line 585

def R.c7
  return 7.0 if @@prec < 1
  R.new("7.0")
end

.c8 ⇒ Object

The constant 8.



590
591
592
593
 
# File 'lib/rnum.rb', line 590

def R.c8
  return 8.0 if @@prec < 1
  R.new("8.0")
end

.c9 ⇒ Object

The constant 9.



595
596
597
598
 
# File 'lib/rnum.rb', line 595

def R.c9
  return 9.0 if @@prec < 1
  R.new("9.0")
end

.e ⇒ Object

The constant e = 2.718281828…



611
612
613
614
 
# File 'lib/rnum.rb', line 611

def R.e
  return Math::E if @@prec < 1
  R.new(@@e)
end

.i(aNumber) ⇒ Object

Returns the inverse of aNumber as a type controled by @@prec.



743
744
745
746
747
748
749
 
# File 'lib/rnum.rb', line 743

def R.i(aNumber)
  if @@prec < 1
    1.0/aNumber.to_f
  else
    R.c1/R.new(aNumber)
  end
end

.i10 ⇒ Object

The constant 1/10.



657
658
659
660
 
# File 'lib/rnum.rb', line 657

def R.i10
  return 0.1 if @@prec < 1
  R.new("0.1")
end

.i2 ⇒ Object

The constant 1/2 (inverse 2)



627
628
629
630
 
# File 'lib/rnum.rb', line 627

def R.i2
  return 0.5 if @@prec < 1
  R.new("0.5")
end

.i3 ⇒ Object

The constant 1/3.



632
633
634
 
# File 'lib/rnum.rb', line 632

def R.i3
  R.c1/R.c3
end

.i4 ⇒ Object

The constant 1/4.



636
637
638
639
 
# File 'lib/rnum.rb', line 636

def R.i4
  return 0.25 if @@prec < 1
  R.new("0.25")
end

.i5 ⇒ Object

The constant 1/5.



641
642
643
644
 
# File 'lib/rnum.rb', line 641

def R.i5
  return 0.2 if @@prec < 1
  R.new("0.2")
end

.i6 ⇒ Object

The constant 1/6.



646
 
# File 'lib/rnum.rb', line 646

def R.i6; R.c1/R.c6; end

.i7 ⇒ Object

The constant 1/7.



648
 
# File 'lib/rnum.rb', line 648

def R.i7; R.c1/R.c7; end

.i8 ⇒ Object

The constant 1/8.



650
651
652
653
 
# File 'lib/rnum.rb', line 650

def R.i8
  return 0.125 if @@prec < 1
  R.new("0.125")
end

.i9 ⇒ Object

The constant 1/9.



655
 
# File 'lib/rnum.rb', line 655

def R.i9; R.c1/R.c9; end

.lge ⇒ Object

The constant log10(e) = 0.43429448…



616
617
618
619
 
# File 'lib/rnum.rb', line 616

def R.lge
  return Math::log10(Math::E) if @@prec < 1
  R.new(@@lge)
end

.ln10 ⇒ Object

The constant log(10) = 2.30258…



621
622
623
624
 
# File 'lib/rnum.rb', line 621

def R.ln10
  return Math::log(10.0) if @@prec < 1
  R.new(@@ln10)
end

.nan ⇒ Object

The constant NaN, not a number.



662
663
664
665
 
# File 'lib/rnum.rb', line 662

def R.nan
  return 0.0/0.0 if @@prec < 1
  R.new(BigDecimal("NaN"))
end

.one ⇒ Object

The constant 1.



535
536
537
538
 
# File 'lib/rnum.rb', line 535

def R.one
  return 1.0 if @@prec < 1
  R.new("1.0")
end

.pi ⇒ Object

The constant pi=3.14159…



606
607
608
609
 
# File 'lib/rnum.rb', line 606

def R.pi
  return Math::PI if @@prec < 1
  R.new(@@pi)
end

.prec ⇒ Object

Returns the number of digits. Usage:

n = R.prec



435
436
437
 
# File 'lib/rnum.rb', line 435

def R.prec
  return @@prec
end

.prec=(anInteger) ⇒ Object

Setting precision, i. e. the number of digits.

If the anInteger.to_i is an integer > 0, this number will be set as the class variable @@prec. All other input will set @@prec to 0. Usage:

R.prec = 100

or

R.prec = "float"

A class state @@prec == 0 lets all generating methods create Float-objects instead of R-objects. Calling methods of R-objects is undefined in this class state. This is behavior is reasonable since all real numbers will be created as Floats for this setting of precision so that all method calls will automatically be invoked by Float-objects instead of R-objects. For this to work one needs

require 'float_ext'

and to use R.c instead of R.new.



414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
 
# File 'lib/rnum.rb', line 414

def R.prec=(anInteger)
  if anInteger.kind_of?(Numeric)
    ai=anInteger.to_i # for robustness
    ai = 0 unless ai.class == Fixnum
  else
    ai = 0
  end
  @@prec = ai > 0 ? ai : 0
  if @@prec > 0
    @@pi = BigMath::PI(@@prec)
    @@e = BigMath::E(@@prec)
    ln5 = BigMath.log(BigDecimal("5"),@@prec)
    ln2 = BigMath.log(BigDecimal("2"),@@prec)
    @@ln10 = ln5.add(ln2,@@prec)
    @@lge = BigDecimal("1").div(@@ln10,@@prec)
  end
end

.ran(anInteger) ⇒ Object

Random value (sine-floor random generator).

Chaotic function from the integers into the unit interval [0,1] of R

Argument condition: R.new(anInteger) defined

Although the main intent is to use the function for integer argments, this is not essential.



696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
 
# File 'lib/rnum.rb', line 696

def R.ran(anInteger)
  if @@prec < 1 # no further usage of R
    x = anInteger
    y = 1e6 * Math::sin(x)
    return y - y.floor
  else
    shift = 6
    fac = R.one.ldexp(shift)
    R.add_prec(shift)
    x = R.c anInteger
    y = fac * x.sin
    res = (y - y.floor).round(@@prec - shift)
    R.add_prec(-shift)
    res
  end
end

.ten ⇒ Object

The constant 10.



545
546
547
548
 
# File 'lib/rnum.rb', line 545

def R.ten
  return 10.0 if @@prec < 1
  R.new("10.0")
end

.test(n0, verbose = false) ⇒ Object

Consistency test for class R This is intended to keep the class consistent despite of modifications. The first argument influences the numbers which are selected for the test. Returned is a sum of numbers each of which should be numerical noise and so the result has to be << 1 if the test is to indicate success. For istance, on my system 

R.prec = 100; R.test(137)

produces The error sum is 0.1654782936431420775338085739363521532600906524358 495733897874347055126769599726793415224324363846749E-29 . Computation time was 0.853 seconds.



1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
 
# File 'lib/rnum.rb', line 1291

def R.test(n0, verbose = false )
  puts "Doing R.test(#{n0}, #{verbose}) for R.prec = #{@@prec}:"
  puts "*************************************************"
  require 'float_ext'
  t1 = Time.now
  small = R.prec <= 0
  s = R.c0
  puts "class of s is " + s.class.to_s
  i = n0
  a = R.tob(i) 
  i += 1
  b = R.tob(i)
  i += 1
  c = R.tob(i)
  i += 1

  r = 2 + a
  l = R.c2 + a
  ds = r.dis(l)
  puts "coerce 2 + a: ds = " + ds.to_s if verbose
  s += ds

  r =  a + 1.234
  l =  a + R.c(1.234)
  ds = r.dis(l)
  puts "coerce a + float: ds = " + ds.to_s if verbose
  s += ds

  ai = a.round
  bool_val = ai.integer?

  ds = bool_val  ? 0 : 1 
  puts "rounding gives integer type: ds = " + ds.to_s if verbose
  s += ds

  diff = (a - ai).abs
  bool_val = diff <= 0.5 
  ds = bool_val ? 0 : 1 
  puts "rounding is accurate: ds = " + ds.to_s if verbose
  s += ds

  r = (a + b) * c
  l = a * c + b * c

  ds = r.dis(l)
  puts "Distributive law for +: ds = " + ds.to_s if verbose
  s += ds

  r = (a - b) * c
  l = a * c - b * c
  ds = r.dis(l)
  puts "Distributive law for -: ds = " + ds.to_s if verbose
  s += ds
  
  r = (a * b) * c
  l = b * (c * a)
  ds = r.dis(l)
  puts "Multiplication: ds = " + ds.to_s if verbose
  s += ds
 
  a = R.tob(i) 
  i += 1
  b = R.tob(i)
  i += 1
  c = R.tob(i)
  i += 1

  r = (a * b) / c
  l = (a / c) * b
  ds = r.dis(l)
  puts "Division: ds = " + ds.to_s if verbose
  s += ds
 
  x = R.c0/R.c0
  y = x.nan?
  ds = y ? 0 : 1
  puts "0/0: ds = " + ds.to_s if verbose
  s += ds

  r = R.c1
  l = a * a.inv
  ds = r.dis(l)
  puts "inv: ds = " + ds.to_s if verbose
  s += ds
 

  r = 1/a
  l = a.inv
  ds = r.dis(l)
  puts "inv and 1/x: ds = " + ds.to_s if verbose
  s += ds

  r = b
  l = -(-b)
  ds = r.dis(l)
  puts "Unary minus is idempotent: ds = " + ds.to_s if verbose
  s += ds

  x = -a
  y = x + a
  r = y
  l = R.c0
  ds = r.dis(l)
  puts "Unary -: ds = " + ds.to_s if verbose
  s += ds

  l = a.sin * b.cos + a.cos * b.sin
  r = (a + b).sin
  ds = r.dis(l)
  puts "Addition theorem for sin: ds = " + ds.to_s if verbose
  s += ds
  
  l = a.sin ** 2 + a.cos ** 2
  r = R.c1
  ds = r.dis(l)
  puts "sin^2 + cos^2: ds = " + ds.to_s if verbose
  s += ds

  x = a.sin
  y = a.cos
  l = x.hypot(y)
  r = R.c1
  ds = r.dis(l)
  puts "hypot: ds = " + ds.to_s if verbose
  s += ds

  phi = (R.ran(i) - R.i2) * R.pi * 2
  x = phi.cos
  y = phi.sin
  r = phi
  l = x.arg(y)
  ds = r.dis(l)
  puts "arg: ds = " + ds.to_s if verbose
  s += ds
 
  l = a.exp * b.exp
  r = (a + b).exp
  ds = r.dis(l)
  puts "Addition theorem for exp: ds = " + ds.to_s if verbose
  s += ds

  l = b
  r = b.exp.log
  ds = r.dis(l)
  puts "exp and log: ds = " + ds.to_s if verbose
  s += ds

  x = c.abs
  l = x
  r = x.log.exp
  ds = r.dis(l)
  puts "log and exp: ds = " + ds.to_s if verbose
  s += ds

  i +=1
  a = R.tob(i)
  l = a.sin
  r = l.asin.sin
  ds = r.dis(l)
  puts "asin and sin: ds = " + ds.to_s if verbose
  s += ds

  i +=1
  a = R.tob(i)
  l = a.cos
  r = l.acos.cos
  ds = r.dis(l)
  puts "acos and cos: ds = " + ds.to_s if verbose
  s += ds

  i +=1
  a = R.tob(i)
  l = a.tan
  r = l.atan.tan
  ds = r.dis(l)
  puts "atan and tan: ds = " + ds.to_s if verbose
  s += ds

  i +=1
  a = R.tob(i)
  l = a.cot
  r = l.acot.cot
  ds = r.dis(l)
  puts "acot and cot: ds = " + ds.to_s if verbose
  s += ds

  i +=1
  a = R.tob(i,true) # smaller version, in order
    # not to overload function exp
  l = a.sinh
  r = l.asinh.sinh
  ds = r.dis(l)
  puts "asinh and sinh: ds = " + ds.to_s if verbose
  s += ds

  i +=1
  a = R.tob(i,true)
  l = a.cosh
  r = l.acosh.cosh
  ds = r.dis(l)
  puts "acosh and cosh: ds = " + ds.to_s if verbose
  s += ds

  i +=1
  a = R.tob(i,true)
  l = a.tanh
  r = l.atanh.tanh
  ds = r.dis(l)
  puts "atanh and tanh: ds = " + ds.to_s if verbose
  s += ds

  i +=1
  a = R.tob(i,true)
  l = a.coth
  r = l.acoth.coth
  ds = r.dis(l)
  puts "acoth and coth: ds = " + ds.to_s if verbose
  s += ds

  i += 1
  a = R.tob(i,true)
  i += 1
  b = R.tob(i,true)
  i += 1
  c = R.tob(i,true)
  i += 1

  ap = a.abs
  l = (ap ** b) ** c
  r = ap ** (b * c)
  ds = r.dis(l)
 "general power: ds = " + ds.to_s if verbose
  s += ds
  
  l = (ap ** b) * (ap ** c)
  r = ap ** (b + c)
  ds = r.dis(l)
  puts "general power, addition theorem: ds = " + ds.to_s if verbose
  s += ds

  x=(a.abs+b.abs+c.abs)
  l = x.sqrt
  r = x ** 0.5
  ds = r.dis(l)
  puts "square root as power: ds = " + ds.to_s if verbose
  s += ds

  bi = i % 11 -6
  ci = i % 7 - 3
  bi = 7 if bi.zero?
  ci = 3 if ci.zero?
  # avoid trivial 0
  l = (a ** bi) ** ci
  r = a ** (bi * ci)
  puts "bi = " + bi.to_s + " ci = " + ci.to_s if verbose
  ds = r.dis(l)
  puts "integer power: ds = " + ds.to_s if verbose
  s += ds

  r = b
  l = b.clone
  ds = r.dis(l)
  puts "cloning: ds = " + ds.to_s if verbose
  s += ds
  ds = (l == r ? 0.0 : 1.0)
  puts "cloning and ==: ds = " + ds.to_s if verbose
  x = a
  y = b
  p, q = x.divmod(y)
  l = x
  r = y * p + q
  ds = r.dis(l)
  puts "divmod 1: ds = " + ds.to_s if verbose
  s += ds

  x = a
  y = -b
  p, q = x.divmod(y)
  l = x
  r = y * p + q
  ds = r.dis(l)
  puts "divmod 2: ds = " + ds.to_s if verbose
  s += ds

  x = b
  y = a
  p, q = x.divmod(y)
  l = x
  r = y * p + q
  ds = r.dis(l)
  puts "divmod 3: ds = " + ds.to_s if verbose
  s += ds

  x = b
  y = -a
  p, q = x.divmod(y)
  l = x
  r = y * p + q
  ds = r.dis(l)
  puts "divmod 4: ds = " + ds.to_s if verbose
  s += ds

  x, y = a.frexp
  l = a
  r = x.ldexp(y)
  ds = r.dis(l)
  puts "frexp and ldexp: ds = " + ds.to_s if verbose
  s += ds

  x1 = R.c 1100000
  x2 = R.c "1100000 with comment which will be ignored"
  x3 = R.c(1200000.12)
  x4 = R.c("1200000.12")
  x5 = R.c("34567.89001953125e2")
  x6 = R.c("345.6789001953125E4")

  l = x1
  r = x2
  ds = r.dis(l)
  puts "input 1: ds = " + ds.to_s if verbose
  s += ds

  l = x3
  r = x4
  ds = r.dis(l)
  puts "input 2: ds = " + ds.to_s if verbose
  s += ds

  l = x5
  r = x6
  ds = r.dis(l)
  puts "input 3: ds = " + ds.to_s if verbose
  s += ds
  
  x = R.c9
  h = R.c(1e-6)
  fp = (x + h).erf
  fm = (x - h).erf
  l = (fp - fm)/(h * 2)
  r = (-x*x).exp*R.c2/R.pi.sqrt
  ds = r.dis(l)
  puts "erf derivative : ds = " + ds.to_s if verbose
  s += ds 
  
  t2 = Time.now
  puts "class of s is " + s.class.to_s + " ."
  puts "The error sum s is " + s.to_s + " ."
  puts "It should be close to 0."
  puts "Computation time was #{t2-t1} seconds."
  s
end

.tob(anInteger, small = false) ⇒ Object

Test object.

Needed for automatic tests of arithmetic relations. Intended to give numbers which rapidly change sign and order of magnitude when the argument grows regularly e.g. as in 1,2,3,… . However, suitibility as a random generator is not the focus. If the second argument is ‘true’, the result is multplied by a number << 1 in order to prevent the result from overloading the exponential function.



723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
 
# File 'lib/rnum.rb', line 723

def R.tob(anInteger, small = false)
  small_a = small || @@prec<=0
  ai=anInteger.to_i
  mag_num1 = 7 # 'magic numbers'
  mag_num2 = 11
  mag_num3 = 17
  mag_num4 = small_a ? 0.0423 : 1.7501
  mag_num5 = small_a ? 0.13 : 0.65432
  r1 = ai % mag_num1
  r2 = ai % mag_num2
  r3 = ai % mag_num3
  y=(-r1 - r2 + r3) * mag_num4 + mag_num5
  if @@prec < 1
    res = Math::exp(y) * (ran(ai) - 0.5)
  else
    res = R.new(y).exp * (ran(ai) - R.i2)
  end
end

.two ⇒ Object

The constant 2.



540
541
542
543
 
# File 'lib/rnum.rb', line 540

def R.two
  return 2.0 if @@prec < 1
  R.new("2.0")
end

.zero ⇒ Object

The constant 0.



530
531
532
533
 
# File 'lib/rnum.rb', line 530

def R.zero
  return 0.0 if @@prec < 1
  R.new
end

Instance Method Details

#%(aR) ⇒ Object

Usual modulo division.



826
827
828
 
# File 'lib/rnum.rb', line 826

def %(aR)
  R.new(@g%R.aa(aR))
end

#*(aR) ⇒ Object

Returns the R-object self * aR.



816
817
818
 
# File 'lib/rnum.rb', line 816

def *(aR)
  R.new(@g.mult(R.aa(aR),@@prec))
end

#**(a) ⇒ Object

Returns the a-th power of self, if a is an integer argument, and the a-th power of self.abs for real, non-integer a. We put 0**a = 0 for all non-integer a. – If a is integer, the functionality is as given by classes Float and BigDecimal. For non-integer a, we take the freedom do make the definition 0 where mathematics suggests complex (e.g. (-1)**0.5) or infinite ( e.g. 0**-1) results.



837
838
839
840
841
842
843
844
845
846
847
848
849
850
 
# File 'lib/rnum.rb', line 837

def **(a)
  return R.nan if nan?
  if a.class == Fixnum || a.class == Bignum
    return R.new(@g.power(a))
  end
  a = R.new(a) if a.class != R
  R.nan if a.nan?
  x = abs
  if x.zero?
    R.c0
  else
    (x.log * a).exp
  end
end

#+(aR) ⇒ Object

Returns the R-object self + aR.



806
807
808
 
# File 'lib/rnum.rb', line 806

def +(aR)
  R.new(@g.add(R.aa(aR),@@prec))
end

#+@ ⇒ Object

Unary plus operator. It returns the R-object self.



755
 
# File 'lib/rnum.rb', line 755

def +@; R.new(@g); end

#-(aR) ⇒ Object

Returns the R-object self - aR.



811
812
813
 
# File 'lib/rnum.rb', line 811

def -(aR)
  R.new(@g.sub(R.aa(aR),@@prec))
end

#-@ ⇒ Object

Unary minus operator. It returns the R-object -self.



752
 
# File 'lib/rnum.rb', line 752

def -@; R.new(-@g); end

#/(aR) ⇒ Object

Returns the R-object self / aR.



821
822
823
 
# File 'lib/rnum.rb', line 821

def /(aR)
  R.new(@g.div(R.aa(aR),@@prec))
end

#<=>(a) ⇒ Object

The basic order relation.



781
782
783
 
# File 'lib/rnum.rb', line 781

def <=>(a)
  @g <=> R.aa(a)
end

#abs ⇒ Object

Returns the absolute value of self.



901
 
# File 'lib/rnum.rb', line 901

def abs; R.new(@g.abs); end

#abs2 ⇒ Object

Returns the square of the absolute value of self.



904
 
# File 'lib/rnum.rb', line 904

def abs2; self * self; end

#acos ⇒ Object

Inverse cosine.



1018
1019
1020
1021
1022
 
# File 'lib/rnum.rb', line 1018

def acos
  x = self
  y = (R.c1 - x * x).sqrt
  x.arg(y)
end

#acosh ⇒ Object

Inverse hyperbolic cosine.



1191
1192
1193
 
# File 'lib/rnum.rb', line 1191

def acosh
  ((self * self - R.c1).sqrt + self).abs.log
end

#acot ⇒ Object

Inverse cotangent.



1032
1033
1034
1035
 
# File 'lib/rnum.rb', line 1032

def acot 
  a = inv
  a.atan
end

#acoth ⇒ Object

Inverse hyperbolic cotangent.



1201
1202
1203
 
# File 'lib/rnum.rb', line 1201

def acoth
  ((self + R.c1)/(self - R.c1)).abs.log * R.i2
end

#arg(y) ⇒ Object

Argument, i.e. polar angle phi of point (x=self,y), -pi < phi <= pi.

This is the basic tool for defining asin, acos, atan, acot. Notice x.arg(y) == y.atan2(x) with function atan2 to be defined next.



994
995
996
997
998
999
 
# File 'lib/rnum.rb', line 994

def arg(y)
  a=(self*self+y*y).sqrt
  res = R.c2*(y/(a+self)).atan_aux
  res -= R.pi * 2 if res > R.pi
  res
end

#asin ⇒ Object

Inverse sine.



1011
1012
1013
1014
1015
 
# File 'lib/rnum.rb', line 1011

def asin
  y = self
  x = (R.c1 - y * y).sqrt
  x.arg(y)
end

#asinh ⇒ Object

Inverse hyperbolic sine.



1186
1187
1188
 
# File 'lib/rnum.rb', line 1186

def asinh
  ((self * self + R.c1).sqrt + self).log
end

#atan ⇒ Object

Inverse tangent.



1025
1026
1027
1028
1029
 
# File 'lib/rnum.rb', line 1025

def atan 
  cosine = (R.c1 + self * self).sqrt.inv
  sine = cosine * self
  cosine.arg(sine)
end

#atan2(x) ⇒ Object

The value y.atan2(x) is the polar angle of point (x,y) and corresponds to Math.atan2(y,x), where the squeere order of arguments is the same as in the (poor) formula atan(y/x).



1004
 
# File 'lib/rnum.rb', line 1004

def atan2(x); x.arg(self); end

#atanh ⇒ Object

Inverse hyperbolic tangent.



1196
1197
1198
 
# File 'lib/rnum.rb', line 1196

def atanh
  ((R.c1 + self)/(R.c1 - self)).abs.log * R.i2
end

#ceil ⇒ Object

Returns the ceil value of self (the smallest integer R >= self).



927
 
# File 'lib/rnum.rb', line 927

def ceil; R.new(@g.ceil); end

#clone ⇒ Object

The clone has the required number of digits.



527
 
# File 'lib/rnum.rb', line 527

def clone; R.new(self); end

#complex? ⇒ Boolean

Supports the unified treatment of real and complex numbers.

Returns:

  • (Boolean) 


803
 
# File 'lib/rnum.rb', line 803

def complex?; false; end

#conj ⇒ Object

(Complex) conjugation, no effect on real numbers. Supports the unified treatment of real and complex numbers.



759
 
# File 'lib/rnum.rb', line 759

def conj; self; end

#cos ⇒ Object

Cosine.



1096
1097
1098
1099
 
# File 'lib/rnum.rb', line 1096

def cos
  piHalf=R.pi*R.new("0.5")
  (self+piHalf).sin
end

#cosh ⇒ Object

Hyperbolic cosine.



1139
 
# File 'lib/rnum.rb', line 1139

def cosh; (self.exp + (-self).exp) * R.i2; end

#cot ⇒ Object

Cotangent.



1107
1108
1109
 
# File 'lib/rnum.rb', line 1107

def cot
  self.cos/self.sin
end

#coth ⇒ Object

Hyperbolic cotangent.



1149
1150
1151
1152
1153
 
# File 'lib/rnum.rb', line 1149

def coth
  s = self.exp - (-self).exp
  c = self.exp + (-self).exp
  c/s
end

#dis(aR) ⇒ Object

Returns a kind of relative distance between self and aR. The return value varies from 0 to 1, where 1 means maximum dissimilarity of the arguments. Such a function is needed for testing the validity of arithmetic laws, which, due to numerical noise, should not be expected to be fulfilled exactly.



912
913
914
915
916
917
918
919
920
921
 
# File 'lib/rnum.rb', line 912

def dis(aR)
  aR = R.aa(aR)
  a = self.abs
  b = aR.abs
  d = (self - aR).abs
  s = a + b
  return R.c0 if s.zero?
  d1 = d/s
  d < d1 ? d : d1
end

#erf ⇒ Object

Returns the error function evaluated at x=self.



1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
 
# File 'lib/rnum.rb', line 1254

def erf
  if self < R.c0
    x = abs
    sig = - R.c1
  else
    x = self
    sig = R.c1
  end
  if x > R.c9
    return (R.c1 - x.erfc_asy)*sig
  else 
    return x.erf_ps * sig
  end
end

#erfc ⇒ Object

Returns the complementary error function erfc evaluated at x=self.



1270
1271
1272
1273
1274
1275
1276
 
# File 'lib/rnum.rb', line 1270

def erfc
  if self > R.c9
    erfc_asy
  else
    R.c1 - erf
  end
end

#exp ⇒ Object

Exponential function.



1131
1132
1133
 
# File 'lib/rnum.rb', line 1131

def exp
  (self * R.lge).exp10
end

#floor ⇒ Object

Returns the floor value of self (the largest integer R <= self).



924
 
# File 'lib/rnum.rb', line 924

def floor; R.new(@g.floor); end

#frexp ⇒ Object

Returns a two-component array containing the normalized fraction (an R) and the exponent (a Fixnum) of self. The exponent refers always to radix 10.



670
671
672
673
674
675
676
 
# File 'lib/rnum.rb', line 670

def frexp # member functions will be called only in situations in
# which @@prec > 0 is guarantied
  n = @g.exponent
  fac = BigDecimal("10").power(-n)
  prel = [@g.mult(fac, @@prec), n]
  [R.new(prel[0]),prel[1]]
end

#hypot(y) ⇒ Object

Returns he hypotenuse of a right-angled triangle with sides x = self and y.



1008
 
# File 'lib/rnum.rb', line 1008

def hypot(y); (self * self + y * y).sqrt; end

#infinite? ⇒ Boolean

Returns:

  • (Boolean) 


791
 
# File 'lib/rnum.rb', line 791

def infinite?; @g.infinite?; end

#integer? ⇒ Boolean

Since R is not Fixnum or Bignum we return ‘false’. In scientific computation there may be the need to use various types of ‘real number types’ but there should be always a clear-cut distinction between integer types and real types.

Returns:

  • (Boolean) 


797
 
# File 'lib/rnum.rb', line 797

def integer?; false; end

#inv ⇒ Object

Returns the inverse 1/self.



859
860
861
862
 
# File 'lib/rnum.rb', line 859

def inv
  return R.nan if nan?
  R.new(BigDecimal("1").div(@g,@@prec))
end

#ldexp(anInteger) ⇒ Object

Adds the argument to the exponent of self. If self is a normalized fraction (and thus has exponent 0) the resulting number has an exponent given by the argument. This is the situation to which the functions name ‘load exponent’ refers.



683
684
685
686
 
# File 'lib/rnum.rb', line 683

def ldexp(anInteger)
  ai=anInteger.to_i
  self * (R.c10 ** ai)
end

#log ⇒ Object

Natural logarithm.



1178
1179
1180
1181
1182
1183
 
# File 'lib/rnum.rb', line 1178

def log
  return R.nan if nan?
  r1=self.log10
  r2=R.e.log10
  r1/r2
end

#log10 ⇒ Object

Logarithm to base 10.

Makes use of knowing the decimal exponent for a reduction to an actually small argument. We understand log10(x) as log10(|x|), which is natural in the real domain.



1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
 
# File 'lib/rnum.rb', line 1162

def log10 
  return R.nan if nan?
  return R.nan if @g.infinite? || @g.nan?
  sp=@g.split
  exponent=sp[3]
  s="0."+sp[1]
  x=BigDecimal(s)
  if x.zero?
    fail "zero argument of log"
  end
  x=BigMath.log(x.abs,@@prec)
  y=BigDecimal(exponent.to_s)
  R.new(x) * R.new(@@lge) + R.new(y)
end

#nan? ⇒ Boolean

Returns ‘true’ if self is ‘not a number’ (NaN).

Returns:

  • (Boolean) 


789
 
# File 'lib/rnum.rb', line 789

def nan?; @g.nan?; end

#prn(name) ⇒ Object

Printing the value together with a label



933
934
935
 
# File 'lib/rnum.rb', line 933

def prn(name)
  puts "#{name} = " + to_s
end

#pseudo_inv ⇒ Object

The pseudo inverse is always defined: the pseudo inverse of 0 is 0.



865
866
867
868
 
# File 'lib/rnum.rb', line 865

def pseudo_inv
  return R.c0 if zero?
  inv
end

#real? ⇒ Boolean

Supports the unified treatment of real and complex numbers.

Returns:

  • (Boolean) 


800
 
# File 'lib/rnum.rb', line 800

def real?; true; end

#round(*arg) ⇒ Object

If the method gets no argument we return the ‘nearest integer’: For the return value res we have 

res.int? == true

and 

(self - res).abs <= 0.5

For an integer argument we return a real number, the significand of which has not more than n digits. Notice that there is also a function.



878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
 
# File 'lib/rnum.rb', line 878

def round(*arg)
  n = arg.size
  case n
  when 0
    (self + 0.5).floor.to_i # here we ask for an integer output
    # notice that R#round maps to R
  when 1
    m = arg[0].to_i
    x = frexp
    y = x[0].ldexp(m)
    (y + 0.5).floor.ldexp(x[1] - m)
  else
    fail "needs 0 or 1 arguments"
  end
end

#sin ⇒ Object

Sine.

This reduces computation of the sine of any angle to computation of sine or cosine of a angle less than pi/4 and thus less than 1. This speeds up the convergence of the sine or cosine power series.



1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
 
# File 'lib/rnum.rb', line 1044

def sin 
  return R.nan if @g.infinite? || @g.nan?
  sign=1
  if @g >= BigDecimal("0") 
    x=@g
  else
    x=-@g
    sign=-1
  end
  twoPi=@@pi*BigDecimal("2")
  piHalf=@@pi*BigDecimal("0.5")
  piQuart=@@pi*BigDecimal("0.25")
  res=x.divmod(twoPi)
  phi=res[1]
  res2=phi.divmod(piHalf)
  qn=res2[0] # number of quadrant 0,1,2,3
  alpha=res2[1] # angle in quadrant
  if qn==0
    # first quadrant
    if alpha<=piQuart
      r=R.new(BigMath.sin(alpha,@@prec))
    else
      r=R.new(BigMath.cos(piHalf-alpha,@@prec)) 
    end
  elsif qn==1
    # second quadrant
    if alpha<=piQuart
      r=R.new(BigMath.cos(alpha,@@prec))
    else
      r=R.new(BigMath.sin(piHalf-alpha,@@prec)) 
    end
  elsif qn==2
    # third quadrant
    if alpha<=piQuart
      r=-R.new(BigMath.sin(alpha,@@prec))
    else
      r=-R.new(BigMath.cos(piHalf-alpha,@@prec)) 
    end
  elsif qn==3
    # fourth quadrant
    if alpha<=piQuart
      r=-R.new(BigMath.cos(alpha,@@prec))
    else
      r=-R.new(BigMath.sin(piHalf-alpha,@@prec)) 
    end
  else
    puts "error in R.sin, not assumed to happen"
  end
  sign == -1 ? -r : r
end

#sinh ⇒ Object

Hyperbolic sine.



1136
 
# File 'lib/rnum.rb', line 1136

def sinh; (self.exp - (-self).exp) * R.i2; end

#sqrt ⇒ Object

Returns the square root of self.



895
896
897
898
 
# File 'lib/rnum.rb', line 895

def sqrt
  return R.nan if nan?
  R.new(@g.sqrt(@@prec))
end

#tan ⇒ Object

Tangent.



1102
1103
1104
 
# File 'lib/rnum.rb', line 1102

def tan
  self.sin/self.cos
end

#tanh ⇒ Object

Hyperbolic tangent.



1142
1143
1144
1145
1146
 
# File 'lib/rnum.rb', line 1142

def tanh
  s = self.exp - (-self).exp
  c = self.exp + (-self).exp
  s/c
end

#to_0 ⇒ Object

Returns the zero-element which belongs to the same class than self



853
 
# File 'lib/rnum.rb', line 853

def to_0; R.c0; end

#to_1 ⇒ Object

Returns the unit-element which belongs to the same class than self



856
 
# File 'lib/rnum.rb', line 856

def to_1; R.c1; end

#to_f ⇒ Object

Conversion to double.



944
945
946
 
# File 'lib/rnum.rb', line 944

def to_f
  @g.to_f
end

#to_i ⇒ Object

Conversion to integer.



938
 
# File 'lib/rnum.rb', line 938

def to_i; @g.to_i; end

#to_int ⇒ Object

Conversion to integer.



941
 
# File 'lib/rnum.rb', line 941

def to_int; @g.to_int; end

#to_s ⇒ Object

Conversion to String.



930
 
# File 'lib/rnum.rb', line 930

def to_s; @g.to_s; end

#zero? ⇒ Boolean

Returns ‘true’ if self equals zero.

Returns:

  • (Boolean) 


786
 
# File 'lib/rnum.rb', line 786

def zero?; @g.zero?; end