laurent.spad

)abbrev category ULSCCAT UnivariateLaurentSeriesConstructorCategory
++ Author: Clifton J. Williamson
++ Date Created: 6 February 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: series, Laurent, Taylor
++ Examples:
++ References:
++ Description:
++   This is a category of univariate Laurent series constructed from
++   univariate Taylor series.  A Laurent series is represented by a pair
++   \spad{[n, f(x)]}, where n is an arbitrary integer and \spad{f(x)}
++   is a Taylor series.  This pair represents the Laurent series
++   \spad{x^n * f(x)}.
UnivariateLaurentSeriesConstructorCategory(Coef, UTS) : _
 Category == Definition where
  Coef : Ring
  UTS : UnivariateTaylorSeriesCategory Coef
  I ==> Integer

  Definition ==> Join(UnivariateLaurentSeriesCategory(Coef), _
                      RetractableTo UTS) with
    laurent : (I, UTS) -> %
      ++ \spad{laurent(n, f(x))} returns \spad{x^n * f(x)}.
    degree : % -> I
      ++ \spad{degree(f(x))} returns the degree of the lowest order term of
      ++ \spad{f(x)}, which may have zero as a coefficient.
    taylorRep : % -> UTS
      ++ \spad{taylorRep(f(x))} returns \spad{g(x)}, where
      ++ \spad{f = x^n * g(x)} is represented by \spad{[n, g(x)]}.
    removeZeroes : % -> %
      ++ \spad{removeZeroes(f(x))} removes leading zeroes from the
      ++ representation of the Laurent series \spad{f(x)}.
      ++ A Laurent series is represented by (1) an exponent and
      ++ (2) a Taylor series which may have leading zero coefficients.
      ++ When the Taylor series has a leading zero coefficient, the
      ++ 'leading zero' is removed from the Laurent series as follows:
      ++ the series is rewritten by increasing the exponent by 1 and
      ++ dividing the Taylor series by its variable.
      ++ Note: \spad{removeZeroes(f)} removes all leading zeroes from f
    removeZeroes : (I, %) -> %
      ++ \spad{removeZeroes(n, f(x))} removes up to n leading zeroes from
      ++ the Laurent series \spad{f(x)}.
      ++ A Laurent series is represented by (1) an exponent and
      ++ (2) a Taylor series which may have leading zero coefficients.
      ++ When the Taylor series has a leading zero coefficient, the
      ++ 'leading zero' is removed from the Laurent series as follows:
      ++ the series is rewritten by increasing the exponent by 1 and
      ++ dividing the Taylor series by its variable.
    coerce : UTS -> %
      ++ \spad{coerce(f(x))} converts the Taylor series \spad{f(x)} to a
      ++ Laurent series.
    taylor : % -> UTS
      ++ taylor(f(x)) converts the Laurent series f(x) to a Taylor series,
      ++ if possible.  Error: if this is not possible.
    taylorIfCan : % -> Union(UTS,"failed")
      ++ \spad{taylorIfCan(f(x))} converts the Laurent series \spad{f(x)}
      ++ to a Taylor series, if possible. If this is not possible,
      ++ "failed" is returned.
    if Coef has Field then QuotientFieldCategory(UTS)
      --++ the quotient field of univariate Taylor series over a field is
      --++ the field of Laurent series

   add

    zero? x == zero? taylorRep x
    retract(x : %) : UTS == taylor x
    retractIfCan(x:%):Union(UTS,"failed") == taylorIfCan x
    laurent(n: I, st: Stream Coef): % == laurent(n, series(st)$UTS)

)abbrev domain ULSCONS UnivariateLaurentSeriesConstructor
++ Authors: Bill Burge, Clifton J. Williamson
++ Date Created: August 1988
++ 14 June 1996: provided missing exquo: (%, %) -> % (Frederic Lehobey)
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: series, Laurent, Taylor
++ Examples:
++ References:
++ Description:
++   This package enables one to construct a univariate Laurent series
++   domain from a univariate Taylor series domain. Univariate
++   Laurent series are represented by a pair \spad{[n, f(x)]}, where n is
++   an arbitrary integer and \spad{f(x)} is a Taylor series.  This pair
++   represents the Laurent series \spad{x^n * f(x)}.
UnivariateLaurentSeriesConstructor(Coef, UTS) : _
 Exports == Implementation where
  Coef    : Ring
  UTS     : UnivariateTaylorSeriesCategory Coef
  I     ==> Integer
  L     ==> List
  NNI   ==> NonNegativeInteger
  OUT   ==> OutputForm
  P     ==> Polynomial Coef
  RF    ==> Fraction Polynomial Coef
  RN    ==> Fraction Integer
  ST    ==> Stream Coef
  TERM  ==> Record(k : I, c : Coef)
  monom ==> monomial$UTS
  EFULS ==> ElementaryFunctionsUnivariateLaurentSeries(Coef, UTS, %)
  STTAYLOR ==> StreamTaylorSeriesOperations Coef

  Exports ==> UnivariateLaurentSeriesConstructorCategory(Coef, UTS)

  Implementation ==> add

