The Maths Library
(v.5
171016)
© W O Riha

1. Introduction The Shen maths library
contains all the functions found in the C  maths
library (C89 standard) and several other
operations and functions which may be available
in one form or another as builtin operations in
other languages, but are not native to Shen (for
example, modulus and integer division). The
entire library has been coded in Shen, and hence,
is portable.
Certain functions
and operations which, mathematically speaking,
only involve ‘whole numbers’ or
‘integers’ as arguments and function
values were originally contained in a separate
‘integer library’, which also provided
a subtype integer. As explained elsewhere, there
now is no type integer, and so these functions
have been included in the maths library. They are
described in Section 4.
A user of the
maths library (or any library, for that matter)
will be familiar with at least some of the
functions that are provided. For this reason, I
have not explained in any detail exactly what
certain functions compute, or how they are
defined. For example, I assume that a potential
user of the library will observe that the
trigonometric functions sin, cos and tan are
available (but not cot), and knows what they are
used for, and so I do not define these functions.
In the same vein, I do not explain the purpose of
atan2 (but provide a reference), because if you
don’t know what it does, you can’t use
it.
On the other
hand, I do give explanations and/or definitions,
when I am aware, or suspect, that there is an
ambiguity, misconception or confusion as regards
the meaning and purpose of a function. (See
Section 3.3 on ‘Rounding’, and the
‘Tutorial on Integer Division’ in
Section 4).
The examples
given in the text are based on SBCLShen.
Numerical results agree with those produced by an
early version of JSShen, but not necessarily,
where large numbers (integers) are involved.
Tests on other platforms have not been carried
out, because easily accessible images are not
available.
Unless stated
otherwise, values are quoted and displayed (by
default) to 16 or 17 places. This will invariably
exhibit rounding errors.
Example:
(
1.37 1)
0.3700000000000001 : number


