home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / manual1.doc < prev    next >
Text File  |  1989-04-19  |  45KB  |  1,651 lines

  1.     M.I.R.A.C.L. - Multiprecision Integer and Rational Arithmetic C Library    _______________________________________________________________________
  2.  
  3.                                  M. Scott
  4.  
  5.  
  6. ABSTRACT:
  7.  
  8.  
  9.  
  10.       A portable C library is described  which  implements  multi-precision
  11.  
  12.       integer and rational data-types, and provides the routines to perform
  13.  
  14.       basic  arithmetic  on them.  Certain useful number-theoretic routines
  15.  
  16.       are also included in the package.
  17.  
  18.  
  19.  
  20. 1. INTRODUCTION:
  21.  
  22.  
  23.  
  24.       Conventional computer  arithmetic  facilities  as  provided  by  most
  25.  
  26.       computer language compilers usually provide one or two floating-point
  27.  
  28.       data  types  (e.g.  single and double precision) to represent all the
  29.  
  30.       real numbers,  together with one or more integer types  to  represent
  31.  
  32.       whole  numbers.  These built-in data-types are closely related to the
  33.  
  34.       underlying computer architecture,  which is sensibly designed to work
  35.  
  36.       quickly with large amounts of small numbers,  rather than slowly with
  37.  
  38.       small amounts of large numbers (given  a  fixed  memory  allocation).
  39.  
  40.       Floating-point allows a relatively small binary number (e.g. 32 bits)
  41.  
  42.       to  represent  real numbers to an adequate precision (e.g.  7 decimal
  43.  
  44.       places) over a large dynamic range.  Integer types allow small  whole
  45.  
  46.       numbers to be represented directly by their binary equivalent,  or in
  47.  
  48.       2's complement  form  if  negative.  Nevertheless  this  conventional
  49.  
  50.       approach to computer arithmetic has several disadvantages.
  51.  
  52.  
  53.  
  54.       (1)   Floating-point  and  Integer  data-types  are  representationly
  55.  
  56.       incompatible.  Note that the set of integers, although infinite, is a
  57.  
  58.  
  59.  
  60.                                       1                                       
  61.  
  62.  
  63.  
  64.  
  65.  
  66.  
  67.       subset of the rationals (i.e.  fractions),  which is in turn a subset
  68.  
  69.       of the reals.  Thus every integer has  an  equivalent  floating-point
  70.  
  71.       representation.  Unfortunately  these  two  representations  will  in
  72.  
  73.       general be different.  For example a small positive whole number will
  74.  
  75.       be  represented  by  its  binary  equivalent  as  an integer,  and as
  76.  
  77.       separated mantissa and exponent as a floating-point. This implies the
  78.  
  79.       need for conversion routines, to convert from one form to the other.
  80.  
  81.  
  82.  
  83.       (2) Most rational numbers can not be expressed  exactly  (e.g.  1/3).
  84.  
  85.       Indeed  the  floating-point  system  can  only  express exactly those
  86.  
  87.       rationals whose denominators are multiples  of  the  factors  of  the
  88.  
  89.       underlying  radix.  For  example our familiar decimal system can only
  90.  
  91.       represent exactly  those  rational  numbers  whose  denominators  are
  92.  
  93.       multiples of 2 and 5; 1/20 is 0.05 exactly, 1/21 is .0476190476190...
  94.  
  95.  
  96.  
  97.       (3) Rounding in floating-point is  base-dependant  and  a  source  of
  98.  
  99.       obscure errors.
  100.  
  101.  
  102.  
  103.       (4)  The  fact that the size of integer and floating-point data types
  104.  
  105.       are dictated by the computer architecture,  defeats  the  efforts  of
  106.  
  107.       language designers to keep their languages portable.
  108.  
  109.  
  110.  
  111.       (5)  Numbers  can  only  be  represented to a fixed machine-dependent
  112.  
  113.       precision. In many applications this can be a crippling disadvantage,
  114.  
  115.       for example in the new and growing field of Public-Key cryptography.
  116.  
  117.  
  118.  
  119.       (6) Base-dependent phenomena cannot easily be studied. For example it
  120.  
  121.       would be difficult to access a particular digit of a decimal  number,
  122.  
  123.       as  represented  by  a  traditional  integer  data-type.
  124.  
  125.  
  126.                                       2                                       
  127.  
  128.  
  129.  
  130.  
  131.  
  132.  
  133.  
  134.  
  135.       Herein is described a set of standard  C  routines  which  manipulate
  136.  
  137.       multi-precision  rational  numbers  directly,   with  multi-precision
  138.  
  139.       integers as a compatible subset. Approximate real arithmetic can also
  140.  
  141.       be performed.
  142.  
  143.  
  144.  
  145.  
  146.  
  147. 2. NUMBER REPRESENTATION
  148.  
  149.  
  150.  
  151.       The two new data-types are called "big" and "flash".  The  former  is
  152.  
  153.       used  to  store  multiprecision  integers,   and  the  latter  stores
  154.  
  155.       multiprecision fractions as numerator and denominator  in  "floating-
  156.  
  157.       slash" form.  Both take the form of a fixed length array of integers,
  158.  
  159.       with sign and length information encoded in the  first  element,  and
  160.  
  161.       the  data  itself  in  subsequent  elements.  Both  new  types can be
  162.  
  163.       introduced into the syntax of the 'C' language by the 'C'  statements
  164.  
  165.                             typedef int* big;
  166.                             typedef int* flash;
  167.  
  168.       Now  'big' and 'flash' variables can be declared just like any built-
  169.  
  170.       in data type, e.g.
  171.  
  172.  
  173.  
  174.                              big x,y[10],z[10][10];
  175.  
  176.  
  177.  
  178.       Note that the user of these data-types is  not  concerned  with  this
  179.  
  180.       internal representation; the library routines allow "big" and "flash"
  181.  
  182.       numbers to be manipulated directly.
  183.  
  184.  
  185.  
  186.       The structure of "big" and "flash" numbers is illustrated in figure (1)
  187.  
  188.  
  189.  
  190.  
  191.  
  192.                                       3                                       
  193.  
  194.  
  195.  
  196.  
  197.  
  198.  
  199.         x[0]    x[1]   x[2]                x[n]                          x[max]
  200.       ┌─────┐ ┌─────┐ ┌─────┐             ┌─────┐ ┌─────┐                ┌─────┐
  201.       │s 0 n│ │ LSW │ │     │ ..........  │ MSW │ │  0  │ .............. │  0  │
  202.       └─────┘ └─────┘ └─────┘             └─────┘ └─────┘                └─────┘
  203.  
  204.  
  205.         x[0]   x[1]    x[2]        x[n]   x[n+1]   x[n+2]     x[n+d]     x[max]
  206.       ┌─────┐ ┌─────┐ ┌─────┐     ┌─────┐ ┌─────┐ ┌─────┐     ┌─────┐    ┌─────┐
  207.       │s d n│ │ LSW │ │     │ ... │ MSW │ │ LSW │ │     │ ... │ MSW │ .. │  0  │
  208.       └─────┘ └─────┘ └─────┘     └─────┘ └─────┘ └─────┘     └─────┘    └─────┘
  209.               <------- numerator ------> / <------ denominator ----->
  210.  
  211.  
  212.        Figure (1):  Structure of "big" and "flash" data-types,  where s  is
  213.                    the  sign of the number,  n and d are the lengths of the
  214.                    numerator and denominator respectively,  and LSW and MSW
  215.                    mean  'Least  significant  word'  and  'Most significant
  216.                    word'
  217.  
  218.  
  219.  
  220.       These   structures   combine   ease   of  use  with  representational
  221.  
  222.       efficiency.  A denominator of length zero (d=0),  implies  an  actual
  223.  
  224.       denominator  of one;  and similarily a numerator of length zero (n=0)
  225.  
  226.       implies  a  numerator of one.  Zero itself is uniquely defined as the
  227.  
  228.       number whose first element is zero (i.e. n=d=0).
  229.  
  230.  
  231.  
  232.       Note that the slash in the  'flash'  data-type  is  not  in  a fixed
  233.  
  234.       position,    and  may  "float"   depending  on the relative  size  of
  235.  
  236.       numerator and denominator.
  237.  
  238.  
  239.  
  240.       A "flash" number is manipulated by  splitting  it  up  into  separate
  241.  
  242.       "big"  numerator  and  denominator  components.  A  "big"  number  is
  243.  
  244.       manipulated by extracting and operating  on  each  of  its  component
  245.  
  246.       integer elements.  To avoid possible confusion with the internal sign
  247.  
  248.       bit,  the  numbers  in each element are limited to a somewhat smaller
  249.  
  250.       range than that of the full word-length, e.g. 0 to 16383  (= 2^14 -1)
  251.  
  252.       on a 16-bit micro-computer.
  253.  
  254.  
  255.  
  256.  
  257.  
  258.                                       4                                       
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265.       When the system is initialised the user specifies the fixed number of
  266.  
  267.       words (or bytes) to be assigned to all "big"  or  "flash"  variables,
  268.  
  269.       and the number base to be used. Any base can be used, up to a maximum
  270.  
  271.       MAXBASE which is dependant on the wordlength of the computer used. If
  272.  
  273.       requested  to  use  a  small  base  b,  the system will,  for optimal
  274.  
  275.       efficiency,  actually use base bⁿ ,  where n is the  largest  integer
  276.  
  277.       such that bⁿ fits in a single computer word.
  278.  
  279.  
  280.  
  281.       The   encoding  of  the  sign  and  numerator  and  denominator  size
  282.  
  283.       information into the first word is possible,  as the C  language  has
  284.  
  285.       standard  constructs  for bit manipulation.  On a 16-bit computer the
  286.  
  287.       maximum number of words in each 'big' or 'flash' variable is  limited
  288.  
  289.       to  128,  the  equivalent  of  about 500 decimal digits.  On a 32-bit
  290.  
  291.       computer there is, effectively, no limit.
  292.  
  293.  
  294.  
  295.       Normally in flash arithmetic full precision is always  used.  However
  296.  
  297.       the  internal  global variable 'workprec' can be used to  temporarily
  298.  
  299.       reduce the working precision.  This is  useful  when  using  Newton's
  300.  
  301.       method for the calculation of roots and certain elementary functions
  302.  
  303.       (Scott 1988).
  304.  
  305.  
  306.  
  307.  
  308.  
  309. 3. AN EXAMPLE
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.                                       5                                       
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332. /*   Program to calculate factorials.  */
  333.  
  334. #include <stdio.h>
  335. #include "miracl.h"   /* include big number system */
  336.  
  337. main()
  338. { /* calculate factorial of number */
  339.     big nf;           /* declare "big" variable nf */
  340.     int n;
  341.     mirsys(5000,10);  /* initialise system to base 10, 5000 digits per "big" */
  342.     nf=mirvar(1);     /* initialise "big" variable nf=1 */
  343.     printf("factorial program\n");
  344.     printf("input number n= \n");
  345.     scanf("%d",&n);
  346.     getchar();
  347.     while (n>1) premult(nf,n--,nf);   /* nf=n!=n*(n-1)*(n-2)*....3*2*1  */
  348.     printf("n!= \n");
  349.     otnum(nf,stdout); /* output result */
  350. }
  351.  
  352.       This program can be used to quickly calculate and print out 1000!  (a
  353.  
  354.       2568 digit number) in about 12 seconds on a VAX11/780,  a task  first
  355.  
  356.       performed  "by  H.S.  Uhler using a desk calculator and much patience
  357.  
  358.       over a period of several years"  (Knuth  1973).  Many  other  example
  359.  
  360.       programs are described in Appendix IV.
  361.  
  362.  
  363.  
  364.  
  365.  
  366. 4. THE USER INTERFACE
  367.  
  368.  
  369.  
  370.       Any  program  that  uses  the  MIRACL  system  must have an '#include
  371.  
  372.       "miracl.h"' statement.  This tells the compiler to  include  the  'C'
  373.  
  374.       header  file  "miracl.h"  with  the  main  program source file before
  375.  
  376.       proceeding with the compilation.  This file  (Appendix  II)  contains
  377.  
  378.       declarations  of  all  the  MIRACL routines available to the user and
  379.  
  380.       permits  access  to  some global variables. The small sub-header file
  381.  
  382.       'mirdef.h' contains hardware-specific details,  other user-tailorable
  383.  
  384.       details, and some useful general purpose definitions.
  385.  
  386.  
  387.  
  388.  
  389.  
  390.                                       6                                       
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.       In the main program the MIRACL system must be initialized by  a  call
  398.  
  399.       to  the  routine  mirsys,  which sets the number base and the maximum                        mirsys                                             
  400.  
  401.       size of the 'big' and 'flash' variables.  It also initializes the
  402.  
  403.       random number system,  and creates several workspace 'big' variables
  404.  
  405.       its own  use.  Also  initialized  is  the  relatively sophisticated
  406.  
  407.       error tracing system which is integrated with the MIRACL package.
  408.  
  409.       Whenever an  error  is  detected  the  sequence  of  routine calls
  410.  
  411.       down to the routine which generated the error is reported, as well as
  412.  
  413.       the  error itself. A typical error message might be
  414.  
  415.                     MIRACL error from routine powmod
  416.                                  called from prime
  417.                                  called from your program
  418.                     Raising integer to a negative power
  419.  
  420.       Such an error report facilitates debugging, and certainly assisted me
  421.  
  422.       during  the  development  of  these routines.  An associated variable
  423.  
  424.       TRACER,  initialised to OFF,  if set by the user to ON,  will cause a
  425.  
  426.       trace  of  the  program's  progress through the MIRACL routines to be
  427.  
  428.       output to the screen.
  429.  
  430.  
  431.  
  432.       Every 'big'  or  'flash'  variable  in  the  users  program  must  be
  433.  
  434.       initialized  by  a call to the routine mirvar,  which also allows the                                             mirvar                        
  435.  
  436.       variable to be given an initial integer value.
  437.  
  438.  
  439.  
  440.       The full set of arithmetic and number-theoretic routines declared  in
  441.  
  442.       "miracl.h"  may  be  used  on  these  variables.  Full flexibility is
  443.  
  444.       (almost always) allowed in parameter usage with these  routines.  For
  445.  
  446.       example the call multiply(x,y,z),  multiplies the 'big' variable x by                       multiply(x,y,z)                                     
  447.  
  448.       the 'big' variable y to give the result as 'big' variable z.  Equally
  449.  
  450.       valid would be multiply(x,y,x),  multiply(y,y,x), or multiply(x,x,x).                     multiply(x,y,x)   multiply(y,y,x)     multiply(x,x,x) 
  451.  
  452.       This last simply squares x.  Note that the first  parameters  are  by
  453.  
  454.  
  455.  
  456.                                       7                                       
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  
  463.       convention  always the inputs to the routines.  Routines are provided
  464.  
  465.       not only to allow arithmetic on 'big' and 'flash' numbers,  but  also
  466.  
  467.       to  allow  these  variables  to  perform arithmetic with the built-in
  468.  
  469.       integer and double precision data-types.
  470.  
  471.  
  472.  
  473.       Conversion routines are also provided to convert  from  one  type  to
  474.  
  475.       another.  For  details of each routine see the relevant documentation
  476.  
  477.       in Appendix III and the example programs of Appendix IV.
  478.  
  479.  
  480.  
  481.       Input and output is handled by the routines innum, otnum,  cinnum and                                                  innum  otnum   cinnum    
  482.  
  483.       cotnum. The first two use the fixed number base specified by the user      cotnum                                                               
  484.  
  485.       in  the  initial call of mirsys.  The latter pair work in conjunction                               mirsys                                      
  486.  
  487.       with the global variable IOBASE which can be assigned dynamically  by
  488.  
  489.       the  user.  A  simple  rule  is that if the program is CPU bound,  or
  490.  
  491.       involves changes of base,  then set the base initially to MAXBASE and
  492.  
  493.       use cinnum and cotnum. If on the other hand the program is I/O bound,          cinnum     cotnum                                                
  494.  
  495.       or needs access to individual digits of numbers (using getdig, putdig                                                             getdig  putdig
  496.  
  497.       and numdig), use innum and otnum.          numdig       innum     otnum 
  498.  
  499.  
  500.  
  501.       Numbers to bases up to 256 can be represented.  Numbers up to base 60
  502.  
  503.       use as many of the symbols 0-9, A-Z, a-x as necessary. If the base is
  504.  
  505.       greater than 60,  the symbols used are the ascii codes 0-255.  A base
  506.  
  507.       of  256 is useful when it is necessary to interpret a line of text as
  508.  
  509.       a large integer,  as is the case  for  the  Public  Key  Cryptography
  510.  
  511.       programs described in Appendix IV.
  512.  
  513.  
  514.  
  515.       Rational numbers may be input using either a radix point (e.g 0.3333)
  516.  
  517.       or  as  a fraction (e.g.  1/3).  Either form can be used on output by
  518.  
  519.       setting the global variable POINT=ON or =OFF.  Wraparound  on  output
  520.  
  521.  
  522.                                       8                                       
  523.  
  524.  
  525.  
  526.  
  527.  
  528.  
  529.       can  similarily  be selected by setting WRAP=ON or OFF.  Input can be
  530.  
  531.       taken from a string instead of from a file or I/O device,  by copying
  532.  
  533.       the  string  to  the  global array IBUFF[] immediately before calling
  534.  
  535.       innum or cinnum.  This method can be used to  set  'big'  or  'flash'      innum    cinnum                                                      
  536.  
  537.       numbers  to  large  constant values.  Similarly by setting the global
  538.  
  539.       variable STROUT to TRUE, output can be redirected to the global array
  540.  
  541.       OBUFF[],  instead  of  going  directly  to  a  file  or  I/O  device.
  542.  
  543.       Subsequent formatting can then take place before actual output.
  544.  
  545.  
  546.  
  547.  
  548.  
  549.       5. THE IMPLEMENTATION
  550.  
  551.  
  552.  
  553.       No  great  originality  is claimed for the routines used to implement
  554.  
  555.       arithmetic on the 'big' data-type.  The algorithms used are  faithful
  556.  
  557.       renditions  of  those described by Knuth (1981 $4.3.1).  However some
  558.  
  559.       effort was made to optimise the  implementation  for  speed.  At  the
  560.  
  561.       heart  of  the  time-consuming multiply and divide routines there is,
  562.  
  563.       typically,  a need need  to  multiply  together  a  digit  from  each
  564.  
  565.       operand,  add  in  a  'carry'  from  a  previous operation,  and then
  566.  
  567.       separate the total into a digit of the result,  and a 'carry' for the
  568.  
  569.       next   operation.   To   illustrate   this   consider  this  base  10
  570.  
  571.       multiplication:
  572.  
  573.                                 8723536221
  574.                                x         9
  575.                                -----------
  576.  
  577.  
  578.       To correctly process the column with  the  '5'  in  it,  we  multiply
  579.  
  580.       5x9=45,  add  in  the 'carry' from the previous column (say a 3),  to
  581.  
  582.       give 48,  keep the '8' as the result for this column,  and carry  the
  583.  
  584.       '4' to the next column.
  585.  
  586.  
  587.  
  588.                                       9                                       
  589.  
  590.  
  591.  
  592.  
  593.  
  594.  
  595.  
  596.  
  597.       This  basic primitive operation is essentially the calculation of the
  598.  
  599.       quotient (a*b+c)/m and its remainder. For the example above a=5, b=9,
  600.  
  601.       c=3 and m=10.  This operation has suprisingly universal  application,
  602.  
  603.       and since it lies at the innermost loop of the arithmetic algorithms,
  604.  
  605.       its efficient implementation is essential.
  606.  
  607.  
  608.  
  609.       There  are three main difficulties with a high-level language general
  610.  
  611.       base implementation of this MAD (Multiply, Add and Divide) operation.
  612.  
  613.  
  614.  
  615.       (1)  It will be slow.
  616.  
  617.  
  618.  
  619.       (2)  Quotient  and  remainder  are  not available simultaneously as a
  620.  
  621.            result of the divide operation.  Therefore the calculation  must
  622.  
  623.            be  essentially done twice,  once to get the quotient,  and once
  624.  
  625.            for the remainder.
  626.  
  627.  
  628.  
  629.        (3) Although the operation results in two single  digit  quantities,
  630.  
  631.            the  intermediate  product (a*b+c) may be double-length.  Indeed
  632.  
  633.            such a Multiply-Add and  Divide  routine  can  be  used  on  all
  634.  
  635.            occasions when a double-length quantity would be required by the
  636.  
  637.            basic  arithmetic  algorithms.  Note  that  the  'C' language is
  638.  
  639.            blessed with a 'long' integer data-type which  may  in  fact  be
  640.  
  641.            capable of temporarily storing this product.
  642.  
  643.  
  644.  
  645.       For  these reasons it is best to implement this critical operation in
  646.  
  647.       the assembly language of the computer used,  although a portable  'C'
  648.  
  649.       version is possible. At machine-code level a transitory double-length
  650.  
  651.       result can often be dealt with, even if the 'C' long data-type is not
  652.  
  653.  
  654.                                       10                                      
  655.  
  656.  
  657.  
  658.  
  659.  
  660.  
  661.       itself double-length (as is the case for the 'C' compiler used on the
  662.  
  663.       VAX11/780, for which ints and longs are both 32-bit quantities).  For
  664.  
  665.       further details see the documentation in the file 'bnmuldv.any'.
  666.  
  667.  
  668.  
  669.       A  criticism  of  the  MIRACL system might be its use of fixed length
  670.  
  671.       arrays for its 'big' and 'flash' data types.  This was done to  avoid
  672.  
  673.       the  difficult  and  time-consuming problems of memory allocation and
  674.  
  675.       garbage  collection  which  would  be  needed  by  a  variable-length
  676.  
  677.       representation. However it does mean that when doing a calculation on
  678.  
  679.       'big'  integers  that  the  results and all intermediate calculations
  680.  
  681.       must be less than or equal to the fixed size initially  specified  to
  682.  
  683.       mirsys.      mirsys 
  684.  
  685.  
  686.  
  687.       In practise most numbers in a stable integer calculation are of  more
  688.  
  689.       or  less  the  same size,  except when two are multiplied together in
  690.  
  691.       which case a double-length intermediate product is created.  This  is
  692.  
  693.       usually  immediately  reduced again by a subsequent divide operation.
  694.  
  695.       A classic example of this would be  in  the  Pollard-Brent  factoring
  696.  
  697.       program (Appendix IV and Scott 1986a).
  698.  
  699.  
  700.  
  701.       Note  that  this  is another manifestation, on a macro level,  of the
  702.  
  703.       problem mentioned in point (3) above.  It would be a pity to have  to
  704.  
  705.       specify each variable to be twice as large as necessary, just to cope
  706.  
  707.       with  these  occasional  intermediate  products.  For  this  reason a
  708.  
  709.       special Multiply, Add and Divide routine mad has been included in the                                               mad                         
  710.  
  711.       MIRACL library.  It has proved very useful  when  implementing  large
  712.  
  713.       programs  (like  the  Pomerance-Silverman-Montgomery factoring
  714.  
  715.       program, Appendix IV) on computers with limited memory.
  716.  
  717.  
  718.  
  719.  
  720.                                       11                                      
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  
  727.       As well  as  the  basic  arithmetic  operations,  routines  are  also
  728.  
  729.       provided:
  730.  
  731.  
  732.  
  733.       (a) to generate and test 'big' prime numbers,  using a  probabilistic
  734.  
  735.           primality test (Knuth 1981 $4.5.4)
  736.  
  737.  
  738.  
  739.       (b) to generate 'big'  and  'flash'  random  numbers,  based  on  the
  740.  
  741.           subtractive method (Knuth 1981 $3.6)
  742.  
  743.  
  744.  
  745.       (c) to calculate powers and roots
  746.  
  747.  
  748.  
  749.       (d) to implement both the normal and extended euclidean GCD algorithm
  750.  
  751.           (Knuth 1981 $4.5.2)
  752.  
  753.  
  754.  
  755.       6. FLOATING-SLASH NUMBERS
  756.  
  757.  
  758.  
  759.       The  straightforward  way to represent rational numbers is as reduced
  760.  
  761.       fractions,  as a numerator and denominator with  all  common  factors
  762.  
  763.       cancelled  out.   These  numbers  can  then  be  added,   subtracted,
  764.  
  765.       multiplied and divided in the obvious way,  and the result reduced by
  766.  
  767.       dividing  both  numerator  and  denominator  by their Greatest Common
  768.  
  769.       Divisor.  An efficient GCD subroutine,  using Lehmers modification of
  770.  
  771.       the  classical euclidean algorithm for multi-precision numbers (Knuth
  772.  
  773.       1981), is included in the MIRACL package.
  774.  
  775.  
  776.  
  777.       An alternative way to  represent  rationals  would  be  as  a  finite
  778.  
  779.       continued fraction (Knuth 1981, $4.5.2).  Every rational number {p/q}
  780.  
  781.       can be written as
  782.  
  783.  
  784.  
  785.  
  786.                                       12                                      
  787.  
  788.  
  789.  
  790.  
  791.  
  792.  
  793.  
  794.                p
  795.                -      =     a0    +      1
  796.                q                 --------------------------------
  797.                                   a1    +    1
  798.                                       ----------------------------
  799.                                         a2    +    1
  800.                                             ----------------------
  801.                                                  +  ..
  802.                                                        +  ..
  803.                                                           +    1
  804.                                                               ----
  805.                                                                an
  806.  
  807.  
  808.       or more elegantly as {p/q} = [a0/a1/a2/..../an],  where  the  ai  are
  809.  
  810.       positive integers, usually quite small.
  811.  
  812.  
  813.  
  814.       For example
  815.  
  816.  
  817.  
  818.                        {277/642}  =  [0/2/3/6/1/3/3]
  819.  
  820.  
  821.  
  822.       Note  that  the  ai  elements  of  the   above   continued   fraction
  823.  
  824.       representation  are  easily found as the quotients generated as a by-
  825.  
  826.       product when the euclidean GCD algorithm is applied to p and q.
  827.  
  828.  
  829.  
  830.       As we are committed to fixed length representation  of  rationals,  a
  831.  
  832.       problem  arises  when the result of some operation exceeds this fixed
  833.  
  834.       length.  There is a necessity  for  some  scheme  of  truncation,  or
  835.  
  836.       rounding. While there is no obvious way to truncate a large fraction,
  837.  
  838.       it   is   a   simple   matter  to  truncate  the  continued  fraction
  839.  
  840.       representation.  The resulting,  smaller,  fraction is called a  best
  841.  
  842.       rational approximation, or a convergent, to the original fraction.
  843.  
  844.  
  845.  
  846.       Consider truncating {277/642} = [0/2/3/6/1/3/3]. Simply drop the last
  847.  
  848.       element from the CF representation,  giving [0/2/3/6/1/3] = {85/197},
  849.  
  850.  
  851.  
  852.                                       13                                      
  853.  
  854.  
  855.  
  856.  
  857.  
  858.  
  859.       which is a very close approximation to {277/642} (error  =  0.0018%).
  860.  
  861.       Chopping  more  terms  from  the  CF  expansion  gives the successive
  862.  
  863.       convergents  as  {22/51},  {19/44},  {3/7},  {1/2},  {0/1}.   As  the
  864.  
  865.       fractions get smaller, the error increases.  Obviously the truncation
  866.  
  867.       rule for a computer implemention  should  be  to  chose  the  biggest
  868.  
  869.       convergent that fits the computer representation.
  870.  
  871.  
  872.  
  873.       The  type  of  rounding  described  above  is  also  called  'Mediant
  874.  
  875.       rounding'.  If {p/q} and {r/s}  are  two  neighbouring  representable
  876.  
  877.       slash   numbers   astride   a   gap,   then   their  mediant  is  the
  878.  
  879.       unrepresentable {(p+r)/(q+s)}.  All larger  fractions  between  {p/q}
  880.  
  881.       and the mediant will round to {p/q},  and those between {r/s} and the
  882.  
  883.       mediant will round  to  {r/s}.  The  mediant  itself  rounds  to  the
  884.  
  885.       'simpler' of {p/q} and {r/s}.
  886.  
  887.  
  888.  
  889.       This is theoretically a very good way to round,  much better than the
  890.  
  891.       rather arbitrary and base-dependent methods  used  in  floating-point
  892.  
  893.       arithmetic,  and is the method used here.  The full theoretical basis
  894.  
  895.       of floating-slash arithmetic is  described  in  detail  by  Matula  &
  896.  
  897.       Kornerup  (1985).  It should be noted that our 'flash' representation
  898.  
  899.       is in fact a cross between  the  fixed-  and  floating-slash  systems
  900.  
  901.       analysed  by  Matula & Kornerup,  as our slash can only float between
  902.  
  903.       words,  and not between bits.  However  the  characteristics  of  the
  904.  
  905.       'flash'  data-type  will  tend  to  those  of floating-slash,  as the
  906.  
  907.       precision is increased.
  908.  
  909.  
  910.  
  911.       The MIRACL routine round implements mediant rounding.  If the  result                         round                                             
  912.  
  913.       of an arithmetic operation is the fraction {p/q},  then the euclidean
  914.  
  915.       GCD algorithm is applied as before to p and q.  However this time the
  916.  
  917.  
  918.                                       14                                      
  919.  
  920.  
  921.  
  922.  
  923.  
  924.  
  925.       objective  is  not  to use the algorithm to calculate the GCD per se,
  926.  
  927.       but to use its quotients to build successive  convergents  to  {p/q}.
  928.  
  929.       This  process is stopped when the next convergent is too large to fit
  930.  
  931.       the 'flash' representation.  The complete algorithm  is  given  below
  932.  
  933.       (Kornerup & Matula 1983)
  934.  
  935.  
  936.  
  937.       Given p>=0 and q>=1
  938.  
  939.  
  940.  
  941.                    b(-2) = p     p(-2)  = 0    q(-2) = 1
  942.  
  943.                    b(-1) = q     p(-1)  = 1    q(-1) = 0
  944.  
  945.  
  946.  
  947.       Now  for  i=0,1,....  and  for  b(i-1)  > 0 ,  find the quotient a(i)
  948.  
  949.       and remainder b(i) when b(i-2) is divided by b(i-1), such that
  950.  
  951.  
  952.  
  953.                    b(i) =-a(i).b(i-1) + b(i-2)
  954.  
  955.  
  956.  
  957.       Then calculate
  958.  
  959.  
  960.  
  961.                    p(i) = a(i).p(i-1) + p(i-2)
  962.  
  963.                    q(i) = a(i).q(i-1) + q(i-2)
  964.  
  965.  
  966.  
  967.       Stop when {p(i)/q(i)} is to big to fit the  'flash'   representation,
  968.  
  969.       and take {p(i-1)/q(i-1)} as the rounded result.
  970.  
  971.  
  972.  
  973.       If applied to {277/642},  this process will give the same sequence of
  974.  
  975.       convergents as stated earlier.
  976.  
  977.  
  978.  
  979.       Since  this  rounding procedure must be applied to the result of each
  980.  
  981.       arithmetic operation, and since it is potentially rather slow,  a lot
  982.  
  983.  
  984.                                       15                                      
  985.  
  986.  
  987.  
  988.  
  989.  
  990.  
  991.       of effort has been made to optimize its implementation. Lehmer's idea
  992.  
  993.       of  operating only with the most significant piece of each number for
  994.  
  995.       as long as possible (Knuth 1981) is used,  so that for  most  of  the
  996.  
  997.       iterations  only single-precision arithmetic is needed.  Special care
  998.  
  999.       is taken to avoid the rounded result overshooting the limits  of  the
  1000.  
  1001.       'flash'  representation  (Scott 1986b).  The application of the basic
  1002.  
  1003.       arithmetic routines to the calculation of elementary  functions  such
  1004.  
  1005.       as log(x), exp(x), sin(x), cos(x), tan(x) etc., is described in Scott
  1006.  
  1007.       (1988).
  1008.  
  1009.  
  1010.  
  1011.       In  many  cases the result given by a program can be guaranteed to be
  1012.  
  1013.       exact.  This can be checked by testing  the  global  variable  EXACT,
  1014.  
  1015.       which  is  initialised  to  TRUE and is only set to FALSE if rounding
  1016.  
  1017.       takes place.
  1018.  
  1019.  
  1020.  
  1021.  
  1022.  
  1023.       7. SOME RESULTS
  1024.  
  1025.  
  1026.  
  1027.       The MIRACL library has been installed  and  tested  on  a  VAX11/780,
  1028.  
  1029.       using Digital's own C compiler, on an IBM PC, using (amongst  others)
  1030.  
  1031.       the Microsoft V5.0, the Mark Williams V2.0,  Zorland V1.1, Aztec V3.4
  1032.  
  1033.       and Borlands Turbo C V1.0 compilers, on an Acorn Archimedes using the
  1034.  
  1035.       Norcroft ARM  C V1.54A, and on an Apple Macintosh, using the  Megamax
  1036.  
  1037.       compiler.  The only necessary changes required are to the  'mirdef.h'
  1038.  
  1039.       header  file.  Special  assembly  language  versions  of the critical
  1040.  
  1041.       muldiv function for these implementations are given in  the  file
  1042.  
  1043.       'bnmuldv.any', together with a portable 'C' version.
  1044.  
  1045.  
  1046.  
  1047.       Many example programs are described in Appendix IV.  Six  different
  1048.  
  1049.  
  1050.                                       16                                      
  1051.  
  1052.  
  1053.  
  1054.  
  1055.  
  1056.  
  1057.       factoring  programs  are  included.  The biggest and most powerful of
  1058.  
  1059.       these   is   the   Pomerance-Silverman-Montgomery algorithm
  1060.  
  1061.       (Pomerance 1985, Silverman 1987) which factored F7 = 2^128 -1 =
  1062.  
  1063.  
  1064.  
  1065.                  340282366920938463463374607431768211457
  1066.  
  1067.  
  1068.  
  1069.       in less than 45 minutes, running on an IBM PC-AT micro-computer. When
  1070.  
  1071.       this number was first factored, it took 90 minutes on an IBM 360
  1072.  
  1073.       mainframe (Morrison & Brillhart 1975), albeit using a somewhat
  1074.  
  1075.       inferior algorithm. Another program implements Lenstra's recently
  1076.  
  1077.       discovered Elliptic Curve method (Montgomery 1987).
  1078.  
  1079.  
  1080.  
  1081.       Two Public-Key cryptography systems, whose strength appears to depend
  1082.  
  1083.       on  the  difficulty  of  factorisation,  are given.  The first is the
  1084.  
  1085.       classic RSA system (Rivest,  Shamir & Adleman 1978).  This is fast to
  1086.  
  1087.       encode a message, but painfully slow at decoding.  A much faster, but
  1088.  
  1089.       less  tried  and  tested,  approach  is  due  to  Okamato  (1986).  A
  1090.  
  1091.       disadvantage  of  this method is  that the  encoded text is  somewhat
  1092.  
  1093.       bigger than  the  input  source.   For  both  methods  the  keys  are
  1094.  
  1095.       constructed from 'strong' primes (Gordon 1984), to improve security.
  1096.  
  1097.  
  1098.  
  1099.       Several programs demonstrate the use of 'flash' variables.  One gives
  1100.  
  1101.       an implementation of Guassian elimination to solve a  set  of  linear
  1102.  
  1103.       equations,  involving  notoriously  ill-conditioned Hilbert matrices.
  1104.  
  1105.       Others show how rational arithmetic can be used to  approximate  real
  1106.  
  1107.       arithmetic,  in,  for  example  the calculation of roots and pi.  The
  1108.  
  1109.       former program detected an error in the value for the square root  of
  1110.  
  1111.       5 given in Knuth (1981) appendix A. The correct value is
  1112.  
  1113.  
  1114.  
  1115.  
  1116.                                       17                                      
  1117.  
  1118.  
  1119.  
  1120.  
  1121.  
  1122.  
  1123.               2.23606 79774 99789 69640 91736 68731 27623 54406                                                    _          
  1124.  
  1125.  
  1126.  
  1127.       The error is in the tenth last digit, which is a 2, and not a 1.
  1128.  
  1129.  
  1130.  
  1131.       The  roots program runs particularly fast when calculating the square
  1132.  
  1133.       roots of single precision integers,  as a simple  form  of  continued
  1134.  
  1135.       fraction generator can be used (Scott 1987a). In one test the  golden
  1136.  
  1137.       ratio was calculated to 100,000 decimal places in 3 hours of CPU time
  1138.  
  1139.       on a VAX11/780.
  1140.  
  1141.  
  1142.  
  1143.       The 'pi' program was used to calculate pi  correct  to  1000  decimal
  1144.  
  1145.       places, taking about 6 minutes of CPU time on the VAX to do so.
  1146.  
  1147.  
  1148.  
  1149.  
  1150.  
  1151.       8. CRITICISM and FUTURE DEVELOPMENTS
  1152.  
  1153.  
  1154.  
  1155.       A  disadvantage  of  using  a 'flash' type of variable to approximate
  1156.  
  1157.       real  arithmetic  is   the   non-uniformity   in   gap-size   between
  1158.  
  1159.       representable  values  (Matula  & Kornerup 1985).  To illustrate this
  1160.  
  1161.       consider a floating-slash system which is  constrained  to  have  the
  1162.  
  1163.       product of numerator and denominator less than 256. All representable
  1164.  
  1165.       values  in  the  range 0-1 are given in Appendix I.  Observe that the
  1166.  
  1167.       first representable fraction less than {1/1}  in  such  a  system  is
  1168.  
  1169.       {15/16},  a  gap  of  {1/16}.  The next fraction larger than {0/1} is
  1170.  
  1171.       {1/255},  a gap of {1/255}.  In general,  for a k-bit  floating-slash
  1172.  
  1173.       system,  the gap size varies from smaller than 2^(-k) to a worst case
  1174.  
  1175.       2^(-k/2). In practise this means that a real value  that  falls  into
  1176.  
  1177.       one of the larger gaps,  will be represented by a fraction which will
  1178.  
  1179.       be accurate to only half its usual precision.  Fortunately such large
  1180.  
  1181.  
  1182.                                       18                                      
  1183.  
  1184.  
  1185.  
  1186.  
  1187.  
  1188.  
  1189.       gaps  are rare,  and increasingly so for higher precision,  occurring
  1190.  
  1191.       only near simple fractions.  However it does mean that  real  results
  1192.  
  1193.       can  only  be completely trusted to half the given decimal places.  A
  1194.  
  1195.       partial solution to this problem  would  be  to  represent  rationals
  1196.  
  1197.       directly as continued fractions.  This gives a much better uniformity
  1198.  
  1199.       of gap-size (Kornerup & Matula 1985),  but would be very difficult to
  1200.  
  1201.       implement using a high level language.
  1202.  
  1203.  
  1204.  
  1205.       Arithmetic  on  'flash'  data-types  is undoubtedly slower than on an
  1206.  
  1207.       equivalent  sized  multiprecision  floating-point  type  (e.g.  Brent
  1208.  
  1209.       1978).  The  advantages  of  the  'flash' approach are its ability to
  1210.  
  1211.       exactly represent rational numbers,  and do exact arithmetic on them.
  1212.  
  1213.       Even  when rounding is needed,  the result often works out correctly,
  1214.  
  1215.       due to the tendency of mediant-rounding to prefer a  simple  fraction
  1216.  
  1217.       over  a  complex  one.  For example the 'roots' program (Appendix IV)
  1218.  
  1219.       when asked to find the square root of 2 and then square  the  result,
  1220.  
  1221.       comes back with the exact answer of 2, despite internal rounding.  So
  1222.  
  1223.       the MIRACL approach should appeal to those who prefer  an  answer  of
  1224.  
  1225.       10, to one of 9.9999997.
  1226.  
  1227.  
  1228.  
  1229.       Many users of the MIRACL package would be disappointed that they have
  1230.  
  1231.       to calculate
  1232.  
  1233.  
  1234.  
  1235.                             t = x.x  + x + 1
  1236.  
  1237.  
  1238.  
  1239.       by the sequence
  1240.  
  1241.                             fmul(x,x,t);
  1242.                             fadd(t,x,t);
  1243.                             fincr(t,1,1,t);
  1244.  
  1245.       rather than by simply t=x*x+x+1;
  1246.  
  1247.  
  1248.                                       19                                      
  1249.  
  1250.  
  1251.  
  1252.  
  1253.  
  1254.  
  1255.  
  1256.  
  1257.       Someone  could  of  course  use the MIRACL library to write a special
  1258.  
  1259.       purpose C compiler which could properly interpret such a command (see
  1260.  
  1261.       Cherry and Morris 1984 for an  example  of  this  approach).  However
  1262.  
  1263.       such  a drastic step is not necessary.  A superset of C,  called C++,
  1264.  
  1265.       has been developed by AT&T (Stroustrup 1986),  which may gain general
  1266.  
  1267.       acceptance  as the natural successor to C.  The enhancements to C are
  1268.  
  1269.       mainly aimed at making it an object-oriented  language.  By  defining
  1270.  
  1271.       'big' and 'flash' variables as 'classes' (in C++ terminology), it
  1272.  
  1273.       will be possible to  'overload'  the usual mathematical operators, so
  1274.  
  1275.       that the compiler will automatically substitute  calls  to  the
  1276.  
  1277.       appropriate  MIRACL  routines  when these operators are used in
  1278.  
  1279.       conjunction with 'big'  or  'flash'  variables. Furthermore  C++
  1280.  
  1281.       will  be able to look after the initialization (and ultimate
  1282.  
  1283.       elimination) of these data-types  automatically,  using  its
  1284.  
  1285.       constructor/destructor constructs,  which are included with the class
  1286.  
  1287.       definition.  Indeed once the classes are properly defined and set up,
  1288.  
  1289.       it  should  be  as simple to work with the new data-types as with the
  1290.  
  1291.       built-in double and int types.  It should then be possible  to  adapt
  1292.  
  1293.       the C source of existing mathematical libraries, with minimal effort,
  1294.  
  1295.       to work with the new data-types.
  1296.  
  1297.  
  1298.  
  1299.       A C++ implementation will be the object of further research.
  1300.  
  1301.  
  1302.  
  1303.  
  1304.  
  1305.  
  1306.  
  1307.  
  1308.  
  1309.  
  1310.  
  1311.  
  1312.  
  1313.  
  1314.                                       20                                      
  1315.  
  1316.  
  1317.  
  1318.  
  1319.  
  1320.  
  1321.   
  1322.  
  1323.       9. REFERENCES
  1324.  
  1325.  
  1326.  
  1327.  
  1328.  
  1329.       BRENT, R.P.  A Fortran Multi-precision Arithmetic Package. ACM Trans.
  1330.  
  1331.       Math. Software 4,1 (March 1978), 57-81.
  1332.  
  1333.  
  1334.  
  1335.       CHERRY, L. AND MORRIS, R. BC - An Arbitrary Precision Desk-Calculator
  1336.  
  1337.       Language. in ULTRIX-32 Supplementary Documents Vol.  1 General Users.
  1338.  
  1339.       Digital Equipment Corporation 1984.
  1340.  
  1341.  
  1342.  
  1343.       GORDON, J. Strong Primes are easy to find. In Advances in Cryptology,
  1344.  
  1345.       Proc. Eurocrypt 84, Lecture Notes in Computer Science, Vol. 209, 216-
  1346.  
  1347.       223, Springer-Verlag, 1985.
  1348.  
  1349.  
  1350.  
  1351.       KNUTH,  D.E.  The  Art  of Computer Programming,  Vol 1:  Fundamental
  1352.  
  1353.       Algorithms. Addison-Wesley, Reading, Mass., 1973.
  1354.  
  1355.  
  1356.  
  1357.       KNUTH,  D.E.  The Art of Computer Programming,  Vol 2:  Seminumerical
  1358.  
  1359.       Algorithms. Addison-Wesley, Reading, Mass., 1981.
  1360.  
  1361.  
  1362.  
  1363.       KORNERUP, P.  AND MATULA, D.W.  Finite Precision Rational Arithmetic:
  1364.  
  1365.       An Arithmetic Unit.  IEEE Trans.  Comput., C-32,4 (April 1983),  378-
  1366.  
  1367.       387.
  1368.  
  1369.  
  1370.  
  1371.       KORNERUP,  P.   AND  MATULA,  D.W.   Finite  Precision  Lexicographic
  1372.  
  1373.       Continued  Fraction  Number  Systems.   Proc.   7th  Sym.   on  Comp.
  1374.  
  1375.       Arithmetic, IEEE Cat.  #85CH2146-9, 1985, 207-214.
  1376.  
  1377.  
  1378.  
  1379.  
  1380.                                       21                                      
  1381.  
  1382.  
  1383.  
  1384.  
  1385.  
  1386.  
  1387.       MATULA, D.W.  AND KORNERUP, P.  Finite Precision Rational Arithmetic:
  1388.  
  1389.       Slash Number Systems. IEEE Trans.  Comput., C-34,1 (January 1985), 3-
  1390.  
  1391.       18.
  1392.  
  1393.  
  1394.  
  1395.       MORRISON,  M.A.  AND  BRILLHART,  J.  A  Method  of Factoring and the
  1396.  
  1397.       Factorization of F7. Math. Comput., 29,129 (January 1975), 183-205.
  1398.  
  1399.  
  1400.  
  1401.       OKAMATO,  T.  Fast Public-Key Cryptosystem using Congruent Polynomial
  1402.  
  1403.       Equations. Electr. Letters, 22,11 (May 1986), 581-582.
  1404.  
  1405.  
  1406.  
  1407.       POMERANCE,  C. The Quadratic Sieve Factoring Algorithm. In Advances
  1408.  
  1409.       in Cryptology, Lecture Notes in Computer Science, Vol. 209, Springer-
  1410.  
  1411.       Verlag, 1985, 169-182
  1412.  
  1413.  
  1414.  
  1415.       RIVEST, R., SHAMIR, A. AND ADLEMAN, L. A Method for obtaining Digital
  1416.  
  1417.       Signatures and Public-Key Cryptosystems.  Comm.  ACM,  21,2 (February
  1418.  
  1419.       1978), 120-126.
  1420.  
  1421.  
  1422.  
  1423.       SCOTT, M.P.J. Review Feedback letter. Byte, 11,6 (June 1986) 289.
  1424.  
  1425.  
  1426.  
  1427.       SCOTT,   M.P.J.   Fast  rounding  in  multiprecision   floating-slash
  1428.  
  1429.       arithmetic.  Accepted for publication IEEE Transactions on Computers.
  1430.  
  1431.  
  1432.  
  1433.       SCOTT, M.P.J.  On the Efficient Generation of Multiprecision Rational
  1434.  
  1435.       Approximations. NIHE Dublin Working Paper CA 0487, 1987.
  1436.  
  1437.  
  1438.  
  1439.       SCOTT, M.P.J. Evaluation of elementary functions using multiprecision
  1440.  
  1441.       rational arithmetic. In preparation. 1988.
  1442.  
  1443.  
  1444.  
  1445.  
  1446.                                       22                                      
  1447.  
  1448.  
  1449.  
  1450.  
  1451.  
  1452.  
  1453.       SILVERMAN, R.D. The Multiple Polynomial Quadratic Sieve, Math. Comp.
  1454.  
  1455.       48, 177, (January 1987), 329-339
  1456.  
  1457.  
  1458.  
  1459.       STROUSTRUP, B. The C++ Programming Language. Addison-Wesley, Reading,
  1460.  
  1461.       Mass., 1986.
  1462.  
  1463.  
  1464.  
  1465.  
  1466.  
  1467.  
  1468.  
  1469.  
  1470.  
  1471.  
  1472.  
  1473.  
  1474.  
  1475.  
  1476.  
  1477.  
  1478.  
  1479.  
  1480.  
  1481.  
  1482.  
  1483.  
  1484.  
  1485.  
  1486.  
  1487.  
  1488.  
  1489.  
  1490.  
  1491.  
  1492.  
  1493.  
  1494.  
  1495.  
  1496.  
  1497.  
  1498.  
  1499.  
  1500.  
  1501.  
  1502.  
  1503.  
  1504.  
  1505.  
  1506.  
  1507.  
  1508.  
  1509.  
  1510.  
  1511.  
  1512.                                       23                                      
  1513.  
  1514.  
  1515.  
  1516.  
  1517.  
  1518.  
  1519.  
  1520.  
  1521.                       Appendix I
  1522.  
  1523.     Representable values  in  range  0-1  for  a  slash  system  for  which
  1524.     numerator*denominator < 256.
  1525.  
  1526.   0/1 1/255 1/254 1/253 1/252 1/251 1/250 1/249 1/248 1/247
  1527. 1/246 1/245 1/244 1/243 1/242 1/241 1/240 1/239 1/238 1/237 1/236
  1528. 1/235 1/234 1/233 1/232 1/231 1/230 1/229 1/228 1/227 1/226 1/225
  1529. 1/224 1/223 1/222 1/221 1/220 1/219 1/218 1/217 1/216 1/215 1/214
  1530. 1/213 1/212 1/211 1/210 1/209 1/208 1/207 1/206 1/205 1/204 1/203
  1531. 1/202 1/201 1/200 1/199 1/198 1/197 1/196 1/195 1/194 1/193 1/192
  1532. 1/191 1/190 1/189 1/188 1/187 1/186 1/185 1/184 1/183 1/182 1/181
  1533. 1/180 1/179 1/178 1/177 1/176 1/175 1/174 1/173 1/172 1/171 1/170
  1534. 1/169 1/168 1/167 1/166 1/165 1/164 1/163 1/162 1/161 1/160 1/159
  1535. 1/158 1/157 1/156 1/155 1/154 1/153 1/152 1/151 1/150 1/149 1/148
  1536. 1/147 1/146 1/145 1/144 1/143 1/142 1/141 1/140 1/139 1/138 1/137
  1537. 1/136 1/135 1/134 1/133 1/132 1/131 1/130 1/129 1/128 1/127 1/126
  1538. 1/125 1/124 1/123 1/122 1/121 1/120 1/119 1/118 1/117 1/116 1/115
  1539. 1/114 1/113 1/112 1/111 1/110 1/109 1/108 1/107 1/106 1/105 1/104
  1540. 1/103 1/102 1/101 1/100 1/99  1/98  1/97  1/96  1/95  1/94  1/93
  1541. 1/92  1/91  1/90  1/89  1/88  1/87  1/86  1/85  1/84  1/83  1/82
  1542. 1/81  1/80  1/79  1/78  1/77  1/76  1/75  1/74  1/73  1/72  1/71
  1543. 1/70  1/69  1/68  1/67  1/66  1/65  1/64  2/127 1/63  2/125 1/62
  1544. 2/123 1/61  2/121 1/60  2/119 1/59  2/117 1/58  2/115 1/57  2/113
  1545. 1/56  2/111 1/55  2/109 1/54  2/107 1/53  2/105 1/52  2/103 1/51
  1546. 2/101 1/50  2/99  1/49  2/97  1/48  2/95  1/47  2/93  1/46  2/91
  1547. 1/45  2/89  1/44  2/87  1/43  2/85  1/42  2/83  1/41  2/81  1/40
  1548. 2/79  1/39  2/77  1/38  2/75  1/37  2/73  1/36  2/71  1/35  2/69
  1549. 1/34  2/67  1/33  2/65  1/32  2/63  1/31  2/61  1/30  2/59  1/29
  1550. 2/57  3/85  1/28  3/83  2/55  3/82  1/27  3/80  2/53  3/79  1/26
  1551. 3/77  2/51  3/76  1/25  3/74  2/49  3/73  1/24  3/71  2/47  3/70
  1552. 1/23  3/68  2/45  3/67  1/22  3/65  2/43  3/64  1/21  3/62  2/41
  1553. 3/61  1/20  3/59  2/39  3/58  1/19  3/56  2/37  3/55  1/18  3/53
  1554. 2/35  3/52  1/17  3/50  2/33  3/49  1/16  4/63  3/47  2/31  3/46
  1555. 4/61  1/15  4/59  3/44  2/29  3/43  4/57  1/14  4/55  3/41  2/27
  1556. 3/40  4/53  1/13  4/51  3/38  2/25  3/37  4/49  1/12  4/47  3/35
  1557. 2/23  3/34  4/45  1/11  4/43  3/32  2/21  3/31  4/41  5/51  1/10
  1558. 5/49  4/39  3/29  5/48  2/19  5/47  3/28  4/37  5/46  1/9   5/44
  1559. 4/35  3/26  5/43  2/17  5/42  3/25  4/33  5/41  1/8   5/39  4/31
  1560. 3/23  5/38  2/15  5/37  3/22  4/29  5/36  1/7   6/41  5/34  4/27
  1561. 3/20  5/33  2/13  5/32  3/19  4/25  5/31  6/37  1/6   6/35  5/29
  1562. 4/23  3/17  5/28  2/11  5/27  3/16  4/21  5/26  6/31  7/36  1/5
  1563. 7/34  6/29  5/24  4/19  7/33  3/14  5/23  7/32  2/9   7/31  5/22
  1564. 3/13  7/30  4/17  5/21  6/25  7/29  1/4   8/31  7/27  6/23  5/19
  1565. 4/15  7/26  3/11  8/29  5/18  7/25  2/7   7/24  5/17  8/27  3/10
  1566. 7/23  4/13  5/16  6/19  7/22  8/25  9/28  1/3   9/26  8/23  7/20
  1567. 6/17  5/14  9/25  4/11  7/19  3/8   8/21  5/13  7/18  9/23  2/5
  1568. 9/22  7/17  5/12  8/19  3/7   10/23 7/16  4/9   9/20  5/11  6/13
  1569. 7/15  8/17  9/19  10/21 11/23 1/2   11/21 10/19 9/17  8/15  7/13
  1570. 6/11  11/20 5/9   9/16  4/7   11/19 7/12  10/17 3/5   11/18 8/13
  1571. 5/8   12/19 7/11  9/14  11/17 2/3   13/19 11/16 9/13  7/10  12/17
  1572. 5/7   13/18 8/11  11/15 3/4   13/17 10/13 7/9   11/14 4/5   13/16
  1573. 9/11  14/17 5/6   11/13 6/7   13/15 7/8   15/17 8/9   9/10  10/11
  1574. 11/12 12/13 13/14 14/15 15/16 1/1
  1575.  
  1576.  
  1577.  
  1578.                                       24                                      
  1579.  
  1580.  
  1581.  
  1582.  
  1583.  
  1584.  
  1585.  
  1586.                       Appendix II - Header files
  1587.  
  1588.   ***** Deleted to save space - insert mirdef.h/miracl.h files here  *******
  1589.  
  1590.  
  1591.  
  1592.  
  1593.  
  1594.  
  1595.  
  1596.  
  1597.  
  1598.  
  1599.  
  1600.  
  1601.  
  1602.  
  1603.  
  1604.  
  1605.  
  1606.  
  1607.  
  1608.  
  1609.  
  1610.  
  1611.  
  1612.  
  1613.  
  1614.  
  1615.  
  1616.  
  1617.  
  1618.  
  1619.  
  1620.  
  1621.  
  1622.  
  1623.  
  1624.  
  1625.  
  1626.  
  1627.  
  1628.  
  1629.  
  1630.  
  1631.  
  1632.  
  1633.  
  1634.  
  1635.  
  1636.  
  1637.  
  1638.  
  1639.  
  1640.  
  1641.  
  1642.  
  1643.  
  1644.                                       25                                      
  1645.  
  1646.  
  1647.  
  1648.  
  1649.  
  1650.  
  1651.