home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / APPROXIM.M next >
Encoding:
Text File  |  1992-07-29  |  38.6 KB  |  1,034 lines
  1.  
  2. (*:Name: NumericalMath`Approximations` *)
  3.  
  4. (*:Title: Approximation of Functions *)
  5.  
  6. (*:Author: Jerry B. Keiper *)
  7.  
  8. (*:Summary:
  9.     This package provides the following tools for approximating
  10.     differentiable functions:
  11.     1. RationalInterpolation - finds a rational approximation to
  12.        a function which has no error at specified abscissae.  The
  13.        abscissae may be specified explicitly or implicitly.
  14.     2. MiniMaxApproximation - finds a rational approximation to
  15.        a function.  The maximum of the relative error in the
  16.        approximation over the interval has the smallest possible value.
  17.     3. GeneralRationalInterpolation - the same as RationalInterpolation 
  18.        except that the function may be specified parametrically.  This
  19.        simplifies the finding of approximations to inverse functions
  20.        that are difficult to specify explicitly (e.g., the inverse
  21.        of Erf[x]).
  22.     4. GeneralMiniMaxApproximation - the same as MiniMaxApproximation
  23.        except that the function may be specified parametrically, and the
  24.        error being minimized need not be the relative error.
  25. *)
  26.  
  27. (*:Context: NumericalMath`Approximations` *)
  28.  
  29. (*:Package Version: 2.0 *)
  30.  
  31. (* :Copyright: Copyright 1990,  Wolfram Research, Inc.
  32.     Permission is hereby granted to modify and/or make copies of
  33.     this file for any purpose other than direct profit, or as part
  34.     of a commercial product, provided this copyright notice is left
  35.     intact.  Sale, other than for the cost of media, is prohibited.
  36.     
  37.     Permission is hereby granted to reproduce part or all of
  38.     this file, provided that the source is acknowledged.
  39. *)
  40.  
  41. (*:History:
  42.     Updated to 2.0 by Jerry B. Keiper, November, 1990.
  43.     Version 1.2 by Jerry B. Keiper, January, 1989.
  44. *)
  45.  
  46. (*:Keywords:
  47.     functional approximation, Chebyshev approximation,
  48.     rational approximation, minimax approximation
  49. *)
  50.  
  51. (*Sources:
  52.     Carl-Erik Froberg, Numerical Mathematics: Theory and Computer
  53.         Applications, Benjamin/Cummings, 1985, pp. 250-266
  54.  
  55.     A. Ralston & P. Rabinowitz, A First Course in Numerical Analysis
  56.         (2nd. ed.), McGraw-Hill, New York, 1978
  57. *)
  58.  
  59. (*:Mathematica Version: 2.0 *)
  60.  
  61. (*:Limitation:
  62.     Each of the tools will fail if you attempt to find an
  63.     approximation to an even or odd function with numerator and
  64.     denominator degrees that are impossible to achieve. For
  65.     example RationalInterpolation[Cos[x], {x, 3, 2}, {x, -1, 1},
  66.     WorkingPrecision -> 30] or RationalInterpolation[Sin[x],
  67.     {x, 3, 3}, {x, -1, 1}].  The denominator of such symmetric,
  68.     rational interpolations should be even, and the numerator
  69.     should be even or odd, as the function is even or odd.
  70.  
  71.     Each of the tools can be very sensitive to the form of the
  72.     function being approximated.  For example, if the function
  73.     has an irrational singularity (algebraic or transcendental),
  74.     MiniMaxApproximation[ ] will be rather ineffective,  However,
  75.     if you either subtract out or divide out such a singularity,
  76.     MiniMaxApproximation[ ] will be much more effective.  Thus
  77.     on the interval [1, 2], Integrate[1/Log[t], {t, x, 2}]
  78.     should be approximated, as the approximation to (Log[x-1] +
  79.     Integrate[1/Log[t], {t, x, 2}]) with Log[x-1] subtracted
  80.     back off.
  81.  
  82.     Can occasionally get lost in its search for a minimax
  83.     approximation.  With care and experience such problems can
  84.     usually be dealt with by solving a sequence of approximations
  85.     with parameters varying from "good" vales to "bad" values.
  86.     *)
  87.  
  88. (*:Discussion:
  89. RationalInterpolation and GeneralRationalInterpolation are for
  90. interpolation.  MiniMaxApproximation and GeneralMiniMaxApproximation
  91. are for approximating a function on an interval.  The "General" in
  92. GeneralRationalInterpolation and GeneralMiniMaxApproximation refers
  93. to the fact that the function need not be specified explicitly
  94. and/or that the error being minimized is not required to be the
  95. relative error.  Each of these functions has the option WorkingPrecision,
  96. which specifies the number of decimal digits to use in calculating
  97. the approximation.
  98.  
  99. MiniMaxApproximation uses RationalInterpolation to get an initial
  100. approximation, which is then improved.  RationalInterpolation (in
  101. the case of implicitly specified abscissae) and MiniMaxApproximation
  102. both use the option Bias, which is a number between -1 and 1 and
  103. which causes the abscissae to be chosen with a bias to the left or
  104. right of the interval.  Bias -> 0 causes the abscissae to be chosen
  105. symmetrically.
  106.  
  107. MiniMaxApproximation has several options in addition to WorkingPrecision and
  108. Bias:
  109.     Brake - {n1, n2} where both n1 and n2 are non-negative integers.
  110.         MiniMaxApproximation tries to follow a "path" from the initial
  111.         approximation provided by RationalInterpolation to the final
  112.         answer.  This "path" is not well marked, and if the approx-
  113.         imation changes too much from one iteration to the next,
  114.         MiniMaxApproximation can lose its way.  Brake provides a way to
  115.         restrict the changes.  n1 specifies how many of the changes are
  116.         to be restricted and n2 specifies how much of a restriction is
  117.         to be placed on the first change.  The restrictions decrease
  118.         quadratically to no restriction at the n1+1-th iteration.
  119.     MaxIterations - the maximum number of iterations allowed after the
  120.         changes are no longer restricted by Brake.
  121.     Derivatives - {func, D[func, x], D[func, x, 2]}.  This option may be
  122.         left Automatic if Mathematica can find the derivatives itself.
  123.         However, it may be more efficient to specify a function that
  124.         returns the required list of values for any x in the interval.
  125.     PrintFlag - if MiniMaxApproximation cannot be made to work (or the
  126.         user would simply like to watch the convergence) setting
  127.         PrintFlag -> True will cause certain information to be printed
  128.         as it iterates.  More precisely, it prints the ordered pairs
  129.         consisting of the abscissae of the extrema of the error and
  130.         the value of the error at those extrema.
  131.     PlotFlag - setting this to True causes a plot of the error to be
  132.         made after each iteration.
  133.  
  134. It should be noted that in minimizing the maximum of the error,
  135. MiniMaxApproximation works with the relative error, i.e., the maximum of
  136. Abs[1 - approx[x]/func[x]] is minimized.  It will not normally be possible to
  137. approximate a function with a zero in the interval, since the relative error
  138. near the zero will approach infinity. More general forms of the error to be
  139. minimized require the use of GeneralMiniMaxApproximation.
  140.  
  141. GeneralRationalInterpolation and GeneralMiniMaxApproximation differ from
  142. RationalInterpolation and MiniMaxApproximation in that the function is
  143. specified parametrically. E.g., to approximate the inverse of the Gamma
  144. function on some interval, we might specify the function as {Gamma[t], t}
  145. or {Log[Gamma[t]], t}.  The latter is probably better, since the derivatives
  146. of Log[Gamma[t]] are simpler to compute and the domain for which the
  147. approximation is to be valid will be much smaller.  Another difference is
  148. that the error that is minimized is 1 - approx[x[t]]/y[t], or as an 
  149. alternative (y[t] - approx[x[t]])/g[t], where g is a third function
  150. specified by the user.
  151. *)
  152.  
  153. (*:Examples:
  154. ri1 = RationalInterpolation[Exp[x],{x,2,2},{-1,-1/2,0,1/2,1}]
  155. Plot[ri1 - Exp[x],{x,-2,2}]
  156.  
  157. ri2 = RationalInterpolation[Exp[x],{x,2,2},{x,-1,1}]
  158. Plot[ri2 - Exp[x],{x,-2,2}]
  159.  
  160. (* Note the effect of changing the option Bias. *)
  161. ri3 = RationalInterpolation[Exp[x],{x,2,2},{x,-1,1}, Bias -> .2]
  162. Plot[ri3 - Exp[x],{x,-2,2}]
  163.  
  164. ri4 = RationalInterpolation[Exp[x],{x,2,2},{x,-1,1}, Bias -> -.7]
  165. Plot[ri4 - Exp[x],{x,-2,2}]
  166.  
  167. mma1 = MiniMaxApproximation[Exp[x],{x,{-1,1},4,1}][[2,1]]
  168. Plot[1 - mma1/Exp[x],{x,-1,1}]
  169.  
  170. mma2 = MiniMaxApproximation[Exp[x],{x,{-1,1},4,1},
  171.     PlotFlag -> True, PrintFlag -> True]
  172.  
  173. mma3 = MiniMaxApproximation[Exp[x],{x,{-1,1},4,1},
  174.     Bias -> -.09, Brake -> {0,0}, PrintFlag -> True]
  175.  
  176. (* If it cannot find the derivative of the function it fails. *)
  177. mma4 = MiniMaxApproximation[Exp[Abs[x]],{x,{1,3},2,2}]
  178.  
  179. (* But we know how to find the derivative and we can tell it how. *)
  180. mma5 = MiniMaxApproximation[Exp[Abs[x]],{x,{1,3},2,2},
  181.     Derivatives -> Module[{exp = Exp[x]}, {exp, exp, exp}]]
  182.  
  183. gr1 = GeneralRationalInterpolation[{Exp[t],t},{t,2,3},x,{-1,-1/2,0,1/2,1,3/2}]
  184. Plot[gr1 - Log[x],{x,Exp[-3/2],Exp[2]}]
  185.  
  186. gr2 = GeneralRationalInterpolation[{Exp[t],t},{t,2,3},x,{t,-1,3/2}]
  187. Plot[gr2 - Log[x],{x,Exp[-3/2],Exp[2]}]
  188.  
  189. gr3 = GeneralRationalInterpolation[{Exp[t],t},{t,2,3},x,{t,-1,3/2},Bias->.4]
  190. Plot[gr3 - Log[x],{x,Exp[-3/2],Exp[2]}]
  191.  
  192. gr4 = GeneralRationalInterpolation[{Exp[t],t},{t,2,3},x,{t,-1,3/2},Bias->-.3]
  193. Plot[gr4 - Log[x],{x,Exp[-3/2],Exp[2]}]
  194.  
  195. (* If the weight function vanishes on the interval, there is a problem. *)
  196. gmma1 = GeneralMiniMaxApproximation[{Exp[t],t},{t,{-1,1},3,0},x]
  197.  
  198. gmma2 = GeneralMiniMaxApproximation[{Exp[t],t},{t,{1,3},3,0},x][[2,1]]
  199. Plot[1 - gmma2/Log[x],{x,Exp[1],Exp[3]}]
  200.  
  201. (* This is how we get it to minimize the maximum absolute error. *)
  202. gmma3 = GeneralMiniMaxApproximation[{Exp[t],t,1},{t,{1,3},3,0},x][[2,1]]
  203. Plot[gmma3 - Log[x],{x,Exp[1],Exp[3]}]
  204.  
  205. Finally we give a more involved use of MiniMaxApproximation.  We desire an
  206. approximation to the complementary error function erfc[x] == Erf[x,Infinity]
  207. on the interval 1 < x < Infinity.
  208.  
  209. We know that 
  210.  
  211.         Sqrt[Pi] x Exp[x^2] erfc[x] ~ 1
  212.  
  213. (cf. e.g., Abramowitz and Stegun, Handbook of Mathematical Functions, Dover,
  214. eqn 7.1.23.)
  215.  
  216. If we can approximate the function x Exp[x^2] erfc[x] we can easily solve
  217. for erfc[x].  (Note that the numerator and the denominator must be of the
  218. same degree, since as x -> Infinity our function approaches neither 0 nor
  219. Infinity.)  We will not be able to go all the way to Infinity; in fact we
  220. can only go to a little more than 25 before MiniMaxApproximation complains
  221. that it has a zero in the denominator of its approximation:
  222.  
  223. prec = 65;
  224. sqrtpi = N[Sqrt[Pi],65];
  225. f =  x  Erf[x,Infinity] Exp[x^2];
  226. ffders = Module[{erf = Erf[x,Infinity], exp = Exp[x^2]},
  227.     {x erf exp,erf exp (1+2x^2) - 2x/sqrtpi,
  228.     -4(x^2+1)/sqrtpi + (4x^3+6x) erf exp}];
  229. era = MiniMaxApproximation[f, {x,{1,30},3,3}, 
  230.     Derivatives->ffders, WorkingPrecision->65, Brake->{40,30}, Bias->-.5];
  231.  
  232. MiniMaxApproximation::zeroden: Denominator of the approximation ... has a
  233.     zero at 1.00456
  234.  
  235. (* The interval is too big; we shorten it somewhat. *)
  236.  
  237. era = MiniMaxApproximation[f, {x,{1,25},3,3}, 
  238.     Derivatives->ffders, WorkingPrecision->65, Brake->{40,30}, Bias->-.5];
  239.  
  240. (* does find an approximation for the interval 1 < x < 25.  We now proceed to
  241. stretch the interval to 1 < x < 70: *)
  242.  
  243. stretch[x1_] := Module[{},
  244.             era[[1,-1]] = x1;
  245.             era = MiniMaxApproximation[f, era, {x,{1,x1},3,3},
  246.                 Derivatives->ffders, WorkingPrecision->25]
  247.             ]
  248.  
  249. stretch[27];
  250. stretch[30];
  251. stretch[33];
  252. stretch[37];
  253. stretch[42];
  254. stretch[58];
  255. stretch[55];
  256. stretch[63];
  257. stretch[70]
  258.  
  259. gives the approximation
  260.  
  261. x Exp[x^2] erfc[x] = (-0.034999362734701385 + 1.18540580001459306*x +
  262.         1.07709037653126903*x^2 + 0.65975310784329749*x^3)/
  263.     (1 + 2.6734238081308433*x + 1.9096936117306481*x^2 +
  264.             1.1693737552523284*x^3)
  265.  
  266. The relative error on the interval 1 < x < 70 is less than .0000016.  On the
  267. interval 70 < x < Infinity the relative error can get as large as .000007.
  268. *)
  269.  
  270. BeginPackage["NumericalMath`Approximations`"]
  271.  
  272. RationalInterpolation::usage = 
  273. "RationalInterpolation[func, {x, m, k}, {x1, x2, ..., xmk1}, (opts)] gives the
  274. rational interpolant to func (a function of the variable x), where m and
  275. k are the degrees of the numerator and denominator, respectively and
  276. {x1, x2, ..., xmk1} is a list of m+k+1 abscissae of the interpolation points.
  277. An alternative form is \n
  278. RationalInterpolation[func, {x, m, k}, {x, x0, x1}, (opts)], which
  279. specifies the list of abscissae implicitly: the abscissae come from the
  280. interval (x0,x1)."
  281.  
  282. Options[RationalInterpolation] = {WorkingPrecision -> Precision[N[1]],
  283.                   Bias -> 0}
  284.  
  285. MiniMaxApproximation::usage = 
  286. "MiniMaxApproximation[func, {x, {x0, x1}, m, k}, (opts)] finds the mini-max
  287. approximation to func (a function of the variable x) on the interval (x0, x1),
  288. where m and k are the degrees of the numerator and denominator, respectively.
  289. The answer returned is {AbscissaList, {Approximation, MaxError}}, where
  290. AbscissaList is a list of the abscissae where the maximum error occurs,
  291. Approximation is the rational approximation desired, and MaxError is the
  292. value of the mini-max error.\n
  293. MiniMaxApproximation[f, approx, {x, {x0, x1}, m, k}, (opts)] is a form that
  294. allows the user to start the iteration from a known approximation.  Here
  295. approx must be in the form of an answer returned by MiniMaxApproximation."
  296.  
  297. Options[MiniMaxApproximation] = {
  298.             Bias -> 0,
  299.             Brake -> {5, 5},
  300.             Derivatives -> Automatic,
  301.             MaxIterations -> 20,
  302.             WorkingPrecision -> Precision[N[1]],
  303.             PrintFlag -> False,
  304.             PlotFlag -> False
  305.             }
  306.  
  307. GeneralRationalInterpolation::usage = 
  308. "GeneralRationalInterpolation[{x, y}, {t, m, k}, xvar, {t1, t2, ..., tmk1}],
  309. (opts)] gives the rational interpolant to the function y[x] on the interval
  310. t0 < t < t1, where x and y are functions of t which parametrically define y[x],
  311. m and k are the degrees of the numerator and denominator, respectively, xvar
  312. is the variable to be used for x in the approximation and {t1, t2, ..., tmk1}
  313. is a list of m+k+1 abscissae (values of t) of the interpolation points.
  314. An alternative form is \n
  315. GeneralRationalInterpolation[{x, y}, {t, m, k}, xvar, {t, t0, t1}, (opts)]
  316. which specifies the list of abscissae implicitly: the abscissae come from the
  317. interval t0 < t < t1."
  318.  
  319. Options[GeneralRationalInterpolation] = {WorkingPrecision -> Precision[N[1]],
  320.                   Bias -> 0}
  321.  
  322. GeneralMiniMaxApproximation::usage = 
  323. "GeneralMiniMaxApproximation[{x, y, (g)}, {t, {t0, t1}, m, k}, xvar, (opts)]
  324. finds the mini-max approximation to the function y[x] on the interval
  325. t0 < t < t1 where x and y are functions of t which parametrically define y[x],
  326. m and k are the degrees of the numerator and denominator, respectively and xvar
  327. is the variable to be used for x in the approximation.  The optional g is also
  328. a function of t, and the error to be minimized is (y[x[t]] - approx[x[t]])/g[t].
  329. The answer returned is {AbscissaList, {Approximation, MaxError}}, where
  330. AbscissaList is a list of the abscissae (values of t) where the maximum error
  331. occurs, Approximation is the rational approximation desired, and MaxError is
  332. the value of the mini-max error.\n
  333. GeneralMiniMaxApproximation[{x, y, (g)}, approx, {t, {t0, t1}, m, k}, xvar,
  334. (opts)] is a form that allows the user to start the iteration from a known
  335. approximation.  Here approx must be in the form of an answer returned by
  336. GeneralMiniMaxApproximation."
  337.  
  338. Options[GeneralMiniMaxApproximation] = {
  339.             Bias -> 0,
  340.             Brake -> {5, 5},
  341.             Derivatives -> Automatic,
  342.             MaxIterations -> 20,
  343.             WorkingPrecision -> Precision[N[1]],
  344.             PrintFlag -> False,
  345.             PlotFlag -> False
  346.             }
  347.  
  348. Bias::usage = "Bias is an option to RationalInterpolation,
  349. GeneralRationalInterpolation, MiniMaxApproximation and
  350. GeneralMiniMaxApproximation. It is a number between -1 and 1
  351. and causes the abscissae to be chosen with a bias to the left
  352. or right of the interval.  Bias -> 0 causes the abscissae to
  353. be chosen symmetrically."
  354.  
  355. Brake::usage = "Brake is an option to MiniMaxApproximation and
  356. GeneralMiniMaxApproximation.  A valid choice for Brake is a list of two
  357. nonnegative integers which control how the changes from one iteration to the
  358. next are to be restricted.  The first integer indicates the number of
  359. iterations that are to be restricted and the second indicates the magnitude
  360. of the first restriction.  The restrictions decrease to 0 as the algorithm 
  361. proceeds." 
  362.  
  363. Derivatives::usage = "Derivatives is an option to MiniMaxApproximation and
  364. GeneralMiniMaxApproximation.  It can be Automatic or an expression that
  365. evaluates to a list containing the function, the first derivative, and the
  366. second derivative evaluated at the variable."
  367.  
  368. Begin["NumericalMath`Approximations`Private`"]
  369.  
  370. RationalInterpolation[func_, {x_, m_Integer, k_Integer}, xlist_, opts___] :=
  371.     Module[{xinfo, bias, answer, biasOK = True,
  372.     prec = WorkingPrecision /.{opts} /.Options[RationalInterpolation]},
  373.     answer /;
  374.         (If[Head[xlist[[1]]] === Symbol,
  375.         xinfo = {x, SetPrecision[N[{xlist[[2]], xlist[[3]]},
  376.                     prec], prec], m, k};
  377.         bias = Bias /. {opts} /. Options[RationalInterpolation];
  378.         bias = SetPrecision[N[bias, prec],prec];
  379.         biasOK = (NumberQ[bias] && bias > -1 && bias < 1),
  380.               (* else *)
  381.         xinfo = {x, m, k};
  382.         bias = SetPrecision[N[xlist,prec], prec]
  383.         ];
  384.         (biasOK &&
  385.         (answer = RI[func, xinfo, prec, bias]; answer =!= Fail)))
  386.     ];
  387.  
  388. MiniMaxApproximation[f_, {x_,{x0_, x1_}, m_Integer, k_Integer}, opts___] :=
  389.     Module[{bias = Bias /. {opts} /. Options[MiniMaxApproximation],
  390.         brake = Brake /. {opts} /. Options[MiniMaxApproximation],
  391.         ffders = Derivatives /. {opts} /. Options[MiniMaxApproximation],
  392.         maxit = MaxIterations /. {opts} /. Options[MiniMaxApproximation],
  393.         prec = WorkingPrecision  /. {opts} /. Options[MiniMaxApproximation],
  394.         sflag = PrintFlag /. {opts} /. Options[MiniMaxApproximation],
  395.         pflag = PlotFlag /. {opts} /. Options[MiniMaxApproximation],
  396.     answer},
  397.     answer /; (If[SameQ[ffders,Automatic],
  398.                     ffders = D[f,x];
  399.                     ffders = {f, ffders, D[ffders,x]}
  400.                     ];
  401.             bias = SetPrecision[N[bias,prec], prec];
  402.             answer = MMA[f, ffders, {x, {x0, x1, bias}, m, k},
  403.             prec, maxit, brake, sflag, pflag];
  404.         Fail =!= answer)
  405.         ];
  406.  
  407. MiniMaxApproximation[f_, {xers_, erars_},
  408.                         {x_,{x0_, x1_}, m_Integer, k_Integer}, opts___] :=
  409.     Module[{brake = Brake /. {opts} /. Options[MiniMaxApproximation],
  410.         ffders = Derivatives /. {opts} /. Options[MiniMaxApproximation],
  411.         maxit = MaxIterations /. {opts} /. Options[MiniMaxApproximation],
  412.         prec = WorkingPrecision  /. {opts} /. Options[MiniMaxApproximation],
  413.         sflag = PrintFlag /. {opts} /. Options[MiniMaxApproximation],
  414.         pflag = PlotFlag /. {opts} /. Options[MiniMaxApproximation],
  415.     answer},
  416.     answer /; (If[SameQ[ffders,Automatic],
  417.                     ffders = D[f,x];
  418.                     ffders = {f, ffders, D[ffders,x]}
  419.                     ];
  420.             answer = MMA[f, ffders, {xers, erars}, {x, {x0, x1}, m, k},
  421.             prec, maxit, brake, sflag, pflag];
  422.         Fail =!= answer)
  423.         ];
  424.  
  425. GeneralRationalInterpolation[flist_, {t_, m_Integer, k_Integer},
  426.                         xvar_, tlist_, opts___] :=
  427.     Module[{tinfo, bias, answer,biasOK = True,
  428.         prec = WorkingPrecision /. {opts} /.
  429.                 Options[GeneralRationalInterpolation]},
  430.     answer /;
  431.         (If[SameQ[Head[tlist[[1]]],Symbol],
  432.         tinfo = {t, SetPrecision[N[{tlist[[2]], tlist[[3]]},
  433.                     prec], prec], m, k};
  434.         bias = Bias /. {opts} /. Options[GeneralRationalInterpolation];
  435.         bias = SetPrecision[N[bias, prec], prec];
  436.         biasOK = (NumberQ[bias] && bias > -1 && bias < 1),
  437.            (* else *)
  438.         tinfo = {t, m, k};
  439.         bias = SetPrecision[N[tlist, prec], prec];
  440.         ];
  441.         (biasOK &&
  442.         (answer = GRI[flist,tinfo,xvar,prec,bias]; Fail =!= answer)))
  443.     ];
  444.  
  445. GeneralMiniMaxApproximation[flist_, {t_, {t0_, t1_}, m_Integer, k_Integer},
  446.                             xvar_, opts___] :=
  447.     Module[{j,flen = Length[flist],
  448.     bias = Bias /. {opts} /. Options[GeneralMiniMaxApproximation],
  449.     brake = Brake /. {opts} /. Options[GeneralMiniMaxApproximation],
  450.     ffders = Derivatives /. {opts} /. Options[GeneralMiniMaxApproximation],
  451.     maxit = MaxIterations /. {opts} /.
  452.                     Options[GeneralMiniMaxApproximation],
  453.     prec = WorkingPrecision  /. {opts} /.
  454.                     Options[GeneralMiniMaxApproximation],
  455.     sflag = PrintFlag /. {opts} /. Options[GeneralMiniMaxApproximation],
  456.     pflag = PlotFlag /. {opts} /. Options[GeneralMiniMaxApproximation],
  457.     answer},
  458.     answer /; (If[SameQ[ffders,Automatic],
  459.             ffders = Array[1,{3,flen}];
  460.             Do[ffders[[1,j]] = flist[[j]], {j,flen}];
  461.             Do[ffders[[2,j]] = D[ffders[[1,j]], t], {j,flen}];
  462.             Do[ffders[[3,j]] = D[ffders[[2,j]], t], {j,flen}];
  463.             ffders = Transpose[ffders]
  464.             ];
  465.         bias = SetPrecision[N[bias,prec], prec];
  466.         answer = GMMA[flist, ffders, {t, {t0, t1, bias}, m, k},
  467.             xvar, prec, maxit, brake, sflag, pflag];
  468.         Fail =!= answer)
  469.     ];
  470.  
  471. GeneralMiniMaxApproximation[flist_, {ters_, erars_},
  472.         {t_, {t0_, t1_}, m_Integer, k_Integer}, xvar_, opts___] :=
  473.     Module[{j,flen = Length[flist],
  474.     bias = Bias /. {opts} /. Options[GeneralMiniMaxApproximation],
  475.     brake = Brake /. {opts} /. Options[GeneralMiniMaxApproximation],
  476.     ffders = Derivatives /. {opts} /. Options[GeneralMiniMaxApproximation],
  477.     maxit = MaxIterations /. {opts} /.
  478.                     Options[GeneralMiniMaxApproximation],
  479.     prec = WorkingPrecision  /. {opts} /.
  480.                     Options[GeneralMiniMaxApproximation],
  481.     sflag = PrintFlag /. {opts} /. Options[GeneralMiniMaxApproximation],
  482.     pflag = PlotFlag /. {opts} /. Options[GeneralMiniMaxApproximation],
  483.     answer},
  484.     answer /; (If[SameQ[ffders,Automatic],
  485.             ffders = Array[1,{3,flen}];
  486.             Do[ffders[[1,j]] = flist[[j]], {j,flen}];
  487.             Do[ffders[[2,j]] = D[ffders[[1,j]], t], {j,flen}];
  488.             Do[ffders[[3,j]] = D[ffders[[2,j]], t], {j,flen}];
  489.             ffders = Transpose[ffders]
  490.             ];
  491.         answer = GMMA[flist, ffders, {ters, erars},
  492.                 {t, {t0, t1}, m, k},
  493.                 xvar, prec, maxit, brake, sflag, pflag];
  494.         Fail =!= answer)
  495.     ];
  496.  
  497. RationalApproximation::precloss = "Drastic loss of precision.  Try starting
  498.     with higher precision."
  499.  
  500. RationalApproximation::notnum = "Function is not numerical at the interpolation
  501.     points: `1`."
  502.  
  503. pa[x_] := If[x == 0, Accuracy[x], Precision[x]];
  504. precacc[l_List] := Min[pa /@ l]
  505.  
  506. RI[f_, xinfo_, prec_, bias_] :=
  507.     Module[{i, mk1, xx, fx, mat, tempvec, x, x0, x1, m, k},
  508.     If[(IntegerQ[prec] && (prec > 0)) =!= True, Return[Fail]];
  509.     x = xinfo[[1]];
  510.     If[SameQ[Length[xinfo],4],
  511.         (x0 = xinfo[[2,1]];
  512.         x1 = xinfo[[2,2]];
  513.         m = xinfo[[3]];
  514.         k = xinfo[[4]];
  515.         mk1 = m+k+1;
  516.         xx = Reverse[Table[Cos[N[Pi (i+1/2)/(mk1),prec]], {i,0,m+k}]];
  517.         xx += bias (1 - xx xx);
  518.         xx = xx*(x1-x0)/2 + (x1+x0)/2),
  519.         (* else *)
  520.         (m = xinfo[[2]];
  521.         k = xinfo[[3]];
  522.         mk1 = m+k+1;
  523.         xx = bias)
  524.     ];
  525.     xx = SetPrecision[N[xx,prec],prec];
  526.     fx = SetPrecision[N[f /. x->xx, prec], prec];
  527.     If[!Apply[And,Map[NumberQ,fx]],
  528.         Message[RationalApproximation::notnum, fx];
  529.         Return[Fail]
  530.     ];
  531.     mat = Table[1,{i,mk1}];
  532.     tempvec = mat;
  533.     mat[[1]] = tempvec;
  534.     Do[tempvec *= xx;
  535.         If[i <= m,mat[[i+1]] = tempvec,,];
  536.         If[i <= k,mat[[i+m+1]] = -tempvec*fx],{i,Max[m,k]}];
  537.     xx = LinearSolve[Transpose[mat],fx];
  538.     If[Head[xx] === LinearSolve, Return[Fail]];
  539.     If[precacc[xx] < 5,
  540.         Message[RationalApproximation::precloss];
  541.         Return[Fail]
  542.     ];
  543.     xx = SetPrecision[xx, prec];
  544.     (xx[[1]]+Sum[xx[[i+1]] x^i,{i,m}])/(1+Sum[xx[[i+m+1]] x^i,{i,k}])
  545.     ];
  546.  
  547. MiniMaxApproximation::dervnotnum = "`1` is not numerical at `2` = `3` : `4`";
  548.  
  549. MiniMaxApproximation::van = "Failed to locate the extrema in `1` iterations.
  550.     The function `2` may be vanishing on the interval `3` or the
  551.     WorkingPrecision may be insufficient to get convergence.";
  552.  
  553. LocateExtrema[ffders_, era_, xinfo_, prec_, xguess_, maxit_,sflag_] :=
  554.     Module[{i, ii=1, mk1, xx, xd, fx, x, x0, x1, m, k, j=0,
  555.         rnum, rnumd, rnumdd, rden, rdend, rdendd, dxmax, oldd,
  556.         rnumx, rnum1, rnum2, rdenx, rden1, rden2, elist, r, repeat},
  557.     If[Head[era] === List,
  558.         r = era[[1]]; repeat = False,
  559.         r = era; repeat = True,
  560.     ];
  561.     x = xinfo[[1]];
  562.     x0 = xinfo[[2,1]];
  563.     x1 = xinfo[[2,2]];
  564.     m = xinfo[[3]];
  565.     k = xinfo[[4]];
  566.     mk1 = m+k+1;
  567.     If[SameQ[Head[xguess],List],xx=xguess,
  568.         xx=Reverse[Table[Cos[N[Pi i/mk1,prec]],{i,0,mk1}]];
  569.         xx += xguess (1 - xx xx);
  570.         xx = xx*(x1-x0)/2 + (x1+x0)/2
  571.     ];
  572.     If[Length[xx] != mk1+1, Return[Fail]];
  573.     xx = SetPrecision[N[xx,prec],prec];
  574.     xd = Take[xx,{2,mk1}];
  575.     dxmax = (rnum1 = xd-Drop[xx,-2];
  576.          rnum2 = Drop[xx,2]-xd;
  577.          Map[Min,Transpose[{rnum1,rnum2}]]/6);
  578.     oldd = Array[0&,{m+k}];
  579.     rnum = Numerator[r];
  580.     rnumd = D[rnum,x];
  581.     rnumdd = D[rnumd,x];
  582.     rden = Denominator[r];
  583.     rdend = D[rden,x];
  584.     rdendd = D[rdend,x];
  585.     While[ii < 5+prec/10,
  586.         If[j++ > 2 maxit,
  587.             Message[MiniMaxApproximation::van, maxit,
  588.                     ffders[[1]], {x0,x1}];
  589.             Return[Fail]
  590.         ];
  591.         fx = SetPrecision[N[ffders /. x->xd,prec],prec];
  592.         If[!Apply[And,Map[NumberQ,Flatten[fx]]],
  593.             Message[MiniMaxApproximation::dervnotnum,
  594.                 ffders, x, xd, fx];
  595.             Return[Fail]
  596.         ];
  597.         rnumx = SetPrecision[N[rnum /. x->xd,prec],prec];
  598.         rnum1 = SetPrecision[N[rnumd /. x->xd,prec],prec];
  599.         rnum2 = SetPrecision[N[rnumdd /. x->xd,prec],prec];
  600.         rdenx = SetPrecision[N[rden /. x->xd,prec],prec];
  601.         rden1 = SetPrecision[N[rdend /. x->xd,prec],prec];
  602.         rden2 = SetPrecision[N[rdendd /. x->xd,prec],prec];
  603.         r = rnum1 rdenx - rnumx rden1;
  604.         rden2 = (rdenx (rnum2 rdenx - rnumx rden2) - 
  605.                             2 r rden1)/(rdenx^3);
  606.         rden1 = r/(rdenx^2);
  607.         rdenx = rnumx/rdenx;
  608.         elist = 1-rdenx/fx[[1]];
  609.         rden1 = (rden1 fx[[1]] - rdenx fx[[2]]);
  610.         rden2 = rden1/(rden2 fx[[1]] - rdenx fx[[3]]);
  611.         rden1 = Sign[rden1] Sign[elist];
  612.         Do[If[Sign[rden1[[i]]] != Sign[rden2[[i]]] ||
  613.                     Abs[rden2[[i]]] > dxmax[[i]],
  614.             ii = 1;
  615.             rden2[[i]] = dxmax[[i]] Sign[rden1[[i]]],,],
  616.             {i,m+k}];
  617.         Do[If[rden2[[i]]+oldd[[i]] == 0,
  618.             rden2[[i]] /= 2;
  619.             dxmax[[i]] /= 2],
  620.             {i,m+k}];
  621.         oldd = rden2;
  622.         xd = xd - rden2;
  623.         xd = SetPrecision[xd,prec];
  624.         rden2 = N[Abs[rden2/xd]];
  625.         If[sflag, Print[rden2]];
  626.         If[repeat, If[Max[rden2] > .01, ii = 1, ii *= 2], ii = prec]
  627.     ];
  628.     xd = Prepend[xd,xx[[1]]];
  629.     xd = Append[xd,xx[[-1]]];
  630.     xd
  631.     ];
  632.  
  633. FindCoefficients[f_, era_, xinfo_, prec_, xlist_, brake_] :=
  634.     Module[{mk2, x, m, k, i, aa, bb, suma, sumb, abe, ff, temp1, temp2, jac,
  635.         r = If[SameQ[Head[era],List],era[[1]],era]},
  636.     x = xinfo[[1]];
  637.     m = xinfo[[3]];
  638.     k = xinfo[[4]];
  639.     mk2 = m+k+2;
  640.     suma = Numerator[r] /. x->xlist;
  641.     sumb = Denominator[r] /. x->xlist;
  642.     abe = Join[CoefficientList[Numerator[r],x],
  643.             Drop[CoefficientList[Denominator[r],x],1]];
  644.     abe = Append[abe,If[SameQ[Head[era],List],era[[2]],0]];
  645.     ff = SetPrecision[N[f /. x->xlist, prec],prec];
  646.     If[!Apply[And,Map[NumberQ,ff]],
  647.         Message[RationalApproximation::notnum, ff];
  648.         Return[Fail]];
  649.     temp1 = -1/(ff sumb);
  650.     temp2 = temp1;
  651.     jac = Range[m+k+2];
  652.     Do[jac[[i]] = temp2; temp2 *= xlist,{i,m+1}];
  653.     temp1 *= -suma;
  654.     temp2 = xlist temp1/sumb;
  655.     Do[jac[[i]] = temp2; temp2 *= xlist,{i,m+2,mk2-1}];
  656.     temp2 = Table[(-1)^i,{i,mk2}];
  657.     jac[[mk2]] = temp2;
  658.     temp2 = Table[1,{i,mk2}] - temp1 + abe[[mk2]] temp2;
  659.     temp2 = SetPrecision[temp2,prec];
  660.     temp2 = LinearSolve[Transpose[jac],temp2]/(brake+1);
  661.     If[Head[temp2] === LinearSolve, Return[Fail]];
  662.     abe = SetPrecision[abe-temp2,prec];
  663.     {Sum[abe[[i]] x^(i-1),{i,m+1}]/
  664.             (1+Sum[abe[[i]] x^(i-m-1),{i,m+2,mk2-1}]),
  665.         abe[[mk2]]}
  666.     ];
  667.  
  668.  
  669. MiniMaxApproximation::zeroden = "Denominator of the approximation `1` has
  670.     a zero at `2`.";
  671.  
  672. MMA[f_, ffders_, {x_,{x0_, x1_, bias_}, m_Integer, k_Integer},
  673.                 prec_, maxit_, brake_, sflag_, pflag_] :=
  674.     Module[{i=0, root, xe, era, xinfo={x, {x0, x1}, m, k}},
  675.     If[(IntegerQ[prec] && (prec > 0) && (bias < 1) && (bias > -1)
  676.             && IntegerQ[maxit] && (maxit > 0)) =!= True,
  677.         Return[Fail]];
  678.     era = RI[f, xinfo, prec, bias];
  679.     If[era === Fail, Return[Fail]];
  680.     If[k > 0,
  681.         xe = NRoots[N[Denominator[era]]==0,x];
  682.         While[i++ < k,
  683.             root = If[k==1, xe[[2]], xe[[i,2]]];
  684.             If[SameQ[Head[root],Real] && root<=x1 && root>=x0,
  685.                 Message[MiniMaxApproximation::zeroden,
  686.                     era, root];
  687.                 Return[Fail]
  688.             ]
  689.         ]
  690.     ];
  691.     xe = LocateExtrema[ffders,era,xinfo,prec,bias,maxit,sflag];
  692.     If[SameQ[xe, Fail], Return[Fail]];
  693.     If[!Apply[And,Map[NumberQ,xe]],
  694.         Message[MiniMaxApproximation::exnotnum, xe];
  695.         Return[Fail]
  696.     ];
  697.     MMA[f, ffders, {xe, era}, xinfo, prec, maxit, brake, sflag, pflag]
  698.     ];
  699.  
  700. MiniMaxApproximation::exnotnum = "The extrema are not numerical: `1`."
  701.  
  702. MiniMaxApproximation::extalt =
  703. "The extrema of the error do not alternate in sign.  It may be that 
  704. MiniMaxApproximation has lost track of the extrema by going too fast.  If so
  705. try increasing the values in the option Brake.  It may be that the
  706. WorkingPrecision is insufficient.  Otherwise there is an extra extreme value of
  707. the error, and MiniMaxApproximation cannot deal with this problem."
  708.  
  709. MiniMaxApproximation::extsort =
  710. "MiniMaxApproximation has lost track of the extrema of the error; they are no
  711. longer in sorted order.  Try increasing the values in the option Brake.";
  712.  
  713. MiniMaxApproximation::conv = "Warning: convergence was not complete.";
  714.  
  715. MMA[f_, ffders_, {xers_, erars_}, {x_,{x0_, x1_}, m_Integer, k_Integer},
  716.                 prec_, maxit_, brake_, sflag_, pflag_] :=
  717.     Module[{xe=SetPrecision[N[xers,prec],prec],
  718.         era=N[erars,prec], olderr=1, count=0, xinfo,
  719.         brvar=brake[[1]], elist, notdone=True, temp, tmp, alpha},
  720.     If[(IntegerQ[prec] && (prec > 0) && IntegerQ[maxit] &&
  721.             (maxit > 0)) =!= True, Return[Fail]
  722.     ];
  723.     xinfo = {x, SetPrecision[N[{x0, x1},prec],prec], m, k};
  724.     alpha = If[brake[[2]] == 0, brvar = 0; 1, brake[[1]]^2/brake[[2]]];
  725.     While[notdone,
  726.         elist = SetPrecision[N[f /. x->xe, prec],prec];
  727.         If[SameQ[List,Head[era]], tmp = era[[1]], tmp = era];
  728.         elist = {xe, 1 - (tmp /. x->xe)/elist};
  729.         temp = Sign[elist[[2]]];
  730.         temp = Drop[temp,1]+Drop[temp,-1];
  731.         If[Max[Abs[temp]]>0,
  732.             Message[MiniMaxApproximation::extalt];
  733.             Return[{xe, tmp, elist[[2]]}]
  734.         ];
  735.         If[sflag, Print[N[Transpose[elist]]]];
  736.         If[brvar != 0, count = 0; brvar--];
  737.         era = FindCoefficients[f,era,xinfo,prec,xe,brvar^2/alpha];
  738.         If[era === Fail, Return[Fail]];
  739.         If[pflag,Plot[1-era[[1]]/f,{x,x0,x1}]];
  740.         xe=LocateExtrema[ffders,era,xinfo,prec,xe,maxit,sflag];
  741.         If[SameQ[xe, Fail], Return[Fail]];
  742.         If[!Apply[And,Map[NumberQ,xe]],
  743.             Message[MiniMaxApproximation::exnotnum, xe];
  744.             Return[Fail]
  745.         ];
  746.         If[!SortedQ[xe],
  747.             Message[MiniMaxApproximation::extsort];
  748.             Return[Fail],,];
  749.         notdone = (brvar != 0) || (Abs[olderr-era[[2]]] >
  750.                 (Abs[olderr]+Abs[era[[2]]])/10000);
  751.         If[notdone && (count == maxit),
  752.             Message[MiniMaxApproximation::conv]
  753.         ];
  754.         notdone = ((count++ < maxit) && notdone);
  755.         olderr = era[[2]]
  756.     ];
  757.     {xe, era}
  758.     ];
  759.  
  760. GRI[f_, tinfo_, x_, prec_, bias_] :=
  761.     Module[{i, mk1, tt, ft, mat, tempvec, t, t0, t1, m, k},
  762.     If[(IntegerQ[prec] && (prec > 0)) =!= True, Return[Fail]];
  763.     t = tinfo[[1]];
  764.     If[SameQ[Length[tinfo],4],
  765.         (t0 = tinfo[[2,1]];
  766.         t1 = tinfo[[2,2]];
  767.         m = tinfo[[3]];
  768.         k = tinfo[[4]];
  769.         mk1 = m+k+1;
  770.         tt = Reverse[Table[Cos[N[Pi (i+1/2)/(mk1),prec]],{i,0,m+k}]];
  771.         tt += bias (1 - tt tt);
  772.         tt = tt*(t1-t0)/2 + (t1+t0)/2),
  773.         (* else *)
  774.         (m = tinfo[[2]];
  775.         k = tinfo[[3]];
  776.         mk1 = m+k+1;
  777.         tt = bias)
  778.     ];
  779.     tt = SetPrecision[N[tt,prec],prec];
  780.     ft = SetPrecision[N[f /. t->tt, prec], prec];
  781.     If[!Apply[And,Map[NumberQ,Flatten[ft]]],
  782.         Message[RationalApproximation::notnum, ft];
  783.         Return[Fail]
  784.     ];
  785.     tt = ft[[1]];    (* tt is now a list of x-values *)
  786.     ft = ft[[2]];    (* ft is now a list of y-values *)
  787.     mat = Table[1,{i,mk1}];
  788.     tempvec = mat;
  789.     mat[[1]] = tempvec;
  790.     Do[tempvec *= tt;
  791.         If[i <= m,mat[[i+1]] = tempvec];
  792.         If[i <= k,mat[[i+m+1]] = -tempvec*ft],
  793.         {i,Max[m,k]}
  794.     ];
  795.     tt = LinearSolve[Transpose[mat],ft];
  796.     If[Head[tt] === LinearSolve, Return[Fail]];
  797.     If[precacc[tt] < 5,
  798.         Message[RationalApproximation::precloss];
  799.         Return[Fail]
  800.     ];
  801.     tt = SetPrecision[tt, prec];
  802.     (tt[[1]]+Sum[tt[[i+1]] x^i,{i,m}])/(1+Sum[tt[[i+m+1]] x^i,{i,k}])
  803.     ];
  804.  
  805. GeneralMiniMaxApproximation::van = "Failed to locate the extrema in `1`
  806.     iterations.  The weight function `2` may be vanishing on the
  807.     interval `3`, or the WorkingPrecision may be insufficient to
  808.     get convergence.";
  809.  
  810. GeneralLocateExtrema[ffders_,era_,tinfo_,x_,prec_,tguess_,maxit_,sflag_] :=
  811.     Module[{i, ii = 1, mk1, tt, td, ft, gt, xd, t, t0, t1, m, k,
  812.         rnum, rnumd, rnumdd, rden, rdend, rdendd, dtmax, oldd, j=0,
  813.         rnumt, rnum1, rnum2, rdent, rden1, rden2, elist, r, repeat},
  814.     If[SameQ[Head[era],List],
  815.         r = era[[1]]; repeat = False,
  816.         r = era; repeat = True,
  817.     ];
  818.     t = tinfo[[1]];
  819.     t0 = tinfo[[2,1]];
  820.     t1 = tinfo[[2,2]];
  821.     m = tinfo[[3]];
  822.     k = tinfo[[4]];
  823.     mk1 = m+k+1;
  824.     If[SameQ[Head[tguess],List],tt=tguess,
  825.         tt=Reverse[Table[Cos[N[Pi i/mk1,prec]],{i,0,mk1}]];
  826.         tt += tguess (1 - tt tt);
  827.         tt = tt*(t1-t0)/2 + (t1+t0)/2
  828.     ];
  829.     If[Length[tt] != mk1+1,Return[Fail]];
  830.     tt = SetPrecision[N[tt,prec],prec];
  831.     td = Take[tt,{2,mk1}];
  832.     dtmax = (rnum1 = td-Drop[tt,-2];
  833.          rnum2 = Drop[tt,2]-td;
  834.          Map[Min,Transpose[{rnum1,rnum2}]]/6);
  835.     oldd = Array[0&,{m+k}];
  836.     rnum = Numerator[r];
  837.     rnumd = D[rnum,x];
  838.     rnumdd = D[rnumd,x];
  839.     rden = Denominator[r];
  840.     rdend = D[rden,x];
  841.     rdendd = D[rdend,x];
  842.     While[ii < 5+prec/10,
  843.         If[j++ > 2 maxit,
  844.             Message[GeneralMiniMaxApproximation::van,
  845.                 maxit,ffders[[-1,1]], {t0,t1}];
  846.             Return[Fail]
  847.         ];
  848.         ft = SetPrecision[N[ffders /. t->td,prec],prec];
  849.         If[!Apply[And,Map[NumberQ,Flatten[ft]]],
  850.             Message[MiniMaxApproximation::dervnotnum,
  851.                 ffders, t, td, ft];
  852.             Return[Fail]
  853.         ];
  854.         xd = ft[[1]];    (* xd is a list of x val's & der's *)
  855.         gt = ft[[-1]];    (* gt is a list of t val's & der's *)
  856.         ft = ft[[2]];    (* ft is a list of y val's & der's *)
  857.         rnumt = SetPrecision[N[rnum /. x->xd[[1]],prec],prec];
  858.         rnum1 = SetPrecision[N[rnumd /. x->xd[[1]],prec],prec];
  859.         rnum2 = SetPrecision[N[rnumdd /. x->xd[[1]],prec],prec];
  860.         rdent = SetPrecision[N[rden /. x->xd[[1]],prec],prec];
  861.         rden1 = SetPrecision[N[rdend /. x->xd[[1]],prec],prec];
  862.         rden2 = SetPrecision[N[rdendd /. x->xd[[1]],prec],prec];
  863.         r = rnum1 rdent - rnumt rden1;
  864.         rden2 = (rdent(rnum2 rdent - rnumt rden2)-2 r rden1)/(rdent^3);
  865.         rden1 = r/(rdent^2);
  866.         (* now we convert derivatives wrt x to der's wrt t *)
  867.         rden2 = rden2 (xd[[2]])^2 + rden1 xd[[3]];
  868.         rden1 *= xd[[2]];
  869.         rdent = rnumt/rdent;
  870.         elist = (ft[[1]]-rdent)/gt[[1]];
  871.         rden1 = (rden1 gt[[1]] - rdent gt[[2]]) -
  872.                 (ft[[2]] gt[[1]] - ft[[1]] gt[[2]]);
  873.         rden2 = rden1/((rden2 gt[[1]] - rdent gt[[3]]) -
  874.                 (ft[[3]] gt[[1]] - ft[[1]] gt[[3]]));
  875.         rden1 = Sign[rden1]  Sign[elist];
  876.         Do[If[Sign[rden1[[i]]] != Sign[rden2[[i]]] ||
  877.                 Abs[rden2[[i]]] > dtmax[[i]],
  878.             ii = 1;
  879.             rden2[[i]] = dtmax[[i]]*Sign[rden1[[i]]]
  880.             ],
  881.             {i,m+k}
  882.         ];
  883.         Do[If[rden2[[i]]+oldd[[i]] == 0,
  884.             rden2[[i]] /= 2;
  885.             dtmax[[i]] /= 2],
  886.             {i,m+k}
  887.         ];
  888.         oldd = rden2;
  889.         td = td - rden2;
  890.         td = SetPrecision[td,prec];
  891.         rden2 = N[Abs[rden2/td]];
  892.         If[sflag, Print[rden2]];
  893.         If[repeat, If[Max[rden2] > .01, ii = 1, ii *= 2], ii = prec]
  894.     ];
  895.     td = Prepend[td,tt[[1]]];
  896.     td = Append[td,tt[[-1]]];
  897.     td
  898.     ];
  899.  
  900. GeneralFindCoefficients[f_, era_, tinfo_, x_, prec_, tlist_, brake_] :=
  901.     Module[{mk2, t, m, k, i, aa, bb, suma, sumb, abe, ff,
  902.         xlist, temp1, temp2, jac, gt,
  903.         r = If[SameQ[Head[era],List],era[[1]],era]},
  904.     t = tinfo[[1]];
  905.     m = tinfo[[3]];
  906.     k = tinfo[[4]];
  907.     mk2 = m+k+2;
  908.     ff = SetPrecision[N[f /. t->tlist, prec],prec];
  909.     If[!Apply[And,Map[NumberQ,Flatten[ff]]],
  910.         Message[RationalApproximation::notnum, ff];
  911.         Return[Fail]
  912.     ];
  913.     xlist = ff[[1]];   (* xlist are the abscissas of the extrema *)
  914.     gt = ff[[-1]];     (* gt is the normalizing function *)
  915.     ff = If[Length[ff] == 2, 1, ff[[2]]/gt];
  916.                 (* ff are the normalized ordinates *)
  917.     abe = Join[CoefficientList[Numerator[r],x],
  918.             Drop[CoefficientList[Denominator[r],x],1]];
  919.     abe = Append[abe,If[SameQ[Head[era],List],era[[2]],0]];
  920.     suma = Numerator[r] /. x->xlist;
  921.     sumb = Denominator[r] /. x->xlist;
  922.     temp1 = -1/(gt sumb);
  923.     If[NumberQ[temp1], temp1 = Table[temp1,{i,m+k+2}]];
  924.     temp2 = temp1;
  925.     jac = Range[m+k+2];
  926.     Do[jac[[i]] = temp2; temp2 *= xlist,{i,m+1}];
  927.     temp1 *= -suma;
  928.     temp2 = xlist temp1/sumb;
  929.     Do[jac[[i]] = temp2; temp2 *= xlist,{i,m+2,mk2-1}];
  930.     temp2 = Table[(-1)^i,{i,mk2}];
  931.     jac[[mk2]] = temp2;
  932.     temp2 = SetPrecision[ff-temp1 + abe[[mk2]] temp2,prec];
  933.     temp2 = LinearSolve[Transpose[jac],temp2]/(brake+1);
  934.     If[Head[temp2] === LinearSolve, Return[Fail]];
  935.     abe = SetPrecision[abe-temp2,prec];
  936.     {Sum[abe[[i]] x^(i-1),{i,m+1}]/
  937.             (1+Sum[abe[[i]] x^(i-m-1),{i,m+2,mk2-1}]),
  938.         abe[[mk2]]}
  939.     ];
  940.  
  941. GMMA[f_, ffders_, {t_,{t0_, t1_, bias_}, m_Integer, k_Integer}, x_,
  942.                 prec_, maxit_, brake_, sflag_, pflag_] :=
  943.     Module[{i=0, root, xe, te, era, tinfo = {t, {t0, t1}, m, k}},
  944.     If[(IntegerQ[prec] && (prec > 0) && (bias < 1) && (bias > -1)
  945.             && IntegerQ[maxit] && (maxit > 0)) =!= True,
  946.         Return[Fail]
  947.     ];
  948.     era = GRI[f, tinfo, x, prec, bias];
  949.     If[era === Fail, Return[Fail]];
  950.     If[k > 0,
  951.         te = NRoots[N[Denominator[era]]==0,x];
  952.         xe = Sort[N[f /. t->{t0,t1}][[1]]];
  953.         While[i++ < k,
  954.             root = If[k==1, te[[2]], te[[i,2]]];
  955.             If[Head[root]===Real && root>=xe[[1]] && root<=xe[[2]],
  956.                 Message[MiniMaxApproximation::zeroden,
  957.                         era, root];
  958.                 Return[Fail]]
  959.             ]
  960.         ];
  961.     te = GeneralLocateExtrema[ffders,era,tinfo,x,prec,bias,maxit,sflag];
  962.     If[SameQ[te, Fail], Return[Fail]];
  963.     If[!Apply[And,Map[NumberQ,xe]],
  964.         Message[MiniMaxApproximation::exnotnum, xe];
  965.         Return[Fail]
  966.     ];
  967.     GMMA[f, ffders, {te, era}, tinfo, x, prec, maxit, brake, sflag, pflag]
  968.     ];
  969.  
  970. GMMA[f_, ffders_, {ters_, erars_}, {t_,{t0_, t1_}, m_Integer, k_Integer}, x_,
  971.                 prec_, maxit_, brake_, sflag_, pflag_] :=
  972.     Module[{te=SetPrecision[N[ters,prec],prec],
  973.         era=N[erars,prec], olderr=1, count=0, tinfo,
  974.         brvar=brake[[1]], elist, notdone=True, temp, tmp, alpha},
  975.     If[(IntegerQ[prec] && (prec > 0) && IntegerQ[maxit] &&
  976.             (maxit > 0)) =!= True, Return[Fail]];
  977.     tinfo = {t, SetPrecision[N[{t0, t1},prec],prec], m, k};
  978.     alpha = If[brake[[2]] == 0, brvar = 0; 1, brake[[1]]^2/brake[[2]]];
  979.     While[notdone,
  980.         elist = SetPrecision[N[f /. t->te,prec],prec];
  981.         If[SameQ[List,Head[era]], tmp = era[[1]], tmp = era];
  982.         elist[[2]] = (elist[[2]] - SetPrecision[
  983.             N[tmp /. x->elist[[1]],prec],prec])/elist[[-1]];
  984.         temp = Sign[elist[[2]]];
  985.         temp = Drop[temp,1]+Drop[temp,-1];
  986.         If[Max[Abs[temp]]>0,
  987.             Message[MiniMaxApproximation::extalt];
  988.             Return[{elist[[1]], tmp, elist[[2]]}]
  989.         ];
  990.         If[sflag, Print[N[Sort[Transpose[Take[elist,{1,2}]]]]]];
  991.         If[brvar !=0,count = 0;brvar--];
  992.         era = GeneralFindCoefficients[ f, era, tinfo, x, prec,
  993.                         te, brvar^2/alpha];
  994.         If[era === Fail, Return[Fail]];
  995.         If[pflag,
  996.             ParametricPlot[Evaluate[
  997.             Module[{ft = f},
  998.                 {ft[[1]],(ft[[2]]-(era[[1]] /. x->ft[[1]]))/
  999.                             ft[[-1]]}
  1000.             ]],
  1001.             {t,t0,t1}]
  1002.         ];
  1003.         te = GeneralLocateExtrema[
  1004.                 ffders, era, tinfo, x, prec, te, maxit, sflag];
  1005.         If[SameQ[te, Fail], Return[Fail]];
  1006.         If[!Apply[And,Map[NumberQ,xe]],
  1007.             Message[MiniMaxApproximation::exnotnum, xe];
  1008.             Return[Fail]
  1009.         ];
  1010.         If[!SortedQ[te],
  1011.             Message[MiniMaxApproximation::extsort];
  1012.             Return[Fail]
  1013.         ];
  1014.         notdone = (brvar != 0) || (Abs[olderr-era[[2]]] >
  1015.                     (Abs[olderr]+Abs[era[[2]]])/10000);
  1016.         If[notdone && (count == maxit),
  1017.                 Message[MiniMaxApproximation::conv]
  1018.         ];
  1019.         notdone = ((count++ < maxit) && notdone);
  1020.         olderr = era[[2]]
  1021.     ];
  1022.     {te, era}
  1023.     ];
  1024.  
  1025. End[]  (* NumericalMath`Approximations`Private` *)
  1026.  
  1027. Protect[RI, LocateExtrema, FindCoefficients, MMA, GRI, GeneralLocateExtrema,
  1028. GeneralFindCoefficients, GMMA, RationalInterpolation, MiniMaxApproximation,
  1029. GeneralRationalInterpolation, GeneralMiniMaxApproximation];
  1030.  
  1031. EndPackage[] (* NumericalMath`Approximations` *)
  1032.  
  1033.  
  1034.