Interval arithmetic and differential equations

Petr Girg, and Petr Hájek, Tomáš Koutný, Jiří Benedikt

Západočeská univerzita v Plzni, Pilsen, Czech Republic

Introduction

Mathematica supports wide range of methods for solving differential equations, both symbolical and numerical.

Symbolical methods are exact. Solutions are closed formulas, that can be later analyzed.  What you calculate with symbolic methods, is mathematicaly correct and can be used in rigorous mathematical proofs. On the other hand, symbolic methods are limited to relatively narrow class of differential equations.

Numerical methods can treat much wider class. On the other hand, it is very difficult to obtain error estimates. Thus usual numerical results can be used to make an idea what is going on, but  their validity in rigorous mathematical sense is questionable.

Example 1

Symbolical :

symbSol[c_]= y /. DSolve[y'[x]==y[x], y,x,GeneratedParameters->(c&)]

Manipulate[Plot[symbSol[c][[1]][x], {x,0,2}, PlotRange->{-1,Exp[2]}], {c,0,2}]

Numerical:

numSol[c_?NumericQ] := y /.NDSolve[{y'[x]==y[x], y[0]==c}, y,{x,0,2}]

numSol[1.]

Manipulate[Plot[Evaluate[numSol[c][[1]][x]], {x,0,2}, PlotRange->{-1,Exp[2]}], {c,0,1}]

Manipulate[Plot[Evaluate[numSol[c][[1]][x]- symbSol[c][[1]][x]], {x,0,2}, PlotRange->All], {c,0,1}]

Example 2

Symbolical

DSolve[y''[x]==-Sin[x]y[x], y,x]

Numerical

numSol2=y /. NDSolve[{y''[x]==-Sin[x] y[x], y[0]==0, y'[0]==1},y,{x,0,10}]

Plot[numSol2[[1]][x], {x,0,10}]

p - Laplacian / eigenvalue problem

1 < p, λ is real

Integrate[1/(1-s^p)^(1/p),{s,0,1}, Assumptions->{p>1}]

We define

In our library, we compute

Computer arithmetic

User classes
Absolutely optimistic user: every result obtained by computer is absolutely correct.
Somehow informed user: some results obtained by computer may be wrong, but such cases are very rare.
Absolutely pesimistic user: ever result obtained by computer is only approximate (therefore mathematicaly wrong), 
because computer uses approximate numbers.
Advanced user: some results are wrong some are correct. Advanced user is aware of the fact that 
arithmetics in computer obeys some rules and the correctness of the results is determined by these rules. Errors are not at random.
Expert user: knows the rules of computer arithmetics and can calculate the result by hand.

Another class (not only a user)
Chip developper: knows how the chip is hardwired and knows the possible exceptions from the rules that nobody else knows.

Here we shortly explain sources of errors in computer arithmetic.

Main source of errors:
Computer has finitely many states. Therefore it can store and calculate with finitely many numbers. 

a) Integers in a 32-bit computer: Integers between  -(2^31) to (2^31) - 1 

These numbers work exactly with respect to addition if no overflow occurs.

b) Reals in a modern computer, finite selection from the real numbers. So called float, double, extended double etc.

These finite sets are used for doing numerical mathematics. The numerical mathematics is used for evaluating continuous mathematical objects where inifinite (and uncountable) reals are needed. 

Concepts of mathematical analysis such as of limit, derivative, integral etc. are based on the fact that we use continuous real numbers (that have the completeness property).

Example of fixed precision arithmetic

In order to explain source of errors let us consider the following selection from the reals.
Let S be a set of numbers that can be written as ±(Subscript[a, 0]).(Subscript[a, -1])(Subscript[a, -2])(Subscript[a, -3])*(10^(Subscript[e, 1])) in a base ten. Thus the number
-12.34=-1.234*(10^1), 1.12*(10^9), 0.1234 = 1.234*(10^(-1))  are in our selection S, while 12345 or 0.12345 are not.

1.000*(10^0)  + 1.000*(10^0)=2.000*(10^0)    Exact result.

1.000*(10^0)+1.000*(10^(4 ))=1.0001*(10^4)  Result not representable in our system S.
Must be rounded to representable number:  1.000*(10^(4 ))   Inexact result !


9.999*(10^9) + 0.001*(10^9) = 10.000*(10^9)  Not representable, too large. Result represented by symbol +(∞^( )).

1.000*(10^0)/ 0.000*(10^0)   Not defined. Result represented by symbol NaN.
 
It is usefull to add symbols -(∞^( )), +∞, NaN

IEEE - 754

TraditionalForm[Hold[Sum[2^i, {i,0,10}]]]

Sum[2^i, {i,0,10}]==2^11-1

2^11-1

Sum[2^i, {i,0,10}]

IEEE - 754  double precision is a computer number format that collects all real numbers that can be written in the form

(((-1))^s) (2^((E-1023)))(1+  (Underoverscript[∑ , (i=1), 52])(Subscript[m, i])(2^(-i)) )

s∈{0,1},  (Subscript[m, i])∈{0,1},E =1, ...,(∑_(i=0)^10)(2^i)-1== (2^11)-2

and the number 0 , which is coded by (Subscript[m, i])=0 and E=0. The system contains also -0.

IEEE - 754 supports also symbols ++∞ and -∞ coded by E=(∑_(i=0)^10)(2^i) =(2^11)-1 and all (Subscript[m, i])=0.

IEEE - 754 supports also symbols NaN coded by E=(∑_(i=0)^10)(2^i) =(2^11)-1 and some (Subscript[m, i])≠ 0.

Storage of the double precision in the memory

Sign bit : 1 bit
Exponent width : 11 bits
Significand precision : 53 bits (52 explicitly stored)

Grid[{Join[{Item[ "",Background->Blue]},Table[Item[ "",Background->Green], {i,1,11}],Table[Item[ "",Background->Red], {i,1,52}]]},Frame->All]

Grid[Partition[Join[{Item[ "",Background->Blue]},Table[Item[ "",Background->Green], {i,1,11}],Table[Item[ "",Background->Red], {i,1,52}]],8], Frame->All]

Smallest positive double precision number:

MinMantisa=1

MinExponent=1-1023

MinMantisa*2^MinExponent

N[%]

MachineNumberQ[%]

MaxMantisa=Sum[2^-i, {i, 0, 52}]

N[MaxMantisa*2^((-1023))]

MachineNumberQ[%]

Largest positive double precision number:

MaxMantisa=Sum[2^-i, {i, 0, 52}]

MaxExponent=(2^11-2-1023)

MaxMantisa*2^MaxExponent

N[%]

MachineNumberQ[%]

N[2^(2^11-1-1023)]

MachineNumberQ[%]

Machine Epsilon

(((-1))^s) (2^((E-1023)))(1+  (Underoverscript[∑ , (i=1), 52])(Subscript[m, i])(2^(-i)) )

s∈{0,1},  (Subscript[m, i])∈{0,1},E =1, ..., (2^11)-2

The number 1 + (2^(-53)) cannot be written in the above form. Thus it must be rounded - source of errors.

Several possible roundings of 1 + (2^(-53)):
To nearest double:    1
Towards +∞: 1 + (2^(-52))
Towards -∞:  1

Smallest representable number ε>0, such that 1 + ε ≠ 1.

Machine Epsilon for the double precision: (2^(-52)).

Effects of Rounding - Machine Precision

1+2^(-52)

1.+2.^(-52) // FullForm

SetPrecision[1.+2.^(-52), Infinity]

1+2^(-53)

SetPrecision[1.+2^(-53), Infinity]

Effects of Rounding - Constant Array

a=ConstantArray[0.1,{10^6}];

10^5-Total[a]

%/10^5

a=ConstantArray[0.1,{10^7}];

10^6-Total[a]

%/10^6

Effects of Rounding - Constant Array

Table[{i, TraditionalForm[HoldForm[Sin[10^j Pi]]] /. j->i,Sin[10^i Pi],Sin[10.^i Pi]}, {i, 1,20}] // TableForm

Interval arithmetic

Interval arithmetic in Mathematica

Interval[{1,1}]+Interval[{2.^(-53), 2.^(-53)}]

SetPrecision[%, Infinity]

IntervalMemberQ[%, 1+2^(-53)]

(1-(9007199254740991/9007199254740992))*2^(53)

((4503599627370499/4503599627370496)-1)*2^(52)

Experiment with constant array.

a=ConstantArray[Interval[{0.1,0.1}],{10^6}];

Total[a] // FullForm

IntervalMemberQ[%, 10^5]

Experiment with the sine function.

Table[{i, Sin[10^i Pi],Sin[Interval[{10.^i Pi, 10.^i Pi}]]}, {i, 1,20}] // TableForm

Simple Interval Integrator

Integrate[Sin[x],{x,0,1}]

Integrate[Sin[x],{x,0,1}]-NIntegrate[Sin[x],{x,0,1}]

Sin[Interval[{0.1,0.1}]]

Total[Sin[#]*Interval[{1/1000., 1/1000.}]&/@Table[Interval[{i/1000,(i+1)/1000}],{i, 0.,999.}]]

IntervalMemberQ[%, Integrate[Sin[x],{x,0,1}]]

The interval results  contain the true value for the result. The Mathematica implementation works on all platforms, but the price being paid is too pesimistic estimate of the endpoint.

Interval[{1.,1.}]+Interval[{0.,0.}]

SetPrecision[%, Infinity]

Acces to Interval data structure is through high level commands only.

We need millions of operations for performing proofs  in differential  equations. Big number of operations requires a lot of execution time and pesimistic estimates may produce useless result. Thus we need to implement interval arithmetics as efficient as possible for our special purposes.

The processor (Intel^®) 8087 and compatible

According to “THE (INTEL^®) 8087 NUMERIC DATA PROCESSOR” by John Palmer (from Intel) the processor Intel^® 8087 was originaly designed to be used for interval computations. 

Thus we decided to implement our library for (INTEL^®) Pentium, the most common processor that contain FPU compatible with the THE (INTEL^®) 8087 NUMERIC DATA PROCESSOR.

Registers of (Intel^®) 8087

Grid[{{Item["Control Word (control bits only)"], SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft},{Item["?", Background-> LightRed],
         Item["?", Background-> LightRed],
         Item["?", Background-> LightRed],
         Item["IC", Background-> LightRed],
         Item["PC", Background-> LightRed], 
         Item["PC", Background-> LightRed],
         Item["RC", Background-> LightRed],
         Item["RC", Background-> LightRed]
}},ItemSize->{3,3}, ItemStyle->17,Frame->All]

IC : infinity control - there are two types of infinity arithmetic provided: 
affine and projective.

PC : precision control - results are rounded
to one of three precisions: 
Extended (64 bits), Double (53 bits), Float (24 bits).

RC : rounding control - results are rounded
in one of four directions: 
unbiased round to nearest, 
round towards +∞ ,
round towards -∞, 
round towards zero.

{Item["Control Word (control bits only)"], SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft}

Grid[{{Item["Register stack"]},{Item["ST0", Background-> LightRed],Item["1.", Background-> LightRed]},
        {Item["ST1", Background-> LightRed],Item["2.", Background-> LightRed]},
        {Item["ST2", Background-> LightRed],Item["?", Background-> LightRed]},
         {Item["ST3", Background-> LightRed],Item["?", Background-> LightRed]},
         {Item["ST4", Background-> LightRed],Item["?", Background-> LightRed]}, 
         {Item["ST5", Background-> LightRed],Item["?", Background-> LightRed]},
         {Item["ST6", Background-> LightRed],Item["?", Background-> LightRed]},
         {Item["ST7", Background-> LightRed],Item["?", Background-> LightRed]
}},ItemSize->{12,3}, ItemStyle->17,Frame->All]

One operand is at the top of the stack, the other operand can be at any register of the stack.  The result of the operation is rounded according the 
current state of the RC register.

Some basic operations are : abs,  addition, subtraction, multiplication, division, (2^x)-1 (-1≤ x≤1), y*(Subscript[Log, 2 ])(x),sine, cosine, arctan.

Intel Extended Precision

Sign bit : 1 bit
Exponent width : 15 bits
Significand precision : 64 bits (all 64 bits are explicitly stored)

Grid[{Join[{Item[ "",Background->Blue]},Table[Item[ "",Background->Green], {i,1,15}],Table[Item[ "",Background->Red], {i,1,64}]]},Frame->All]

Grid[Partition[Join[{Item[ "",Background->Blue]},Table[Item[ "",Background->Green], {i,1,15}],Table[Item[ "",Background->Red], {i,1,64}]],8], Frame->All]

Normals

(((-1))^s) (2^((E-16383)))(1+  (Underoverscript[∑ , (i=1), 63])(Subscript[m, i])(2^(-i)) )

s∈{0,1},  (Subscript[m, i])∈{0,1},E =1, ..., (2^15)-2

Subnormals

(((-1))^s) (2^(-16382))( (Underoverscript[∑ , (i=0), 64])(Subscript[m, i])(2^(-i)) )

s∈{0,1},  (Subscript[m, i])∈{0,1}, (Subscript[m, (Subscript[i, 0])])≠0 for some (Subscript[i, 0]),and  E =0.

The number 1 + (2^(-64)) cannot be written in the above form. Thus it must be rounded - source of errors.

Several possible roundings of 1 + (2^(-64)):
To nearest double:    1
Towards +∞: 1 + (2^(-63))
Towards -∞:  1

Smallest representable number ε>0, such that 1 + ε ≠ 1.

Machine Epsilon for the double precision: (2^(-52)).

offsetExp=2^14-1

Sum[2^i, {i,0,14}]-1-offsetExp

N[2^(1-offsetExp)]

maxExp=Sum[2^i, {i,0,14}]-1-offsetExp

maxMantisa=Sum[2^(-i), {i, 0, 64}]

N[maxMantisa*2^maxExp]

MachineNumberQ[%]

Machine Epsilon for Extended

N[2^(-63)]

N[2^(-52)/2^(-63)]

Computing in Extended = much more precision with almost no extra cost (8087 works in Extended internaly anyway).

Development of Library with Mathematica

Bad and good treatment of interval extensions of functions - Example

Normal[Series[ArcTan[x],{x,0,10}]]

Normal[Series[ArcTan[x],{x,0,12}]] /.{x->  Interval[{0.1, 0.8}], θ->Interval[{0., 2.}]}

Result is correct but the error is too overestimated.

Plot[ArcTan[x],{x, 0.,1}]

Hull[Interval[{a_, b_}], Interval[{c_, d_}]]:=Interval[{a,d}]

Hull[Interval[{0.,0.}], Interval[{1.,1.}]]

Normal[Series[ArcTan[x],{x,0,10}]] /.{x->  Interval[{0.1, 0.1}]}

Sum[q^(2k+1), {k, m, Infinity}]

% /. {q-> 1/2, m->10.}

error[q_, m_]:=-((q^(1+2 m))/(-1+(q^2)))/m

error[0.5, 10.]

lower=Normal[Series[ArcTan[x],{x,0,12}]]+Hull[error[x, 12.],-error[x, 12.]] /.{x->  Interval[{0.1, 0.1}]}

FullForm[%]

upper=Normal[Series[ArcTan[x],{x,0,12}]]+Hull[error[x, 12.],-error[x, 12.]] /.{x->  Interval[{0.9, 0.9}]}

Hull[lower, upper]

ArcTan[0.9]

Mathematica 8 Library Link

From the Mathematica Help: “Wolfram LibraryLink provides a powerful way to connect external code to Mathematica, enabling high-speed and memory-efficient execution. It does this by allowing dynamic libraries to be directly loaded into the Mathematica kernel so that functions in the libraries can be immediately called from Mathematica. 
     You can exchange not only C-like data types such as integers, reals, packed arrays, and strings, but also arbitrary Mathematica expressions. In addition, there are useful functions such as sending errors and calling back to Mathematica.”

Our Library

Needs["CCompilerDriver`"]

createTensor="
#include \"WolframLibrary.h\"

extern void funkcicka(int lengthD, double* D);

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
	return 0;
}

DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
	return;
}

DLLEXPORT int constantzero(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){
   MArgument_setInteger(Res, 0);
   return LIBRARY_NO_ERROR;
}


DLLEXPORT int integratorEulerIntervalyC(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
	MTensor T0;
	mint i, I0, rank, dims[1];
	int err = LIBRARY_NO_ERROR;

	I0 = MArgument_getInteger(Args[0]);
	dims[0] = 2;
    dims[1] = I0;
	
	err = libData->MTensor_new(MType_Real, 2, dims, &T0);
    int len = libData->MTensor_getFlattenedLength( T0); 
    mreal* elems = libData->MTensor_getRealData(T0);
    funkcicka(len, elems);
	MArgument_setMTensor(Res, T0);
	return err;
}
";

createTensLib=CreateLibrary[createTensor,"FUN53Tensor", "Compiler"->CCompilerDriver`MinGWCompiler`MinGWCompiler, "ExtraObjectFiles"->
{"C:\\POKUSYC\\objektove\\rounding_control.o", 
"C:\\POKUSYC\\objektove\\operations.o", 
"C:\\POKUSYC\\objektove\\interval.o",
"C:\\POKUSYC\\objektove\\controls.o", 
"C:\\POKUSYC\\objektove\\aux_ops.o", 
"C:\\POKUSYC\\pokus\\funkcickaObdelnik.o"
},
"ShellOutputFunction"->Print]

integrateEuler=LibraryFunctionLoad[createTensLib, "integratorEulerIntervalyC",{Integer},{Real, 2}]

a=integrateEuler[2^16];

FullForm[a[[1, -100;;-1]]-a[[2, -100;;-1]]]

ListPlot[a[[1]], Joined->True]

ListPlot[a[[2]], Joined->True]

Show[%, %%]

Exp[1]

IntervalMemberQ[Interval[{a[[1, -1]],a[[2,-1]]}], E]

ListPlot[a[[2]]-a[[1]], Joined->True]

Quit[]

Needs["CCompilerDriver`"]

createTensor="
#include \"WolframLibrary.h\"

extern void funkceSinus(int lengthD, double* D);

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
	return 0;
}

DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
	return;
}

DLLEXPORT int constantzero(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){
   MArgument_setInteger(Res, 0);
   return LIBRARY_NO_ERROR;
}



DLLEXPORT int integratorSinusIntervalyC(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
	MTensor T0;
	mint i, I0, rank, dims[1];
	int err = LIBRARY_NO_ERROR;

	I0 = MArgument_getInteger(Args[0]);
	dims[0] = 4;
    dims[1] = I0;
	
	err = libData->MTensor_new(MType_Real, 2, dims, &T0);
    int len = libData->MTensor_getFlattenedLength( T0); 
    mreal* elems = libData->MTensor_getRealData(T0);
    funkceSinus(len, elems);
	MArgument_setMTensor(Res, T0);
	return err;
}
";

createTensLib=CreateLibrary[createTensor,"FUN53Tensor", "Compiler"->CCompilerDriver`MinGWCompiler`MinGWCompiler, "ExtraObjectFiles"->
{"C:\\POKUSYC\\objektove\\rounding_control.o", 
"C:\\POKUSYC\\objektove\\operations.o", 
"C:\\POKUSYC\\objektove\\interval.o",
"C:\\POKUSYC\\objektove\\controls.o", 
"C:\\POKUSYC\\objektove\\aux_ops.o", 
"C:\\POKUSYC\\pokus\\funkceSinus.o"
},
"ShellOutputFunction"->Print]

integrateSinus=LibraryFunctionLoad[createTensLib, "integratorSinusIntervalyC",{Integer},{Real, 2}]

a=integrateSinus[7*2^16];

7*2^16

a[[1,2^16]]

a[[2,2^16]]

Sin[1.]

ListPlot[a[[1]], Joined->True]

ListPlot[a[[2]]-a[[1]], Joined->True]

Quit[]

Needs["CCompilerDriver`"]

createTensor="
#include \"WolframLibrary.h\"

extern void funkceSinusP(int lengthD, double* D, double pInf, double pSup);

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
	return 0;
}

DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
	return;
}

DLLEXPORT int constantzero(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){
   MArgument_setInteger(Res, 0);
   return LIBRARY_NO_ERROR;
}



DLLEXPORT int integratorSinusIntervalyC(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
	MTensor T0;
	mint i, I0, rank, dims[1];
	int err = LIBRARY_NO_ERROR;

	I0 = MArgument_getInteger(Args[0]);
    double pInf = MArgument_getReal(Args[1]);
    double pSup = MArgument_getReal(Args[2]);
	dims[0] = 4;
    dims[1] = I0;
	
	err = libData->MTensor_new(MType_Real, 2, dims, &T0);
    int len = libData->MTensor_getFlattenedLength( T0); 
    mreal* elems = libData->MTensor_getRealData(T0);
    funkceSinusP(len, elems, pInf, pSup);
	MArgument_setMTensor(Res, T0);
	return err;
}
";

createTensLib=CreateLibrary[createTensor,"FUN68Tensor", "Compiler"->CCompilerDriver`MinGWCompiler`MinGWCompiler, "ExtraObjectFiles"->
{"C:\\POKUSYC\\objektove\\rounding_control.o", 
"C:\\POKUSYC\\objektove\\operations.o", 
"C:\\POKUSYC\\objektove\\interval.o",
"C:\\POKUSYC\\objektove\\controls.o", 
"C:\\POKUSYC\\objektove\\aux_ops.o", 
"C:\\POKUSYC\\pokus\\funkceSinusP.o"
},
"ShellOutputFunction"->Print]

