Chapter 22: First encounters of an Aldor novice
22.2 First vignette: a modified factorial function
22.3 Second vignette: the Lambert W function
22.4 Third vignette: binary doubling
22.5 Fourth vignette: compact finite differences
22.6 Fifth vignette: rounded rationals
22.1 : Overview
This chapter is a report of my first experiences with Aldor. The programs that I wrote are disjointed and (mostly) unrelated. They will be presented as a series of vignettes which can be read independently, though they will be given in historical order.
I first learned to program in ALGOL-W, which was rather like Pascal (only better). Subsequently, I have used lots of ``blockish'' languages like Fortran. My most usual languages for programming at the moment are Maple, MATLAB, Derive, and occasionally RPL, the language for the Hewlett-Packard 48 calculators. I wrote my last LISP program more than 15 years ago. I have never programmed in C, ML, Scheme, SmallTalk, or C++, and before learning Aldor I had never really used AXIOM. So my background is really ``orthogonal'' to the roots of Aldor. Nonetheless, I found Aldor easy to ``get into''. Some comments refer to the environment in which Aldor is used, not Aldor per se. However, as a typical user will come across these as well, they remain part of these Aldor experiences.
[Editor's note. Readers may find problems described here
which are no longer encountered in the current version of Aldor.
However, the problem-solving approach should still be found useful
and, since the status of a ``novice'' cannot be re-created, no attempt
has been made to recast the author's descriptions to match the
behaviour of the current system.]
22.2 : First vignette: a modified factorial function
The first thing I did when I wanted to learn how to write an Aldor program was to type in and modify a model program. I chose the factorial program (from section23.2) first.
The final version of my modified file, complete with
comments, is
in figure 22.1.
22.2.1 : Lessons learned in writing this program
Figure 22.1: Modified factorial program
#include "aldor.as"
fact(n:Integer):Integer == {
if n < 0 then error "Boom!"
else if n < 2 then 1
else n*fact(n-1) ;
}
import from Integer;
print << fact 7 << newline;
print << fact 78 << newline;
print << fact(-3) << newline;
-- The following call never occurs, because the error above stops
-- the file from being read further.
print << fact 0 << newline;
fact(-3)
.
Note that parentheses are not always necessary for function
application, in the absence of ambiguity: fact 7
is the same
as fact(7)
.
"fact.as", line 25: print << fact 7 << newline; ..............^ [L25 C15] #1 (Error) No meaning for integer-style literal `7'. "fact.as", line 26: print << fact 78 << newline; ..............^ [L26 C15] #2 (Error) No meaning for integer-style literal `78'. "fact.as", line 27: print << fact(-3) << newline; ..............^^ [L27 C15] #3 (Error) There are no suitable meanings for the operator `-'. [L27 C16] #4 (Error) No meaning for integer-style literal `3'. "fact.as", line 30: print << fact 0 << newline; ..............^ [L30 C15] #5 (Error) No meaning for identifier `0'. |
This is a ``sanitized'' version of my first program, in the sense that it works now. Of course, the program went through several iterations. The output of the program is now
5040 1132428117820629783145752115873204622873174957948825199004 8962825668835325234200766245086213177344000000000000000000 Boom!(except that I manually split the big integer answer for factorial 78 over two lines). The answers are correct, as can be verified by comparing them with Stirling's formula and counting the zeros (there are fifteen 5's in 78 and three 25's, making eighteen factors of 5 in factorial 78, which gives eighteen trailing zeros).
It would, perhaps, be more useful to have presented seven variations
of this program, removing an error each time (which is more or less
how it was created) but I think that just reading the ``lessons learned''
in the comments above gives the flavour. It is useful to see two of
the error messages from the compiler, though.
Suppose we comment
out the line import from Integer. Then the error messages
we get are shown in figure 22.2;.
tmoderr1Leaving out import from Integerexamples/moderr1.out
I found them very puzzling to begin with.
The compiler does not use its knowledge about
Let's put the import from Integer statement back in and investigate
another error message.
Suppose that we
leave the semicolon off the first print statement.
Then compilation
produces the error messages in figure 22.3;
In fact, the error messages have improved since I first wrote that
program: now, it is clear that
a missing semicolon is causing the problem. What about
the first error though? What does There are no suitable meanings
for the operator `newline' mean? The compiler interpreted
lines 25 and 26 as a one-line statement, juxtaposing
newline and print. If there had been an operator named
newline which took as arguments something of the same type as
print, then the compiler would have thought we meant
newline(print) --
that is, apply the operator newline to the object print. That is
what the compiler is complaining about -- it couldn't find any suitable
operators called newline (and a good thing, too, because that isn't
what I meant). 22.3 : Second vignette: the Lambert W functionInteger
s
unless you tell it that it can --
but this is a feature: after all, we might have wanted to use
the single precision kind, SingleInteger,
instead of the infinite precision kind.
Figure 22.3: Leaving off a semicolon
"fact.as", line 25: print << fact 7 << newline
......................^
[L25 C23] #2 (Error) There are no suitable meanings for
the operator `newline'.
"fact.as", line 26: print << fact 78 << newline;
^
[L26 C1] #1 (Warning) Suspicious juxtaposition. Check for missing `;'.
Check indentation if you are using `#pile'.
22.3.3 Lessons learned in writing this program
Let's now look at a simple program to evaluate a mathematical
function known as the Lambert W function;(see the paper ``The Lambert W
Function'', by Corless, Gonnet, Hare, Jeffrey and Knuth, or the article
``The Lambert W Function in Maple'', in the Maple Technical Newsletter no. 9).
The mathematics for this programming example is encapsulated below.
You may skip directly to section if you wish: you
need only know that this procedure computes the values of a
mathematical function.
22.3.1 : Mathematics of W
The program uses a third-order accurate variant of Newton's method known as
Halley's method to solve the equation defining the Lambert W function,
namely
Originally, the program just used the initial guess w0 = 0, but this
is a very poor starng pointt and we can do better. The program now
uses three different initial guesses, depending on the value of x.
If -1/e <= x < -.26, we use the initial guess
where p2 = 2(ex + 1). This uses the first four terms of the series
expansion of W(x) about the branch point x=-1/e. It is exact if x = -1/e, and out by about 0.2 at the origin.
which uses the first four terms of the Taylor series expansion for W(x) about x=0. The magic number -0.26 was chosen because the two initial guesses agree near there.
The third initial guess is from an asymptotic series. If x >= 2 we use
where ,
, and
.
22.3.2 : The program
The program appears below. The line numbers are added for this book -- they have no meaning for Aldor (and indeed if you type them in you will get a syntax error).
#include "aldor.as" import from Float; import from Integer; +++ The Lambert W function by Halley's iteration +++ with initial guess w given by a polyalgorithm. LambertW(x:Float):Float == { local w:Float := 0.0; local residual:Float := 1.0; local ew,p,sigma,tau,Ltau:Float; import from Integer; { -- A nifty way to do cases. Thx BT. x < -0.26 => { p := sqrt(2.0*(exp(1.0)*x+1.0)); w := -1.0 + p - p^2/3.0 + (11.0)/(72.0)*p^3; } x < 2.0 => w := x - x^2 + 1.5*x^3 - 8.0/3.0*x^4; sigma := 1.0/log(x); tau := log(log(x))*sigma; Ltau := -log(1.0-tau); w := (1.0 - tau)/sigma + Ltau - Ltau*sigma/(1.0 - tau + sigma) -- 4.27 } ew := exp(w); residual := w*ew - x; print << "x = " << x << newline; print << "w = " << w << newline; print << "residual = " << residual << newline; while abs residual > 10.0^(1-digits()) repeat { w := w - residual/( (w+1.0)*ew - 0.5*(w+2.0)*residual/(w+1.0) ); ew := exp(w); residual := w*ew - x; print << "w = " << w << ", residual = " << residual << newline; }; return w } digits 30; print << LambertW(-0.2) << newline; print << LambertW(7.0) << newline; -- This following one causes an error message and stops execution, -- because LambertW is not defined for x < -1/e. print << LambertW(-1.0) << newline;
22.3.3 : Lessons learned in writing this program
.0
after floats to distinguish them from integers.
local
(of the same type) in
one statement.
The output is as below:
x = -0.2 w = -0.2562666666 6666666666 6666666667 residual = 0.0016661142 1954631044 636425262 w = -0.2591710831 8621608732 8396776583 residual = 0.1065223677 6279665202 97 E -7 w = -0.2591711018 1907374505 664700903 residual = 0.2824831 E -23 w = -0.2591711018 1907374505 6651950215 residual = 0.1 E -30 -0.2591711018 1907374505 6651950215 x = 7.0 w = 1.5152707108 8819613933 992903568 residual = -0.1045289885 5515935357 16815778 w = 1.5243450637 7087054754 909586332 residual = -0.0000016369 6326817813 3919098 w = 1.5243452049 8414436912 194500079 residual = -0.615611821 E -20 w = 1.5243452049 8414436912 247606068 residual = 0.3 E -29 1.5243452049 8414436912 247606068 negative sqrt
We can see that Aldor breaks floating-point numbers up into groups of ten when it prints, and that it doesn't print trailing zeros. We also see that we get to choose how many digits to use for floating-point arithmetic.
Remarks and Questions. It occurs to me now, as I write this, that I may sometimes want to evaluate this function over SingleFloat or DoubleFloat (for efficiency reasons, because that arithmetic can be performed in hardware), or in interval arithmetic (for guaranteed accuracy), or over complex forms of any of those. So I should think about writing a ``generic'' version of this program which is ``parameterized'' by the domain I wish to work in.
If this were AXIOM and I were a pure mathematician, I might try to make
this program work for an arbitrary field -- but that is too general an
object, because I also need the function ex as well as the field
operations. And, thinking a bit more, we see that strictly speaking
a field is too restrictive anyway: for Float
s, associativity of +
does not hold (try adding a very large number plus its negative plus 1
with different groupings; in one case you will get zero and in the
other you will get 1). So we will need to define a new class of object,
an ``arithmetic system'' if you will, that includes any object which has
the field operations (but may not satisfy the field axioms). A useful
subclass of this would consist of sets which allowed elementary
functions, such as ex, sin x,
to be defined on them.
Clearly, this is something to be mulled over a bit more. On with the
Aldor vignettes.
As a parting shot, I note that the program for the Lambert W function can
be compiled successfully against the AXIOM library and used from within AXIOM,
with only a trivial change.
Since AXIOM used **
instead of ^
for raising to a power,
line 31 of the program had to be changed. In addition, AXIOM
refused to compute 1 - digits()
because the answer might
be a negative integer, whereas digits()
must return a positive
integer. To avoid this over-niceness, we rewrite line 31 with
10.0/10.0**digits()
instead. After that change, the program works inside of AXIOM.
Some other programs from this chapter are not yet able to be compiled
against the AXIOM library.
22.4 : Third vignette: binary doubling
22.4.1 Binary doubling: theory
22.4.2 The program and results
22.4.3 Lessons learned in writing this program
This next example looks at binary doubling as a method of
embedding integers into any arithmetic system. The original motivation
for this program was that in an earlier version of the program compact.as
(which is covered in the next vignette) I needed to be able to embed integers
into an arbitrary field. The program turned into a nice Aldor exercise all
on its own, and I probably learned more from it than I did from the compact
finite difference program that motivated it.
While I was writing this program, the developers went ahead and included
the functionality in the library, so this program is not, in fact, necessary.
Further, there is a package for binary powering (see section 26.2;)
which implements these same ideas still more efficiently. Nonetheless, this
program taught me a lot and hence is included here.
22.4.1 : Binary doubling: theory
Skip this paragraph if you know how binary doubling
or binary powering works, and go straight to section 22.4.2.
Binary doubling
is a reasonably efficient way of computing n x when all you can
do is add things like x together. We are going to use it to compute n x
1 where 1 is the unit of the field under consideration.
The basic idea is that we put q = c + vs
where c =0, s = 1, and v = n initially, and test the bits of
n. (Binary powering uses q = c sv, with
appropriate changes in the initial values.)
If the lowest-order bit of n is 1, then v
is odd, and we rearrange q = c + vs as (c + s) + (v - 1)s. Thus, we can
replace c by c+s and v by v-1. If instead
the lowest-order bit of n is zero, then v is even, and we can
rearrange q = c + vs as c + (v/2)(s + s). Thus we can replace s by s+s and
v by v/2.
Of course, we needn't actually carry out the division of v
by 2 because that's just shifting, and we can implicitly do that
just by looking at the next order bit of n. Note that our
answer is always equal to q, and that each operation decreases the
integer v.
When v = 0, we know q = c, and we are finished.
22.4.2 : The program and results
The code, which is found below, contains an example of a ``generic''
function. The function accepts as input three things: an integer of one
type or another (e.g. Integer or SingleInteger), its
type, and the type of the output that is desired, which must be one
ArithmeticSystem or another (for example, SingleFloat
is
an ArithmeticSystem
because it has +, -, and
*, though these operations need not (and do not) satisfy the
axioms defining a ring, which imply, for example, that
a + (b + c) = ((a + b) +
c).
1 #include "aldor.as" 2 3 SI ==> SingleInteger; import from SI; FI ==> DoubleFloat; import from FI; 4 5 -- "embed then into this " 6 embed(I:IntegerNumberSystem, n:I, F:ArithmeticSystem):F == { 7 local m:I; local s,c:F; local len:SI; 8 m := abs n; len := length(m); 9 s := 1; c := 0; 10 11 -- case statement: simple cases handled efficiently. 12 13 { m=0 => c := 0; 14 m=1 => c := s; 15 16 -- otherwise, (loop invariant `answer' = c + v*s, v = m initially) 17 18 for i:SI in 1..(len-1) repeat { 19 if bit(m,i-1) then c := c + s; -- v <== v - 1 20 s := s + s; -- v <== v/2 21 } 22 c := c + s; -- bit(m,length(m)-1) = 1 always. Save 1 add by putting 23 } -- this statement out here. 24 25 -- return the value of c with the appropriate sign. 26 if n < 0 then -c else c 27 } 28 -- Now a simple set of tests to try the above program out. 29 30 print << embed(SI, 1, FI) << newline; 31 print << embed(SI, 0, FI) << newline; 32 print << embed(SI, 2, FI) << newline; 33 print << embed(SI, -2, FI) << newline; 34 print << embed(SI, 3, FI) << newline; 35 print << embed(SI, 16, FI) << newline; 36 print << embed(SI, 5, FI) << newline; 37 print << embed(SI, 420, FI) << newline; 38 pow:SI := 30; big:SI := 2^pow; -- 2^30 = 1,073,741,824 39 print << embed(SI, big, FI) << newline; 40 41 import from Integer; 42 bignum:Integer := 2^64; -- 2^64 = 18446 74407 37095 51616 43 print << embed(Integer, bignum, FI) << newline;
The fact that we can ``parameterize'' algorithms by the type of input and output is very powerful, and a major reason to learn Aldor. It allows us to write code once for a very wide variety of applications.
The output from this program was:
1 0 2 -2 3 16 5 420 1073741824 1.8446744073709552e+19and this is correct. Note that trailing zeros are trimmed from floating point output.
22.4.3 : Lessons learned in writing this program
foo(n:I, I:IntegerNumberSystem, ...)
doesn't work, because, when the signature is being compiled, n:I can't be figured out. Instead, you have to put I:IntegerNumberSystem first. Perhaps this could be changed in a later release.
local
s to be declared explicitly,
and in fact the program can be shorter and more readable if they
aren't. However, if you don't declare a local explicitly, say in a
procedure which is contained in a file, say file1
, and if
this file is #include
d into another file, say file2
which declares a global variable with the same name, then the
compiler, when it sees file2
, will produce
a warning message telling you about the duplication.
Explicit declaration of locals will prevent this warning
message from being issued.
Addendum to the lesson:
Using #include
is possible, but rarely done if the included
files were written by someone else. Separate compilation is by far
to be preferred, and if file1
and file2
are compiled
separately, then it is safe to use implicit declarations.
Integer
is so long that its length doesn't
fit in a SingleInteger
(and I hope never to find out).
The number would be slightly more than a billion digits long, and
would require about 1.5 billion field + operations to embed in the
field.
Addendum: An ``engineering decision'' was
made here: if you have an integer so long that its length doesn't fit
into a SingleInteger
, then that integer will take up an eighth
of your address space; it was decided that length
will just
return the wrong answer in that case.
SingleInteger
example, since 231
is too big to fit in a SingleInteger
and so looks like a
negative number. ``That's not a bug, that's a feature''.
Caveat emptor.
SingleInteger
to the SingleInteger
or are
Integer
s involved?)
import
ArithmeticSystem
,
F
. So this program was a fun exercise but, ultimately, a
reinvention.
BinaryPowering
package in section 26.2.
length(0)
is 0 and of length(1)
is 1.
The program bit(m,n)
returns the n bit of m, not
the other way around. You can find this out from the manual,
but it is easier to write an Aldor program which prints out the
length
of various numbers and the various bits of a number,
say 23.
22.5 : Fourth vignette: compact finite differences
22.5.1 Mathematics of compact finite differences
22.5.2 Lessons learned in writing this program
22.5.4 Lessons learned in writing this program
Now we look at some code which uses ``compact finite differences'' to find fourth-order accurate approximations to the derivative of a function given on a grid of equally-spaced values. The technique of compact finite differences is more general than this, but for this application has an interesting advantage: we can choose the formulas so that the tridiagonal matrix that occurs in the evaluation of the finite differences can be factored exactly. This leads to a very efficient algorithm.
The program is as below.
My original version of this program used the ``embed'' program
discussed in the previous vignette. However, that was superseded
by the library coercion of integers into arbitrary fields: as a
result, the code is now much shorter and easier to read. There
remains a section which explicitly tells the compiler how to multiply
elements of the field by integers, and similarly for other operations.
Perhaps that, too, will soon be unnecessary.
22.5.1 : Mathematics of compact finite differences
Skip this section if you like.
Compact finite differences are a relatively new
method of efficiently approximating derivatives. The basic trick is that
instead of the usual explicit formulas for finite differences, for
example:
which we use at every interior grid point, we write down a (banded)
linear system of equations for the unknown derivatives, for
example:
for k = 2, 3, ..., n - 1.
Together with equations for differences at the boundaries, this gives us a tridiagonal linear system to solve, to find the derivatives at all the grid points. This approach has several advantages:
I implemented an Aldor program to find the fourth-order accurate approximation to the first derivatives on a uniform grid using the scheme
for k = 2, 3, ..., nn -1, and here I chose slightly different formulas for the endpoints so the linear system is
a = 1/c and D = diag(c,c,...,c). (It turns out
that c = 2 - would also work , here.) Then solving A
{f '} = {
b} can be reduced to the two trivial problems L{y} = {b} followed by U{x} = {y}. The solution of the system can therefore be expressed by the following recurrence relations:
22.5.2 : Lessons learned in writing this program
103::F * x
is no fun when you want 103*x
to
be automatic. With the current compiler, it's cleaner
to add eight one-line programs which define integer*field
etc. These are the lines marked by -*-*-*-*
in the program
on.
(These may be subsumed by library code at some future release).
PrimitiveArray
s'' don't know how big they are but ``Array
s'' do.
The compiler automatically figures out that 3*c
is an element
of an ArithmeticSystem and evaluates this by a call to the
particular multiplication routine which takes ``SingleInteger
s''
as its first argument and ArithmeticSystem
elements
as its second argument, however...
It turns out that there are two possible
routines that could be used to subtract 1
from an
ArithmeticSystem
element: namely, one that takes an integer
as its second argument and one that takes an ArithmeticSystem
element there - this is because 1 is already an
ArithmeticSystem
element (both 1 and 0 are members
of any ArithmeticSystem
).
Now, the routine that takes (x:F) - (y:F)
to F
will give
the same results as the one that takes (x:F) - (n:SI)
to F
,
because the second actually does a coercion and calls the first -- but the
Aldor compiler isn't (yet) smart enough to figure out that hierarchy.
So we have to help Aldor out by telling it which routine to use, and
the simplest solution is to force the choice by using 1@F.
22.5.3 : Testing
We need to test this program. The program below gives a simple test, to ensure that no blunders were committed.
1 -- compile with aldor -Fx -Q3 -M no-mactext comtest.as compact.ao 2 3 #include "aldor.as" 4 #library compactlib "compact.ao" 5 import from compactlib ; 6 7 SI==> SingleInteger; F ==> DoubleFloat; AF==> Array F; 8 9 import from SI, F, AF, InFile, StandardIO, NumberScanPackage Integer, String; 10 11 print << "Enter n : "; 12 line: String := readline! stdin; 13 n:SI := (retract$Integer) (scanNumber line); 14 print << "n is " << n << newline; -- Good people echo their input. 15 16 h:F := 1@F/(n::F); 17 18 x:AF := new(n, 0); -- x values 19 f:AF := new(n, 0); -- function values 20 xdf:AF := new(n, 0); -- exact derivatives 21 adf:AF := new(n, 0); -- approximate derivatives 22 23 fun(t:F):F == { 1@F/(1@F + t*t) }; 24 dfun(t:F):F == { -2::F*t/(1@F + t*t)/(1@F + t*t) }; 25 26 for i:SI in 1..n repeat { 27 x(i) := (i::F)*h ; 28 f(i) := fun(x(i)); 29 xdf(i) := dfun(x(i)); 30 } 31 32 rt3:F := sqrt(3::F); 33 adf := derivatives(F, f, h, rt3); 34 35 errors : AF := new(n,0); 36 for i in 1..n repeat errors(i) := abs (xdf(i) - adf(i))/(h*h*h*h) ; 37 38 -- The following bit of code is modelled on the list routine "reduce". 39 maxerror(e:AF):F == { 40 ans:F := 0; 41 for er in e repeat ans := max(ans,er); 42 ans 43 } 44 print << "The errors divided by h^4 are less than or equal to " << maxerror(errors) << newline;The results are
[1.296] % comtest Enter n : 10 n is 10 The errors divided by h^4 are less than or equal to 18.440306446613757 [1.297] % comtest Enter n : 100 n is 100 The errors divided by h^4 are less than or equal to 3.7878196627871832 [1.298] % comtest Enter n : 1000 n is 1000 The errors divided by h^4 are less than or equal to 2.9709568138969185
We can see that the errors are O(h4) as promised.
22.5.4 : Lessons learned in writing this program
^2
is not defined.
You can rewrite x^2
as x*x
to avoid importing
from Integer. This saves confusion if
``SingleInteger
s are present.
.as
files just like you can
the library. You can nest #include
#includes as well.
DoubleFloat
to have a square root function in it.
#include
statements.
22.6 : Fifth vignette: rounded rationals
This last vignette demonstrates the most powerful features of Aldor I have used to date. Paradoxically, however, it was one of the easiest programs to write. It appears to be generally true that writing ``parameterized types'' is even easier than writing ordinary functions.
What are ``rounded rationals''? They are ordinary rational numbers, restricted so that their denominators do not grow ``too big''. They are like floating- point numbers in many ways: for example, they are not exact, and they can be implemented efficiently. Suppose we are working with rationals rounded so that the denominators cannot exceed 9 (that is, only one-digit denominators are allowed). Then if we add 1/2 + 1/3, we can get 5/6 exactly right in this system, but if we add 1/4 + 1/7 we get 11/28, which cannot be represented in this system, and we round to the nearest, namely 2/5. How do we know 2/5 is the closest one-digit denominator rational to 11/28? The answer lies in the classical theory of simple continued fractions.
We expect that all the usual rational arithmetic operations will be available
for rounded rationals, with the understanding that results of the operations
+, -, *, and / will be rounded to nearest.
22.6.1 : Simple continued fractions: theory
22.6.1 Simple continued fractions: theory
22.6.2 Lessons learned in writing the program
22.6.4 Lessons learned in writing the program
This section
can be skipped, once you realize that it defines a ``magic algorithm'' that
finds the nearest rational number with denominator within the limits.
A simple continued fraction (see C. D. Olds, Continued Fractions,
Dover) is a fraction of the form
The next paragraph explains how we can compute a continued fraction for a rational number by (essentially) the Euclidean algorithm for finding the greatest common divisor of two integers.
Recall that the Euclidean algorithm runs as follows: given two integers
a and b as input, we put e0 = a > e1
= b (interchanging a and b if necessary so that a > b),
and the numbers ak and ek in the followig sequence:
that is to say, defining ek+1 as the (integer) remainder on dividing ek-1 by ek, and defining the integers ak as the quotients that arise. We terminate the computation when some ek = 0, which must eventually happen since 0 <= ek+1 by the properties of the division algorithm.
Now suppose that we wish to compute a continued fraction expansion for a/b.
Then
The main trick is now to rewrite e2/e1 as the
reciprocal of its reciprocal:
and at last
This process is continued in the obvious way, which is to replace e3/e2 recursively by its continued fraction.
Therefore, given a rational number x we can find its continued fraction
by the Euclidean
algorithm, so the cost of finding the continued fraction is no more than the cost
of computing a GCD (in the naive way). What happens if instead we are given the
continued fraction first? It turns out that, given the partial quotients
ak of
the continued fraction, we can compute a sequence of approximants xk =
pk/qk by the following recurrences, with
p0 = a0, q0 = 1, p-1 = 1 and
q-1 = 0.
These approximants have a ``best approximation'' property: pk/qk is the best approximation to x with denominator smaller than qk2.
The algorithm we will use should be obvious now. Given two rounded rationals and an operation to perform on them, we first perform the operation exactly and then round the result to the nearest acceptable rational by computing the partial quotients of our result, until the denominator of xk exceeds the maximum allowed. The code below does exactly that. I did not write this code from scratch: I looked at the source code for Ratio, which I found in the file
$ALDORROOT/samples/libaldor/ratio.as
1 -- 2 -- First new type: Rounded Rationals. RMC 25/5/1994 3 -- 4 -- Basic idea: use continued fractions to keep rationals short. 5 -- 6 7 -- Lessons. 8 -- 6) Compile with aldor -Fo RoundedRatio.as for later use. 9 10 11 #include "aldor.as" 12 13 14 -- The parameter MAXDENOM ensures that denominators 15 -- are not too big. Currently numerators are allowed to 16 -- be MUCH larger---probably this is a mistake. 17 18 RoundedRatio( I: IntegerNumberSystem, 19 MAXDENOM: I): 20 Join(Field, OrderedRing) with { 21 -- the following extra operations: 22 23 +: (I, %) -> %; -- These eight operations allow easy use of 24 -: (I, %) -> %; -- embedded integers. 25 *: (I, %) -> %; 26 /: (I, %) -> %; 27 +: (%, I) -> %; 28 -: (%, I) -> %; 29 *: (%, I) -> %; 30 /: (%, I) -> %; 31 32 /: (%, %) -> %; 33 \: (%, %) -> %; 34 /: (I, I) -> %; 35 inv: % -> %; 36 numer: % -> I; 37 denom: % -> I; 38 -- coerce: I -> %; 39 quo: (%, %) -> %; -- needed if we want to 40 rem: (%, %) -> %; -- pretend we're a EuclideanDomain 41 gcd: (%, %) -> %; 42 divide: (%, %) -> (%, %); 43 sqrt: (%) -> %; 44 45 export from I; 46 47 } 48 49 == add { 50 Rep ==> Record(numer: I, denom: I); 51 52 import from Rep; 53 54 default a, b: %; 55 56 ratio(n: I, d: I): % == per [n, d]; 57 58 -- Use continued fractions to keep denominators small. 59 60 round( x: Rep) : % == { 61 local a0,a1,q,r,p0,p1,p2,q0,q1,q2,s : I; 62 a0 := abs x.numer; 63 a1 := abs x.denom; 64 s := sign(x.numer)*sign(x.denom); 65 66 (q,r) := divide(a0,a1); 67 68 a0 := a1; a1 := r; 69 70 p0 := 1; p1 := q; 71 q0 := 0; q1 := 1; 72 73 74 while r ~= 0 repeat { 75 76 (q,r) := divide(a0,a1); 77 78 a0 := a1; a1 := r; 79 80 p2 := q*p1 + p0; 81 q2 := q*q1 + q0; 82 83 q2 > MAXDENOM => break; 84 85 p0 := p1; p1 := p2; 86 q0 := q1; q1 := q2; 87 } 88 89 per [s* abs p1, abs q1] 90 } 91 92 93 -- Public Part -- 94 95 -- 1) Needed for all types (from BasicType) 96 97 (stream: TextWriter) << (z: %) : TextWriter == 98 stream << numer z << "/" << denom z ; 99 100 -- 2) Operations necessary to qualify as an OrderedRing 101 -- (though we don't guarantee the ring axioms hold) 102 103 (a: %) = (b: %): Boolean == numer a = numer b and denom a = denom b; 104 (a: %) ~= (b: %): Boolean == ~(a = b); 105 (a: %) < (b: %): Boolean == denom b * numer a < denom a * numer b; 106 (a: %) <= (b: %): Boolean == denom b * numer a <= denom a * numer b; 107 (a: %) > (b: %): Boolean == denom b * numer a > denom a * numer b; 108 (a: %) >= (b: %): Boolean == denom b * numer a >= denom a * numer b; 109 max(a: %, b: %): % == if a > b then a else b; 110 min(a: %, b: %): % == if a < b then a else b; 111 112 sign (a: %): % == ratio(sign numer a, 1); 113 abs (a: %): % == if negative? a then -a else a; 114 zero? (a: %): Boolean == zero? numer a; 115 negative? (a: %): Boolean == negative? numer a; 116 positive? (a: %): Boolean == positive? numer a; 117 118 -- 3) Operations necessary to qualify as an ArithmeticSystem 119 120 -- Default coercion used for coerce(n: Integer) : % 121 -- unless I=Integer. Decide that at compile time. 122 123 -- coerce(n: I) :% == ratio(n, 1); 124 coerce(n: Integer) :% == ratio(n::I, 1); -- craps out if I=SingleInteger 125 coerce(n: SingleInteger) :% == ratio(n::I, 1); 126 127 0: % == ratio(0,1); 128 1: % == ratio(1,1); 129 130 +(a: %): % == a; 131 -(a: %): % == ratio(-numer a, denom a); 132 133 (a: %) + (b: %): % == { 134 z := [numer a * denom b + numer b * denom a, 135 denom a * denom b ]; 136 round z; 137 } 138 139 (a: %) - (b: %): % == a + (-b); 140 141 (a: %) * (b: %): % == round( [ numer a * numer b, denom a * denom b] ); 142 143 144 -- Use binary powering, with rounding at each step. 145 146 (alpha: %) ^ (n: Integer): % == { 147 s: % := alpha; 148 m: Integer := abs n; 149 len: SingleInteger := length m; 150 c: % := 1; 151 { m=0 => c := 1; 152 m=1 => c := s; 153 -- otherwise (loop invariant `answer' = c*s^v, v=m initially) 154 { for i:SingleInteger in 1..(len-1) repeat { 155 if bit(m,i-1) then c := c*s; -- v <== v - 1 156 s := s*s; -- v <== v/2 157 } 158 c := c*s; -- bit(m, length(m) - 1) = 1 always. 159 } 160 161 } 162 if n >= 0 then c else round [denom c, numer c] 163 } 164 165 -- 4) Extra operations. 166 167 numer(a: %): I == rep(a).numer; 168 denom(a: %): I == rep(a).denom; 169 170 (a: %) quo (b: %): % == a/b; 171 (a: %) rem (b: %): % == 0; 172 gcd(a:%, b:%): % == 1; 173 divide(a:%, b:%): (%, %) == (a/b, 0); 174 175 -- We experiment here with allowing denormalized numbers to appear 176 -- temporarily as a result of "inv". 177 178 inv(a: %): % == per( [denom a, numer a] ); 179 180 -- The results of other operations, however, are rounded. 181 182 (n: I) / (d: I): % == round( [n, d] ) ; 183 (a: %) / (b: %): % == round( [ denom b * numer a, numer b * denom a] ); 184 (a: %) \ (b: %): % == round( [ denom a * numer b, denom b * numer a] ); 185 (n: I) + (b: %): % == ratio(n, 1) + b; 186 (a: %) + (n: I): % == a + ratio(n, 1); 187 (n: I) - (b: %): % == ratio(n, 1) - b; 188 (a: %) - (n: I): % == a - ratio(n, 1); 189 (n: I) * (b: %): % == ratio(n, 1) * b; 190 (a: %) * (n: I): % == a * ratio(n, 1); 191 (n: I) / (b: %): % == ratio(n, 1) / b; 192 (a: %) / (n: I): % == a / ratio(n, 1); 193 194 sqrt(x:%) :% == { 195 x < 0 => error "negative square root"; 196 guess: % := 1; 197 residual: % := guess * guess - x; 198 while abs residual > 100/(MAXDENOM * MAXDENOM) repeat { 199 guess := guess - residual/(2 * guess); 200 residual := guess * guess - x; 201 } 202 guess 203 } 204 205 }
22.6.2 : Lessons learned in writing the program
%^
Integer have to be made.
(Though you may make them give an error message, as Ratio used to do.)
% aldor -gloop %1 >> #include "aldor" %2 >> ArithmeticSystem () @ Third( with 0: % +: % -> % +: (%, %) -> % -: % -> % -: (%, %) -> % 1: % *: (%, %) -> % ^: (%, Integer) -> % coerce: Integer -> % coerce: SingleInteger -> %) == with 0: % +: % -> % +: (%, %) -> % -: % -> % -: (%, %) -> % 1: % *: (%, %) -> % ^: (%, Integer) -> % coerce: Integer -> % coerce: SingleInteger -> %
So we must, necessarily, provide ten operations (those listed after
the first ``with''). This ArithmeticSystem
will
provide more functionality, so we will have to add more
operations.
22.6.3 : Testing
This implementation of a type needs to be tested. Evaluating a polynomial
at 1000 points on an interval and comparing the results with the exact results
might be one way to do it. Another is to do the compact finite difference problem
with ``Rounded Rationals'' instead of DoubleFloat
s, and see if the
errors are of similar size.
1 #include "aldor.as" 2 --#assert broken 3 4 -- This is a double test. 5 -- 1) You need to say .ao in the #library command (confusion with .lib?) 6 -- 2) compile with aldor -Fx -M no-mactext testround.as RoundedRatio.ao 7 -- no-mactext means don't point at the macro definition when you find 8 -- an error in the text 9 -- 3) 1@R means `I know this is ambiguous but pick the 1 from R'. It's 10 -- better than 1::R because that might mean coerce the SingleInteger 1 11 -- to an R or coerce the Integer 1 to an R. 12 13 #library rrlib "rratio.ao" 14 #library compactlib "compact.ao" 15 16 17 SI ==> SingleInteger; 18 import from CommandLine, Array String; 19 20 strarg:String := if #arguments = 1 then arguments(1) else "10000"; 21 k:SI == (retract) (scanNumber strarg); 22 23 F ==> RoundedRatio (SI, k); -- Errors should be less than 10^(-6) 24 SI==> SingleInteger; 25 AF==> Array F; 26 27 import from rrlib, compactlib, Integer, SI, F, AF, InFile, StandardIO, 28 NumberScanPackage Integer, String; 29 30 print << "Enter n : "; 31 line: String := readline! stdin; 32 n:SI := (retract$Integer) (scanNumber line); 33 print << "n is " << n << newline; -- Good people echo their input. 34 35 h:F := 1@F/(n::F); 36 37 x:AF := new(n, 0); -- x values 38 f:AF := new(n, 0); -- function values 39 xdf:AF := new(n, 0); -- exact derivatives 40 adf:AF := new(n, 0); -- approximate derivatives 41 42 fun(t:F):F == { 1@F/(1@F + t*t) }; 43 dfun(t:F):F == { -2*t/(1@F + t*t)/(1@F + t*t) }; 44 45 for i:SI in 1..n repeat { 46 x(i) := (i::F)*h ; 47 f(i) := fun(x(i)); 48 xdf(i) := dfun(x(i)); 49 } 50 51 rt3 := sqrt(3@SI::F); 52 adf := derivatives(F, f, h, rt3); 53 54 errors : AF := new(n,0); 55 for i:SI in 1..n repeat errors(i) := abs (xdf(i)/(h*h*h*h) - adf(i)/(h*h*h*h)) ; 56 57 -- The following bit of code is modelled on the list routine "reduce". 58 maxerror(e:AF):F == { 59 ans:F := 0; 60 for er in e repeat ans := max(ans,er); 61 ans 62 } 63 print << "The errors divided by h^4 are less than or equal to " << maxerror(errors) << newline;
Typical output is shown below.
[1.402] % testroun Enter n : 5 n is 5 Using <978122/564719> for the square root of 3 The errors divided by h^4 are less than or equal to 3.6281588780112313 [1.403] % testroun Enter n : 10 n is 10 Using <978122/564719> for the square root of 3 The errors divided by h^4 are less than or equal to 18.440609676946824 [1.404] % testroun Enter n : 20 n is 20 Using <978122/564719> for the square root of 3 The errors divided by h^4 are less than or equal to 15.856048015270291 [1.405] % testroun Enter n : 100 n is 100 Using <978122/564719> for the square root of 3 The errors divided by h^4 are less than or equal to 38.027868204740237 [1.406] % testroun Enter n : 400 n is 400 Using <978122/564719> for the square root of 3 The errors divided by h^4 are less than or equal to 453313.35558103083
The similarity of behaviour is a good indicator that things are (now) working
correctly. I had to add a routine root3
which calculated the
square root of 3 as a rounded rational. I could have written a general
square root program, but it seemed easiest to use the known simple continued
fraction, = 1 + [1,2,1,2,1,2,...], to do the job here.
22.6.4 : Lessons learned in writing the program
SingleIntegers
which are too large
become negative (that's not a bug, I'm told, it's a feature).
This leads to incorrect results if your intermediate results are large.
Since I don't want incorrect results, however quickly they are obtained,
I would like to see some feature in Aldor to allow me to generate
an error message if a SingleInteger
overflows. It has been
suggested that what I wanted was CheckedSingleInteger
(which will
presumably be implemented in a later release), and I concur.
I have fixed the program so that choosing MAXINTSIZE
smaller
than 32768 will avoid any of these overflows with SingleInteger
.
This provides roughly the same precision as ``single precision''
floating point arithmetic, namely
that numbers can be represented accurately to 1 part in 230.
In searching for the cause of that problem,
I was forced to go over the code for
RoundedRatio
very carefully, and I found several potential bugs:
in particular I found some places where I was not careful to keep the
denominators positive. (I also discovered that you can't name your
program file what, apparently because it interferes with a Unix
command of the same name.)
.ao
in the #library
command, and to
include the .ao
file in the compile command.
-M no-mactext
means that errors involving
macros are pointed to where they occur, and not where the macro is defined.
This made several error messages much clearer.
1
is ambiguous,
but pick the 1
from R
''.
It's different from 1::R, which might mean coerce the
SingleInteger
1
to type R
or coerce the
Integer
1
to type R
.
CommandLine
, Array(String)
, InFile
and StandardIO
types, together with the NumberScanPackage
.
DoubleFloats
is easy. It ought to be only
slightly harder to write a coercion the other way. (Continued fraction
algorithms are helpful here too.)
1/24/3
instead of (1/2)(4/3)
-- obviously spaces might help
too). Ratio Integer
s are printed with parentheses (( )).
I chose to print Rounded Rationals with angle brackets (< >).
22.7 : Retrospective
Though I have really only been programming in Aldor for about three weeks (not counting time spent doing other things), I don't consider myself an Aldor novice anymore. Of course, I still have an enormous amount to learn (I now know that I haven't even peeled off more than one or two layers of the ``onion'') but I can now make Aldor do nontrivial tasks, efficiently and easily. I want to use Aldor to investigate Gröbner basis algorithms for polynomials with rounded rational coefficients, to solve delay differential equations and to implement an efficient and useful interval arithmetic. I have some confidence that all these will be possible -- and interesting.
I also have the feeling that Aldor is an efficient, powerful and elegant language. It is still in its infancy, so I am really talking about ``potential'' here. But if this is an infant, what will the adult be like?