2. Data Type number
number is one of the
basic data types in Shen (see Shen Document Numbers). For
the sake of simplicity, this type is equipped with only a
few functions (see Shen Document The Primitive Functions
of KLambda).
The arithmetic operations
all have the signature number > number > number
The following functions
are available (the first four are builtin – unary
minus has been defined for convenience)
a. + addition
b. – subtraction
c. * multiplication
d. / division (if Y is 0 , (/ X Y) raises the error
“division by zero”
e. ~ unary – (~ X) is equivalent to (– 0 X)
Notice that + and * are
polyadic functions.
(+ 7 8 9
10)
34 : number
(* 5.9 8.2 9.4)
454.772 : number
The comparison operations
all have the signature number > number >
boolean
The following functions
are provided:
a. = equal (equality
is defined for all native types).
b. != not equal (!= is defined for all types as
shorthand for (not (= X Y)))
c. < less than
d. > greater than
e. <= less than or equal
f. >= greater than or equal
Note: The comparison
operators are affected by (floatingpoint) rounding
errors, in particular = and !=, however, this is
perfectly normal.
(= 12 (*
0.006 2000))
true : boolean
(= 12 (* 0.0006 20000))
false : boolean
The system predicate integer? : A >
boolean
tests if an object is an ‘integer’. Shen
integers are different from integers in other programming
languages!
(integer?
12.0)
true : boolean
(integer? (+ 1.2 2.8))
true : boolean
(integer? 15.67e3)
true : boolean
(integer? (* 0.42 100))
false : boolean
(integer? (* 0.41 100))
true : boolean
3. Maths Library Functions
3.1
Predicates
In addition to the
predicate integer?, an equivalent predicate
int? :
number > boolean
has been provided. This
function runs in time O (log N) as opposed to integer? which runs in O (log^{2}
N). For smaller values of N, and/or a small number of
invocations there is no significant difference in run
time. However, for larger values of N the difference is
quite noticeable.
(time
(for I = 1 to 100000 (int? 123456789012.1)))
run time: 8.20 secs
(time (for I = 1 to 100000 (integer?
123456789012.1)))
run time: 27 secs
natural? :
number > boolean
Input: A number X.
Output: true if X is a nonnegative integer, otherwise
false.
(natural?
12.5)
false : boolean
(natural? 100)
true : boolean
but also
(natural?
1.54e5)
true : boolean
positive? :
number > boolean
Input: A
number X.
Output: true if X > 0, otherwise false.
(positive?
12.5)
true : boolean
negative? :
number > boolean
Input: A number X.
Output: true if X < 0, otherwise false.
(negative?
12.5)
false : boolean
even? :
number > boolean
Input: A number X.
Output: true if X is an even integer, otherwise false.
(even?
2.5)
false : boolean
(even? 12)
true : boolean
odd? :
number > boolean
Input: A
number X.
Output: true if X is an odd integer, otherwise false.
(odd?
12)
false : boolean
(odd? 13)
true : boolean
Note: not even is not
equivalent to odd:
(odd?
13.3)
false : boolean
(not (even? 13.3))
true : boolean
3.2 Sign and Absolute Value
sign :
number > number
Input: A
number X.
Output: The sign of X, which is defined as follows
sign (X) = 1 if X > 0
0 if X = 0
1 if X < 0
(sign
123.5)
1 : number
abs : number
> number
Input: A
number X.
Output: The absolute value, X, i.e. the value X × sign
(X).
(abs
123.5)
123.5 : number
(abs 0.0)
0 : number
3.3 Rounding Functions
Rounding means replacing
a numerical value by another value that is approximately
equal, but is shorter or simpler. For example, rounding
pi = 3.1415926535… to three decimal places gives the
value “3.142”; rounding 29.15 to the nearest
integer gives the value “29”.
The most fundamental case
is that of rounding a (real) number to an integer. There
are five methods that are in common use (see, for
example, Rounding).
3.3.1 Rounding Up or Down
Definition: The largest
integer not greater than a, called floor of a, is denoted
by _{}a_{}
or floor(a). Similarly, the smallest integer not
less than a, called ceiling of a, is denoted by ^{}a^{}
or ceiling (a).
Here, floor (a) rounds
down to the nearest integer (or towards minus infinity),
whilst ceiling (a) rounds up (or towards plus infinity).
The two functions have many interesting mathematical
properties (see Floor and Ceiling Functions). We shall
later make use of three of them
P1. _{}a_{}
= a iff a is an integer
P2. ^{}a^{}
= _{}a_{}
+ 1 if a > 0
P3. _{}a_{}
=  ^{}a^{}
for any a
Note: Property P1 is the
basis of the system predicate integer? (and int?).
floor :
number > number
Input: A
number X.
Output: The integer value ^{}X^{}.
(see Table 1)
ceiling :
number > number
Input: A
number X.
Output: The integer value _{}X_{}
(see Table 1)
Table 1
X 
(floor X) 
(ceiling X) 
11.7 
11 
12 
12 
12 
12 
11.7 
12 
11 
12 
12 
12 
3.3.2 Rounding Towards Zero/Away
from Zero
A real number such as
123.89 is composed of an ‘integer part’, 123,
and a ‘fractional part’ 0.89. If the number is
negative, e.g. 123.89 then, obviously, the integer part
is 123. But what is the fractional part? It could be
0.89 or 0.89. In order to preserve the identity integer
part (a)+ fractional part (a)= a one usually opts for
0.89. This motivates the following.
Definition: The number
sign (a) × floor (abs (a)) is called the integer part of
a, denoted by intpart (a). The number a – intpart
(a), is called the fractional part of a, and is denoted
by fracpart (a).
It is not difficult to
see that
intpart(a) = floor(a)
where a = 0 or a > 0
= ceiling(a) where a
< 0
which means that intpart
(a) rounds towards zero. [intpart is often called
truncate (or trunc), because the value is obtained by
ignoring (throwing away) the fractional part].
trunc :
number > number
intpart : number > number
Input: A number X.
Output: The integer part of X. (see Table 2)
fracpart :
number > number
Input: A
number X.
Output: The fractional part of X.
Table 2
(Note: fractional values have been rounded)
X 
(intpart X) 
(fracpart X) 
11.7 
11 
0.7 
12 
12 
0.0 
11.7 
11 
0.7 
There is also a function
(with a rather unfortunate name) which combines the
preceding two functions
modf :
number > (number * number)
Input: A
number X.
Output: The pair consisting of the integer and
fractional parts of X.
(modf
11.5)
(@p 11 0.5) : (number * number)
The counterpart to
truncation (‘rounding towards zero’) is
‘rounding away from zero’. This operation has
no specific name and has not been implemented in the
library. It is given by
For examples, see Table 3
below.
3.3.4 Rounding to the Nearest
Integer
This way of rounding is
used most frequently. The idea is quite simple: either
take the floor or the ceiling of the number, whichever is
the closer. There is a problem, however, what to do if
the number is precisely in the middle, i.e. when its
fractional part is exactly equal to 0.5. There are
various ways to resolve this dilemma by means of a
‘tiebreak’. The main concern is to avoid or to
minimise a bias (for a discussion, see Rounding).
Some libraries may
provide a choice of different tiebreaks. The method that
is used in the Shen mathslibrary is the default rounding
mode of the IEEE 754 standard, which is also widely
employed in bookkeeping. The rule, known as ‘round
half to even’, is quite simple.
If fracpart (a) = 0.5,
then round to the even integer nearest to a.
For example, 13.5 is
rounded to 14, 12.5 to 12, 12.5 to 12, and 13.5 to
14.
round0 :
number > number (see below for a more convenient
alternative)
Input: A number X.
Output: X rounded to the nearest integer.
Example: (also see Table 3)
(round0
123.78)
124 : number
(round 1234.5)
1234 : number
Table 3
(from Rounding)
Comparison of the five integer rounding methods
N 
round
down (towards infinity) 
round
up (towards +infinity) 
round
towards zero 
round
away from zero 
round
to nearest 

floor 
ceiling 
trunc 


+23.67 
+23

+24

+23

+24

+24 
+23.50 
+23 
+24 
+23 
+24 
+24 
+23.35 
+23 
+24

+23

+24

+23 
+23.00 
+23 
+23

+23

+23

+23 
0 
0 
0 
0 
0 
0 
23.00 
23

23

23

23 
23 
23.35 
24

23

23

24 
23 
23.50 
24

23

23 
24 
24 
23.67 
24

23

23

24 
24 
3.3.5 Rounding to a Specified
Number of Decimal Places
round’
: number > number > number (see below for a more convenient
alternative)
Input: A number X and an integer N.
Output: X rounded to N decimal places, or an error
message, if N is not an integer.
(round'
123.7823457 5)
123.78235 : number
Note: The function also
works for negative integers N.
(round'
87451.167 2)
87500.0
Functions round0 and round’ are available as a single
function round with either one or two arguments
round :
number > number (corresponds to round0)
round :
number > number > number (corresponds to round')
(round
12.789)
13 : number
(round 12.789 2)
12.79 : number
Notice that (round X 0) is equivalent to (round X).
3.4 Mathematical Constants
Constants are defined as
0ary functions. Thus, to multiply sqrt(2) by pi/4 we can
write
(*
(sqrt2) (pi/4))
1.1107207345395917 : number
(e) =
2.71828182845904523536
(log10) = 2.30258509299404568402
(log2) = 0.69314718055994530942
(pi) = 3.14159265358979323846
(pi/2) = 1.57079632679489661923
(pi/4) = 0.78539816339744830962
(twopi) = 6.28318530717958647692
(one/pi) = 0.31830988618379067154
(two/pi) = 0.63661977236758134308
(log2e) = 1.44269504088896340736
(log10e) = 0.43429448190325182765
(sqrt2) = 1.41421356237309504880
(one/sqrt2) = 0.70710678118654752440
(two/sqrtpi) = 1.12837916709551257390
3.5 Degree Conversions
Angles can be measured or
specified either in degrees or in radians. Measured in
degrees, a full circle is 360°, a quarter is 90° (right
angle), whereas in radians a circle corresponds to 2 pi,
and a right angle to pi/2. Furthermore, 1° is equal to
60’ (60 minutes), and 1 minute is 60’’ (60
seconds). Specifying an angle in degrees, minutes and
seconds, or in fractional degrees (minutes or seconds) is
very common in astronomy and geodesy. In mathematics and
in numerical computations, radian is assumed by default.
The syntax of a degree
specification – what is allowed and what is not.
38°
15’ 23.78’’ 
ok 
38°
15.45’ 23’’ 
not ok
(should be 38° 15’ 50’’,
since 0.45’ equals 27’’) 
–38°
15.45’ 
ok 
38°
15.45’ 
not ok
(minutes must be nonnegative, if degrees are
positive) 
0°
15.45’ 
ok 
38.5°
15’ 23’’ 
not ok
(should be 38° 45’ 23’’) 
38°
75’ 23.78’’ 
not ok
(minutes must be < 60) 
38°
23.78’’ 
not ok
(should be 38° 0’ 23.78’’) 
38° 0’
23.78’’ 
ok 
0°46’
10.5’’ 
ok 
An angle is represented
as a list of up three numbers based on the
“grammar”:
<int59> := 0 1
 2 … 59 <int–59> := –59 
–58  …  0 1  2 … 59
<num60> := <number : 0 = number < 60 >
<numabs60> := <number : number < 60
>
<dlist> := [<number>]
<dmlist> := [<integer> <num60>]
<mslist> := [0 <int–59>
<num60>]
<mlist> := [0 <numabs60>]
<slist> := [0 0 <numabs60>]
<dmslist> := [<integer> <int59>
<num60>]
<angle> := <dlist>  <dmlist> 
<mslist>  <mlist>  <slist> 
<dmslist>
The following functions
are available:
dms>degs
: (list number) > number
Input: A
list of three numbers, intended to specify an angle in
degs, mins and secs.
Output: The corresponding angle in fractional degrees (or
an error message).
(dms>degs
[38 0 23.8])
38.00661111111111 : number
(dms>degs [38 10.4 23.8])
dms>degs  type error [38 10.4 23.8]
(dms>degs [0 10 30.88])
0.17524444444444445 : number
(dms>degs [0 0 30.88])
0.008577777777777778 : number
(dms>degs [10 20 30 40.8])
dms>degs  too many arguments [10 20 30 40.8]
The inverse of dms>degs is
degs>dms
: number > (list number)
Input: A number representing an angle in fractional
degrees.
Output: The angle as a list of degrees, minutes and
seconds.
(degs>dms
38.175)
[38 10 30] : (list number)
(degs>dms (dms>degs [30 39 59.9]))
[30 39 59.9] : (list number)
degs>rad
: number > number
Input: A
number representing an angle in fractional degrees
Output: The angle expressed in radians.
(degs>rad
30.0)
0.5235987755982988 : number \\ this is pi/6
(degs>rad 123.456)
2.1547136813421197 : number
And its inverse
rad>degs
: number > number
Input: A
number representing an angle in radians.
Output: The angle expressed in fractional degrees.
(rad>degs
(pi/4))
45.00000000000001 : number
(rad>degs (degs>rad 123.456))
123.45600000000002 : number
By combining these
functions other conversions are possible. To convert an
angle from radians to degs, mins and secs, we compose
rad>degs and degs>dms. For example, 1 radian has
the following degree equivalent
(degs>dms
(rad>degs 1))
[57 17 44.8062470964] : (list number)
3.6 Trigonometric Functions
The trigonometric
functions can be defined in terms of convergent infinite
series
sin x = x/1! – x^{3}/3!
+ x^{5}/5! – x^{7}/7! + …
for x < inf (3.6.1)
cos x = 1 – x^{2}/2! + x^{4}/4!
– x^{6}/6! + … for x < inf (3.6.2)
These series can be used
to compute sin and cos, provided that x is
“reasonably small”. There is also a series
expansion for tan, which is more complicated – it
involves the Bernouilli numbers – and does not lend
itself, easily, to numerical computation.
All trigonometric
functions expect arguments specified in radians (see
degree conversion). They are periodic with periods 2pi
(sin and cos) and pi (tan). Thus, sin x = sin (2x + 2pi)
= sin (x + 4pi) = …
There are many elementary
relationships between the trigonometric functions, for
example
sin (x) = sin (x) (3.6.3)
cos (x) = sin (x) (3.6.4)
cos x = sin (pi/2  x) (3.6.5)
tan x = sin x / cos x (3.6.6)
These can be used to
compute cos and tan in terms of sin, which is evaluated
via its series expansion.
3.6.1 Computing sin x :
In view of (3.6.3), we
may assume that x = 0. Furthermore, because of the
periodicity of sin, we only need to consider arguments x
in [0 2pi); x = fmod (y, 2pi) will do this, but if y is
large the rounding error incurred can be quite
significant (see below). If x is not already in the
interval [0 pi/2) then we make use of the following
properties:
sin x = sin (pi
– x) for x in [pi/2 pi)
sin x = –sin (pi – x) for x in [pi 3pi/2)
sin x = –sin (2pi – x) for x in [3pi/2 2pi)
We then proceed to
evaluate the sum (3.6.1), stopping when the next term is
smaller in magnitude than a prescribed tolerance (set as
10^{–15}). Since (pi/2)^{21}/21!
< 2.5810^{–16}, at most 11 terms will be
required.
Before returning the
result, the correct sign has to be chosen!
3.6.2 Computing cos x :
In view of (3.6.4), we
may assume that x >= 0. Furthermore, because of the
periodicity of cos, we only need to consider arguments x
in [0 2pi); x = fmod (y, 2pi) will do this, but if y is
large the rounding error incurred can be quite
significant. If x is not already in the interval [0 pi/2)
then we make use of the following properties:
cos x = sin (2pi
– x) for x in [pi 2pi)
cos x = – cos (pi – x) for x in [pi/2 pi)
We then proceed to
evaluate the sum (3.6.2), stopping when the next term is
smaller in magnitude than a prescribed tolerance (set as
10^{–15}). Since (pi/2)^{22}/22!
< 10^{–17}, at most 11 terms will be
required.
Before returning the
result, the correct sign has to be chosen! This is easier
than for sin.
sin : number
> number
Input: A
number X.
Output: sine of X – a value between 1 and +1,
inclusive.
(sin
(degs>rad 45))
0.7071067811865475 : number \\ The sine of 45° is
1/sqrt(2) – check!
(sin (pi))
2.4821693935912637e16 : number \\ should be 0
(sin 2.59)
0.5240443416872763 : number
(sin (+ 2.59 (* 5 (pi*2))))
0.5240443416872789 : number \\ should be same as (sin
2.59)
cos : number
> number
Input: A number X.
Output: cosine of X – a value between 1 and +1,
inclusive.
(cos (/
(pi) 3))
0.5000000000000001 : number \\ The cosine of 60° is
½.
tan : number > number
Input: A number X.
Output: tangent of X.
(tan (degs>rad 60))
1.7320508075688765 : number
Note: The tangent of 60° is sqrt(3).
(sqrt 3)
1.7320508075688774 : number
3.7 Inverse Trigonometric
Functions
Since the trigonometric
functions are periodic, their inverses are multivalued
functions, that is, for any given value, y, there exist
infinitely many values of x for which sin x = y, cos x =
y and tan x = y. For example, since sin has period 2pi,
we have, for its inverse arcsin
arcsin y = x + (2pi *
k), k is in Z, where sin x = y.
We say that arcsin has an
infinite number of branches. In order to define a unique
inverse, one selects a particular branch as the socalled
principal branch. For the inverse trigonometric functions
the principal branch corresponds to k = 0. The principal
branch of any multivalued function is usually denoted by
an initial capital letter, for example, Arcsin, Arccos,
Arctan. One can then express the multivalued functions
as
arcsin = Arcsin x +
(2pi * k), k is in Z etc
3.7.2 Computing arctan x
The inverse tan function,
arctan is defined for all real numbers, and assumes
values in the interval (pi/2, pi/2). It can be
represented by the MacLaurin series
arctan x = x –
x3/3 + x5/5 – x7/7 … x= 1 (3.7.1)
This expansion can be
used as a basis for computing arctan, at least if x is
“not too large” (not too close to 1). The
(alternating) series converges very slowly, when x is
near 1. In principle, it could be used to compute p,
since arctan 1 = 1 – 1/3 + 1/5 – 1/7 … =
pi/4.
But in order get an
accuracy of about 6 decimal places, one would need to
consider some 500000 terms.
For arguments x, greater
than 1 in magnitude (3.7.1) is not applicable. But we can
make use of the formula
arctan x = pi/2
– arctan (1/ x) x > 1 (3.7.2)
As to the actual
computation of arctan, we first note that
arctan (x) = 
arctan x.
Therefore we may assume
that x >= 0.
Now, since (3.7.1) is no
good when x is near 1, we need to transform x to some
smaller value (and still be able to work out arctan).
The addition
(subtraction) formula for tan is
tan (x – y) =
(tan x – tan y) / (1 + (tan x)(tan y))
With x = pi/4 and y = x,
we obtain, noting that tan pi/4 = 1
tan (pi/4 – x) =
(1 – tan x)/(1 + tan x)
Replacing tan x by x, and
applying arctan on both sides, yields
pi/4 – arctan x
= arctan (1– x)/(1 + x) (3.7.3)
The function x >
(1– x)/(1 + x) is monotone decreasing from 1 to 0.
So, the closer x is to 1 the faster the RHS of (3.7.3)
will converge, when expressed as the series (3.7.1). The
question is, when should we use (3.7.1) and, exactly
when, should we switch to (3.7.3)?
Since the function x
> x increases in [0 1] and x > (1– x)/(1 +
x) decreases, the point where the two functions
intersect, is the deciding factor. So we must solve the
equation
x = (1– x) / (1 + x)
This leads to the
quadratic equation x^{2} + 2x – 1 = 0, which
has the positive solution sqrt (2) – 1 =
0.4142…
To obtain a 16digit
accuracy, we never have to work out more than 21 terms.
asin :
number > number
Input: A
number X, where –1 = X = +1.
Output: arcsine of X – a value between –pi/2
and +pi/2, inclusive, or an error message.
(asin
(sin 0.89))
0.8899999999999999 : number
but
(asin
(sin 2.89))
0.2515926535897928 : number
which is equal to
( (pi)
2.89)
0.251592653589793 : number
(asin 1.5)
asin(x) for x > 1!
acos :
number > number
Input: A number X, where –1 = X = +1.
Output: arccosine of X – a value between 0 and pi,
inclusive, or an error message.
(acos
1)
3.141592653589793 : number \* should be pi =
3.14159265358979323846 *\
(acos 1)
0.0 : number
(acos 1.3)
acos(x) for x > 1!
(cos (acos 0.789))
0.7890000000000001 : number
atan :
number > number
Input: A
number X.
Output: arctangent of X – a value strictly between
–pi/2 and +pi/2.
(atan
1.5e15)
1.570796326794896 : number \* close to –pi/2 *\
(atan (tan 0.45))
0.4500000000000001 : number
(tan (atan 2.45))
2.449999999999999 : number
atan2 :
number > number > number
Input: Two numbers X and Y.
Output: The twoargument arctangent of Y and X – a
value between –pi (not included) and +pi
(inclusive).
(atan2
3.5 0)
1.5707963267948966 : number
(atan2 3.5 2.3)
2.152176510599796 : number
(atan2 0 0)
atan2 – undefined
(atan2 0 4.7)
3.141592653589793 : number
3.8 Exponential and Hyperbolic
Functions
exp : number
> number
Input: A
number X.
Output: The value of e^{X}.
(exp
1.5)
4.481689070338065 : number
(exp 1.0)
0.3678794411714423 : number
(exp 1.0)
2.7182818284590455 : number
(exp 100.0)
2.688117141816169e43 : number
sinh :
number > number
Input: A number X.
Output: The value of the hyperbolic sine of X.
(sinh
10.12)
12417.385397399632 : number
(sinh 8.23)
1875.9167427768757 : number
(sinh 0)
0 : number
cosh :
number > number
Input: A number X.
Output: The value of the hyperbolic cosine of X.
(cosh
10.12)
12417.385437665756 : number
(cosh 8.23)
1875.917009313206 : number
(cosh 0)
1 : number
tanh :
number > number
Input: A
number X.
Output: The value of the hyperbolic tangent of X.
(tanh
8.23)
0.9999998579167795 : number
(tanh 1)
0.761594155955765 : number
3.9
Logarithms
Mathematically most
important is the inverse of exp, called the natural
logarithm, or logarithm to base e. It is
commonly denoted by log_{e }(or simply log or
ln). By definition
log_{e} e^{X}
= e^{log}e^{X} = X
Sometimes, logarithms to
other bases, especially 2 and 10 are needed. If the base
is b (> 0) then
log_{b} b^{X}
= b^{log}b^{X} = X
The library provides the
following log functions
log : number
> number
Input: A positive number X.
Output: The value of the natural logarithm of X, or an
error message if X <= 0.
(log
(e))
1.0000000000000004 : number
(log (exp 12.45))
12.45 : number
(log 5.1)
(log 5.1) is not a real number!
log2 :
number > number
Input: A positive number X.
Output: The value of the base2 logarithm of X, or an
error message if X <= 0.
(log2
1024)
9.999999999999998 : number
(log2 10)
3.321928094887361 : number
log10 :
number > number
Input: A
positive number X.
Output: The value of the base10 logarithm of X, or an
error message if X <= 0.
(log10
(power 10 17))
16.999999999999996
(log10 2)
0.3010299956639811
Logarithms for general
base b are computed by the formula log_{b} x =
log x/log b. In Shen, by calling log with two arguments.
For negative b, the logarithms are therefore always
complex numbers, for example,
log_{5} 25 =
0.41578048 – 0.8115957i
notwithstanding that log_{5}
25 = log_{5} (5)^{2} = 2. Log is
multivalued for negative arguments and bases.
log : number
> number > number
Input: Two
positive numbers X and B.
Output: The value of the baseB logarithm of X, or an
error message if X <= 0 or B <= 0.
(log
234.89 9) \\ base 9
2.4845513634614647 : number
(log (expt 17 3.456) 17) \\ base 17
3.456 : number
(log 25 5) \\ base 5
log with base 5 is not a real number!
13.10 Inverse Hyperbolic Functions
asinh :
number > number
Input: A
number X.
Output: The inverse of (sinh X).
(asinh
3.788)
2.0419697670129047 : number
(asinh (sinh 12.3))
12.3 : number
acosh :
number > number
Input: A
number X = 1
Output: The inverse of (cosh X), or an error message if X
< 1
(acosh
(cosh 12.567))
12.567 : number
(acosh 0.5)
(acosh 0.5) is not a real number!
atanh :
number > number
Input: A
number X where X < 1
Output: The inverse of (tanh X), or an error message if
X >= 1
(atanh
0.5)
0.5493061443340549 : number
(atanh 1.5)
(atanh 1.5) is not a real number!
13.11 Power Functions
square :
number > number
Input: A
number X.
Output: The square of X.
(square
12.89)
166.15210000000002
(square 500)
250000
A fast integer power
power :
number > number > number
Input: A
number X and an integer N.
Output: The value of X^{N} or an error message,
if N is not an integer.
Note: The computation runs in log N  time. We have
defined 0^{0} = 1. (See note under expt below)
(power
1.5 10)
57.6650390625 : number
(power 1.5 10.9)
power  exponent must be an integer
(power 1.5 11)
0.011561019943888409 : number
The general power
function X^{Y}
expt :
number > number > number
Input: Two
numbers X and Y.
Output: The value of X^{Y} or an error message,
if X^{Y} is not defined or not a real number.
Note: We have defined that 0^{0} = 1. There is
considerable controversy whether it should be 1, 0 or
undefined. (click here for a discussion).
(expt
(pi) (/ 1 (pi))) \* p1/p *\
1.4396194958475907 : number \* the figures quoted are
all correct! *\
(expt 1.5 1.2)
(expt 1.5 1.2) is not a real number!
(expt 1.5 3)
3.375 : number
(expt 1.5 0)
1 : number
(expt 0 3)
(expt 0 3) is undefined!
Function expt can of course be used to compute
square roots (and other roots). A special function is
provided which computes square roots using an iterative
algorithm. The accuracies of the two methods are very
similar.
sqrt :
number > number
Input: A
nonnegative number X.
Output: The value of sqrt(X), or an error message, if X
is negative.
(sqrt
1000000)
1000.0 : number
(expt 1000000 0.5)
999.9999999999994 : number
The square root of the
golden ratio: ((sqrt(5) + 1)/2)^{½} = 1.272 019
649 514 068 962 … (see Mathematical Constants)
(sqrt (/
(+ (sqrt 5) 1) 2))
1.272019649514069 : number
(expt (/ (+ (expt 5 0.5) 1) 2) 0.5)
1.2720196495140688 : number
13.12 Various Functions
The following function,
called ‘free the exponent’ splits a number into
its mantissa (significand) and exponent.
frexp :
number > (number * number)
Input: A
number X.
Output: A pair (@p M E) such that the relationship M × 2^{E}
holds, where E is an integer and 0.5 <= M < 1, or
M = 0.
(frexp
0.1)
(@p 0.8 3) : (number * number)
(frexp 256.0)
(@p 0.5 9) : (number * number)
(frexp 0.0)
(@p 0 0) : (number * number)
An ‘inverse’ of
frexp, called ‘load exponent’
is
ldexp :
number > number > number
Input: A number M and an integer E.
Output: The value M × 2^{E} , or an error
message if E is not an integer.
(ldexp
0.8 3)
0.1
(ldexp 1.5 9.1)
power  exponent must be an integer
(ldexp 0.95 4)
15.200000000000001
The following function
called ‘floatingpoint modulus’ is a
generalisation of the remainder function rem for integers (see Section 4.1.3)
fmod :
number > number > number
Input: Two numbers A and B.
Output: The value A – K × B, where K is an integer
such that the result is 0 or has the same sign as A, and
magnitude less than the magnitude of B.
(fmod
123.45 17.4)
1.65000000000002 : number
(fmod 417.2 29.8)
0.0 : number
factorial :
number > number
Input: A
nonnegative integer N.
Output: The value of N!, or if N < 0, the value 1. For
nonintegral positive N, returns I^{}I
(N  K).
(factorial
10)
3628800 : number
(factorial 4)
1 : number
(factorial 4.5)
59.0625 : number \\ = 4.5 x 3.5 x 2.5 x 1.5
Two polyadic functions
are available
max :
number > . . . > number
Input: Any
number of numerical values A, B, C . . .
Output: The maximum value of the given arguments.
(max
12.3 1 33.4567 2.45)
33.4567 : number
(max 3 2)
3 : number
(max 12)
12 : number
min : number
> . . . > number
Input: Any number of numerical values A, B, C .
. .
Output: The minimum value of the given arguments.
(min
12.3 1 33.4567 2.45)
1 : number
(min 3 2)
2 : number
(min 12)
12 : number
4. Integer Functions
4.1 Integer Division and Remainder
– A Short Tutorial
4.1.1 Mathematical Background
In general, dividing a
whole number by another, does not result in a whole
number. For example, 19/5 = 3.8, which is not an integer.
In many practical situations it may be necessary to
pretend that the result is in fact an integer
‘close’ to the actual result, as in the
situation when 19 eggs have to be distributed equally
amongst 5 people. The mathematical answer is to give each
person 3.8 eggs, which is clearly impossible and
preposterous. In practice, one will probably give each
recipient 3 eggs with 4 eggs left over. Thus 19 = 5 × 3
+ 4, where 3 is the ‘integer quotient 19/5’ and
4 is the ‘remainder’.
Programming languages
often provide operations for finding
‘remainders’, variously called mod, rem, % (see
below), but not always operators for ‘integer
division’. Examples are LISP and JavaScript. Other
programming languages provide both: C/C++, Pascal,
Python, and often two different remainder operations (see
table in Modulo operation). The semantics of these
operations varies from language to language, and % in one
can have a different meaning in another. This, of course,
is confusing and can cause enormous problems when
translating programs between languages. The different
semantics will be explained below.
The mathematical basis of
integer division rests on a fundamental theorem in
elementary number theory:
Theorem:
Let a, b in Z be two integers, with b > 0. There exist
q, r in Z, such that
a = qb + r, 0 <= r
< b (4.1)
q and r are unique. q is
said to be the (integer) quotient, and r the remainder,
on dividing a by b.
Notice the condition b
> 0. The result can be generalised as follows (see Division algorithm)
Theorem:
Let a, b in Z be two integers, with ~(b = 0). There exist
q, r in Z, such that
a = qb + r, 0 <= r
< b (4.2)
Note that (here) the
remainder, r, is always nonnegative and satisfies
–b < r < b. This is known as the Euclidean
remainder (according to Raymond T. Boute, see Modulo operation).
4.1.2 Integer Division and
Remainder in Programming Languages
How do programming
languages deal with integer division and remainder? The
simple answer is: in various ways! All programming
languages tend to agree that the remainder should not be
greater than the (absolute value of the) divisor, i.e.
that r < b. There are a number of possibilities,
besides the above, Euclidean definition of quotient and
remainder (which is only used in Algol68, some versions
of Pascal and in Stata). Note that the positivity of the
remainder (4.2) is useful in certain applications, e.g.
calendar computations.
To satisfy the condition
r < b, it is possible to choose the sign of r.
This will then define the value of the quotient q, and
hence, the meaning of “integer division”.
Name
1. sign(r) = sign(a)
remainder has the same sign as the dividend rem
2. sign(r) = sign(b) remainder has the same sign as the
divisor mod
3. r = 0 remainder is nonnegative (Euclidean) % (not in
common use!)
Some programming
languages provide both (1) and (2) [see Modulo operation]
Here is an example that
illustrates the three choices
Example: Let a = ±100, b
= ±14
1. Truncated integer
division – remainder has same sign as the
dividend
Used in versions of C/C++, C#, Java, JavaScript, Lisp
Here q = trunc (a/b)
r = rem (a, b)
a 
b 
a/b 
q 
r 
a = q ×
b + r 
100 
14 
7.142… 
7 
2 
100 = 7
× 14 + 2 
100 
14 
7.142… 
7 
2 
100 =
(7) × (14) + 2 
100 
14 
7.142… 
7 
2 
100 =
(7) × 14 –2 
100 
14 
7.142…

7 
2 
100 = 7
× (14) –2 
2. Floored integer
division – remainder has same sign as the
divisor
Used in Clojure, Mathematica, Python, Lisp.
Here q = floor (a/b)
r = mod (a, b) This is perhaps the most useful
definition.
a 
b 
a/b 
q 
r 
a = q ×
b + r 
100 
14 
7.142… 
7 
2 
100 = 7
× 14 + 2 
100 
14 
7.142… 
8 
12 
100 =
(8) × (14)  12 
100 
14 
7.142… 
8 
12 
100 =
(8) × 14 + 12 
100 
14 
7.142… 
7 
2 
100 = 7
× (14) – 2 
3. Euclidean integer
division – remainder is always nonnegative
Used rarely, because it is awkward to implement
(Algol68, some versions of Pascal, Stata)
Here
r = % (a, b)
a 
b 
a/b 
q 
r 
a = q ×
b + r 
100 
14 
7.142… 
7 
2 
100 = 7
× 14 + 2 
100 
14 
7.142… 
7 
2 
100 =
(7) × (14) + 2 
100 
14 
7.142… 
8 
12 
100 =
(8) × 14 + 12 
100 
14 
7.142… 
8 
12 
100 = 8
× (14) + 12 
Observation 1
for a > 0, b
> 0 (1), (2) and (3) agree (4.3)
for a > 0, b < 0, (1) and (3) agree
for a < 0, b > 0, (2) and (3) agree
for a < 0, b < 0, (1) and (2) agree
Observation 2:
for b > a = 0
all quotients are equal to 0, and the remainders
are equal to a. (4.4)
4.1.3 Integer Division and
Remainder in Shen
The maths library
provides all three types of integer division with
corresponding remainder. As one frequently requires both
the integer quotient and the remainder, functions are
available that return both. This avoids unnecessary
recomputation of the quotient.
Note that all references to number in the signatures
below are references to integers. Also note that all
functions in this section are evaluated using integer
arithmetic only (provided, of course, that the platform
supports it)!
/rem :
number > number > (number * number)
Input: Two integers A and B.
Output: The pair consisting of the truncated integer
quotient and remainder, or an error message if B is 0 or
not an integer
(/rem
1000 33)
(@p 30 10) : (number * number)
(/rem 1000 33.3)
divisor must be an integer!
truncdiv :
number > number > number
Input: Two
integers A and B.
Output: The truncated integer quotient of A and B, or an
error message if B is 0 or not an integer.
Note: (truncdiv A B) is mathematically equivalent to
(trunc (/ A B)), but using this would involve
floatingpoint arithmetic, and might, therefore,
introduce inaccuracies.
(truncdiv
100 6)
16 : number
(truncdiv 10000000000000000000001 33) \\ correct
303030303030303030303 : number
(trunc (/ 10000000000000000000001 33))
303030303030302998528 : number \\ not correct
rem : number
> number > number
Input: Two integers A and B.
Output: The value of A – B × (truncdiv A B), or an
error message if B is 0 or not an integer.
Note: The given expression is not the basis for computing
rem !
(rem 100
6)
4 : number
/mod :
number > number > (number * number)
Input: Two
integers A and B.
Output: The pair consisting of the floored integer
quotient and remainder,
or an error message if B is 0 or not an integer
(/mod
1000 33)
(@p 31 23) : (number * number)
(/mod 1000 33)
(@p 30 10) : (number * number)
(/mod 1000 33.3)
divisor must be an integer!
div : number
> number > number
Input: Two integers A and B.
Output: The floored integer quotient of A and B, or an
error message if B is 0 or not an integer.
Note: (div A B) could be computed as (floor (/ A B)), but
this would involve floatingpoint arithmetic. (See note
under truncdiv)
(div 100
6)
17 : number
(div 100 6.2)
divisor must be an integer!
mod : number
> number > number
Input: Two
integers A and B.
Output: The value of A – B × (div A B), or an error
message if B is 0 or not an integer.
Note: The given expression is not the basis for computing
mod !
(mod 100
6)
2 : number
/% : number
> number > (number * number)
Input: Two
integers A and B.
Output: The pair consisting of the Euclidean integer
quotient and remainder,
or an error message if B is 0 or not an integer
(/% 1000
33)
(@p 30 10) : (number * number)
(/% 1000 33)
(@p 31 23) : (number * number)
(/% 1000 33.3)
divisor must be an integer!
diveucl :
number > number > number
Input: Two integers A and B.
Output: The Euclidean integer quotient of A and B, or an
error message if B is 0 or not an integer.
(diveucl
100 6)
16 : number
% : number
> number > number
Input: Two integers A
and B.
Output: The value of A – B × (diveucl A B), or an
error message if B is 0 or not an integer.
Note: The given expression is not the basis for computing
% !
(% 100
6)
4 : number
Note: In the above
functions, it is stated that both inputs, A and B, should
be integers, and if B is not, then an error will result.
If A is not an integer then the functions will produce
results, which are actually meaningful. Here is the
explanation for /rem. The behaviour of /mod and /% can be
explained similarly.
Example:
(/rem
101.45 3)
(@p 33 2.450000000000003) : (number * number)
We observe that 101.45 =
3 × (33)  2.45. In general, the quotient q is trunc
(A/B) and the remainder is A – B × q.
4.2 Basic Number Theory Functions
Closely related to the
computation of remainders is the notion of divisibility:
an integer N is divisible by an integer M, precisely when
the remainder on dividing N by M is 0. It is clearly
irrelevant which of the three remainders discussed above
one takes. For reasons of efficiency a special function
is provided
divisibleby?
: number > number > boolean
Input: Two integers A
= 0 and B > 0.
Output: true if B divides A, otherwise false.
Note: For the sake of efficiency, the inputs are not
tested for validity. Illegal inputs will cause the
program to crash!
(divisibleby?
10175 25)
true : boolean
(divisibleby? 101 13)
false : boolean
Also included is a simple
primality test (based on trial division), which makes use
of the preceding function.
prime? :
number > boolean
Input: A positive integer N.
Output: true if N is a prime number, otherwise false.
(prime?
16127)
true : boolean
(prime? 11111)
false : boolean
This simple
primalitytest function is reasonably fast. On platforms
supporting only 32bit integers, the largest prime number
that can be represented is the Mersenne prime 2^{31} – 1 =
2147483647
(time
(prime? 2147483647))
run time: 0.09 secs
true : boolean
For larger primes the run
time may become prohibitive
(time
(prime? 274876858367))
run time: 2.48 secs
true : boolean
(time (prime? 4398050705407))
run time: 11.81 secs
true : boolean
(time (prime? 333130634930268667))
run time: 55.68 secs
false : boolean
A closely related
function factorises a given number
factorise :
number > (list number)
Input: A
positive integer N.
Output: The list of prime factors of N.
(factorise
(factorial 13))
[2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 5 5 7 11 13] : (list
number)
(factorise 3888472005563)
[197807 19657909] : (list number)
Two polyadic functions
for the ‘greatest common divisor’ and the
‘least common multiple’ are available
gcd : number
> . . . > number
Input: Any
number of integer values A, B, C . . .
Output: The greatest common divisor of the given
arguments. (Always >= 0).
(gcd 123
78 888888)
3 : number
(gcd 2262 2726 1682 1414562)
58 : number
lcm : number
> . . . > number
Input: Any
number of integer values A, B, C . . .
Output: The least common multiple of the given arguments.
(Always >= 0).
(lcm
2262 2726 1682)
3083106 : number
(lcm 12 34)
204 : number