integrateSinusP=LibraryFunctionLoad[createTensLib, "integratorSinusIntervalyC",{Integer, Real, Real},{Real, 2}]

$HistoryLength=2

a=integrateSinusP[3*2^16, 1.5,1.51];

12*2^16

a[[1,-100;;-1]]

a[[1,-1]]

Sin[1.]

ListPlot[{a[[1]],a[[2]]}, Joined->True]

a=integrateSinusP[3*2^16, 3.5,3.51];

12*2^16

a[[1,-100;;-1]]

ListPlot[{a[[1]],a[[2]]}, Joined->True]

ListPlot[a[[2]]-a[[1]], Joined->True]

Quit[]

Needs["CCompilerDriver`"]

createTensor="
#include \"WolframLibrary.h\"

extern void arcsp(double aInf, double aSup, double dInf, double dSup, double* D);

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
	return 0;
}

DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
	return;
}

DLLEXPORT int constantzero(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){
   MArgument_setInteger(Res, 0);
   return LIBRARY_NO_ERROR;
}



DLLEXPORT int ArcSinpC(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
	MTensor T0;
	mint i, I0, rank, dims[1];
	int err = LIBRARY_NO_ERROR;

	double aInf = MArgument_getReal(Args[0]);
	double aSup = MArgument_getReal(Args[1]);
	double dInf = MArgument_getReal(Args[2]);
	double dSup = MArgument_getReal(Args[3]);

    rank = 1;
	dims[0] = 2;
	
	err = libData->MTensor_new(MType_Real, rank, dims, &T0);
    int len = libData->MTensor_getFlattenedLength( T0); 
    mreal* elems = libData->MTensor_getRealData(T0);
    arcsp(aInf, aSup, dInf, dSup, elems);
	MArgument_setMTensor(Res, T0);
	return err;
}
";

