NUMBERS MODULE version 0.35


(C) Copyright Nick Craig-Wood 1992-3 & Gavin Wraith 2003-05

The module Numbers and this documentation NumModTxt are freeware.  
They remain copyright Nick Craig-Wood at all times.  You may use and 
distribute them freely provided that the following conditions are 
adhered to
 1) both files are always kept and distributed together
 2) the files are not altered in any way
 3) no profit is made from any distribution with any of the files in
If you wish to break any of these conditions get in touch with me at the 
address below first.

Under no circumstances shall the author be liable for any damage, loss 
of profits, time or data or any indirect or consequential loss arising 
out of the use of this software or inability to use this software or 
documentation.

The software and documentation was first published in BBC Acorn User 
Magazine, September 1992, and for more information see that issue.

INTRODUCTION

These routines started life in 1989 written in C on the university 
mainframe in response to a  Numbers Count puzzle in PCW.  The problem 
(or part of it) was to calculate the biggest prime you could find of the 
form N!+1.   However the routines did not get finished in time  for the 
puzzle deadline.   In the quest for speed I converted the routines to 
68000 code for my  Atari ST, and then a year later into ARM code.  With 
an ARM3 the routines run at the same speed as they did on the mainframe!

Numbers are held in application memory area in a heap, and maintained by 
OS_Heap.  Before the module can be used, it must be told about this heap 
with Num_HeapInit.  Numbers come in two pieces, the head and the tail.   
All numbers have a head, but need not have a tail.  The head of the 
number is a short block in the heap, which points to the tail and keeps 
the length and sign of the number, and various other things.  Numbers 
are passed by reference only, as pointers to the head of the number.  
These are referred to as NUMs.  For example NUM r0 indicates that r0 is 
a pointer to the head of a number.

Since NUMs are passed by reference only, if you pass them into 
procedures and then alter the values of the input parameters, you are 
actually altering the NUM.  So to avoid problems, if you write a 
procedure which takes NUMs as parameters for both input and output, it 
is best to define some temporary NUMs using Num_Init for the output, and 
then at the end Num_Swap the results into the correct place, and 
Num_Remove the temporary variables.  All the internal routines of the 
module work like this, so there is never any problem with passing the 
same NUM as input and output.

Some of the routines take either NUMs or scalars as arguments.  Scalars 
are normal signed or unsigned 32 bit numbers.  All routines expect NUMs 
to be tidy (except Num_Tidy).  (A num is tidy if its most significant 
digit is non-zero, and if it is zero, then it is positive, and of length 
1.  Num_Tidy accomplishes this.  It is unlikely that a user will need 
this function.)

The numbers module checks the type of any parameters passed to any of 
its SWIs quite carefully.  The NUM parameters are checked by looking for 
a special word in the NUM header block.  The heap pointers are checked 
in a similar way.  The pointers to memory are checked to see whether 
they point to logical memory above &8000 (the system workspace).  This 
checking takes very little time, and increases the reliability of the 
module enormously, especially when developing new programs.

Any error from the numbers module will be prefixed with Numbers: .  
For a full list see the ERRORS section.  For a short example program see 
the EXAMPLE section.

If you find any bugs, have any comments or queries or wish to use the 
module in commercial software please write to

Nick Craig-Wood
26 Wodeland Avenue
Guildford
Surrey
GU2 5JZ

e-mail: ncw@axis.demon.co.uk


NUMBERS SWI

On Entry:       r0-r3 are parameters
On Exit:        r0-r3 are returned as results or corrupted
                r4-r14 are preserved

Interrupts:     Interrupts are enabled
                Fast interrupts are enabled

Processor mode: Processor is in SVC mode

Re-entrancy:    SWIs are not re-entrant

Use:            See QUICK REFERENCE, and DEFINITIONS for details
                All SWIs that are capable of looping respond to 
ESCAPE
                See ERRORS section for details of errors returned

