Creating a conversion of indefinite integrals (antiderivatives) of real functions

[binghe]

Following some recent discussions on Slack, I'm considering writing a HOL conversion for computing indefinite integrals (antiderivatives) of real functions, which may be called INTEGRATE_CONV.

After some investigations on the underlying algorithms used by CAS software, I found a direct porting of these algorithm into HOL (as ML functions, with intermediate values as HOL theorems) is infeasible: not only these algorithms are extremely complex but also sometime they requires theories beyond HOL's current status.

A feasible approach which can be implemented in relatively short time is to use CAS software as external oracles. For an input function f(x), the CAS software will return another function F(x). Then, on the HOL side, all we need to do is to verify if D(F) - f = 0 is true, and this can be done by simplifying D(F) - f using existing conversions provided by realLib. If transcendental functions are involved and the existing simplifier is too weak, we can improve the simplifier.

(Note that the existing Diff.DIFF_CONV does not do any simplifications on its outputs. But users of DIFF_CONV should have been used to further simplify the outputs on their own.)

If we use the following FTC theorem as the working basis:

   [FUNDAMENTAL_THEOREM_OF_CALCULUS]  Theorem      
      ⊢ ∀f f' a b.
          a ≤ b ∧
          (∀x. x ∈ interval [(a,b)] ⇒
               (f has_vector_derivative f' x)
                 (at x within interval [(a,b)])) ⇒
          (f' has_integral f b − f a) (interval [(a,b)])

Let f be the input function and F be the desired antiderivative (from external oracle), a call of INTEGRATE_CONV ``\x. f`` should output a theorem like this:

      ⊢ ∀a b.
          a ≤ b ∧ (∀x. x ∈ interval [(a,b)] ⇒
          (f has_integral F b − F a) (interval [(a,b)])

or just (maybe)

      ⊢  ∀x. (F has_vector_derivative f x) (at x)

Below I describe how a CAS software called "Axiom" [1] is used as an external oracle. This software is quite strong in doing symbolic integration as it implements a decision procedure: for any so-called "elementary function" (NOTE: all transcendental functions are elementary), if the software cannot give the result, then it means the antiderivative function does NOT exist in finite form.

To use Axiom as an oracle from HOL, on HOL side the input function should be written into a temporarily file, then Axiom command is called, then the output is expected in another temporarily file, in an acceptable format.

For example, once Axiom is installed, I can put the following file as my $HOME/.axiom.input: (which comes from Axiom distribution)

--Copyright The Numerical Algorithms Group Limited 1994.

-- This file is read on start-up by the AXIOM system.  It is
-- most commonly used to customize the user's environment via the
-- )set system commands.
 
-- do not enter lisp break on error
)set break nobreak
 
-- use the simple character set
)set output char plain

-- everyone uses this
)set userlevel development

-- make the highlighting characters go away
)set message high off

)set message prompt frame

-- read the user's synonym file for system commands, if it exists
)read "mysyns.input"  )ifthere )quiet

The idea is that, Axiom will try to read mysyns.input in current directory as the input, if exists. Suppose the input function is x pow 2, on HOL side this function should be converted into Axiom's format, x ^ 2 or x ** 2, and be written into a file mysyns.input together with other Axiom commands:

)set output algebra off
)set output fortran on
)set output fortran output
)set fortran optlevel 0
)set fortran segment off
)set fortran ints2floats off
)set fortran precision single
)set quit unprotected
)set message type off
integrate(x^2,x)
)quit

The commands in the above file set the output format to FORTRAN and output target to a file named output.sfort (the part .sfort is added by Axiom), and then call the integrate function on the input function. If I start Axiom in the same directory, with the following command:

$ axiom -noht -nogr -noclef -nox

I will see the following outputs on screen, which indicates that the result has been written into output.sfort:

                        AXIOM Computer Algebra System
                          Version: Axiom (May 2017)
               Timestamp: Thursday June 30, 2022 at 10:59:57 
-----------------------------------------------------------------------------
   Issue )copyright to view copyright notices.
   Issue )summary for a summary of useful system commands.
   Issue )quit to leave AXIOM and return to shell.
   Visit http://axiom-developer.org for more information
-----------------------------------------------------------------------------
 
   Re-reading interp.daase
   Re-reading operation.daase
   Re-reading category.daase
   Re-reading browse.daase
   FORTRAN output will be written to file 
      /Users/binghe/Lisp/axiom/output.sfort.

The contents of output.sfort is this: (note that there are six whitespaces before R1, this whole file is a valid FORTRAN (77?) program)

      R1=(1/3)*x**3

For a more complicated input function like tan(atan(x)/3), the output is the following: (note that the 2nd line continues the 1st line by having & after five whitespaces. This seems also FORTRAN-derived.)

      R1=(8*ALOG(3*TAN(ATAN(x)/3)**2-1)+(-3*TAN(ATAN(x)/3)**2)+18*x*TAN(
     &ATAN(x)/3))/18

which means

                  atan(x) 2             atan(x) 2           atan(x)
        8log(3tan(-------)  - 1) - 3tan(-------)  + 18x tan(-------)
                     3                     3                   3
   (1)  ------------------------------------------------------------
                                     18

At this step, either the FORTRAN-compatible expressions are parsed by some ML code to be converted to HOL terms, or at Axiom side we can inject some user library code to output function expressions in HOL's directly acceptable format.

Of course, another open source CAS software can be used instead, e.g. Maxima (what else could be here?). But from the above notes I guess there's no easy way to support a wide range of different CAS software because each of them has different input/output formats.

[1] http://www.axiom-developer.org

[binghe]

Alternative to FORTRAN output (based on unparse):

Nasser wrote: 
> 
> newbie here. 
> I can't figure how to tell fricas to output result of command in 1D plain 
> text, instead of the default 2D. 
> 
> I tried all combinations of the )set output commands, but nothing is 
> working for me. All I want, is that instead of 
> 
> integrate(sin(x)*cos(x),x) 
> 
> 2 
> cos(x) 
> - ----------- 
> 2 
> 
> is to get 
> 
> 
> - cos(x)^2/2 
> 
> Spend 30 minutes, reading docs, googling, and I can't find one example how 
> to do this. Closest I found is 
> 
> (7) -> )set output algebra off 
> (7) -> )set output fortran on 
> (7) -> integrate(sin(x)*cos(x),x) 
> 
> R7=-DCOS(x)**2/2 
> 
> But I do not want fortran format. I just want 1D format instead of 2D for 
> standard algebra. 
> 
> Please show me the complete commands to use, as I am very new to axiom and 
> fricas 
> and assume I know nothing about the system. 

Unfortunately, normal output routines always produce 2D output. 
Linear output can be obtained via 'unparse': 

(1) -> ii := integrate(sin(x)*cos(x),x) 

2 
cos(x) 
(1) - ------- 
2 
Type: Union(Expression(Integer),...) 
(2) -> unparse(ii::InputForm) 

(2) "((-1)*cos(x)^2)/2" 
Type: String 

(3) -> ((-1)*cos(x)^2)/2 

2 
cos(x) 
(3) - ------- 
2 
Type: Expression(Integer) 

Note: result of 'unparse' when given back to FriCAS is supposed 
to produce orignal value. To make sure that result can be 
parsed correctly 'unparse' may add extra parenthesis (sometimes 
unnecessary). 

-- 
Waldek Hebisch 

[binghe]

Updates: )set message type off can be used to disable the type printing of each output terms.