createTensLib=CreateLibrary[createTensor,"FUN59Tensor", "Compiler"->CCompilerDriver`MinGWCompiler`MinGWCompiler, "ExtraObjectFiles"->
{"C:\\POKUSYC\\objektove\\rounding_control.o", 
"C:\\POKUSYC\\pokus\\operations.o", 
"C:\\POKUSYC\\objektove\\interval.o",
"C:\\POKUSYC\\objektove\\controls.o", 
"C:\\POKUSYC\\objektove\\aux_ops.o", 
"C:\\POKUSYC\\pokus\\arcsinpOptimalTest.o"
},
"ShellOutputFunction"->Print]

ArcSinp=LibraryFunctionLoad[createTensLib, "ArcSinpC",{Real, Real, Real, Real},{Real, 1}]

ArcSinp[1.4, 1.5,0.8, 0.81]

ArcSinp[1.4, 1.4,0.8, 0.81]

ArcSinp[1.5, 1.5,0.8, 0.81]

ArcSinpL[p_,z_]:=First[ArcSinp[p, p,z, z]]

ArcSinpH[p_,z_]:=Last[ArcSinp[p, p,z, z]]

Manipulate[Plot[ArcSinpL[p,z], {z,0,1}], {p,2,10}]

Quit[]

Needs["CCompilerDriver`"]

