Professional

Home
Learn Shen
Videos
Community Wiki
Community
OS Kernel
OS Library
Shen Professional

 

The Maths Library

(v.5 17-10-16)
W O Riha

1. Introduction

The Shen maths library contains all the functions found in the C - maths library (C-89 standard) and several other operations and functions which may be available in one form or another as built-in 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 SBCL-Shen. Numerical results agree with those produced by an early version of JS-Shen, 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

1. Introduction

2. Data Type number

3. Maths Library Functions

3.1 Predicates
3.2 Sign and Absolute Value
3.3 Rounding Functions
3.3.1 Rounding Up or Down
3.3.2 Rounding Towards Zero/Away from Zero
3.3.3 Rounding to the Nearest Integer
3.3.4 Rounding to a Specified Number of Decimal Places
3.4 Mathematical Constants
3.5 Degree Conversions
3.6 Trigonometric Functions
3.6.1 Computing sin x
3.6.2 Computing cos x
3.7 Inverse Trigonometric Functions
3.8 Exponential and Hyperbolic Functions
3.9 Logarithms
3.10 Inverse Hyperbolic Functions
3.11 Power Functions
3.12 Various Functions

4. Integer Functions

4.1 Integer Division and Remainder – A Short Tutorial
4.1.1 Mathematical Background
4.1.2 Integer Division and Remainder in Programming Languages
4.1.3 Integer Division and Remainder in Shen
4.2 Basic Number Theory Functions

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 built-in – 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 (floating-point) 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 (log2 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 non-negative 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 int-part (a). The number a – int-part (a), is called the fractional part of a, and is denoted by frac-part (a).

It is not difficult to see that

int-part(a) = floor(a) where a = 0 or a > 0

= ceiling(a) where a < 0

which means that int-part (a) rounds towards zero. [int-part is often called truncate (or trunc), because the value is obtained by ignoring (throwing away) the fractional part].

trunc : number --> number
int-part : number --> number

Input: A number X.
Output: The integer part of X. (see Table 2)


frac-part : number --> number
Input: A number X.
Output: The fractional part of X.


Table 2 (Note: fractional values have been rounded)

X (int-part X) (frac-part 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 ‘tie-break’. The main concern is to avoid or to minimise a bias (for a discussion, see Rounding).

Some libraries may provide a choice of different tie-breaks. The method that is used in the Shen maths-library 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 frac-part (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 0-ary 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/sqrt-pi) = 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 non-negative, 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 > <num-abs60> := <number : |number| < 60 >
<d-list> := [<number>]
<dm-list> := [<integer> <num60>]
<ms-list> := [0 <int–59> <num60>]
<m-list> := [0 <num-abs60>]
<s-list> := [0 0 <num-abs60>]
<dms-list> := [<integer> <int59> <num60>]
<angle> := <d-list> | <dm-list> | <ms-list> | <m-list> | <s-list> | <dms-list>

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! – x3/3! + x5/5! – x7/7! + … for |x| < inf (3.6.1)
cos x = 1 – x2/2! + x4/4! – x6/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.4821693935912637e-16 : 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 multi-valued 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 so-called principal branch. For the inverse trigonometric functions the principal branch corresponds to k = 0. The principal branch of any multi-valued function is usually denoted by an initial capital letter, for example, Arcsin, Arccos, Arctan. One can then express the multi-valued 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 x2 + 2x – 1 = 0, which has the positive solution sqrt (2) – 1 = 0.4142…

To obtain a 16-digit 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 two-argument 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 eX.

(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 loge (or simply log or ln). By definition

loge eX = elogeX = X

Sometimes, logarithms to other bases, especially 2 and 10 are needed. If the base is b (> 0) then

logb bX = blogbX = 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 base-2 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 base-10 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 logb 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 multi-valued for negative arguments and bases.


log : number --> number --> number
Input: Two positive numbers X and B.
Output: The value of the base-B 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 XN or an error message, if N is not an integer.
Note: The computation runs in log N - time. We have defined 00 = 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 XY

expt : number --> number --> number
Input: Two numbers X and Y.
Output: The value of XY or an error message, if XY is not defined or not a real number.
Note: We have defined that 00 = 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 non-negative 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 2E 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 2E , 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 ‘floating-point 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 non-negative integer N.
Output: The value of N!, or if N < 0, the value 1. For non-integral 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 : numbe
r

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 non-negative 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 non-negative (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 non-negative
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 re-computation 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!


trunc-div : 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: (trunc-div A B) is mathematically equivalent to (trunc (/ A B)), but using this would involve floating-point arithmetic, and might, therefore, introduce inaccuracies.

(trunc-div 100 -6)
-16 : number

(trunc-div 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 (trunc-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 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 floating-point arithmetic. (See note under trunc-div)

(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!


div-eucl : 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.

(div-eucl 100 -6)
-16 : number


% : number --> number --> number
Input: Two integers A and B.
Output: The value of A – B (div-eucl 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


divisible-by? : 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!

(divisible-by? 10175 25)
true : boolean

(divisible-by? 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 primality-test function is reasonably fast. On platforms supporting only 32-bit integers, the largest prime number that can be represented is the Mersenne prime 231 – 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