home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / COMPUTER.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  22.8 KB  |  753 lines
  1.  
  2. (* :Name: NumericalMath`ComputerArithmetic` *)
  3.  
  4. (* :Title: Fixed Precision, Correctly Rounded Computer Arithmetic *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Implements fixed precision, rounded arithmetic.  The arithmetic
  10.     can be in any base from 2 to 16, and and any of several rounding
  11.     schemes can be used.  The range of the exponent can also be
  12.     varied within limits.  This package is not suitable for
  13.     computational purposes; it is much too slow.  Its real use is
  14.     educational, but it can be used in conjunction with other
  15.     packages such as GaussianElimination.m.
  16. *)
  17.  
  18. (* :Context: NumericalMath`ComputerArithmetic` *)
  19.  
  20. (* :Package Version: Mathematica 2.0 *)
  21.  
  22. (* :Copyright: Copyright 1990  Wolfram Research, Inc.
  23.     Permission is hereby granted to modify and/or make copies of
  24.     this file for any purpose other than direct profit, or as part
  25.     of a commercial product, provided this copyright notice is left
  26.     intact.  Sale, other than for the cost of media, is prohibited.
  27.  
  28.     Permission is hereby granted to reproduce part or all of
  29.     this file, provided that the source is acknowledged.
  30. *)
  31.  
  32. (* :History:
  33.     Revised by Jerry B. Keiper, December, 1990.
  34.     Originally by Jerry B. Keiper, May, 1990.
  35. *)
  36.  
  37. (* :Keywords: arithmetic, rounding, fixed precision *)
  38.  
  39. (* :Source:
  40.     Any elementarty numerical analysis textbook.
  41. *)
  42.  
  43. (* :Mathematica Version: 2.0 *)
  44.  
  45. (* :Warning:
  46.     The function Normal[ ] is extended to convert ComputerNumber[ ]s
  47.     to their exact rational equivalent.
  48. *)
  49.  
  50. (* :Limitations:
  51.     Changing the arithmetic can result in previous ComputerNumber[ ]s
  52.     becoming invalid.
  53.  
  54.     Only one type of arithmetic is allowed at a time.  It would not
  55.     be difficult to extend this package to support different data
  56.     types simultaneously, e.g., ComputerInteger[ ], ComputerNumber[ ],
  57.     DoubleComputerNumber[ ], etc.
  58. *)
  59.  
  60. (* :Discussion:
  61.     ComputerNumber[sign, mantissa, exp, value, x] is a data object that
  62.     represents a computer number, but its default print format is a
  63.     simple number in the base specified by SetArithmetic[ ].  The fifth
  64.     element of a ComputerNumber is the value that the number would have
  65.     if the arithmetic had been done with high-precision arithmetic,
  66.     i.e., by comparing the value of the ComputerNumber with this value
  67.     you get the total accumulated roundoff error in the number.
  68.     ComputerNumber[x] gives the complete data object representing
  69.     the number x.  ComputerNumber[sign, mantissa, exp] gives the
  70.     complete data object ComputerNumber[sign, mantissa, exp, value, x],
  71.     where value and x have the value sign * mantissa * base^exp, where
  72.     base is the the base specified by SetArithmetic[ ].  The sign must
  73.     be +1 or -1, the integer mantissa must be between base^(digits-1)
  74.     and base^digits - 1, and the exponent must be an integer within a
  75.     range specified by SetArithmetic[ ].  The fourth element, value,
  76.     is used only for efficiency.
  77.  
  78.     Basic aritmetic with these objects is automatic using any of
  79.     several rounding schemes.  Although the rounding within the
  80.     basic operations is correct, numerically converting a complicated
  81.     expression involving transcendental functions to a ComputerNumber
  82.     can result in incorrect rounding.  (Such errors are unavoidable.
  83.     Although the package could be designed to correctly round any
  84.     particular expression, there would always be some expressions
  85.     for which the rounding would be incorrect.)  For typical expressions
  86.     the rounding error will be less than .50000000000000000001 ulps
  87.     for rounding and less than 1.00000000000000000001 ulps for
  88.     truncation.  The basic aritmetic is easily extensible to the
  89.     elementary functions and even to special functions.
  90.  
  91.     Note that the default division is really two operations:
  92.     multiplication by the reciprocal.  True division is implemented
  93.     by the function IdealDivide[x, y].  These two forms of division
  94.     can give different results.
  95.     
  96.     The arithmetic implemented in this package is slightly better
  97.     than most computer arithmetic: integers and rational numbers
  98.     used in multiplication and division, and integers used as
  99.     exponents are NOT first converted to ComputerNumber[ ]s, but
  100.     rather they are used in their given form, and the final result is
  101.     then converted to a ComputerNumber.
  102. *)
  103.  
  104. (* :Examples:
  105.  
  106. In[1]:= << NumericalMath`ComputerArithmetic`    (* read in the package *)
  107.  
  108. Out[1]= NumericalMath`ComputerArithmetic`
  109.  
  110. In[2]:= Arithmetic[ ]
  111.  
  112. Out[2]= {4, 10, RoundingRule -> RoundToEven, ExponentRange -> {-50, 50},
  113.  
  114. >    MixedMode -> False}
  115.  
  116.         (* the default arithmetic is 4 digits in base 10 with a
  117.         rounding rule of RoundToEven and numbers between
  118.         10^-50 and .9999 10^50 allowed.  Mixed-mode arithmetic
  119.         is not allowed. *)
  120.  
  121. In[3]:= ComputerNumber[Pi]    (* expressions that can be interpretted as
  122.                 numbers are converted to ComputerNumber[ ]s. *)
  123.  
  124. Out[3]= 3.142
  125.  
  126. In[4]:= FullForm[%]
  127.  
  128. Out[4]//FullForm= ComputerNumber[1, 3142, -3, Rational[1571, 500],
  129.  
  130. >        3.14159265358979323846264]
  131.  
  132. In[5]:= {ComputerNumber[-1, 1234, -6], ComputerNumber[-1, 123, 7]}
  133.  
  134. Out[5]= {-0.001234, NaN}    (* you also can enter ComputerNumber[ ]s in
  135.             terms of a sign, mantissa, and an exponent, but only
  136.             if it forms a valid ComputerNumber[ ].  (in this
  137.             example the one mantissa was not 4 digits.) *)
  138.  
  139. In[6]:= ComputerNumber[Pi - 22/7]
  140.         (* expressions are evaluated numerically before
  141.                 they become "ComputerNumbers". *)
  142.  
  143. Out[6]= -0.001264
  144.  
  145. In[7]:= ComputerNumber[Pi] - ComputerNumber[22/7] (* this result is different. *)
  146.  
  147. Out[7]= -0.001
  148.  
  149. In[8]:= sum = 0 
  150.  
  151. Out[8]= 0
  152.  
  153. In[9]:= Do[sum += ComputerNumber[i]^(-2), {i, 200}]; FullForm[sum]
  154.  
  155. Out[9]= FullForm=
  156.  
  157. >   ComputerNumber[1, 1625, -3, Rational[13, 8], 1.63994654601499726794569]
  158.  
  159. In[10]:= (sum = 0; Do[sum += ComputerNumber[i]^(-2), {i, 200, 1, -1}];
  160.     FullForm[sum])
  161.  
  162. Out[10]//FullForm=
  163.  
  164. >   ComputerNumber[1, 1640, -3, Rational[41, 25], 1.63994654601499726794569]
  165.  
  166.     (* As a general rule it is better to sum the smaller terms first,
  167.     but it does not guarantee a better result: *)
  168.  
  169. In[11]:= sum = 0; Do[sum += 1/ComputerNumber[i], {i, 300}]; FullForm[sum]
  170.  
  171. Out[11]//FullForm=
  172.  
  173. >   ComputerNumber[1, 6281, -3, Rational[6281, 1000], 6.28266388029950346191949]
  174.  
  175. In[12]:= sum = 0; Do[sum += 1/ComputerNumber[i], {i,300,1,-1}]; FullForm[sum]
  176.  
  177. Out[12]//FullForm=
  178.  
  179. >   ComputerNumber[1, 6280, -3, Rational[157, 25], 6.28266388029950346191949]
  180.  
  181.     (* The difference is slight, and such examples are rare. *)
  182.  
  183. In[13]:= ComputerNumber[Sin[Pi/7]]
  184.  
  185. Out[13]= 0.4339
  186.  
  187. In[14]:= FullForm[Sin[ComputerNumber[N[Pi]/7]]]
  188.  
  189. Out[14]//FullForm= Sin[ComputerNumber[1, 4488, -4, Rational[561, 1250],
  190.  
  191. >        0.448798950512827587999709]]
  192.  
  193.     (* Basic arithmetic is all that is implemented in the package.
  194.        We could easily extend things to include elementary functions. *)
  195.  
  196. In[15]:= sq = ComputerNumber[Sqrt[47]]
  197.  
  198. Out[15]= 6.856
  199.  
  200. In[16]:= sq sq
  201.  
  202. Out[16]= 47.
  203.  
  204.     (* It is a theorem that correctly rounded square roots of small
  205.     integers will always square back to the original integer if
  206.     the arithmetic is correct.  Such is not the case for cube roots: *)
  207.  
  208. In[17]:= cr = ComputerNumber[3^(1/3)]
  209.  
  210. Out[17]= 1.442
  211.  
  212. In[18]:= cr cr cr
  213.  
  214. Out[18]= 2.998
  215.  
  216.     (* We now want to work with 7 significant digits: *)
  217.  
  218. In[19]:= SetArithmetic[7]
  219.  
  220. Out[19]= {7, 10, RoundingRule -> RoundToEven, ExponentRange -> {-50, 50}, 
  221.  
  222. >    MixedMode -> False}
  223.  
  224.     (* the arithmetic is now 7 digits in base 10 with a rounding rule
  225.     of RoundToEven and an exponent range of -50 to 50 *)
  226.  
  227. In[20]:= ComputerNumber[.9999999499999999999999999]
  228.  
  229. Out[20]= 0.9999999
  230.  
  231. In[21]:= ComputerNumber[.9999999500000000000000001]
  232.  
  233. Out[21]= 1.
  234.  
  235.     (* Note that correct rounding is performed even near the discontinuity
  236.     in the exponent. *)
  237.  
  238.  
  239.     (* The reciprocal of the reciprocal is not the original number; in
  240.     fact it may be quite different: *)
  241.  
  242. In[22]:= x = ComputerNumber[9010004]
  243.  
  244.                     6
  245. Out[22]= 9.010004 10
  246.  
  247. In[23]:= y = 1/x
  248.  
  249.             -7
  250. Out[23]= 1.109877 10
  251.  
  252. In[24]:= z = 1/y
  253.  
  254.             6
  255. Out[24]= 9.010007 10
  256.  
  257.     (* Likewise division (as a single operation) is different from
  258.     multiplication by the reciprocal: *)
  259.  
  260. In[25]:= ComputerNumber[2]/x    (* this is multiplication by the reciprocal. *)
  261.  
  262.             -7
  263. Out[25]= 2.219754 10
  264.  
  265. In[26]:= IdealDivide[ComputerNumber[2], x]    (* this is true division. *)
  266.  
  267.             -7
  268. Out[26]= 2.219755 10
  269.  
  270.     (* Note: you can set the arithmetic to use IdealDivide[ ]
  271.     automatically whenever it encounters the `/' divide symbol.
  272.     This is an option to SetArithmetic[ ], but it uses $PreRead
  273.     and may interfere with other behavior that also uses $PreRead. *)
  274.  
  275.     (* Division by 0 and overflow conditions result in NaN: *)
  276.  
  277. In[27]:= IdealDivide[ComputerNumber[1], ComputerNumber[0]]
  278.  
  279. ComputerNumber::divzer: Division by 0 occurred in computation.
  280.  
  281. Out[27]= NaN
  282.  
  283. In[28]:= x = ComputerNumber[1000000000]
  284.  
  285.           9
  286. Out[28]= 1. 10
  287.  
  288. In[29]:= x^7
  289.  
  290. Out[29]= NaN
  291.  
  292. In[30]:= Pi < 22/7
  293.  
  294.           22
  295. Out[30]= Pi < --
  296.           7
  297.  
  298. In[31]:= ComputerNumber[Pi] < ComputerNumber[22/7]
  299.  
  300. Out[31]= True
  301.  
  302. In[32]:= SetArithmetic[3, 2, RoundingRule -> Truncation,
  303.     ExponentRange -> {-3,3}]
  304.  
  305.     (* set the arithmetic to be 3 digits in base 2 using the rounding
  306.         rule of Truncation and allowed numbers between -7 and -1/8
  307.         and 1/8 and 7. *)
  308.  
  309. Out[32] = {3, 2, RoundingRule -> Truncation, ExponentRange -> {-3, 3},
  310.  
  311. >    MixedMode -> False}
  312.  
  313. In[33]:= ComputerNumber[Pi]
  314.  
  315. Out[33]= 11.
  316.         2
  317.  
  318.     (* Pi is simply 3 in this very limited arithmetic. *)
  319.  
  320. In[34]:= Plot[Normal[ComputerNumber[x]] - x, {x, -10, 10}, PlotPoints -> 47]
  321.  
  322.     (* plot the representation error of the numbers between -10 and 10. *)
  323.  
  324. In[35]:= Plot[Normal[ComputerNumber[x]] - x, {x, -1, 1}, PlotPoints -> 47]
  325.  
  326.     (* resolve the wiggles near the hole at zero. *)
  327.  
  328. In[36]:= 2 ComputerNumber[Pi] - 4
  329.  
  330. Out[36]= -4 + 2 11.
  331.            2
  332.  
  333.     (* the default is to disallow mixed-mode arithmetic, i.e.,
  334.         arithmetic between ComputerNumbers and ordinary numbers.
  335.         (Integer and Rational exponents are allowed however.) *)
  336.  
  337. In[37]:= SetArithmetic[3, 2, RoundingRule -> Truncation, 
  338.     ExponentRange -> {-3,3}, MixedMode -> True];
  339.  
  340. In[38]:= 2 ComputerNumber[Pi] - 4
  341.  
  342. Out[38]= 10.        (* now mixed-mode arithmetic is allowed. *)
  343.         2
  344. *)
  345.  
  346. BeginPackage["NumericalMath`ComputerArithmetic`"]
  347.  
  348. MixedMode::usage =
  349. "MixedMode is an option to SetArithmetic[ ] that can be either True or False.
  350. It specifies whether mixed-mode arithmetic is to be allowed.  The default is
  351. False."
  352.  
  353. RoundingRule::usage =
  354. "RoundingRule is an option to SetArithmetic[ ] and specifies the rounding
  355. scheme to use.  The choices are RoundToEven, RoundToInfinity, and Truncation.
  356. The default is RoundToEven."
  357.  
  358. RoundToEven::usage =
  359. "RoundToEven is a choice for the option RoundingRule of SetArithmetic[ ].
  360. It specifies that the rounding is to be to the nearest representable number
  361. and, in the case of a tie, round to the one represented by an even mantissa."
  362.  
  363. RoundToInfinity::usage =
  364. "RoundToInfinity is a choice for the option RoundingRule of SetArithmetic[ ].
  365. It specifies that the rounding is to be to the nearest representable number
  366. and, in the case of a tie, round away from 0."
  367.  
  368. Truncation::usage =
  369. "Truncation is a choice for the option RoundingRule of SetArithmetic[ ].
  370. It specifies that the ``rounding '' is to simply discard excess digits, as
  371. Floor[ ] does for positive numbers."
  372.  
  373. ExponentRange::usage =
  374. "ExponentRange is an option to SetArithmetic[ ] and specifies the range of
  375. exponents that are to be allowed.  The exponent range must be of the form
  376. {minexp, maxexp} where -1000 <= minexp < maxexp <= 1000.  The default
  377. ExponentRange is {-50, 50}."
  378.  
  379. SetArithmetic::usage =
  380. "SetArithmetic[dig] evaluates certain global constants used in the package
  381. ComputerArithmetic.m to make the arithmetic work properly with dig digits 
  382. precision.  The value of dig must be an integer between 1 and 10, inclusive,
  383. and the default value is 4.  SetArithmetic[dig, base] causes the arithmetic
  384. to be dig digits in base base.  The value of base must be an integer between
  385. 2 and 16, and the default is 10.  Changing the arithmetic and attempting
  386. to refer to ComputerNumbers that were defined prior to the change can
  387. lead to unpredictable results."
  388.  
  389. Arithmetic::usage =
  390. "Arithmetic[ ] gives a list containing the number of digits, the base,
  391. the rounding rule and the exponent range that are currently in effect."
  392.  
  393. ComputerNumber::usage =
  394. "ComputerNumber[sign, mantissa, exp, value, x] is a data object that represents
  395. a computer number, but its default print format is a simple number in the base
  396. specified by SetArithmetic[ ].  ComputerNumber[x] gives the complete data object
  397. representing the number x.  ComputerNumber[sign, mantissa, exp] likewise gives
  398. the complete data object ComputerNumber[sign, mantissa, exp, value, x] where
  399. value and x have the value sign * mantissa * base^exp.  The arithmetic with
  400. value is computer arithmetic; the arithmetic with x is ordinary high-precision
  401. arithmetic.  Normal[ComputerNumber[...]] gives value."
  402.  
  403. NaN::usage =
  404. "NaN is the symbol used in ComputerArithmetic.m to represent a nonrepresentable
  405. number.  NaN stands for Not-a-Number."
  406.  
  407. IdealDivide::usage =
  408. "IdealDivide[x, y] is gives the correctly rounded result of x divided by
  409. y.  The default `/' division operator in Mathematica is, in fact, multiplication
  410. by the reciprocal, involves two rounding errors, and can result in an
  411. incorrectly rounded quotient.  IdealDivide is also an option to SetArithmetic
  412. that can be set to True or False indicating whether $PreRead should be
  413. used to translate the default `/' division operator to use IdealDivide."
  414.  
  415. Options[SetArithmetic] =
  416.     {RoundingRule -> RoundToEven,
  417.     ExponentRange -> {-50, 50},
  418.     MixedMode -> False,
  419.     IdealDivide -> False}
  420.  
  421. Begin["NumericalMath`ComputerArithmetic`Private`"]
  422.  
  423. Arithmetic[ ] := {$digits, $base, RoundingRule -> $roundrule,
  424.     ExponentRange -> {$minexp - 1, $maxexp} + $digits,
  425.     MixedMode -> $mixedmode, IdealDivide -> $idealdivide}
  426.  
  427. SetArithmetic::digs =
  428. "The number of digits requested is not an integer between 1 and 10."
  429.  
  430. SetArithmetic::base = "The base requested is not an integer between 2 and 16."
  431.  
  432. SetArithmetic::rr =
  433. "The rounding rule `1` is not RoundToEven, RoundToInfinity, or Truncation."
  434.  
  435. SetArithmetic::er =
  436. "The exponent range `1` is not a pair of integers between -1000 and 1000."
  437.  
  438. SetArithmetic::mm =
  439. "The option MixedMode -> `1` is neither True nor False."
  440.  
  441. SetArithmetic::id =
  442. "The option IdealDivide -> `1` is neither True nor False."
  443.  
  444. SetArithmetic[n_Integer:4, base_Integer:10, opts___] :=
  445.     Module[{tmp, er, mm, id},
  446.     If[!(1 <= n <= 10), Message[SetArithmetic::digs]; Return[$Failed]];
  447.     If[!(2 <= base <= 16), Message[SetArithmetic::base]; Return[$Failed]];
  448.     tmp = (RoundingRule /. {opts} /. Options[SetArithmetic]);
  449.     If[!MemberQ[{RoundToEven, RoundToInfinity, Truncation}, tmp],
  450.         Message[SetArithmetic::rr, tmp];
  451.         Return[$Failed]];
  452.     er = (ExponentRange /. {opts} /. Options[SetArithmetic]);
  453.     If[!ListQ[er] || (Length[er] != 2) || !IntegerQ[er[[1]]] ||
  454.             !IntegerQ[er[[2]]] || (er[[1]] >= er[[2]]) ||
  455.             (er[[1]] < -1000) || (er[[2]] > 1000),
  456.         Message[SetArithmetic::er, er];
  457.         Return[$Failed]];
  458.     mm = (MixedMode /. {opts} /. Options[SetArithmetic]);
  459.     If[(mm =!= True) && (mm =!= False),
  460.         Message[SetArithmetic::mm, mm];
  461.         Return[$Failed]];
  462.     id = (IdealDivide /. {opts} /. Options[SetArithmetic]);
  463.     If[(id =!= True) && (id =!= False),
  464.         Message[SetArithmetic::id, id];
  465.         Return[$Failed]];
  466.     $mixedmode = mm;
  467.     $idealdivide = id;
  468.     $roundrule = tmp;
  469.     $digits = n;
  470.     $base = base;
  471.     $minman = base^(n-1);
  472.     $maxman = base $minman - 1;
  473.     $minexp = er[[1]] + 1 - n;
  474.     $maxexp = er[[2]] - n;
  475.     $prec = Round[n Log[10., base]];
  476.     $nfdigits = Max[$prec+1, $MachinePrecision+3];
  477.     $prec += 20;
  478.     If[$idealdivide,
  479.         $PreRead = DivideReplace,
  480.         If[$PreRead === DivideReplace, $PreRead = .]
  481.         ];
  482.     Update[ ];
  483.     Arithmetic[ ]
  484.     ]
  485.  
  486. DivideReplace[s_String] :=
  487.     Module[{ss = StringReplace[s , "/" -> "~IdealDivide~"]},
  488.         StringReplace[ss, {"~IdealDivide~~IdealDivide~" -> "//",
  489.         "~IdealDivide~;" -> "/;", "~IdealDivide~:" -> "/:",
  490.         "~IdealDivide~@" -> "/@", "~IdealDivide~." -> "/.",
  491.         "~IdealDivide~=" -> "/="}]];
  492.  
  493. If[!NumberQ[$digits], SetArithmetic[]];    (* initialization *)
  494.  
  495. ComputerNumber::undflw = "Underflow occured in computation.  Exponent is ``."
  496.  
  497. ComputerNumber::ovrflw = "Overflow occured in computation.  Exponent is ``."
  498.  
  499. ComputerNumber[sign_Integer, mantissa_Integer, exp_Integer] := 
  500.     Module[{tmp, value},
  501.         tmp = ((sign == 1) || (sign == -1));
  502.         If[tmp,
  503.             If[$minexp > exp,
  504.                 Message[ComputerNumber::undflw, exp + $digits];
  505.                 tmp = False];
  506.             If[exp > $maxexp,
  507.                 Message[ComputerNumber::ovrflw, exp + $digits];
  508.                 tmp = False];
  509.             ];
  510.         If[tmp && ($minman <= mantissa <= $maxman),
  511.             value = sign mantissa $base^exp;
  512.             tmp = SetPrecision[value, $prec];
  513.             tmp = ComputerNumber[sign, mantissa, exp, value, tmp],
  514.             (* else *)
  515.             If[tmp && (mantissa == 0),
  516.                 tmp = ComputerNumber[1,0,0,0,0],
  517.                 (* else *)
  518.                 tmp = NaN
  519.             ]
  520.         ];
  521.         tmp
  522.     ]
  523.  
  524. round[x_] :=
  525.     Module[{rndx1, rndx2, d},
  526.         If[$roundrule === Truncation, Return[Floor[x]]];
  527.         rndx1 = Round[x];
  528.         rndx2 = rndx1 + If[rndx1 < x, 1, -1];
  529.         d = Abs[x - rndx1] - Abs[x - rndx2];
  530.         Which[  d < 0, rndx1,
  531.             d > 0, rndx2,
  532.             True, If[$roundrule === RoundToEven,
  533.                     If[ EvenQ[rndx1], rndx1, rndx2],
  534.                     Max[{rndx1, rndx2}]]
  535.         ]
  536.     ]
  537.  
  538. ComputerNumber[x_] := x /; !NumberQ[N[x]]
  539.  
  540. ComputerNumber[x_] :=
  541.     Module[{mantissa, exp, nx = x, absnx, tmp, ok},
  542.         If[!NumberQ[nx], nx = SetPrecision[N[nx,$prec],$prec]];
  543.         If[!NumberQ[nx] || (Head[nx] === Complex), Return[NaN]];
  544.         If[nx == 0, Return[ComputerNumber[1,0,0,0,0]]];
  545.         absnx = Abs[nx];
  546.         exp = Floor[Log[$base,N[absnx]]]-$digits+1;
  547.         tmp = $base^exp;
  548.         absnx /= tmp;
  549.         mantissa = round[absnx];
  550.         If[mantissa > $maxman,
  551.             exp++;
  552.             mantissa = round[absnx/$base]
  553.         ];
  554.         ok = ($minman <= mantissa <= $maxman);
  555.         If[ok,
  556.             If[$minexp > exp,
  557.                 Message[ComputerNumber::undflw, exp + $digits];
  558.                 ok = False];
  559.             If[exp > $maxexp,
  560.                 Message[ComputerNumber::ovrflw, exp + $digits];
  561.                 ok = False];
  562.             ];
  563.         If[ok,
  564.             nx = SetPrecision[nx, $prec];
  565.             tmp = Sign[nx] mantissa $base^exp;
  566.             ComputerNumber[Sign[nx], mantissa, exp, tmp, nx],
  567.             (* else *)
  568.             NaN
  569.         ]
  570.     ]
  571.  
  572. Format[ComputerNumber[s_, m_, e_, v_, x_]] ^:= BaseForm[N[v, $nfdigits], $base]
  573.  
  574. Normal[ComputerNumber[s_, m_, e_, v_, x_]] ^:= v
  575.  
  576. NaN /: Abs[NaN] = NaN;
  577. NaN /: Less[NaN, _] := False
  578. NaN /: Less[_, NaN] := False
  579. NaN /: LessEqual[NaN, _] := False
  580. NaN /: LessEqual[_, NaN] := False
  581. NaN /: Greater[NaN, _] := False
  582. NaN /: GreaterEqual[_, NaN] := False
  583. NaN /: GreaterEqual[NaN, _] := False
  584. NaN /: Greater[_, NaN] := False
  585. NaN /: Equal[NaN, _] := False
  586. NaN /: Equal[_, NaN] := False
  587. NaN /: Unequal[NaN, _] := True
  588. NaN /: Unequal[_, NaN] := True
  589. NaN /: Plus[NaN, ___] := NaN
  590. NaN /: Times[NaN, ___] := NaN
  591. NaN /: Power[NaN, _] := NaN
  592. NaN /: Power[_, NaN] := NaN
  593. NaN /: IdealDivide[NaN, _] := NaN
  594. NaN /: IdealDivide[_, NaN] := NaN
  595.  
  596. ComputerNumber/:
  597.     Abs[x_ComputerNumber] :=
  598.     ComputerNumber[1, x[[2]], x[[3]], Abs[x[[4]]], Abs[x[[5]]]];
  599.  
  600. ComputerNumber/:
  601.     Less[x_ComputerNumber, y_ComputerNumber] := x[[4]] < y[[4]]
  602. ComputerNumber/:
  603.     LessEqual[x_ComputerNumber, y_ComputerNumber] := x[[4]] <= y[[4]]
  604. ComputerNumber/:
  605.     Greater[x_ComputerNumber, y_ComputerNumber] := x[[4]] > y[[4]]
  606. ComputerNumber/:
  607.     GreaterEqual[x_ComputerNumber, y_ComputerNumber] := x[[4]] >= y[[4]]
  608. ComputerNumber/:
  609.     Equal[x_ComputerNumber, y_ComputerNumber] := x[[4]] == y[[4]]
  610. ComputerNumber/:
  611.     Unequal[x_ComputerNumber, y_ComputerNumber] := x[[4]] != y[[4]]
  612.  
  613. ComputerNumber/:
  614. Plus[x_ComputerNumber, y_ComputerNumber] :=
  615.     Module[{tmp},
  616.         tmp = ComputerNumber[x[[4]] + y[[4]]];
  617.         If[Head[tmp] === ComputerNumber,
  618.             tmp[[5]] = x[[5]] + y[[5]],
  619.             (* else *)
  620.             tmp = NaN
  621.         ];
  622.         tmp
  623.     ]
  624.  
  625. ComputerNumber/:
  626. Times[-1, ComputerNumber[s_, m_, e_, v_, x_]] :=
  627.     ComputerNumber[-s, m, e, -v, -x];
  628.  
  629. ComputerNumber/:
  630. Times[x_ComputerNumber, y_ComputerNumber] :=
  631.     Module[{tmp},
  632.         tmp = ComputerNumber[x[[4]] y[[4]]];
  633.         If[Head[tmp] === ComputerNumber,
  634.             tmp[[5]] = x[[5]] y[[5]],
  635.             (* else *)
  636.             tmp = NaN
  637.         ];
  638.         tmp
  639.     ]
  640.  
  641. ComputerNumber::divzer = "Division by 0 occured in computation."
  642.  
  643. IdealDivide[x_ComputerNumber, y_ComputerNumber] :=
  644.     Module[{tmp},
  645.         If[y[[2]] == 0,
  646.             Message[ComputerNumber::divzer];
  647.             Return[NaN]];
  648.         tmp = ComputerNumber[x[[4]]/y[[4]]];
  649.         If[Head[tmp] === ComputerNumber,
  650.             tmp[[5]] = x[[5]]/y[[5]],
  651.             (* else *)
  652.             tmp = NaN
  653.         ];
  654.         tmp
  655.     ]
  656.  
  657. ComputerNumber/:
  658. Power[x_ComputerNumber, (n_Integer | n_Rational)] :=
  659.     Module[{tmp},
  660.         If[(x[[2]] == 0) && (n <= 0),
  661.             Message[ComputerNumber::divzer];
  662.             Return[NaN]];
  663.         tmp = ComputerNumber[x[[4]]^n];
  664.         If[Head[tmp] === ComputerNumber,
  665.             tmp[[5]] = x[[5]]^n,
  666.             (* else *)
  667.             tmp = NaN
  668.         ];
  669.         tmp
  670.     ]
  671.  
  672. ComputerNumber/:
  673. Power[x_ComputerNumber, y_ComputerNumber] :=
  674.     Module[{tmp},
  675.         If[y[[2]] == 0, Return[ComputerNumber[1]]];
  676.         tmp = ComputerNumber[x[[4]]^y[[4]]];
  677.         If[Head[tmp] === ComputerNumber,
  678.             tmp[[5]] = x[[5]]^y[[5]],
  679.             (* else *)
  680.             tmp = NaN
  681.         ];
  682.         tmp
  683.     ]
  684.  
  685. (* the following are only active under mixed-mode arithmetic. *)
  686.  
  687. ComputerNumber/:
  688. Times[(n_Integer | n_Rational), y_ComputerNumber] :=
  689.     Module[{tmp},
  690.         tmp = ComputerNumber[n y[[4]]];
  691.         If[Head[tmp] === ComputerNumber,
  692.             tmp[[5]] = n y[[5]],
  693.             (* else *)
  694.             tmp = NaN
  695.         ];
  696.         tmp
  697.     ] /; $mixedmode
  698.  
  699. IdealDivide[(n_Integer | n_Rational), y_ComputerNumber] :=
  700.     Module[{tmp},
  701.         If[y[[2]] == 0,
  702.             Message[ComputerNumber::divzer];
  703.             Return[NaN]];
  704.         tmp = ComputerNumber[n/y[[4]]];
  705.         If[Head[tmp] === ComputerNumber,
  706.             tmp[[5]] = n/y[[5]],
  707.             (* else *)
  708.             tmp = NaN
  709.         ];
  710.         tmp
  711.     ] /; $mixedmode
  712.  
  713. IdealDivide[x_ComputerNumber, (n_Integer | n_Rational)] :=
  714.     Module[{tmp},
  715.         If[n == 0,
  716.             Message[ComputerNumber::divzer];
  717.             Return[NaN]];
  718.         tmp = ComputerNumber[x[[4]]/n];
  719.         If[Head[tmp] === ComputerNumber,
  720.             tmp[[5]] = x[[5]]/n,
  721.             (* else *)
  722.             tmp = NaN
  723.         ];
  724.         tmp
  725.     ] /; $mixedmode
  726.  
  727. ComputerNumber/:
  728. Plus[x_ComputerNumber, y_?NumberQ] := x + ComputerNumber[y] /; $mixedmode
  729.  
  730. ComputerNumber/:
  731. Times[x_ComputerNumber, y_?NumberQ] := x ComputerNumber[y] /; $mixedmode
  732.  
  733. IdealDivide[x_ComputerNumber, y_?NumberQ] :=
  734.     IdealDivide[x,ComputerNumber[y]] /; $mixedmode
  735.  
  736. IdealDivide[x_?NumberQ, y_ComputerNumber] :=
  737.     IdealDivide[ComputerNumber[x],y] /; $mixedmode
  738.  
  739. IdealDivide[x_, y_] := Divide[x, y];    (* all other cases *)
  740.  
  741. ComputerNumber/:
  742. Power[x_ComputerNumber, y_?NumberQ] := x^ComputerNumber[y] /; $mixedmode
  743.  
  744. ComputerNumber/:
  745. Power[x_?NumberQ, y_ComputerNumber] := ComputerNumber[x]^y /; $mixedmode
  746.  
  747. End[ ] (* "NumericalMath`ComputerArithmetic`Private`" *)
  748.  
  749. Protect[SetArithmetic, Arithmetic, ComputerNumber, NaN, IdealDivide,
  750.     RoundToEven, RoundToInfinity, Truncation, RoundingRule, ExponentRange];
  751.  
  752. EndPackage[ ] (* "NumericalMath`ComputerArithmetic`" *)
  753.