createTensor="
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include \"WolframLibrary.h\"

extern void sp(double aInf, double aSup, double dInf, double dSup, double* D);

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
	return 0;
}

DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
	return;
}

DLLEXPORT int constantzero(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){
   MArgument_setInteger(Res, 0);
   return LIBRARY_NO_ERROR;
}



DLLEXPORT int SinpC(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
	MTensor T0;
	mint i, I0, rank, dims[1];
	int err = LIBRARY_NO_ERROR;

	double aInf = MArgument_getReal(Args[0]);
	double aSup = MArgument_getReal(Args[1]);
	double dInf = MArgument_getReal(Args[2]);
	double dSup = MArgument_getReal(Args[3]);

    rank = 1;
	dims[0] = 2;
	
	err = libData->MTensor_new(MType_Real, rank, dims, &T0);
    int len = libData->MTensor_getFlattenedLength( T0); 
    mreal* elems = libData->MTensor_getRealData(T0);
    sp(aInf, aSup, dInf, dSup, elems);
	MArgument_setMTensor(Res, T0);
	return err;
}
";

createTensLib=CreateLibrary[createTensor,"FUN61Tensor", "Compiler"->CCompilerDriver`MinGWCompiler`MinGWCompiler, "ExtraObjectFiles"->
{"C:\\POKUSYC\\objektove\\rounding_control.o", 
"C:\\POKUSYC\\objektove\\operations.o", 
"C:\\POKUSYC\\objektove\\interval.o",
"C:\\POKUSYC\\objektove\\controls.o", 
"C:\\POKUSYC\\objektove\\aux_ops.o",
"C:\\POKUSYC\\pokus\\extraOps.o",
"C:\\POKUSYC\\pokus\\sinpOpt.o"
},
"ShellOutputFunction"->Print]