--% representation

    Rep := Record(expon : I, ps : UTS)

    getExpon : % -> I
    getUTS   : % -> UTS

    getExpon x == x.expon
    getUTS   x == x.ps

--% creation and destruction

    laurent(n: I, psr: UTS): % == [n, psr]
    taylorRep x    == getUTS x
    degree x       == getExpon x

    0 == laurent(0, 0)
    1 == laurent(0, 1)

    monomial(s, e) == laurent(e, s::UTS)

    coerce(uts : UTS) : % == laurent(0, uts)
    coerce(r : Coef) : %  == r :: UTS  :: %
    coerce(i : I) : %     == i :: Coef :: %

    taylorIfCan uls ==
      n := getExpon uls
      n < 0 =>
        uls := removeZeroes(-n, uls)
        getExpon(uls) < 0 => "failed"
        getUTS uls
      n = 0 => getUTS uls
      getUTS(uls) * monom(1, n :: NNI)

    taylor uls ==
      (uts := taylorIfCan uls) case "failed" =>
        error "taylor: Laurent series has a pole"
      uts

    termExpon : TERM -> I
    termExpon term == term.k
    termCoef : TERM -> Coef
    termCoef term == term.c
    rec : (I, Coef) -> TERM
    rec(exponent, coef) == [exponent, coef]

    recs : (ST, I) -> Stream TERM
    recs(st, n) == delay
      empty? st => empty()
      zero? (coef := frst st) => recs(rst st, n + 1)
      concat(rec(n, coef), recs(rst st, n + 1))

    terms x == recs(coefficients getUTS x, getExpon x)

    recsToCoefs : (Stream TERM, I) -> ST
    recsToCoefs(st, n) == delay
      empty? st => empty()
      term := frst st; ex := termExpon term
      n = ex => concat(termCoef term, recsToCoefs(rst st, n + 1))
      concat(0, recsToCoefs(st, n + 1))

    series st ==
      empty? st => 0
      ex := termExpon frst st
      laurent(ex, series recsToCoefs(st, ex))

--% normalizations

    removeZeroes x ==
      empty? coefficients(xUTS := getUTS x) => 0
      coefficient(xUTS, 0) = 0 =>
        removeZeroes laurent(getExpon(x) + 1, quoByVar xUTS)
      x

    removeZeroes(n, x) ==
      n <= 0 => x
      empty? coefficients(xUTS := getUTS x) => 0
      coefficient(xUTS, 0) = 0 =>
        removeZeroes(n - 1, laurent(getExpon(x) + 1, quoByVar xUTS))
      x

--% predicates

    x = y ==
      EQ(x, y)$Lisp => true
      (expDiff := getExpon(x) - getExpon(y)) = 0 =>
        getUTS(x) = getUTS(y)
      abs(expDiff) > _$streamCount$Lisp => false
      expDiff > 0 =>
        getUTS(x) * monom(1, expDiff :: NNI) = getUTS(y)
      getUTS(y) * monom(1, (- expDiff) :: NNI) = getUTS(x)

    pole? x ==
      (n := degree x) >= 0 => false
      x := removeZeroes(-n, x)
      degree x < 0