QUICK REFERENCE

NUMBERS HEAP:   area in application workspace where all the numbers 
                are kept
HEAD:           16 byte parameter block in the numbers heap pointing 
                to
TAIL:           variable length block containing the value of the 
                number
NUM:            meaning a 32 bit integer pointing to a HEAD
INT:            meaning a 32 bit integer
SCALAR:         meaning a 32 bit integer
FLAG:           is either 0 for false, or 1 for true
STRING:         a pointer to a 0 terminated ASCII string

a,b,c,d         represent pointers to NUMs
i,j,k           represent SCALARS
p,q             represent a pointer to memory
f               represents a FLAG
h               represents a HeapPointer as returned by Num_HeapInit 
                in r0
          
SWI             Params  Results Comments

Num_Author      -       -       Prints out some info
Num_HeapInit    p,i     h,a,b,c h<-HeapPointer, a<-0, b<-1, c<-2
Num_MakeSmallPrimes i   j       Makes primes up to i, # found to j

        (All the SWIs below need a valid NUM or HeapPointer in r0)

Num_Allocate    a,i     -       Allocates i words to TAIL of a
Num_Deallocate  a       -       Removes the TAIL of a
Num_Set         a,i     -       Sets a to signed i
Num_USet        a,i     -       Sets a to unsigned i
Num_Init        h       a       Makes a HEAD and TAIL and pointer a
Num_Remove      a       -       Removes the HEAD and TAIL of a
Num_Equals      a,i     f       Returns truth of a = i
Num_Swap        a,b     -       Swaps a and b.  This is quick.
Num_Move        a,b     -       Sets b to a.  Not so quick.
Num_Clear       a       -       Sets TAIL of a to all zeros
Num_Tidy        a       -       Reduces a to shortest length
Num_UCmp        a,b     i       Unsigned compare of a-b to i
Num_Cmp         a,b     i       Returns the sign of a-b to i
Num_ScalarCmp   a,i     j       Returns the sign of a-i to j
Num_ScalarAdd   a,i,c   -       c <- a + i
Num_ScalarSub   a,i,c   -       c <- a - i
Num_ScalarMul   a,i,c   -       c <- a * i
Num_ScalarDiv   a,i,c   j       c <- a / i, j <- a MOD i
Num_ScalarMod   a,i     j       j <- a MOD i
Num_SmallFactorN a,i    k       k=smallest factor of a or 0, try i
Num_SmallFactor a       k       k=smallest factor of a or 0, try all
Num_Add         a,b,c   -       c <- a + b
Num_Sub         a,b,c   -       c <- a - b
Num_Mul         a,b,c   -       c <- a * b
Num_Div         a,b,c,d -       c <- a / b, d <- a MOD b
Num_Mod         a,b,c   -       c <- a MOD b
Num_Dump        a       -       Prints a in hex, and other info
Num_ToString    a       p,q,b   Makes STRING of a to p, end=q
                                After use do SYS Num_Remove,b
Num_Print       a       -       Prints a in the current base
Num_FromString  a,p     f       Converts STRING at p to a, f=ok
Num_Input       a       f       Gets a input in current base to a, f=ok
Num_Info        h       -       Prints memory usage info
Num_RndScalar   h       i       Makes a random number 0-&FFFFFFFF
Num_SetSeed     h,i     -       Sets seed of rnd generator
Num_Rnd         a,b     -       Makes random number to b, 0 <= b < a
Num_Gcd         a,b,c   -       c <- greatest common divisor a,b
Num_Sqr         a,b     -       b <- square root of a
Num_Pow         a,b,c   -       c <- a ^ b (a to the power b)
Num_PowMod      a,b,c,d -       d <- (a^b) MOD c
Num_Factorial   a,b     -       b <- a! = a*(a-1)*(a-1)*..*3*2*1
Num_Inv         a,b,c,d -       Finds c,d: a*c MOD b=d & d=GCD(a,b)
Num_FermatTest  a,i     j       Computes truth of i^(a-1) MOD a = 1 to j
Num_ProbablyPrime a     j       j <- 0 if not prime, 1 probablyprime
Num_Base        h,i     j       Sets base to i, old base in j
Num_ToMem       a,p,i,f -       Converts Nums to Memory
Num_FromMem     a,p,i,f -       Converts Memory to Nums


