Class: AppMath::R

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.

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)

  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

Methods

% * ** + +@ - -@ / <=> aa abs abs2 acos acosh acot acoth add_prec arg asin asin_aux asinh atan atan2 atan_aux atanh c c0 c1 c10 c2 c3 c4 c5 c6 c7 c8 c9 ceil clone coerce complex? conj cos cosh cot coth dis e erf erf_ps erfc erfc_asy exp exp10 floor frexp hypot i i10 i2 i3 i4 i5 i6 i7 i8 i9 infinite? integer? inv ldexp lge ln10 log log10 nan nan? new one pi prec prec= prn pseudo_inv ran real? round sin sinh sqrt tan tanh ten test to_0 to_1 to_f to_i to_int to_s tob two zero zero?

Attributes

g [R]

Public Class methods

aa(a)

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

Methods

Attributes

Public Class methods

Public Instance methods

Protected Instance methods