--% arithmetic

    x + y  ==
      n := getExpon(x) - getExpon(y)
      n >= 0 =>
        laurent(getExpon y, getUTS(y) + getUTS(x) * monom(1, n::NNI))
      laurent(getExpon x, getUTS(x) + getUTS(y) * monom(1, (-n)::NNI))

    x - y  ==
      n := getExpon(x) - getExpon(y)
      n >= 0 =>
        laurent(getExpon y, getUTS(x) * monom(1, n::NNI) - getUTS(y))
      laurent(getExpon x, getUTS(x) - getUTS(y) * monom(1, (-n)::NNI))

    x : % * y : % == laurent(getExpon x + getExpon y, getUTS x * getUTS y)

    x : % ^ n : NNI ==
      zero? n => 1
      laurent(n * getExpon(x), getUTS(x) ^ n)

    -- FIXME: should be lazy
    recip x ==
      x := removeZeroes(1000, x)
      zero? coefficient(x,d := degree x) => "failed"
      (uts := recip getUTS x) case "failed" => "failed"
      laurent(-d, uts)

    elt(uls1 : %, uls2 : %) ==
      (uts := taylorIfCan uls2) case "failed" =>
        error "elt: second argument must have positive order"
      uts2 := uts@UTS
      not zero? coefficient(uts2, 0) =>
        error "elt: second argument must have positive order"
      if (deg := getExpon uls1) < 0 then uls1 := removeZeroes(-deg, uls1)
      (deg := getExpon uls1) < 0 =>
        (recipr := recip(uts2 :: %)) case "failed" =>
          error "elt: second argument not invertible"
        uts1 := taylor(uls1 * monomial(1, -deg))
        (elt(uts1, uts2) :: %) * (recipr@%) ^ ((-deg) :: NNI)
      elt(taylor uls1, uts2) :: %

    eval(uls : %, r : Coef) ==
      if (n := getExpon uls) < 0 then uls := removeZeroes(-n, uls)
      uts := getUTS uls
      (n := getExpon uls) < 0 =>
        zero? r => error "eval: 0 raised to negative power"
        (recipr := recip r) case "failed" =>
          error "eval: non-unit raised to negative power"
        (recipr@Coef) ^ ((-n) :: NNI) *$STTAYLOR eval(uts, r)
      zero? n => eval(uts, r)
      r ^ (n :: NNI) *$STTAYLOR eval(uts, r)

--% values

    variable x == variable getUTS x
    center   x == center   getUTS x

    coefficient(x, n) ==
      a := n - getExpon(x)
      a >= 0 => coefficient(getUTS x, a :: NNI)
      0

    elt(x : %, n : I) == coefficient(x, n)