Sinp=LibraryFunctionLoad[createTensLib, "SinpC",{Real, Real, Real, Real},{Real, 1}]

Sinp[2., 2.,1., 1.]

Sin[1.]

IntervalMemberQ[Interval[Sinp[2., 2.,1., 1.]],Sin[1]]

Solve[1/p+1/q==1,q]

Integrate[1/(1-z^p)^(1/p), {z,0,1}]

SinpH[p_,z_]:=Last[Sinp[p, p,z, z]]

SinpL[p_,z_]:=First[Sinp[p, p,z, z]]

piP[p_]:=((π Csc[(π/p)])/p)

Plot[SinpL[2.3, z], {z,0,piP[2.3]-0.000001}, PlotRange->All, PlotStyle->Thick]

Manipulate[Plot[{SinpL[p, z*piP[p]],SinpH[p, z*piP[p]]}, {z,0,1-0.000001}, PlotRange->All, PlotStyle->Thick, PlotPoints->1000], {p, 2, 20}]

TraditionalForm[Integrate[(Subscript[Sin, p])[x], {x,0.001,0.999}]]

Total[Interval[Sinp[2.1, 2.1,Min[#], Max[#]]]*Interval[{1/1000., 1/1000.}]&/@Table[Interval[{i/1000,(i+1)/1000}],{i, 0.001,999.}]]

Conlusion

We have implemented arithmetic operations and functions (Subscript[arcsin, p]), (Subscript[sin, p]) in GAS for Intel Pentium. This library can be easily used with the new Mathematica 8 technology to perform sofisticated computations with these functions. There is still a lot of work to be done on our library (periodic extension of (Subscript[sin, p]), rewritting the library for 64-bit platform etc.)