DEFINITIONS

Num_HeapInit

This expects r0 as a pointer to a block of memory to be used as a heap, 
and r1 as the length of this memory in bytes. This returns a pointer to 
the heap in r0 and some NUMs which are useful constants.  NUM r1 <- 
zero, NUM r2 <- one, NUM r3 <- two.  The HeapPointer returned in r0 is 
required by some of the other routines.

Num_Allocate

This updates NUM r0 to point to a new area of memory of length r1 (in 
32-bit words).  This memory is not initialised.

Num_Deallocate

This releases the memory used by a NUM (tail), and sets its pointers to 
NULL.

Num_Set

This sets the NUM supplied in r0 to the signed scalar supplied in r1.

Num_Uset

This sets the NUM supplied in r0 to the unsigned scalar supplied in r1.

Num_Init

Expects r0=HeapPointer.  This returns the address of a new NUM in r0.   
For creating new or LOCAL variables.

Num_Remove

This removes the allocation for NUM r0 and removes NUM r0.  For removing 
LOCAL variables.

Num_Equals

This function returns TRUE in r0 if NUM r0 is equal to the SCALAR 
supplied in r1

Num_Swap

This swaps the two NUMs in r0 and r1.  It is a fast way of getting 
information out of temporary variables, before destroying them

Num_Move

This moves NUM r0 to NUM r1.  It does this by actually copying the NUM.  
It is slow compared to Num_Swap.

Num_Clear

This clears the tail of NUM r0 to all zeros.  This might be useful for 
security.

Num_Tidy

This makes the NUM r0 use as few digits as possible by removing all the 
leading zeros.  All numbers should be tidied.

Num_UCmp

This returns an unsigned comparison between NUM r0 and NUM r1 in r0. This 
is SGN(ABS(NUM r0) - ABS(NUM r1))

Num_Cmp

This function returns a signed comparison between NUM r0 and NUM r1. It 
returns SGN(NUM r0 - NUM r1)

Num_ScalarCmp

This function returns a signed comparison between NUM r0 and SCALAR r1 
It returns SGN(NUM r0 - NUM r1)

Num_ScalarAdd

This adds NUM r0 to (signed int) r1 into NUM r2.

Num_ScalarSub

This subtracts NUM r0 - (signed int) r1 into NUM r2.

Num_ScalarMul

This multiplies NUM r0 by (unsigned int) r1 into NUM r2.

Num_ScalarDiv

This divides NUM r0 by (unsigned int) r1 into NUM r2, remainder 
(unsigned int) r0.

Num_ScalarMod

This divides NUM r0 by (unsigned int) r1 returning the (unsigned int) 
remainder in r0.

Num_Mul

NUM r0 * NUM r1 -> NUM r2

Num_Div

This divides NUM r0 by NUM r1 to give a quotient NUM r2 remainder NUM 
r3.  It does this using a base 2^32 method for speed.  For divide and 
correct algorithm see Knuth Seminumerical Algorithms section 4.3.1.

Num_Mod

This divides NUM r0 by NUM r1 to remainder NUM r2.

Num_MakeSmallPrimes

This makes all the primes up to r0 (using a sieve method), and stores 
them in the RMA pointed to by SmallPrimes.  Thus the SmallPrimes array 
is common to all the users of the numbers module.  NSmallPrimes is set 
to the number of primes found.  It deallocates any old SmallPrimes array 
so this may be called more than once.  It returns the number of 
smallprimes found in r0, and a pointer to the array in r1.  If there are 
already more than enough SmallPrimes in the RMA then this will not 
recalculate.