--% other functions

    order x == getExpon x + order getUTS x
    order(x, n) ==
      (m := n - (e := getExpon x)) < 0 => n
      e + order(getUTS x, m :: NNI)

    truncate(x, n) ==
      (m := n - (e := getExpon x)) < 0 => 0
      laurent(e, truncate(getUTS x, m :: NNI))

    truncate(x, n1, n2) ==
      if n2 < n1 then (n1, n2) := (n2, n1)
      (m1 := n1 - (e := getExpon x)) < 0 => truncate(x, n2)
      laurent(e, truncate(getUTS x, m1 :: NNI, (n2 - e) :: NNI))

    if Coef has IntegralDomain then
      rationalFunction(x, n) ==
        (m := n - (e := getExpon x)) < 0 => 0
        poly := polynomial(getUTS x, m :: NNI) :: RF
        zero? e => poly
        v := variable(x) :: RF; c := center(x) :: P :: RF
        positive? e => poly * (v - c) ^ (e :: NNI)
        poly / (v - c) ^ ((-e) :: NNI)

      rationalFunction(x, n1, n2) ==
        if n2 < n1 then (n1, n2) := (n2, n1)
        (m1 := n1 - (e := getExpon x)) < 0 => rationalFunction(x, n2)
        poly := polynomial(getUTS x, m1 :: NNI, (n2 - e) :: NNI) :: RF
        zero? e => poly
        v := variable(x) :: RF; c := center(x) :: P :: RF
        positive? e => poly * (v - c) ^ (e :: NNI)
        poly / (v - c) ^ ((-e) :: NNI)

      -- FIXME: This looks bogus due to eager search for zeros
      --
      x exquo y ==
        x := removeZeroes(1000, x)
        y := removeZeroes(1000, y)
        zero? coefficient(y, d := degree y) => "failed"
        (uts := (getUTS x) exquo (getUTS y)) case "failed" => "failed"
        laurent(degree x-d, uts)

    if Coef has coerce : Symbol -> Coef then
      if Coef has "^": (Coef,I) -> Coef then

        approximate(x, n) ==
          (m := n - (e := getExpon x)) < 0 => 0
          app := approximate(getUTS x, m :: NNI)
          zero? e => app
          app * ((variable(x) :: Coef) - center(x)) ^ e

    complete x == laurent(getExpon x, complete getUTS x)
    extend(x, n) ==
      e := getExpon x
      (m := n - e) < 0 => x
      laurent(e, extend(getUTS x, m :: NNI))

    map(f : Coef -> Coef, x : %) == laurent(getExpon x, map(f, getUTS x))

    multiplyCoefficients(f, x) ==
      e := getExpon x
      laurent(e, multiplyCoefficients((z1 : I) : Coef +-> f(e + z1), getUTS x))

    multiplyExponents(x, n) ==
      laurent(n * getExpon x, multiplyExponents(getUTS x, n))

    differentiate x ==
      e := getExpon x
      laurent(e - 1, multiplyCoefficients((z1 : I) : Coef +-> (e+z1)::Coef, getUTS x))

    if Coef has PartialDifferentialRing(Symbol) then
      differentiate(x : %, s : Symbol) ==
        (s = variable(x)) => differentiate x
        map((z1 : Coef) : Coef +-> differentiate(z1, s), x) - differentiate(center x, s)*differentiate(x)

    characteristic() == characteristic()$Coef

    if Coef has Field then

      retract(x : %) : UTS                      == taylor x
      retractIfCan(x:%):Union(UTS,"failed") == taylorIfCan x

      (x : %) ^ (n : I) ==
        zero? n => 1
        n > 0 => laurent(n * getExpon(x), getUTS(x) ^ (n :: NNI))
        xInv := inv x; minusN := (-n) :: NNI
        laurent(minusN * getExpon(xInv), getUTS(xInv) ^ minusN)

      (x : UTS) * (y : %) == (x :: %) * y
      (x : %) * (y : UTS) == x * (y :: %)

      inv x ==
        (xInv := recip x) case "failed" =>
          error "multiplicative inverse does not exist"
        xInv

      (x : %) / (y : %) ==
        (yInv := recip y) case "failed" =>
          error "inv: multiplicative inverse does not exist"
        x * (yInv@%)

      (x : UTS) / (y : UTS) == (x :: %) / (y :: %)

      numer x ==
        (n := degree x) >= 0 => taylor x
        x := removeZeroes(-n, x)
        (n := degree x) = 0 => taylor x
        getUTS x

      denom x ==
        (n := degree x) >= 0 => 1
        x := removeZeroes(-n, x)
        (n := degree x) = 0 => 1
        monom(1, (-n) :: NNI)

