home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / NUTUG11.ZIP / REALFAST.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-11-22  |  18.5 KB  |  476 lines
  1. Unit RealFast;
  2.  
  3. Interface
  4. Uses    Crt;
  5.  
  6.      TYPE  RealRec = record case byte of
  7.                        1 : (R : real);
  8.                        2 : (E : byte);
  9.                        3 : (F : array [1..5] of byte;
  10.                             S : byte);
  11.                      end; {record}
  12.  
  13.        DoubleRec = Record
  14.                      case byte of
  15.                        1 : (R  : real);
  16.                        2 : (F1 : array [1..6] of byte;
  17.                             E  : integer);
  18.                        3 : (F2 : array [1..7] of byte;
  19.                             S  : byte);
  20.                    end;
  21.  
  22. FUNCTION Intel8087Present : Boolean;
  23.          {Returns TRUE if machine has a 80x87 Processor           }
  24.  
  25. PROCEDURE MulDiv (VAR Input : RealRec; N : integer);
  26.          {Multiplies or divides (depending on the sign of 'N') the R
  27.          component of 'Input' by two to the power 'N'. Generates a
  28.          run-time error if a multiplication would cause an overflow.
  29.          Protection is also provided against underflow.            }
  30.  
  31. PROCEDURE MulDiv87 (VAR Input : DoubleRec; N : integer);
  32.           {This is an 8087 version of MulDiv.               }
  33.  
  34.  
  35. (*     Note: Version 4.0 changes the name of the 8087 version 3.0 real to
  36.        Double; real now stands the 6 byte real regardless of whether the
  37.        program was compiled to work with the 8087.
  38.              {NORTHWESTERN UNIVERSITY TURBO USERS GROUP UTILITIES}
  39.  
  40.                 (** NUtility Fast Floating Point Operations **)
  41.  
  42.                            {(C) J. E. Hilliard 1986}
  43.  
  44.  
  45.    {/TEN MICROSECOND MULTIPLICATION, DIVISION AND NEGATION OF REAL NUMBERS
  46.      ---------------------------------------------------------------------
  47.  
  48.      A ten microsecond floating-point multiplication, division and nega-
  49.      tion on a plain PC?  Yes, but (surprise) there is a catch. As Henry
  50.      Ford might have put it: 'You can multiply or divide by any number you
  51.      want so long as it is a power of two'. This condition is not quite so
  52.      restrictive as it might first appear, because several commonly used
  53.      numerical techniques such as the Simplex procedure for curve fitting
  54.      and fast Fourier transforms require just that operation. With respect
  55.      to negation (changing the sign) there are no restrictions.
  56.  
  57.      The trick is an old one in assembly language programming and in-
  58.      volves, as you will probably have guessed, manipulating the exponent
  59.      of the number. Floating point numbers are stored as two fields: a
  60.      mantissa and an exponent. In addition, one bit is used to define the
  61.      sign of the mantissa. Since the number is in binary format, adding an
  62.      integer N to the exponent corresponds to a multiplication of the
  63.      number by two to the power N. Similarly, a substraction yields a
  64.      division. Negation is performed by flipping the sign bit.
  65.  
  66.      Because only integer arithmetic is involved, these operations are
  67.      very fast as compared with the built in floating-point commands. The
  68.      implementation of these operations depends on the number storage
  69.      format which differs for the 8087 and non-8087 versions of Turbo
  70.      Pascal. We will therefore consider the two versions separately.
  71.  
  72.  
  73.      Non-8087 Version
  74.      ----------------
  75.      In this version real numbers occupy six bytes. The first holds the
  76.      exponent and the remainder the mantissa which, following the usual
  77.      topsy-turvy convention, is stored with the most significant byte
  78.      last. The first bit of this byte is the sign bit. A convenient way of
  79.      accessing the exponent byte and the sign bit is to define a variant
  80.      record type:                                                      /}
  81.  
  82. (*     TYPE  RealRec = record case byte of
  83.                        1 : (R : real);
  84.                        2 : (E : byte);
  85.                        3 : (F : array [1..5] of byte;
  86.                             S : byte);
  87.                      end; {record}  *)
  88.  
  89.      {/Being a variant record, the E byte overlays the first byte (ie. the
  90.      exponent) of the number R. The only purpose of component F is to
  91.      place component S in the byte holding the sign bit. It should be
  92.      noted that a variable of type RealRec occupies no more space than a
  93.      variable defined as real and, somewhat surprisingly, the executions
  94.      times for normal operations on the R component are identical with
  95.      those for a real variable. Thus nothing is sacrificed by using this
  96.      declaration.
  97.  
  98.      The following illustrates a multiplication by two:
  99.  
  100.      VAR DataPoint : RealRec;
  101.  
  102.      begin
  103.        DataPoint.R := 1.11111;
  104.        DataPoint.E := succ (DataPoint.E);
  105.        write (DataPoint.R);
  106.      end.
  107.  
  108.      [Note: The succ () command takes up less memory and executes faster
  109.      than the numerically equivalent: DataPoint.E + 1.] For division by
  110.      two the succ () command is replaced by pred (). Obviously, adding 2,
  111.      3, 4 . . to the exponent is equivalent to multiplications by 4, 16,
  112.      64 . . and similarly for divisions.
  113.  
  114.      The above coding provides no protection against overflow (ie. incre-
  115.      menting the exponent beyond the allowable maximum). If such protec-
  116.      tion is required the following should be substituted for the second
  117.      line:
  118.  
  119.        if DataPoint.E = $FF
  120.          then
  121.            DataPoint.R := Exp (1E3)
  122.          else
  123.            DataPoint.E := succ (DataPoint.E);
  124.  
  125.      If executed, the third line will generate a run-time error 01.
  126.  
  127.      Protection against underflow when dividing by two is simpler because
  128.      no error message is required:
  129.  
  130.        if DataPoint.E <> 0 then
  131.          DataPoint.E := pred (DataPoint.E);
  132.  
  133.      To change the sign it is only necessary to OR the sign bit; thus:
  134.  
  135.        DataPoint.S := DataPoint.S or $80;
  136.  
  137.  
  138.      8087 Version
  139.      ------------
  140.  
  141.      This version uses an 8087 format known as Real8 which is 8 bytes
  142.      long. The first 6 bytes are devoted to the mantissa. The exponent and
  143.      the sign bit occupy parts of the last two bytes and are stored as
  144.      follows:
  145.  
  146.                           Byte 7              Byte 8
  147.                       E E E E M M M M    S E E E E E E E
  148.  
  149.      in which E, M and S denote exponent, mantissa and sign bits respec-
  150.      tively. (The reason for this odd packing is for compatibility with a
  151.      shorter 8087 format, Real4). The appropriate variant record for this
  152.      format is:                                                         /}
  153.  
  154. (*       DoubleRec = Record
  155.                      case byte of
  156.                        1 : (R  : real);
  157.                        2 : (F1 : array [1..6] of byte;
  158.                             E  : integer);
  159.                        3 : (F2 : array [1..7] of byte;
  160.                             S  : byte);
  161.                    end; *)
  162.  
  163.      {/Again, the only purpose of F1 and F2 is to position the other compo-
  164.      nents. Unlike the non-8087 case, the E component is not actually the
  165.      exponent but is related to it by:
  166.  
  167.        Exponent = (DataPoint.E shr 4) and $7FF,
  168.  
  169.      in which the $7FF mask is required to eliminate the sign bit. It is
  170.      important to note that since E is declared as an integer, Bytes 7 and
  171.      8 will be reversed. It therefore follows that a multiplication by two
  172.      requires adding 10 hex to E as shown in the following example:
  173.  
  174.      VAR DataPoint : DoubleRec;
  175.  
  176.      Begin
  177.        DataPoint.R := 1.1111;
  178.        DataPoint.E := DataPoint.E + $10;
  179.        write (DataPoint.R);
  180.      End.
  181.  
  182.      For protection against overflow when multiplying by two the coding
  183.      should be changed to:
  184.  
  185.        if (DataPoint.E shr 4) and $7FF = $7FF
  186.          then
  187.            DataPoint.R := Exp (1E3)
  188.          else
  189.            DataPoint.E := DataPoint.E + $10;
  190.  
  191.      and, for protection against underflow when dividing by two:
  192.  
  193.        if (DataPoint.E shr 4) and $7FF <> 0 then
  194.          DataPoint.E := DataPoint.E - $10;
  195.  
  196.      The code for a change in sign is identical to that for the non-8087
  197.      version.
  198.  
  199.      At a considerable sacrifice in execution speed, the multiplication
  200.      and division routines can be replaced by a call to one of the MulDiv
  201.      procedures listed at the end. A function Compiled87 is also included
  202.      for checking which version of Turbo was used for compiling the
  203.      source.
  204.  
  205.      TIMINGS
  206.      -------
  207.      The execution times listed in the Table are for an Intel 8088
  208.      processor running at 4.77 Mhz. All times are in microsec. and were
  209.      determined by executing the code 10,000 times. The time required was
  210.      measured with a resolution of one ms using the JTimer function.
  211.      Corrections were applied for the time required for executing a null
  212.      loop and for accessing the clock. The times are rounded to three
  213.      significant figures.
  214.  
  215.      It will be seen that without range checking a 10 to 13 microsec.
  216.      multiplication or division by two can be achieved. This is some 125
  217.      times faster than normal in the plain version and 25 times faster in
  218.      the 8087 version. Although using one of the MulDiv procedures is more
  219.      convenient and saves space, there is a large penalty to be paid in
  220.      speed. Most of the extra time required is spent calling the procedure
  221.      and in parameter passing. In fact, if the 8087 Turbo version were as
  222.      fast as it should and could be, MulDiv87 would probably be slower
  223.      than a regular multiplication.
  224.  
  225.      With respect to changing the sign, flipping the sign bit is 10 to 13
  226.      times faster than the conventional instruction R := - R.
  227.  
  228.  
  229.                                    TABLE
  230.                                    -----
  231.                          Execution Times in Microsec.
  232.  
  233.            ------------------------------------------------------
  234.            Turbo V3                           Plain         8087
  235.            ------------------------------------------------------
  236.  
  237.            In-code multiply and
  238.               divide by 2 (NRC)                10.5         12.5
  239.  
  240.            In-Code multiply by 2 (WRC)         25.9         35.9
  241.  
  242.            In-Code divide by 2 (WRC)           23.1         32.1
  243.  
  244.            Procedure MulDiv(87)               177.0        200.0
  245.  
  246.            Regular multiplication            1320.0        283.0
  247.  
  248.            Regular division                  1890.0        306.0
  249.  
  250.            Fast negation (flip bit)            13.0         13.0
  251.  
  252.            Regular negation (R := - R)        127.0        180.0
  253.  
  254.            ------------------------------------------------------
  255.            Notes : NRC - No Range Check for underflow or overflow.
  256.                    WRC - With Range Check.
  257.            ------------------------------------------------------
  258.  
  259. LISTINGS
  260. --------   /}
  261.  
  262. Implementation
  263.  
  264. FUNCTION Intel8087Present : Boolean;
  265.  
  266.          {Returns TRUE if machine has a 80x87 Processor           }
  267.  
  268. Begin
  269.  
  270.   Intel8087Present := False;
  271.  
  272.  {$IFDEF CPU87}
  273.   Intel8087Present := True;
  274.  {$ENDIF}
  275.  
  276. End; {Intel8087Present : Boolean}
  277.  
  278.  
  279. PROCEDURE MulDiv (VAR Input : RealRec; N : integer);
  280.  
  281.          {Multiplies or divides (depending on the sign of 'N') the R
  282.          component of 'Input' by two to the power 'N'. Generates a
  283.          run-time error if a multiplication would cause an overflow.
  284.          Protection is also provided against underflow.            }
  285.  
  286. Begin
  287.  
  288.   if N > 0                             {Multiplication.                     }
  289.     then
  290.       if Input.E + N >= $FF
  291.         then
  292.           Input.R := Exp (1E3)         {Generate run-time error 01.         }
  293.         else
  294.           Input.E := Input.E + N
  295.     else                               {Division.                           }
  296.       if Input.E + N < 0
  297.         then                           {Reduce Input.E to zero.             }
  298.           while Input.E > 0 do
  299.             Input.E := pred (Input.E)
  300.         else
  301.           Input.E := Input.E + N;
  302.  
  303. End; {MulDiv (VAR Input : RealRec; N : integer)}
  304.  
  305.  
  306. PROCEDURE MulDiv87 (VAR Input : DoubleRec; N : integer);
  307.  
  308.           {This is an 8087 version of MulDiv.               }
  309.  
  310. Begin
  311.  
  312.   if N > 0                             {Multiplication.                     }
  313.     then
  314.       if (Input.E shr 4) and $7FF + N >= $7FF
  315.         then
  316.           Input.R := Exp (1E3)         {Generate run-time error 01.         }
  317.         else
  318.           Input.E := Input.E + (N shl 4)
  319.     else                               {Division.                           }
  320.       if (Input.E shr 4) and $7FF + N < 0
  321.         then                            {Reduce exponent to zero.           }
  322.           while (Input.E shr 4) and $7FF > 0 do
  323.               Input.E := Input.E - $10
  324.         else
  325.           Input.E := Input.E + (N shl 4);
  326.  
  327. End; {MulDiv87 (VAR Input : DoubleRec; N : integer)}
  328.  
  329.  
  330.  
  331.  
  332. PROCEDURE RealFastDemo;
  333.  
  334.          {This demo serves the following purposes: (1) Shows the coding
  335.          required for the two versions of TURBO. (2) Gives confirmation
  336.          that the operations function as they should. (3) Provides a
  337.          convincing, if only qualitative, demonstration of the speed of
  338.          the fast procedures. NOTE: This procedure can be compiled with
  339.          either the non-8087 or 8087 versions without any changes.      }
  340.  
  341. VAR
  342.  
  343.   Data   : RealRec;
  344.   Data87 : DoubleRec;
  345.   R      : Real;
  346.   J      : integer;
  347.   Ch     : Char;
  348.  
  349. Begin
  350.  
  351.   ClrScr;
  352.   GoToXY (29, 3); write ('REALFAST DEMONSTRATION');
  353.   GoToXY (29, 4); write ('----------------------');
  354.   GoToXY (1, 7);
  355.   if not Intel8087Present
  356.     then                               {Non-8087 version.                   }
  357.       begin
  358.         writeln ('Non-8087 Version':30);
  359.         writeln;
  360.         Data.R := 1.111111111;
  361.         writeln (Data.R:30, 'Starting number':30);
  362.         Data.E := succ (Data.E);          {X 2.                                }
  363.         writeln (Data.R:30, 'Fast multiplication by two':30);
  364.         Data.S := Data.S or $80;          {Change sign.                        }
  365.         writeln (Data.R:30, 'Fast change of sign':30);
  366.         Data.E := Data.E - 2;             {Divide by 4.                      }
  367.         writeln (Data.R:30, 'Fast division by four':30);
  368.         writeln;
  369.         R := 1.234567;
  370.         writeln;
  371.         writeln ('        Press any key to start 3000 REGULAR multiplications');
  372.         write ('and divisions by two ':35);
  373.         repeat {Null} until KeyPressed;
  374.         write ('Running':10);
  375.         for J := 1 to 300 do           {Coding repeated to reduce the effect}
  376.           begin                        {of the time required to execute the }
  377.             R := 2 * R;  R := R / 2;   {the loop.                           }
  378.             R := 2 * R;  R := R / 2;
  379.             R := 2 * R;  R := R / 2;
  380.             R := 2 * R;  R := R / 2;
  381.             R := 2 * R;  R := R / 2;
  382.             R := 2 * R;  R := R / 2;
  383.             R := 2 * R;  R := R / 2;
  384.             R := 2 * R;  R := R / 2;
  385.             R := 2 * R;  R := R / 2;
  386.             R := 2 * R;  R := R / 2;
  387.           end;
  388.         write ('    Finished');
  389.         Data.R := 1.234567;
  390.         writeln; writeln;
  391.         writeln ('        Press any key to start 3000 FAST multitplications');
  392.         write ('and divisions by two ':35);
  393.         repeat {Null} until KeyPressed;
  394.         write ('Running':10);
  395.         for J := 1 to 300 do
  396.           begin
  397.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  398.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  399.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  400.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  401.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  402.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  403.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  404.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  405.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  406.             Data.E := succ (Data.E);  Data.E := pred (Data.E);
  407.           end;
  408.         write ('    Finished');
  409.  
  410.       end
  411.     else                               {8087 Version.                       }
  412.       begin
  413.         writeln ('8087 Version':30);
  414.         writeln;
  415.         Data87.R := 1.1111111111111;
  416.         writeln (Data87.R:30, 'Starting number':30);
  417.         Data87.E := Data87.E + $10;          {X 2.                           }
  418.         writeln (Data87.R:30, 'Fast multiplication by two':30);
  419.         Data87.S := Data87.S or $80;          {Change sign.                  }
  420.         writeln (Data87.R:30, 'Fast change of sign':30);
  421.         Data87.E := Data87.E - $20;           {Divided by 4.                 }
  422.         writeln (Data87.R:30, 'Fast division by four':30);
  423.         R := 1.234567;
  424.         writeln;
  425.         writeln ('        Press any key to start 5000 regular multitplications');
  426.         write ('and divisions by two ':35);
  427.         repeat {Null} until KeyPressed;
  428.         Ch := ReadKey;
  429.         write ('Running':10);
  430.         for J := 1 to 500 do
  431.           begin
  432.             R := 2 * R;  R := R / 2;
  433.             R := 2 * R;  R := R / 2;
  434.             R := 2 * R;  R := R / 2;
  435.             R := 2 * R;  R := R / 2;
  436.             R := 2 * R;  R := R / 2;
  437.             R := 2 * R;  R := R / 2;
  438.             R := 2 * R;  R := R / 2;
  439.             R := 2 * R;  R := R / 2;
  440.             R := 2 * R;  R := R / 2;
  441.             R := 2 * R;  R := R / 2;
  442.           end;
  443.         write ('    Finished');
  444.         Data87.R := 1.234567;
  445.         writeln; writeln;
  446.         writeln (
  447.         '        Press any key to start 50000 FAST multitplications');
  448.         write ('and divisions by two ':35);
  449.         repeat {Null} until KeyPressed;
  450.         write ('Running':10);
  451.         for J := 1 to 5000 do
  452.           begin
  453.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  454.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  455.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  456.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  457.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  458.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  459.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  460.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  461.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  462.             Data87.E := Data87.E + $10;  Data87.E := Data87.E - $10;
  463.           end;
  464.         write ('    Finished');
  465.       end;
  466.   writeln; writeln;
  467.  
  468. End; {RealFastDemo}
  469.  
  470. BEGIN (* of Initialization for Unit RealFast  *)
  471.  
  472.  
  473.   RealFastDemo;  (* Delete this for actual use of unit *)
  474.  
  475. END.
  476.