Num_Dump

This dumps the number supplied in r0 in Hexadecimal, for debugging, 
along with some other information.

Num_ToString

This converts the number pointed to by r0 into a string in the current 
base pointed to by r0.  It returns the position of the null at the end 
of the string in r1 and the pointer to the num in r2.  When finished, 
this memory should be released with SYS Num_Remove,r2

Num_Print

This prints out the number supplied in r0 (in the current base).  It 
prints no spaces or newlines.

Num_Info

Expects r0 = HeapPointer.  This prints information on the numbers 
module, and current heap.

Num_FromString

This converts a ASCII number in the current base pointed to by r1 into 
the NUM pointed to by r0, until a character which is not valid in the 
current base is reached.  If this is a control char then true is 
returned otherwise false is returned.  It deals with leading +,-,
 .  r1 is updated to point to the last unconverted character.

Num_Input

Inputs a number in the current base to NUM r0.  Returns r0 flagging ok 
conversion (1 => good conversion, 0 => bad conversion).

Num_RndScalar

Expects r0 = HeapPointer.  Produces a random number into (unsigned int) 
r0 from Seed, and updates Seed, in the range 0-&FFFFFFFF.  It works 
using the algorithm x(n+1)=(1664525 * x(n) + 907633393) MOD 2^32 It is 
known that the least significant bits are not as random as the most 
significant bits, so two consecutive random numbers are generated, and 
combined with EOR after having been rotated by 16 bits with respect to 
each other.  This halves the period.  Reference: Knuth, Seminumerical 
Algorithms 3.3 p102 p84 This is Line 26, with the best spectral co-
efficients for a MOD 2^32 generator, listed in Knuth.  This form is 
especially easy to calculate with the ARMs mul instruction. The A term 
is calculated as the nearest prime to (1/2-1/3*SQR(3))*2^32.

Num_SetSeed

Expects r0 = HeapPointer.  This sets the internal 32-bit random number 
generator to the seed supplied in (unsigned int) r1.

Num_Rnd

This generates a random number 0 <= rnd < NUM r0 to NUM r1.

Num_Gcd

This finds the GCD(NUM r0, NUM r1) -> r2 using Euclids algorithm.  GCD 
is the greatest common divisor, that is the largest number which divides 
exactly into NUM r0 and NUM r1.

Num_Sqr

This takes the square root of NUM r0 to NUM r1.  It uses a second order 
Newton-Raphson convergence.  This doubles the number of significant 
figures on each iteration.  Square root of a negative number is returned 
as 0.  It returns the number of iterations in r0.

Num_Pow

This finds NUM r0 to the power of NUM r1 to NUM r2.

Num_PowMod

This finds NUM r0 to the power of NUM r1 MOD NUM r2 to NUM r3.

Num_Factorial

This finds NUM r0 factorial to NUM r1 if NUM r0 <=1 then NUM r1=1.

Num_SmallFactorN

This sees whether NUM r0 has any small factors.  It requires 
Num_MakeSmallPrimes to have been called.  It tests r1 small primes.  It 
returns either the factor found, or 0 to show none found in r0.

Num_SmallFactor

This sees whether NUM r0 has any small factors.  It requires 
Num_MakeSmallPrimes to have been called.  It tests all the small primes 
that have been made.  It returns either the factor found, or 0 to show 
none found in r0.

Num_Inv

This finds NUM r2 such that r0*r2 MOD r1=r3 AND r3=GCD(r0,r1) so if 
GCD(r0,r1)=1 (for instance if r1 is prime) then this will compute the 
inverse MOD r1. This is adapted from Knuth: Seminumerical Algorithms 
4.5.2 pp325

Num_FermatTest