--% algebraic and transcendental functions

    if Coef has Algebra Fraction Integer then

      coerce(r : RN) == r :: Coef :: %

      if Coef has Field then
         (x : %) ^ (r : RN) == x ^$EFULS r

      exp x   == exp(x)$EFULS
      log x   == log(x)$EFULS
      sin x   == sin(x)$EFULS
      cos x   == cos(x)$EFULS
      tan x   == tan(x)$EFULS
      cot x   == cot(x)$EFULS
      sec x   == sec(x)$EFULS
      csc x   == csc(x)$EFULS
      asin x  == asin(x)$EFULS
      acos x  == acos(x)$EFULS
      atan x  == atan(x)$EFULS
      acot x  == acot(x)$EFULS
      asec x  == asec(x)$EFULS
      acsc x  == acsc(x)$EFULS
      sinh x  == sinh(x)$EFULS
      cosh x  == cosh(x)$EFULS
      tanh x  == tanh(x)$EFULS
      coth x  == coth(x)$EFULS
      sech x  == sech(x)$EFULS
      csch x  == csch(x)$EFULS
      asinh x == asinh(x)$EFULS
      acosh x == acosh(x)$EFULS
      atanh x == atanh(x)$EFULS
      acoth x == acoth(x)$EFULS
      asech x == asech(x)$EFULS
      acsch x == acsch(x)$EFULS

      ratInv : I -> Coef
      ratInv n ==
        zero? n => 1
        inv(n :: RN) :: Coef

      integrate x ==
        not zero? coefficient(x, -1) =>
          error "integrate: series has term of order -1"
        e := getExpon x
        laurent(e + 1, multiplyCoefficients((z1 : I) : Coef +-> ratInv(e + 1 + z1), getUTS x))

      if Coef has integrate : (Coef, Symbol) -> Coef and _
         Coef has variables : Coef -> List Symbol then
        integrate(x : %, s : Symbol) ==
          (s = variable(x)) => integrate x
          not entry?(s, variables center x) => map((z1 : Coef) : Coef +-> integrate(z1, s), x)
          error "integrate: center is a function of variable of integration"

      if Coef has TranscendentalFunctionCategory and _
         Coef has PrimitiveFunctionCategory and _
         Coef has AlgebraicallyClosedFunctionSpace Integer then

        integrateWithOneAnswer : (Coef, Symbol) -> Coef
        integrateWithOneAnswer(f, s) ==
          res := integrate(f, s)$FunctionSpaceIntegration(I, Coef)
          res case Coef => res
          first(res@List Coef)

        integrate(x : %, s : Symbol) ==
          (s = variable(x)) => integrate x
          not entry?(s, variables center x) =>
            map((z1 : Coef) : Coef +-> integrateWithOneAnswer(z1, s), x)
          error "integrate: center is a function of variable of integration"

    termOutput : (I, Coef, OUT) -> OUT
    termOutput(k, c, vv) ==
    -- creates a term c * vv ^ k
      k = 0 => c :: OUT
      mon :=
        k = 1 => vv
        vv ^ (k :: OUT)
      c = 1 => mon
      c = -1 => -mon
      (c :: OUT) * mon

    showAll? : () -> Boolean
    -- check a global Lisp variable
    showAll?() == _$streamsShowAll$Lisp

    termsToOutputForm : (I, ST, OUT) -> OUT
    termsToOutputForm(m, uu, xxx) ==
      l : L OUT := empty()
      empty? uu => (0$Coef) :: OUT
      n : NNI ; count : NNI := _$streamCount$Lisp
      for n in 0..count while not empty? uu repeat
        if frst(uu) ~= 0 then
          l := concat(termOutput((n :: I) + m, frst(uu), xxx), l)
        uu := rst uu
      if showAll?() then
        uu1 := uu
        for n in (count + 1).. while explicitEntries? uu and _
               not eq?(uu1, rst uu) repeat
          if frst(uu) ~= 0 then
            l := concat(termOutput((n::I) + m, frst(uu), xxx), l)
          if odd?(n) then uu1 := rst(uu1)
          uu := rst uu
      l :=
        explicitlyEmpty? uu => l
        eq?(uu, rst uu) and frst uu = 0 => l
        concat(prefix('O :: OUT, [xxx^((n::I) + m) :: OUT]), l)
      empty? l => (0$Coef) :: OUT
      reduce("+",reverse! l)

    coerce(x : %) : OUT ==
      x := removeZeroes(_$streamCount$Lisp, x)
      m := degree x
      uts := getUTS x
      p := coefficients uts
      var := variable uts; cen := center uts
      xxx :=
        zero? cen => var :: OUT
        paren(var :: OUT - cen :: OUT)
      termsToOutputForm(m, p, xxx)

)abbrev domain ULS UnivariateLaurentSeries
++ Author: Clifton J. Williamson
++ Date Created: 18 January 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: series, Laurent
++ Examples:
++ References:
++ Description: Dense Laurent series in one variable
++   \spadtype{UnivariateLaurentSeries} is a domain representing Laurent
++   series in one variable with coefficients in an arbitrary ring.  The
++   parameters of the type specify the coefficient ring, the power series
++   variable, and the center of the power series expansion.  For example,
++   \spad{UnivariateLaurentSeries(Integer, x, 3)} represents Laurent series in
++   \spad{(x - 3)} with integer coefficients.
UnivariateLaurentSeries(Coef, var, cen) : Exports == Implementation where
  Coef : Ring
  var  : Symbol
  cen  : Coef
  I   ==> Integer
  UTS ==> UnivariateTaylorSeries(Coef, var, cen)

  Exports ==> UnivariateLaurentSeriesConstructorCategory(Coef, UTS) with
    coerce : Variable(var) -> %
      ++ \spad{coerce(var)} converts the series variable \spad{var} into a
      ++ Laurent series.
    differentiate : (%, Variable(var)) -> %
      ++ \spad{differentiate(f(x), x)} returns the derivative of
      ++ \spad{f(x)} with respect to \spad{x}.
    if Coef has Algebra Fraction Integer then
      integrate : (%, Variable(var)) -> %
        ++ \spad{integrate(f(x))} returns an anti-derivative of the power
        ++ series \spad{f(x)} with constant coefficient 0.
        ++ We may integrate a series when we can divide coefficients
        ++ by integers.

  Implementation ==> UnivariateLaurentSeriesConstructor(Coef, UTS) add

    variable x == var
    center   x == cen

    coerce(v : Variable(var)) ==
      zero? cen => monomial(1, 1)
      monomial(1, 1) + monomial(cen, 0)

    differentiate(x : %, v : Variable(var)) == differentiate x

    if Coef has Algebra Fraction Integer then
      integrate(x : %, v : Variable(var)) == integrate x

)abbrev package ULS2 UnivariateLaurentSeriesFunctions2
++ Author: Clifton J. Williamson
++ Date Created: 5 March 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: Laurent series, map
++ Examples:
++ References:
++ Description: Mapping package for univariate Laurent series
++   This package allows one to apply a function to the coefficients of
++   a univariate Laurent series.
UnivariateLaurentSeriesFunctions2(Coef1, Coef2, var1, var2, cen1, cen2) : _
 Exports == Implementation where
  Coef1 : Ring
  Coef2 : Ring
  var1 : Symbol
  var2 : Symbol
  cen1 : Coef1
  cen2 : Coef2
  ULS1  ==> UnivariateLaurentSeries(Coef1, var1, cen1)
  ULS2  ==> UnivariateLaurentSeries(Coef2, var2, cen2)
  UTS1  ==> UnivariateTaylorSeries(Coef1, var1, cen1)
  UTS2  ==> UnivariateTaylorSeries(Coef2, var2, cen2)
  UTSF2 ==> UnivariateTaylorSeriesFunctions2(Coef1, Coef2, UTS1, UTS2)

  Exports ==> with
    map : (Coef1 -> Coef2, ULS1) -> ULS2
      ++ \spad{map(f, g(x))} applies the map f to the coefficients of the Laurent
      ++ series \spad{g(x)}.

  Implementation ==> add

    map(f, ups) == laurent(degree ups, map(f, taylorRep ups)$UTSF2)

--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.