home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 247_01 / manual1.doc < prev    next >
Encoding:
Text File  |  1989-04-19  |  43.5 KB  |  1,651 lines

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