This finds whether NUM r0 is a prime with respect to SCALAR r1, using 
the fermat test.  That is r0 is not a prime if r1^(r0-1) MOD r0 <> 1 IF 
r1^(r0-1) MOD r0=1 THEN r0 may be a prime. r1 must be prime for the test 
to work. IE this test returns false if r0 is not prime, or true if r0 
might be prime.  This test is less reliable than Num_ProbablyPrime and 
takes about the same time to execute.

Num_ProbablyPrime

This tests whether r0 is prime or not.  The function returns (in r0) 
false if the number is not prime, and true if the number is probably 
prime. The algorithm will return true with r0 composite with a 
probability of less than 0.25.  This algorithm is as in Knuth - 
Seminumerical algorithms 4.5.4P

Num_Base

This sets the current base to that supplied in r1.  r0 should be 
supplied as the HeapPointer or a valid NUM.  That base will be used for 
Num_ToString, Num_FromString, Num_Input and Num_Print.  r1 can be from 2 
to 32.  If it is outside that range then this SWI will not change the 
base.  The old base will be returned in r0.  If you want to know the 
current base then pass a 0 to Num_Base in r1.

Num_ToMem 

This converts variable length NUMs into fixed length integers in memory.  
The numbers in memory are assumed to unsigned.  If the NUM is too long 
for the fixed space the most significant digits will be truncated.
r0 is pointer to NUM to be converted
r1 is the pointer to the memory space for the fixed integer
r2 is the length in words of the memory integer
r3 = 0 means lsw first, 1 = msw first

Num_FromMem

This converts fixed length integers in memory into variable length NUMs.  
The numbers in memory are assumed to unsigned.
r0 is pointer to NUM to be converted
r1 is the pointer to the memory space for the fixed integer
r2 is the length in words of the fixed space
r3 = 0 means lsw first, 1 = msw first


EXAMPLE

In the spirit of an illustration is worth 1000 words, here is an example 
BASIC program using the numbers module.  It tries to find any primes of 
the form 111...111 (these are called rep-units).

REM Check to see we have the numbers module
*RMEnsure Numbers 0.0 RMLoad Numbers
*RMEnsure Numbers 0.0 Error 1 Numbers module not found

REM put aside some BASIC memory for the Numbers Heap
HeapSize=64*1024
DIM Numbers HeapSize

REM Initialise the Heap
SYS Num_HeapInit,Numbers,HeapSize TO hp%,zero%,one%,two%

REM Make some small primes for finding small factors
SYS Num_MakeSmallPrimes,5000 TO a%
PRINT ;a%; small primes found

REM This is the number of times we will run Num_ProbablyPrime to be sure
REM a repunit is prime
timestotest=20

REM Make the NUM repunit% and start it off at 1
SYS Num_Init,hp% TO repunit%
SYS Num_Set,repunit%,1

FOR n%=2 TO 100000
  REM repunit=10*repunit+1, ie add 1 digit to the repunit
  SYS Num_ScalarMul,repunit%,10,repunit%
  SYS Num_ScalarAdd,repunit%,1,repunit%

  PRINT .;
  REM try to find a small factor of repunit%
  SYS Num_SmallFactor,repunit% TO composite
  IF composite=0 THEN
   composite=timestotest
   WHILE composite>0
    PRINT |;
    REM test the primality of repunit% timestotest times
    SYS Num_ProbablyPrime,repunit% TO pprime
    IF pprime THEN
     composite-=1
    ELSE
     REM if the number is composite but passes a test, this is unusual
     IF composite<>timestotest THEN PRINT Unusual prime:;FNprint(repunit%)
     composite=-1
    ENDIF
   ENDWHILE
  ENDIF
  REM print out a prime if we have found one
  IF composite=0 THEN PRINT R(;n%;) = ;FNprint(repunit%); is prime
NEXT n%
END

REM This prints NUM a% returning a null string

DEF FNprint(a%): SYS Num_Print,a%: =

ERRORS

The following errors are unique to the Numbers module.  However the 
numbers Module will return any OS errors (such as OS_Heap errors) 
unchanged, except for the addition of a prefix of Numbers: 

Numbers module has become corrupted
Returned on initialisation if the Numbers code has been changed (say by 
a virus or by a disk error).

Numbers: Unknown Operation
Returned when an unimplemented numbers SWI is called.

Numbers: Escape
Returned when the user pressed Escape in a Looping SWI.

Numbers: Invalid Heap Pointer
Returned when a bad heap pointer is passed to a SWI or possibly if a NUM 
header has become corrupt.

Numbers: Invalid Number Pointer
Returned if a bad NUM pointer was passed to a SWI.

Numbers: Invalid Pointer
Returned if a pointer was passed to a SWI which was not in logical 
memory above &8000 (system workspace).

Numbers: Divide by 0 in Num_ScalarDiv

Numbers: Divide by 0 in Num_Div

Numbers: N > NSmallPrimes in Num_SmallFactorN
Returned when you ask for more SmallPrimes than have been calculated by 
Num_MakeSmallPrimes.

Numbers: Internal: Bad Validation Character
Numbers: Internal: Addback Failed in Num_Div
Numbers: Internal: Trial quotient wrong in Num_Div
These errors should never occur.  If they do, I would be interested to 
know!

HISTORY

Version 0.08, 01 Aug 92

First release, published in Acorn User magazine, September 1992.

Version 0.10, 15 Oct 92

SWI argument validation was added.  This means all arguments to the SWIs 
are checked as rigorously as possible to catch mistakes (for example 
passing an integer rather that a pointer to a num) when writing programs 
with the numbers module.  It makes debugging a lot easier, takes very 
little time and ensures the numbers module wont crash if you pass the 
wrong thing to it.  Special words were put in the various blocks in the 
heap to enable checking.  The errors Bad Validation Character, 
Invalid Heap Pointer, Invalid Number Pointer and Invalid Pointer 
were added for this purpose.

Through popular request the following SWIs were added:

Num_Base: To allow number IO to be in any base 2-36.  Required a major 
re-write of the conversion routines.

Num_ToMem: To allow numbers to be put into the heap in a user friendly 
way for putting to a file or other manipulation.  Also capability for 
LSW or MSW first added for compatibility with other programs.

Num_FromMem: To allow numbers to be got from the heap in a user friendly 
way.

Version 0.11, 07 Sep 93

A bug in Num_Div normalisation routines removed.  It only showed itself 
when the divisor was longer than 1 word and the top word was &FFFFFFFF.  
A new error Internal: Trial quotient wrong in Num_Div and checking 
(very quick) added to Num_Div to catch any further bugs in Num_Div 
(should there be any!). 

Version 0.12, 23 Sep 93

Conditions of use changed slightly.
Bug in documentation corrected.

Version 0.12x, 05 Nov 2003

A 32-bit compatible version, tweaked by Gavin Wraith.

Version 0.32, 06 Nov 2003

A bug in Num_Sqr in v. 0.12x corrected.  It showed itself when trying to
find the root of a number that was one less than a perfect square.
(Gavin Wraith)

Version 0.33, 06 Jun 2005

Another bug in Num_Sqr corrected.  It only showed itself when the input
number was 2^32 or greater.     (Gavin Wraith)
Some more documentation bugs found and corrected.     (Gavin Wraith)

Version 0.34, 06 Jun 2005
Old bug replaced in 0.33. Corrected. (Gavin Wraith)

Version 0.35, 25 July 2005
Wrong algorithm being used in 0.34.  Corrected. (Gavin Wraith)

REFERENCES

Books I have found useful whilst writing this module

SemiNumerical Algorithms (2nd Edition) D.E.Knuth, Addison-Wesley
Elementary Number Theory D.M.Burton, Allyn and Bacon
Cryptography: an introduction to computer security J.Seberry, Prentice 
Hall

