Go back to the main page of Sollya

- abs
- absolute
- accurateinfnorm
- acos
- acosh
- &&
- :.
- ~
- asciiplot
- asin
- asinh
- atan
- atanh
- autodiff
- autosimplify
- bashevaluate
- bashexecute
- binary
- boolean
- canonical
- ceil
- checkinfnorm
- coeff
- @
- constant
- cos
- cosh
- D
- DD
- DE
- decimal
- default
- degree
- denominator
- diam
- dieonerrormode
- diff
- dirtyfindzeros
- dirtyinfnorm
- dirtyintegral
- display
- /
- double
- doubledouble
- doubleextended
- dyadic
- ==
- erfc
- erf
- error
- evaluate
- execute
- expand
- exp
- expm1
- exponent
- externalplot
- externalproc
- false
- file
- findzeros
- fixed
- floating
- floor
- fpminimax
- fullparentheses
- function
- >=
- >
- guessdegree
- halfprecision
- head
- hexadecimal
- honorcoeffprec
- hopitalrecursions
- horner
- HP
- implementconstant
- implementpoly
- inf
- infnorm
- in
- integer
- integral
- isbound
- isevaluable
- <=
- length
- libraryconstant
- library
- list of
- log10
- log1p
- log2
- log
- <
- mantissa
- max
- mid
- midpointmode
- min
- -
- *
- nearestint
- !=
- nop
- !
- numberroots
- numerator
- off
- on
- ||
- parse
- perturb
- pi
- plot
- +
- points
- postscriptfile
- postscript
- ^
- powers
- precision
- prec
- .:
- printdouble
- printexpansion
- printsingle
- printxml
- procedure
- proc
- QD
- quad
- quit
- range
- rationalapprox
- rationalmode
- RD
- readfile
- readxml
- relative
- remez
- rename
- restart
- return
- revert
- RN
- roundcoefficients
- roundcorrectly
- roundingwarnings
- round
- RU
- RZ
- searchgal
- SG
- simplifysafe
- simplify
- single
- sinh
- sin
- sort
- sqrt
- string
- subpoly
- substitute
- supnorm
- sup
- tail
- tanh
- tan
- taylorform
- taylorrecursions
- taylor
- TD
- time
- timing
- tripledouble
- true
- var
- verbosity
- void
- worstcase
- write

- The use of absolute in the command externalplot indicates that during
plotting in externalplot an absolute error is to be considered.

See externalplot for details. - Used with fpminimax, absolute indicates that fpminimax must try to minimize
the absolute error.

See fpminimax for details. - When given in argument to supnorm, absolute indicates that an absolute error
is to be considered for supremum norm computation.

See supnorm for details.

> bashexecute("gcc -shared -o externalplotexample externalplotexample.o -lgmp -lmpfr");

> externalplot("./externalplotexample",absolute,exp(x),[-1/2;1/2],12,perturb);

- function represents the function whose infinity norm is to be computed
- range represents the infinity norm is to be considered on
- constant represents the number of bits in the significant of the result
- exclusion range 1 through exclusion range n represent ranges to be excluded

- The command accurateinfnorm computes an upper bound to the infinity norm of
function function in range. This upper bound is the least
floating-point number greater than the value of the infinity norm that
lies in the set of dyadic floating point numbers having constant
significant mantissa bits. This means the value accurateinfnorm evaluates to
is at the time an upper bound and a faithful rounding to constant
bits of the infinity norm of function function on range range.

If given, the fourth and further arguments of the command accurateinfnorm, exclusion range 1 through exclusion range n the infinity norm of the function function is not to be considered on. - Users should be aware about the fact that the algorithm behind accurateinfnorm is highly inefficient and that other, better suited algorithms, such as supnorm, are available inside Sollya. As a matter of fact, while accurateinfnorm is maintained for compatibility reasons with legacy Sollya codes, users are advised to avoid using accurateinfnorm in new Sollya scripts and to replace it, where possible, by the supnorm command.

> accurateinfnorm(p - exp(x), [-1;1], 20);

4.52055246569216251373291015625e-5

> accurateinfnorm(p - exp(x), [-1;1], 30);

4.520552107578623690642416477203369140625e-5

> accurateinfnorm(p - exp(x), [-1;1], 40);

4.5205521043867324948450914234854280948638916015625e-5

> midpointmode = on!;

> infnorm(p - exp(x), [-1;1]);

0.45205~5/7~e-4

> accurateinfnorm(p - exp(x), [-1;1], 40);

4.5205521043867324948450914234854280948638916015625e-5

- expr1 and expr2 represent boolean expressions

false

true

- L is a list (possibly empty).
- x is an object of any type.

- :. adds the element x at the end of the list L.
- Note that since x may be of any type, it can in particular be a list.

[|2, 3, 4, 5|]

[|1, 2, 3, [|4, 5, 6|]|]

[|1|]

- expression stands for an expression that is a constant
- something stands for some language element that is not a constant expression

- ~ expression evaluates the expression that is a constant term to a floating-point constant. The evaluation may involve a rounding. If expression is not a constant, the evaluated constant is a faithful rounding of expression with precision bits, unless the expression is exactly 0 as a result of cancellation. In the latter case, a floating-point approximation of some (unknown) accuracy is returned.
- ~ does not do anything on all language elements that are not a constant expression. In other words, it behaves like the identity function on any type that is not a constant expression. It can hence be used in any place where one wants to be sure that expressions are simplified using floating-point computations to constants of a known precision, regardless of the type of actual language elements.
- ~ error evaluates to error and provokes a warning.
- ~ is a prefix operator not requiring parentheses. Its precedence is the same as for the unary + and - operators. It cannot be repeatedly used without brackets.

exp(5)

> print(~ exp(5));

1.48413159102576603421115580040552279623487667593878e2

0

exp(x)

> print(~ "Hello");

Hello

exp((pi) * 5 * x)

> print(exp(x* ~(5*Pi)));

exp(x * 1.57079632679489661923132169163975144209858469968757e1)

1.48413159102576603421115580040552279623487667593878e2 * x

> print( (~exp(5))*x);

1.48413159102576603421115580040552279623487667593878e2 * x

> print(~(exp(5)*x));

exp(5) * x

- function represents a function to be plotted
- range represents a range the function is to be plotted in

- asciiplot plots the function function in range range using ASCII characters. On systems that provide the necessary TIOCGWINSZ ioctl, Sollya determines the size of the terminal for the plot size if connected to a terminal. If it is not connected to a terminal or if the test is not possible, the plot is of fixed size 77x25 characters. The function is evaluated on a number of points equal to the number of columns available. Its value is rounded to the next integer in the range of lines available. A letter "x" is written at this place. If zero is in the hull of the image domain of the function, an x-axis is displayed. If zero is in range, a y-axis is displayed. If the function is constant or if the range is reduced to one point, the function is evaluated to a constant and the constant is displayed instead of a plot.

x

xx

xxx

xx

xx

xxx

xx

xxx

xx

xxx

xxx

xxx

xxx

xxx

xxx

xxxx

xxxx

xxx

xxxx

xxxxx

xxxxx

xxxx

xxxxx

xxxx

| x

| x

| x

| x

| x

| xx

| x

| x

| xx

| xx

| xx

| xx

| x

| xxx

| xx

| xxx

| xxx

| xxxx

| xxxx

| xxxx

| xxxxxx

---------------------xxxxxxxx-----------------------------------------------

xxxxxxxxxxxx |

xxxxxxxxx |

5

2.71828182845904523536028747135266249775724709369998

- f is the function to be differentiated.
- n is the order of differentiation.
- x0 is the point at which the function is differentiated.
- I is the interval over which the function is differentiated.

- autodiff computes the first n derivatives of f at point x0. The computation is performed numerically, without symbolically differentiating the expression of f. Yet, the computation is safe since small interval enclosures are produced. More precisely, autodiff returns a list [f_0, ..., f_n] such that, for each i, f_i is a small interval enclosing the exact value of f^(i)(x0).
- Since it does not perform any symbolic differentiation, autodiff is much more efficient than diff and should be prefered when only numerical values are necessary.
- If an interval I is provided instead of a point x0, the list returned by autodiff satisfies: for each i, f^(i)(I) is included in f_i. A particular use is when one wants to know the successive derivatives of a function at a non representable point such as pi. In this case, it suffices to call autodiff with the (almost) point interval I = [pi].
- When I is almost a point interval, the returned enclosures f_i are also almost point intervals. However, when the interval I begins to be fairly large, the enclosures can be deeply overestimated due to the dependecy phenomenon present with interval arithmetic.
- As a particular case, f_0 is an enclosure of the image of f over I. However, since the algorithm is not specially designed for this purpose it is not very efficient for this particular task. In particular, it is not able to return a finite enclosure for functions with removable singularities (e.g. sin(x)/x at 0). The command evaluate is much more efficient for computing an accurate enclosure of the image of a function over an interval.

> midpointmode = on!;

> for f_i in L do f_i;

0.3559752813266941742012789792982961497379810154498~2/4~e1

0.5403023058681397174009366074429766037323104206179~0/3~

-0.3019450507398802024611853185539984893647499733880~6/2~e1

-0.252441295442368951995750696489089699886768918239~6/4~e1

0.31227898756481033145214529184139729746320579069~1/3~e1

-0.16634307959006696033484053579339956883955954978~3/1~e2

> L = autodiff(log(cos(x)+x), 5, [2,4]);

> L[0];

[0;1.27643852425465597132446653114905059102580436018893]

> evaluate(f, [2,4]);

[0.45986058925497069206106494332976097408234056912429;1.20787210589964169595901037621103012113048821362855]

> fprime = diff(f);

> L[1];

[2.53086745013099407167484456656211083053393118778677e-2;1.75680249530792825137263909451182909413591288733649]

> evaluate(fprime,[2,4]);

[2.71048755415961996452136364304380881763456815673085e-2;1.10919530663943290837397225788623531405558431279949]

> L[0];

[-@Inf@;@Inf@]

> evaluate(sin(x)/x, [-1,1]);

[0.5403023058681397174009366074429766037323104206179;1]

- An assignment autosimplify = activation value, where activation value
is one of on or off, activates respectively deactivates the
automatic safe simplification of expressions of functions generated by
the evaluation of commands or in argument of other commands.

Sollya commands like remez, taylor or rationalapprox sometimes produce expressions that can be simplified. Constant subexpressions can be evaluated to dyadic floating-point numbers, monomials with coefficients 0 can be eliminated. Further, expressions indicated by the user perform better in many commands when simplified before being passed in argument to a command. When the automatic simplification of expressions is activated, Sollya automatically performs a safe (not value changing) simplification process on such expressions.

The automatic generation of subexpressions can be annoying, in particular if it takes too much time for not enough benefit. Further the user might want to inspect the structure of the expression tree returned by a command. In this case, the automatic simplification should be deactivated.

If the assignment autosimplify = activation value is followed by an exclamation mark, no message indicating the new state is displayed. Otherwise the user is informed of the new state of the global mode by an indication.

> print(x - x);

0

> autosimplify = off ;

Automatic pure tree simplification has been deactivated.

> print(x - x);

x - x

> print(rationalapprox(sin(pi/5.9),7));

0.5

> autosimplify = off !;

> print(rationalapprox(sin(pi/5.9),7));

1 / 2

See also: print, prec, points, diam, display, verbosity, canonical, taylorrecursions, timing, fullparentheses, midpointmode, hopitalrecursions, remez, rationalapprox, taylor

Go back to the list of commands- command is a command to be interpreted by the shell.
- input is an optional character sequence to be fed to the command.

- bashevaluate(command) will execute the shell command command in a shell. All output on the command's standard output is collected and returned as a character sequence.
- If an additional argument input is given in a call to bashevaluate(command,input), this character sequence is written to the standard input of the command command that gets executed.
- All characters output by command are included in the character sequence to which bashevaluate evaluates but two exceptions. Every NULL character ('\0') in the output is replaced with '?' as Sollya is unable to handle character sequences containing that character. Additionally, if the output ends in a newline character ('\n'), this character is stripped off. Other newline characters which are not at the end of the output are left as such.

Mon Oct 22 12:02:46 CEST 2012

[|"Hello"|]

> a;

Hello

See also: bashexecute

Go back to the list of commands- command is a command to be interpreted by the shell.

- bashexecute(command) lets the shell interpret command. It is useful to execute some external code within Sollya.
- bashexecute does not return anything. It just executes its argument. However, if command produces an output in a file, this result can be imported in Sollya with help of commands like execute, readfile and parse.

Mon Oct 22 12:02:51 CEST 2012

- boolean represents the boolean type for declarations
of external procedures by means of externalproc.

Remark that in contrast to other indicators, type indicators like boolean cannot be handled outside the externalproc context. In particular, they cannot be assigned to variables.

- The command canonical rewrites the expression representing the function function in a way such that all polynomial subexpressions (or the whole expression itself, if it is a polynomial) are written in canonical form, i.e. as a sum of monomials in the canonical base. The canonical base is the base of the integer powers of the global free variable. The command canonical does not endanger the safety of computations even in Sollya's floating-point environment: the function returned is mathematically equal to the function function.
- An assignment canonical = activation value, where activation value
is one of on or off, activates respectively deactivates the
automatic printing of polynomial expressions in canonical form,
i.e. as a sum of monomials in the canonical base. If automatic
printing in canonical form is deactivated, automatic printing yields to
displaying polynomial subexpressions in Horner form.

If the assignment canonical = activation value is followed by an exclamation mark, no message indicating the new state is displayed. Otherwise the user is informed of the new state of the global mode by an indication.

1 + x^2 + 3 * x^3

> print(canonical((x + 1)^7));

1 + 7 * x + 21 * x^2 + 35 * x^3 + 35 * x^4 + 21 * x^5 + 7 * x^6 + x^7

exp(1 + 5 * x + 10 * x^2 + 10 * x^3 + 5 * x^4 + x^5) - log(asin(16 + 80 * x + 160 * x^2 + 160 * x^3 + 80 * x^4 + 16 * x^5) + x)

off

> (x + 2)^9;

512 + x * (2304 + x * (4608 + x * (5376 + x * (4032 + x * (2016 + x * (672 + x * (144 + x * (18 + x))))))))

> canonical = on;

Canonical automatic printing output has been activated.

> (x + 2)^9;

512 + 2304 * x + 4608 * x^2 + 5376 * x^3 + 4032 * x^4 + 2016 * x^5 + 672 * x^6 + 144 * x^7 + 18 * x^8 + x^9

> canonical;

on

> canonical = off!;

> (x + 2)^9;

512 + x * (2304 + x * (4608 + x * (5376 + x * (4032 + x * (2016 + x * (672 + x * (144 + x * (18 + x))))))))

- function represents the function whose infinity norm is to be checked
- range represents the infinity norm is to be considered on
- constant represents the upper bound the infinity norm is to be checked to

- The command checkinfnorm checks whether the infinity norm of the given
function function in the range range can be proven (by Sollya) to
be less than the given bound bound. This means, if checkinfnorm
evaluates to true, the infinity norm has been proven (by Sollya's
interval arithmetic) to be less than the bound. If checkinfnorm evaluates
to false, there are two possibilities: either the bound is less than
or equal to the infinity norm of the function or the bound is greater
than the infinity norm but Sollya could not conclude using its
internal interval arithmetic.

checkinfnorm is sensitive to the global variable diam. The smaller diam, the more time Sollya will spend on the evaluation of checkinfnorm in order to prove the bound before returning false although the infinity norm is bounded by the bound. If diam is equal to 0, Sollya will eventually spend infinite time on instances where the given bound bound is less or equal to the infinity norm of the function function in range range. In contrast, with diam being zero, checkinfnorm evaluates to true iff the infinity norm of the function in the range is bounded by the given bound.

true

> checkinfnorm(sin(x),[0;1.75], 1/2); checkinfnorm(sin(x),[0;20/39],1/2);

false

true

> b = dirtyinfnorm(p - exp(x), [-1;1]);

> checkinfnorm(p - exp(x), [-1;1], b);

false

> b1 = round(b, 15, RU);

> checkinfnorm(p - exp(x), [-1;1], b1);

true

> b2 = round(b, 25, RU);

> checkinfnorm(p - exp(x), [-1;1], b2);

false

> diam = 1b-20!;

> checkinfnorm(p - exp(x), [-1;1], b2);

true

- f is a function (usually a polynomial).
- n is an integer

10

0

- L1 and L2 are two lists.
- string1 and string2 are two strings.
- proc is a procedure.

- In its first usage form, @ concatenates two lists or strings.
- In its second usage form, @ applies the elements of a list as arguments to a procedure. In the case when proc is a procedure with a fixed number of arguments, a check is done if the number of elements in the list corresponds to the number of formal parameters of the procedure. An empty list can therefore applied only to a procedure that does not take any argument. In the case of a procedure with an arbitrary number of arguments, no such check is performed.

[|1, 2, 3, 7, 8, 9|]

Hello World!

write(a,", ", b," and ",c," are cool guys.\n");

};

> cool @ [| "Christoph", "Mioara", "Sylvain" |];

Christoph, Mioara and Sylvain are cool guys.

"Hello! how are you?";

};

> sayhello();

Hello! how are you?

> sayhello @ [||];

Hello! how are you?

var acc, i;

acc = 0;

for i in L do acc = i + acc;

return acc;

};

> add(1,2);

3

> add(1,2,3);

6

> add @ [|1, 2|];

3

> add @ [|1, 2, 3|];

6

> add @ [||];

0

- constant represents the constant type for declarations
of external procedures externalproc.

Remark that in contrast to other indicators, type indicators like constant cannot be handled outside the externalproc context. In particular, they cannot be assigned to variables.

- cos is the usual cosine function.
- It is defined for every real number x.

- cosh is the usual hyperbolic function: cosh(x) = (exp(x)+exp(-x))/2.
- It is defined for every real number x.

- default is a special value and is replaced by something depending on the context where it is used. It can often be used as a joker, when you want to specify one of the optional parameters of a command and not the others: set the value of uninteresting parameters to default.
- Global variables can be reset by affecting them the special value default.

> q = remez(exp(x),5,[0;1],1,1e-5);

> p==q;

true

165

> prec=200;

The precision has been set to 200 bits.

- f is a function (usually a polynomial).

3

> degree(0);

0

-1

- expr represents an expression

- If expr represents a fraction expr1/expr2, denominator(expr)
returns the denominator of this fraction, i.e. expr2.

If expr represents something else, denominator(expr) returns 1.

Note that for all expressions expr, numerator(expr) / denominator(expr) is equal to expr.

3

1

> b = numerator(a)/denominator(a);

> print(a);

5 / 3

> print(b);

5 / 3

> b = numerator(a)/denominator(a);

> print(a);

exp(x / 3)

> print(b);

exp(x / 3)

- width represents the maximal relative width of the intervals used

- diam is a global variable. Its value represents the maximal width allowed for intervals involved in safe algorithms of Sollya (namely infnorm, checkinfnorm, accurateinfnorm, integral, findzeros, supnorm).
- More precisely, diam is relative to the width of the input interval of the command. For instance, suppose that diam=1e-5: if infnorm is called on interval [0;1], the maximal width of an interval will be 1e-5. But if it is called on interval [0;1e-3], the maximal width will be 1e-8.

- activation value controls if Sollya is exited on an error or not.

- dieonerrormode is a global variable. When its value is off, which is the default, Sollya will not exit on any syntax, typing, side-effect errors. These errors will be caught by the tool, even if a memory might be leaked at that point. On evaluation, the error special value will be produced.
- When the value of the dieonerrormode variable is on, Sollya will exit on any syntax, typing, side-effect errors. A warning message will be printed in these cases at appropriate verbosity levels.

> dieonerrormode = off;

Die-on-error mode has been deactivated.

> for i from true to false do i + "Salut";

Warning: one of the arguments of the for loop does not evaluate to a constant.

The for loop will not be executed.

> exp(17);

Warning: rounding has happened. The value displayed is a faithful rounding of the true result.

2.41549527535752982147754351803858238798675673527224e7

> dieonerrormode = off!;

> 5 */ 4;

Warning: syntax error, unexpected DIVTOKEN.

The last symbol read has been "/".

Will skip input until next semicolon after the unexpected token. May leak memory.

exp(17);

Warning: rounding has happened. The value displayed is a faithful rounding of the true result.

2.41549527535752982147754351803858238798675673527224e7

> dieonerrormode;

off

> dieonerrormode = on!;

> dieonerrormode;

on

> for i from true to false do i + "Salut";

Warning: one of the arguments of the for loop does not evaluate to a constant.

The for loop will not be executed.

Warning: some syntax, typing or side-effect error has occurred.

As the die-on-error mode is activated, the tool will be exited.

> dieonerrormode = on!;

> 5 */ 4;

Warning: syntax error, unexpected DIVTOKEN.

The last symbol read has been "/".

Will skip input until next semicolon after the unexpected token. May leak memory.

Warning: some syntax, typing or side-effect error has occurred.

As the die-on-error mode is activated, the tool will be exited.

> dieonerrormode = on!;

> 5 */ 4;

- function represents a function

cos(x)

1

x^x * (1 + log(x))

- f is a function.
- I is an interval.

- dirtyfindzeros(f,I) returns a list containing some zeros of f in the interval I. The values in the list are numerical approximation of the exact zeros. The precision of these approximations is approximately the precision stored in prec. If f does not have two zeros very close to each other, it can be expected that all zeros are listed. However, some zeros may be forgotten. This command should be considered as a numerical algorithm and should not be used if safety is critical.
- More precisely, the algorithm relies on global variables prec and points and it performs the following steps:
let n be the value of variable points and t be the value
of variable prec.
- Evaluate |f| at n evenly distributed points in the interval I. The working precision to be used is automatically chosen in order to ensure that the sign is correct.
- Whenever f changes its sign for two consecutive points, find an approximation x of its zero with precision t using Newton's algorithm. The number of steps in Newton's iteration depends on t: the precision of the approximation is supposed to be doubled at each step.
- Add this value to the list.

[|-3.14159265358979323846264338327950288419716939937508, 0, 3.14159265358979323846264338327950288419716939937508|]

> points=1000!;

> L2=dirtyfindzeros(x^2*sin(1/x),[0;1]);

> length(L1); length(L2);

18

25

- f is a function.
- I is an interval.

- dirtyinfnorm(f,I) computes an approximation of the infinity norm of the given function f on the interval I, e.g. max {|f(x)|, x in I}.
- The interval must be bound. If the interval contains one of -Inf or +Inf, the result of dirtyinfnorm is NaN.
- The result of this command depends on the global variables prec and points. Therefore, the returned result is generally a good approximation of the exact infinity norm, with precision prec. However, the result is generally underestimated and should not be used when safety is critical. Use infnorm instead.
- The following algorithm is used: let n be the value of variable points
and t be the value of variable prec.
- Evaluate |f| at n evenly distributed points in the interval I. The evaluation are faithful roundings of the exact results at precision t.
- Whenever the derivative of f changes its sign for two consecutive points, find an approximation x of its zero with precision t. Then compute a faithful rounding of |f(x)| at precision t.
- Return the maximum of all computed values.

1

> dirtyinfnorm(exp(cos(x))*sin(x),[0;5]);

1.45856

> prec=40!;

> dirtyinfnorm(exp(cos(x))*sin(x),[0;5]);

1.458528537135

> prec=100!;

> dirtyinfnorm(exp(cos(x))*sin(x),[0;5]);

1.458528537136237644438147455024

> prec=200!;

> dirtyinfnorm(exp(cos(x))*sin(x),[0;5]);

1.458528537136237644438147455023841718299214087993682374094153

@NaN@

- f is a function.
- I is an interval.

- dirtyintegral(f,I) computes an approximation of the integral of f on I.
- The interval must be bound. If the interval contains one of -Inf or +Inf, the result of dirtyintegral is NaN, even if the integral has a meaning.
- The result of this command depends on the global variables prec and points. The method used is the trapezium rule applied at n evenly distributed points in the interval, where n is the value of global variable points.
- This command computes a numerical approximation of the exact value of the integral. It should not be used if safety is critical. In this case, use command integral instead.
- Warning: this command is currently known to be unsatisfactory. If you really need to compute integrals, think of using an other tool or report a feature request to sylvain.chevillard@ens-lyon.org.

-0.54402111088936981340474766185137728168364301291621

> dirtyintegral(cos(x),[0;10]);

-0.54400304905152629822448058882475382036536298356281

> points=2000!;

> dirtyintegral(cos(x),[0;10]);

-0.54401997751158321972222697312583199035995837926892

- notation value represents a variable of type decimal|binary|dyadic|powers|hexadecimal

- An assignment display = notation value, where notation value is
one of decimal, dyadic, powers, binary or hexadecimal, activates
the corresponding notation for output of values in print, write or
at the Sollya prompt.

If the global notation variable display is decimal, all numbers will be output in scientific decimal notation. If the global notation variable display is dyadic, all numbers will be output as dyadic numbers with Gappa notation. If the global notation variable display is powers, all numbers will be output as dyadic numbers with a notation compatible with Maple and PARI/GP. If the global notation variable display is binary, all numbers will be output in binary notation. If the global notation variable display is hexadecimal, all numbers will be output in C99/ IEEE754-2008 notation. All output notations can be parsed back by Sollya, inducing no error if the input and output precisions are the same (see prec).

If the assignment display = notation value is followed by an exclamation mark, no message indicating the new state is displayed. Otherwise the user is informed of the new state of the global mode by an indication.

Display mode is decimal numbers.

> a = evaluate(sin(pi * x), 0.25);

> a;

0.70710678118654752440084436210484903928483593768847

> display = binary;

Display mode is binary numbers.

> a;

1.01101010000010011110011001100111111100111011110011001001000010001011001011111011000100110110011011101010100101010111110100111110001110101101111011000001011101010001_2 * 2^(-1)

> display = hexadecimal;

Display mode is hexadecimal numbers.

> a;

0xb.504f333f9de6484597d89b3754abe9f1d6f60ba88p-4

> display = dyadic;

Display mode is dyadic numbers.

> a;

33070006991101558613323983488220944360067107133265b-165

> display = powers;

Display mode is dyadic numbers in integer-power-of-2 notation.

> a;

33070006991101558613323983488220944360067107133265 * 2^(-165)

- function1 and function2 represent functions
- interval1 and interval2 represent intervals (ranges)
- constant represents a constant or constant expression

- / represents the division (function) on reals. The expression function1 / function2 stands for the function composed of the division function and the two functions function1 and function2, where function1 is the numerator and function2 the denominator.
- / can be used for interval arithmetic on intervals (ranges). / will evaluate to an interval that safely encompasses all images of the division function with arguments varying in the given intervals. If the intervals given contain points where the division function is not defined, infinities and NaNs will be produced in the output interval. Any combination of intervals with intervals or constants (resp. constant expressions) is supported. However, it is not possible to represent families of functions using an interval as one argument and a function (varying in the free variable) as the other one.

2.5

x * 0.5

1

@NaN@

(exp(x) * cos(x) - sin(x) * exp(x)) / exp(x)^2

[0.25;0.66666666666666666666666666666666666666666666666668]

> [1;2] / 17;

[5.8823529411764705882352941176470588235294117647058e-2;0.11764705882352941176470588235294117647058823529412]

> -13 / [4;17];

[-3.25;-0.76470588235294117647058823529411764705882352941175]

- double is both a function and a constant.
- As a function, it rounds its argument to the nearest IEEE 754 double precision (i.e. IEEE754-2008 binary64) number. Subnormal numbers are supported as well as standard numbers: it is the real rounding described in the standard.
- As a constant, it symbolizes the double precision format. It is used in contexts when a precision format is necessary, e.g. in the commands round, roundcoefficients and implementpoly. See the corresponding help pages for examples.

> D(0.1);

1.100110011001100110011001100110011001100110011001101_2 * 2^(-4)

> D(4.17);

1.000010101110000101000111101011100001010001111010111_2 * 2^(2)

> D(1.011_2 * 2^(-1073));

1.1_2 * 2^(-1073)

See also: halfprecision, single, doubleextended, doubledouble, quad, tripledouble, roundcoefficients, implementpoly, round, printdouble

Go back to the list of commands- doubledouble is both a function and a constant.
- As a function, it rounds its argument to the nearest number that can be written as the sum of two double precision numbers.
- The algorithm used to compute doubledouble(x) is the following: let xh = double(x) and let xl = double(x - xh). Return the number xh + xl. Note that if the current precision is not sufficient to exactly represent xh + xl, a rounding will occur and the result of doubledouble(x) will be useless.
- As a constant, it symbolizes the double-double precision format. It is used in contexts when a precision format is necessary, e.g. in the commands round, roundcoefficients and implementpoly. See the corresponding help pages for examples.

> a = 1+ 2^(-100);

> DD(a);

1.0000000000000000000000000000007888609052210118054

> prec=50!;

> DD(a);

Warning: double rounding occurred on invoking the double-double rounding operator.

Try to increase the working precision.

1

See also: halfprecision, single, double, doubleextended, quad, tripledouble, roundcoefficients, implementpoly, round

Go back to the list of commands- doubleextended is a function that computes the nearest floating-point number with 64 bits of mantissa to a given number. Since it is a function, it can be composed with other Sollya functions such as exp, sin, etc.
- doubleextended now does handle subnormal numbers for a presumed exponent width of the double-extended format of 15 bits. This means, with respect to rounding, doubleextended behaves as a IEEE 754-2008 binary79 with a 64 bit significand (with a hidden bit normal range), one sign bit and a 15 bit exponent field would behave. This behavior may be different from the one observed on Intel-based IA32/Intel64 processors (or compatible versions from other vendors). However it is the one seen on HP/Intel Itanium when the precision specifier is double-extended and pseudo-denormals are activated.
- Since it is a function and not a command, its behavior is a bit different from the behavior of round(x,64,RN) even if the result is exactly the same. round(x,64,RN) is immediately evaluated whereas doubleextended(x) can be composed with other functions (and thus be plotted and so on).

> DE(0.1);

1.100110011001100110011001100110011001100110011001100110011001101_2 * 2^(-4)

> round(0.1,64,RN);

1.100110011001100110011001100110011001100110011001100110011001101_2 * 2^(-4)

0

> DE(2^(-20000));

0

> f = sin(DE(x));

> f(pi);

Warning: rounding has happened. The value displayed is a faithful rounding of the true result.

-5.0165576126683320235573270803307570138315616702549e-20

> g = sin(round(x,64,RN));

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

- expr1 and expr2 represent expressions

- The operator == evaluates to true iff its operands expr1 and expr2 are syntactically equal and different from error or constant expressions that are not constants and that evaluate to the same floating-point number with the global precision prec. The user should be aware of the fact that because of floating-point evaluation, the operator == is not exactly the same as the mathematical equality. Further remark that according to IEEE 754-2008 floating-point rules, which Sollya emulates, floating-point data which are NaN do not compare equal to any other floating-point datum, including NaN.

true

> "Hello" == "Salut";

false

> "Hello" == 5;

false

> 5 + x == 5 + x;

true

true

> asin(1) * 2 == pi;

true

> exp(5) == log(4);

false

false

The precision has been set to 12 bits.

> 16384.1 == 16385.1;

true

false

> b = NaN;

> a == a;

true

> b == b;

false

- erfc is the complementary error function defined by erfc(x) = 1 - erf(x).
- It is defined for every real number x.

See also: erf

Go back to the list of commands- erf is the error function defined by: erf(x) = (2/sqrt(Pi)) * integral(exp(-t^2), [0;x])
- It is defined for every real number x.

- The variable error represents an input during the evaluation of
which a type or execution error has been detected or is to be
detected. Inputs that are syntactically correct but wrongly typed
evaluate to error at some stage. Inputs that are correctly typed
but containing commands that depend on side-effects that cannot be
performed or inputs that are wrongly typed at meta-level (cf. parse),
evaluate to error.

Remark that in contrast to all other elements of the Sollya language, error compares neither equal nor unequal to itself. This provides a means of detecting syntax errors inside the Sollya language itself without introducing issues of two different wrongly typed inputs being equal.

error

error

false

> error != error;

false

> incorrect = 5 + "foo";

> correct == correct;

true

> incorrect == incorrect;

false

> errorhappened = !(incorrect == incorrect);

> errorhappened;

true

- function represents a function
- constant represents a constant point
- range represents a range
- function2 represents a function that is not constant

- If its second argument is a constant constant, evaluate evaluates its first argument function at the point indicated by constant. This evaluation is performed in a way that the result is a faithful rounding of the real value of the function at constant to the current global precision. If such a faithful rounding is not possible, evaluate returns a range surely encompassing the real value of the function function at constant. If even interval evaluation is not possible because the expression is undefined or numerically unstable, NaN will be produced.
- If its second argument is a range range, evaluate evaluates its first argument function by interval evaluation on this range range. This ensures that the image domain of the function function on the preimage domain range is surely enclosed in the returned range.
- In the case when the second argument is a range that is reduced to a single point (such that [1;1] for instance), the evaluation is performed in the same way as when the second argument is a constant but it produces a range as a result: evaluate automatically adjusts the precision of the intern computations and returns a range that contains at most three floating-point consecutive numbers in precision prec. This correponds to the same accuracy as a faithful rounding of the actual result. If such a faithful rounding is not possible, evaluate has the same behavior as in the case when the second argument is a constant.
- If its second argument is a function function2 that is not a constant, evaluate replaces all occurrences of the free variable in function function by function function2.

> print(evaluate(sin(pi * x), 2.25));

0.70710678118654752440084436210484903928483593768847

> print(evaluate(sin(pi * x), [2.25; 2.25]));

0.707106781186547524400844362104849039284835937688~4/5~

[-1.72986452514381269516508615031098129542836767991679e-12715;7.5941198201187963145069564314525661706039084390067e-12716]

[-5.143390272677254630046998919961912407349224165421e-50;0.70710678118654752440084436210484903928483593768866]

sin((pi) * (2 + 0.25 * x))

[@NaN@;@NaN@]

See also: isevaluable

Go back to the list of commands- filename is a string representing a file name

- execute opens the file indicated by filename, and executes the sequence of commands it contains. This command is evaluated at execution time: this way you can modify the file filename (for instance using bashexecute) and execute it just after.
- If filename contains a command execute, it will be executed recursively.
- If filename contains a call to restart, it will be neglected.
- If filename contains a call to quit, the commands following quit in filename will be neglected.

> a;

2

> print("a=1;") > "example.sollya";

> execute("example.sollya");

> a;

1

> print("a=1; restart; a=2;") > "example.sollya";

> execute("example.sollya");

Warning: a restart command has been used in a file read into another.

This restart command will be neglected.

> a;

2

> print("a=1; quit; a=2;") > "example.sollya";

> execute("example.sollya");

Warning: the execution of a file read by execute demanded stopping the interpretation but it is not stopped.

> a;

1

- function represents a function

- expand(function) expands all polynomial subexpressions in function function as far as possible. Factors of sums are multiplied out, power operators with constant positive integer exponents are replaced by multiplications and divisions are multiplied out, i.e. denomiators are brought at the most interior point of expressions.

x * x * x

8 + 12 * x + 6 * x * x + x * x * x + 2 * x

exp(243 + 405 * x + 270 * x * x + 90 * x * x * x + 15 * x * x * x * x + x * x * x * x * x + x * 405 + 108 * x * 5 * x + 54 * x * x * 5 * x + 12 * x * x * x * 5 * x + x * x * x * x * 5 * x + x * x * 270 + 27 * x * x * x * 10 + 9 * x * x * x * x * 10 + x * x * x * x * x * 10 + x * x * x * 90 + 6 * x * x * x * x * 10 + x * x * x * x * x * 10 + x * x * x * x * 5 * x + 15 * x * x * x * x + x * x * x * x * x)

- expm1 is defined by expm1(x) = exp(x)-1.
- It is defined for every real number x.

See also: exp

Go back to the list of commands- x is a dyadic number.

> e=exponent(a);

> e;

-17

> m=mantissa(a);

> a-m*2^e;

0

- The command externalplot plots the error of an external function
evaluation code sequence implemented in the object file named
filename with regard to the function function. If mode
evaluates to absolute, the difference of both functions is
considered as an error function; if mode evaluates to relative,
the difference is divided by the function function. The resulting
error function is plotted on all floating-point numbers with
precision significant mantissa bits in the range range.

If the sixth argument of the command externalplot is given and evaluates to perturb, each of the floating-point numbers the function is evaluated at gets perturbed by a random value that is uniformly distributed in +/-1 ulp around the original precision bit floating-point variable.

If a sixth and seventh argument, respectively a seventh and eighth argument in the presence of perturb as a sixth argument, are given that evaluate to a variable of type file|postscript|postscriptfile respectively to a character sequence of type string, externalplot will plot (additionally) to a file in the same way as the command plot does. See plot for details.

The external function evaluation code given in the object file name filename is supposed to define a function name f as follows (here in C syntax): void f(mpfr_t rop, mpfr_t op). This function is supposed to evaluate op with an accuracy corresponding to the precision of rop and assign this value to rop.

> bashexecute("gcc -shared -o externalplotexample externalplotexample.o -lgmp -lmpfr");

> externalplot("./externalplotexample",relative,exp(x),[-1/2;1/2],12,perturb);

See also: plot, asciiplot, perturb, absolute, relative, file, postscript, postscriptfile, bashexecute, externalproc, library

Go back to the list of commands- identifier represents the identifier the code is to be bound to
- filename of type string represents the name of the object file where the code of procedure can be found
- argumenttype represents a definition of the types of the arguments of the Sollya procedure and the external code
- resulttype represents a definition of the result type of the external code

- externalproc allows for binding the Sollya identifier
identifier to an external code. After this binding, when Sollya
encounters identifier applied to a list of actual parameters, it
will evaluate these parameters and call the external code with these
parameters. If the external code indicated success, it will receive
the result produced by the external code, transform it to Sollya's
internal representation and return it.

In order to allow correct evaluation and typing of the data in parameter and in result to be passed to and received from the external code, externalproc has a third parameter argumenttype -> resulttype. Both argumenttype and resulttype are one of void, constant, function, range, integer, string, boolean, list of constant, list of function, list of range, list of integer, list of string, list of boolean.

If upon a usage of a procedure bound to an external procedure the type of the actual parameters given or its number is not correct, Sollya produces a type error. An external function not applied to arguments represents itself and prints out with its argument and result types.

The external function is supposed to return an integer indicating success. It returns its result depending on its Sollya result type as follows. Here, the external procedure is assumed to be implemented as a C function.

If the Sollya result type is void, the C function has no pointer argument for the result. If the Sollya result type is constant, the first argument of the C function is of C type mpfr_t *, the result is returned by affecting the MPFR variable. If the Sollya result type is function, the first argument of the C function is of C type node**, the result is returned by the node * pointed with a new node *. If the Sollya result type is range, the first argument of the C function is of C type sollya_mpfi_t *, the result is returned by affecting the interval-arithmetic variable. If the Sollya result type is integer, the first argument of the C function is of C type int *, the result is returned by affecting the int variable. If the Sollya result type is string, the first argument of the C function is of C type char **, the result is returned by the char * pointed with a new char *. If the Sollya result type is boolean, the first argument of the C function is of C type int *, the result is returned by affecting the int variable with a boolean value. If the Sollya result type is list of type, the first argument of the C function is of C type chain **, the result is returned by the chain * pointed with a new chain *. This chain contains for Sollya type constant pointers mpfr_t * to new MPFR variables, for Sollya type function pointers node * to new nodes, for Sollya type range pointers sollya_mpfi_t * to new interval-arithmetic variables, for Sollya type integer pointers int * to new int variables for Sollya type string pointers char * to new char * variables and for Sollya type boolean pointers int * to new int variables representing boolean values.

The external procedure affects its possible pointer argument if and only if it succeeds. This means, if the function returns an integer indicating failure, it does not leak any memory to the encompassing environment.

The external procedure receives its arguments as follows: If the Sollya argument type is void, no argument array is given. Otherwise the C function receives a C void ** argument representing an array of size equal to the arity of the function where each entry (of C type void *) represents a value with a C type depending on the corresponding Sollya type. If the Sollya type is constant, the C type the void * is to be casted to is mpfr_t *. If the Sollya type is function, the C type the void * is to be casted to is node *. If the Sollya type is range, the C type the void * is to be casted to is sollya_mpfi_t *. If the Sollya type is integer, the C type the void * is to be casted to is int *. If the Sollya type is string, the C type the void * is to be casted to is char *. If the Sollya type is boolean, the C type the void * is to be casted to is int *. If the Sollya type is list of type, the C type the void * is to be casted to is chain *. Here depending on type, the values in the chain are to be casted to mpfr_t * for Sollya type constant, node * for Sollya type function, sollya_mpfi_t * for Sollya type range, int * for Sollya type integer, char * for Sollya type string and int * for Sollya type boolean.

The external procedure is not supposed to alter the memory pointed by its array argument void **.

In both directions (argument and result values), empty lists are represented by chain * NULL pointers.

In contrast to internal procedures, externally bounded procedures can be considered to be objects inside Sollya that can be assigned to other variables, stored in list etc.

> bashexecute("gcc -fPIC -shared -o externalprocexample externalprocexample.o");

> externalproc(foo, "./externalprocexample", (integer, integer) -> integer);

> foo;

foo(integer, integer) -> integer

> foo(5, 6);

11

> verbosity = 1!;

> foo();

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

error

> a = foo;

> a(5,6);

11

See also: library, libraryconstant, externalplot, bashexecute, void, constant, function, range, integer, string, boolean, list of

Go back to the list of commands- false is the usual boolean value.

false

> 2<1;

false

- file is a special value used in commands plot and externalplot to save the result of the command in a data file.
- As any value it can be affected to a variable and stored in lists.

> name="plotSinCos";

> plot(sin(x),0,cos(x),[-Pi,Pi],savemode, name);

- f is a function.
- I is an interval.

- findzeros(f,I) returns a list of intervals I1, ... ,In such that, for every zero z of f, there exists some k such that z is in Ik.
- The list may contain intervals Ik that do not contain any zero of f. An interval Ik may contain many zeros of f.
- This command is meant for cases when safety is critical. If you want to be sure not to forget any zero, use findzeros. However, if you just want to know numerical values for the zeros of f, dirtyfindzeros should be quite satisfactory and a lot faster.
- If d denotes the value of global variable diam, the algorithm ensures that for each k, |Ik| < d*|I|.
- The algorithm used is basically a bisection algorithm. It is the same algorithm that the one used for infnorm. See the help page of this command for more details. In short, the behavior of the algorithm depends on global variables prec, diam, taylorrecursions and hopitalrecursions.

[|[-3.14208984375;-3.140869140625], [-1.220703125e-3;1.220703125e-3], [3.140869140625;3.14208984375]|]

> diam=1e-10!;

> findzeros(sin(x),[-5;5]);

[|[-3.14159265370108187198638916015625;-3.141592652536928653717041015625], [-1.16415321826934814453125e-9;1.16415321826934814453125e-9], [3.141592652536928653717041015625;3.14159265370108187198638916015625]|]

0.9999997480772435665130615234375 + x^2 * (-0.4999928693287074565887451171875 + x^2 * (4.163351492024958133697509765625e-2 + x^2 * (-1.338223926723003387451171875e-3)))

0.99999974816012948686250183527590706944465637207031 + x * (5.5210044061222495131782035802443168321913900126185e-14 + x * (-0.4999928698019768802396356477402150630950927734375 + x * (-3.95371609372064761555136192612768146546591008227978e-13 + x * (4.16335155285858099505347240665287245064973831176758e-2 + x * (5.2492670395835122748014980938834327670386437070249e-13 + x * (-1.33822408807599468535953768366653093835338950157166e-3))))))

- f is the function to be approximated
- n is the degree of the polynomial that must approximate f
- monomials is the list of monomials that must be used to represent the polynomial that approximates f
- formats is a list indicating the formats that the coefficients of the polynomial must have
- range is the interval where the function must be approximated
- L is a list of interpolation points used by the method
- indic1 (optional) is one of the optional indication parameters. See the detailed description below.
- indic2 (optional) is one of the optional indication parameters. See the detailed description below.
- indic3 (optional) is one of the optional indication parameters. See the detailed description below.
- P (optional) is the minimax polynomial to be considered for solving the problem.

- fpminimax uses a heuristic (but practically efficient) method to find a
good polynomial approximation of a function f on an interval range. It
implements the method published in the article:

Efficient polynomial L^\infty - approximations Nicolas Brisebarre and Sylvain Chevillard

Proceedings of the 18th IEEE Symposium on Computer Arithmetic (ARITH 18)

pp. 169-176 - The basic usage of this command is fpminimax(f, n, formats, range). It computes a polynomial approximation of f with degree at most n on the interval range. formats is a list of integers or format types (such as double, doubledouble, etc.). The polynomial returned by the command has its coefficients that fit the formats indications. For instance, if formats[0] is 35, the coefficient of degree 0 of the polynomial will fit a floating-point format of 35 bits. If formats[1] is D, the coefficient of degree 1 will be representable by a floating-point number with a precision of 53 bits (which is not necessarily an IEEE 754 double precision number. See the remark below), etc.
- The second argument may be either an integer or a list of integers
interpreted as the list of desired monomials. For instance, the list
[|0,2,4,6|] indicates that the polynomial must be even and of
degree at most 6. Giving an integer n as second argument is equivalent
as giving [|0,...,n|].

The list of formats is interpreted with respect to the list of monomials. For instance, if the list of monomials is [|0,2,4,6|] and the list of formats is [|161,107,53,24|], the coefficients of degree 0 is searched as a floating-point number with precision 161, the coefficient of degree 2 is searched as a number of precision 107, and so on. - The list of formats may contain either integers or format types (halfprecision, single, double, doubledouble, tripledouble, doubleextended and quad). The list may be too large or even infinite. Only the first indications will be considered. For instance, for a degree n polynomial, formats[n+1] and above will be discarded. This lets one use elliptical indications for the last coefficients.
- The floating-point coefficients considered by fpminimax do not have an exponent range. In particular, in the format list, double is an exact synonym for 53. Currently, fpminimax only ensures that the corresponding coefficient has at most 53 bits of mantissa. It does not imply that it is an IEEE-754 double.
- By default, the list of formats is interpreted as a list of floating-point
formats. This may be changed by passing fixed as an optional argument (see
below). Let us take an example: fpminimax(f, 2, [107, DD, 53], [0;1]).
Here the optional argument is missing (we could have set it to floating).
Thus, fpminimax will search for a polynomial of degree 2 with a constant
coefficient that is a 107 bits floating-point number, etc.

Currently, doubledouble is just a synonym for 107 and tripledouble a synonym for 161. This behavior may change in the future (taking into account the fact that some double-doubles are not representable with 107 bits).

Second example: fpminimax(f, 2, [25, 18, 30], [0;1], fixed). In this case, fpminimax will search for a polynomial of degree 2 with a constant coefficient of the form m/2^25 where m is an integer. In other words, it is a fixed-point number with 25 bits after the point. Note that even with argument fixed, the formats list is allowed to contain halfprecision, single, double, doubleextended, doubledouble, quad or tripledouble. In this this case, it is just a synonym for 11, 24, 53, 64, 107, 113 or 161. This is deprecated and may change in the future. - The fourth argument may be a range or a list. Lists are for advanced users that know what they are doing. The core of the method is a kind of approximated interpolation. The list given here is a list of points that must be considered for the interpolation. It must contain at least as many points as unknown coefficients. If you give a list, it is also recommended that you provide the minimax polynomial as last argument. If you give a range, the list of points will be automatically computed.
- The fifth, sixth and seventh arguments are optional. By default, fpminimax
will approximate f while optimizing the relative error, and interpreting
the list of formats as a list of floating-point formats.

This default behavior may be changed with these optional arguments. You may provide zero, one, two or three of the arguments in any order. This lets the user indicate only the non-default arguments.

The three possible arguments are:- relative or absolute: the error to be optimized;
- floating or fixed: formats of the coefficients;
- a constrained part q.

- The last argument is for advanced users. It is the minimax polynomial that
approximates the function f in the monomial basis. If it is not given
this polynomial will be automatically computed by fpminimax.

This minimax polynomial is used to compute the list of interpolation points required by the method. In general, you do not have to provide this argument. But if you want to obtain several polynomials of the same degree that approximate the same function on the same range, just changing the formats, you should probably consider computing only once the minimax polynomial and the list of points instead of letting fpminimax recompute them each time.

Note that in the case when a constrained part is given, the minimax polynomial must take that into account. For instance, in the previous example, the minimax would be obtained by the following command: P = remez(1-(1+x+x^2/2)/exp(x), [|3,4|], range, 1/exp(x)); Note that the constrained part is not to be added to P. - Note that fpminimax internally computes a minimax polynomial (using the same algorithm as remez command). Thus fpminimax may encounter the same problems as remez. In particular, it may be very slow when Haar condition is not fulfilled. Another consequence is that currently fpminimax has to be run with a sufficiently high working precision.

> printexpansion(P);

(0x3ff0000000000000 + 0xbc09fda20235c100) + x * ((0x3b29ecd485d34781 + 0xb7c1cbc97120359a) + x * (0xbfdfffffffffff98 + x * (0xbbfa6e0b3183cb0d + x * (0x3fa5555555145337 + x * (0x3ca3540480618939 + x * 0xbf56c138142d8c3b)))))

> display = powers!;

> P;

x * (1 + x^2 * (-357913941 * 2^(-31) + x^2 * (35789873 * 2^(-32))))

> display = powers!;

> P;

1 + x * (1 + x * (1 * 2^(-1) + x * (375300225001191 * 2^(-51) + x * (5592621 * 2^(-27)))))

> pstar = remez(f, 5, [-1b-7;1b-7]);

> listpoints = dirtyfindzeros(f-pstar, [-1b-7; 1b-7]);

> P1 = fpminimax(f, 5, [|DD...|], listpoints, absolute, default, default, pstar);

> P2 = fpminimax(f, 5, [|D...|], listpoints, absolute, default, default, pstar);

> P3 = fpminimax(f, 5, [|D, D, D, 24...|], listpoints, absolute, default, default, pstar);

> print("Error of pstar: ", dirtyinfnorm(f-pstar, [-1b-7; 1b-7]));

Error of pstar: 7.9048441305459735102879831325718745399379329453102e-16

> print("Error of P1: ", dirtyinfnorm(f-P1, [-1b-7; 1b-7]));

Error of P1: 7.9048441305459735159848647089192667442047469404883e-16

> print("Error of P2: ", dirtyinfnorm(f-P2, [-1b-7; 1b-7]));

Error of P2: 8.2477144579950871061147021597406077993657714575238e-16

> print("Error of P3: ", dirtyinfnorm(f-P3, [-1b-7; 1b-7]));

Error of P3: 1.08454277156993282593701156841863009789063333951055e-15

See also: remez, dirtyfindzeros, absolute, relative, fixed, floating, default, halfprecision, single, double, doubleextended, doubledouble, quad, tripledouble, implementpoly, coeff, degree, roundcoefficients, guessdegree

Go back to the list of commands- An assignment fullparentheses = activation value, where activation value
is one of on or off, activates respectively deactivates the output
of expressions with full parenthesising. In full parenthesising mode,
Sollya commands like print, write and the implicit command when an
expression is given at the prompt will output expressions with
parenthesising at all places where it is necessary for expressions
containing infix operators to be parsed back with the same
result. Otherwise parentheses around associative operators are
omitted.

If the assignment fullparentheses = activation value is followed by an exclamation mark, no message indicating the new state is displayed. Otherwise the user is informed of the new state of the global mode by an indication.

> fullparentheses = off;

Full parentheses mode has been deactivated.

> print(1 + 2 + 3);

1 + 2 + 3

> fullparentheses = on;

Full parentheses mode has been activated.

> print(1 + 2 + 3);

(1 + 2) + 3

- procedure is a procedure of type (range, integer, integer) -> range

- For the sake of safety and mathematical consistency, Sollya distinguishes clearly between functions, seen in the mathematical sense of the term, i.e. mappings, and procedures, seen in the sense Computer Science gives to functions, i.e. pieces of code that compute results for arguments following an algorithm. In some cases however, it is interesting to use such Computer Science procedures as realisations of mathematical functions, e.g. in order to plot them or even to perform polynomial approximation on them. The function keyword allows for such a transformation of a Sollya procedure into a Sollya function.
- The procedure to be used as a function through function(procedure) must be of type (range, integer, integer) -> range. This means it must take in argument an interval X, a degree of differentiation n and a working precision p. It must return in result an interval Y encompassing the image f^(n)(X) of the n-th derivative of the implemented function f, i.e. f^(n)(X) c Y. In order to allow Sollya's algorithms to work properly, the procedure must ensure that, whenever (p, diam(X)) tends to (infinity, 0), the computed over-estimated bounding Y tends to the actual image f^(n)(X).
- The user must be aware that they are responsible of the correctness of the procedure. If, for some n and X, procedure returns an interval Y such that f^n(X) is not included in Y, function will successfully return a function without any complain, but this function might behave inconsistently in further computations.
- For cases when the procedure does not have the correct signature or does not return a finite interval as a result function(procedure) evaluates to Not-A-Number (resp. to an interval of Not-A-Numbers for interval evaluation).
- function also represents the function type for declarations
of external procedures by means of externalproc.

Remark that in contrast to other indicators, type indicators like function cannot be handled outside the externalproc context. In particular, they cannot be assigned to variables.

var res, oldPrec;

oldPrec = prec;

prec = p!;

res = exp(X);

prec = oldPrec!;

return res;

};

> f = function(EXP);

> f(1);

2.71828182845904523536028747135266249775724709369998

> exp(1);

2.71828182845904523536028747135266249775724709369998

> f(x + 3);

(function(proc(X, n, p)

{

var res, oldPrec;

oldPrec = prec;

prec = p!;

res = exp(X);

prec = oldPrec!;

return res;

}))(3 + x)

> diff(f);

diff(function(proc(X, n, p)

{

var res, oldPrec;

oldPrec = prec;

prec = p!;

res = exp(X);

prec = oldPrec!;

return res;

}))

> (diff(f))(0);

1

> g = f(sin(x));

> g(17);

0.382358169993866834026905546416556413595734583420876

> diff(g);

(diff(function(proc(X, n, p)

{

var res, oldPrec;

oldPrec = prec;

prec = p!;

res = exp(X);

prec = oldPrec!;

return res;

})))(sin(x)) * cos(x)

> (diff(g))(1);

1.25338076749344683697237458088447611474812675164344

> p = remez(f,3,[-1/2;1/2]);

> p;

0.9996712090142519365811043588840936667986880903378 + x * (0.99973702983570053280233869785694438940067223265505 + x * (0.51049729360256555535800202052281444451304355667385 + x * 0.1698143246071767617700502198641549152447429302716))

See also: proc, library, procedure, externalproc, boolean, constant, integer, list of, range, string

Go back to the list of commands- expr1 and expr2 represent constant expressions

- The operator >= evaluates to true iff its operands expr1 and expr2 evaluate to two floating-point numbers a1 respectively a2 with the global precision prec and a1 is greater than or equal to a2. The user should be aware of the fact that because of floating-point evaluation, the operator >= is not exactly the same as the mathematical operation "greater-than-or-equal-to".

true

> 5 >= 5;

true

> 5 >= 6;

false

> exp(2) >= exp(1);

true

> log(1) >= exp(2);

false

The precision has been set to 12 bits.

> 16384.1 >= 16385.1;

true

- expr1 and expr2 represent constant expressions

- The operator > evaluates to true iff its operands expr1 and expr2 evaluate to two floating-point numbers a1 respectively a2 with the global precision prec and a1 is greater than a2. The user should be aware of the fact that because of floating-point evaluation, the operator > is not exactly the same as the mathematical operation "greater-than".

true

> 5 > 5;

false

> 5 > 6;

false

> exp(2) > exp(1);

true

> log(1) > exp(2);

false

The precision has been set to 12 bits.

> 16385.1 > 16384.1;

false

- f is the function to be approximated.
- I is the interval where the function must be approximated.
- eps is the maximal acceptable error.
- w (optional) is a weight function. Default is 1.
- bound (optional) is a bound on the degree. Default is currently 128.

- guessdegree tries to find the minimal degree needed to approximate f on I by a polynomial with an error err=p*w-f whose infinity norm not greater than eps. More precisely, it finds n minimal such that there exists a polynomial p of degree n such that ||p*w-f|| < eps.
- guessdegree returns an interval: for common cases, this interval is reduced
to a single number (i.e. the minimal degree). But in certain cases,
guessdegree does not succeed in finding the minimal degree. In such cases the
returned interval is of the form [n,p] such that:
- no polynomial of degree n-1 gives an error less than eps.
- there exists a polynomial of degree p giving an error less than eps.

- The fifth optional argument bound is used to prevent guessdegree from trying to find too large degrees. If guessdegree does not manage to find a degree n satisfying the error and such that n<=bound, an interval of the form [*, +Inf] is returned. Note that bound must be a positive integer.

[10;10]

[10;@Inf@]

[8;9]

- halfprecision is both a function and a constant.
- As a function, it rounds its argument to the nearest IEEE 754 half-precision (i.e. IEEE754-2008 binary16) number. Subnormal numbers are supported as well as standard numbers: it is the real rounding described in the standard.
- As a constant, it symbolizes the half-precision format. It is used in contexts when a precision format is necessary, e.g. in the commands round and roundcoefficients. It is not supported for implementpoly. See the corresponding help pages for examples.

> HP(0.1);

1.100110011_2 * 2^(-4)

> HP(4.17);

1.00001011_2 * 2^(2)

> HP(1.011_2 * 2^(-23));

1.1_2 * 2^(-23)

See also: single, double, doubleextended, doubledouble, quad, tripledouble, roundcoefficients, fpminimax, implementpoly, round, printsingle

Go back to the list of commands- L is a list.

- head(L) returns the first element of the list L. It is equivalent to L[0].
- If L is empty, the command will fail with an error.

1

> head([|1,2...|]);

1

See also: tail

Go back to the list of commands- hexadecimal is a special value used for the global state display. If
the global state display is equal to hexadecimal, all data will be
output in hexadecimal C99/ IEEE 754-2008 notation.

As any value it can be affected to a variable and stored in lists.

- Used with command implementpoly, honorcoeffprec makes implementpoly honor the precision of the given polynomial. This means if a coefficient needs a double-double or a triple-double to be exactly stored, implementpoly will allocate appropriate space and use a double-double or triple-double operation even if the automatic (heuristic) determination implemented in command implementpoly indicates that the coefficient could be stored on less precision or, respectively, the operation could be performed with less precision. See implementpoly for details.

> q = implementpoly(1 - simplify(TD(1/6)) * x^2,[-1b-10;1b-10],1b-60,DD,"p","implementation.c");

Warning: at least one of the coefficients of the given polynomial has been rounded in a way

that the target precision can be achieved at lower cost. Nevertheless, the implemented polynomial

is different from the given one.

> printexpansion(q);

0x3ff0000000000000 + x^2 * 0xbfc5555555555555

> r = implementpoly(1 - simplify(TD(1/6)) * x^2,[-1b-10;1b-10],1b-60,DD,"p","implementation.c",honorcoeffprec);

Warning: the infered precision of the 2th coefficient of the polynomial is greater than

the necessary precision computed for this step. This may make the automatic determination

of precisions useless.

> printexpansion(r);

0x3ff0000000000000 + x^2 * (0xbfc5555555555555 + 0xbc65555555555555 + 0xb905555555555555)

- n represents the number of recursions

- hopitalrecursions is a global variable. Its value represents the number of steps of recursion that are tried when applying L'Hopital's rule. This rule is applied by the interval evaluator present in the core of Sollya (and particularly visible in commands like infnorm).
- If an expression of the form f/g has to be evaluated by interval arithmetic on an interval I and if f and g have a common zero in I, a direct evaluation leads to NaN. Sollya implements a safe heuristic to avoid this, based on L'Hopital's rule: in such a case, it can be shown that (f/g)(I) C (f'/g')(I). Since the same problem may exist for f'/g', the rule is applied recursively. The number of step in this recursion process is controlled by hopitalrecursions.
- Setting hopitalrecursions to 0 makes Sollya use this rule only once; setting it to 1 makes Sollya use the rule twice, and so on. In particular: the rule is always applied at least once, if necessary.

The number of recursions for Hopital's rule has been set to 0.

> evaluate(log(1+x)^2/x^2,[-1/2; 1]);

[-@Inf@;@Inf@]

> hopitalrecursions=1;

The number of recursions for Hopital's rule has been set to 1.

> evaluate(log(1+x)^2/x^2,[-1/2; 1]);

[-2.52258872223978123766892848583270627230200053744108;6.7725887222397812376689284858327062723020005374411]

- function represents the expression to be rewritten in Horner form

- The command horner rewrites the expression representing the function function in a way such that all polynomial subexpressions (or the whole expression itself, if it is a polynomial) are written in Horner form. The command horner does not endanger the safety of computations even in Sollya's floating-point environment: the function returned is mathematically equal to the function function.

1 + x * (2 + x * 3)

> print(horner((x + 1)^7));

1 + x * (7 + x * (21 + x * (35 + x * (35 + x * (21 + x * (7 + x))))))

exp(1 + x * (5 + x * (10 + x * (10 + x * (5 + x))))) - log(asin(x * (1 + x^2)) + x)

- The command implementconstant implements the constant expression expr in arbitrary precision. More precisely, it generates the source code (written in C, and using MPFR) of a C function const_something with the following signature: void const_something (mpfr_ptr y, mp_prec_t prec) Let us denote by c the exact mathematical value of the constant defined by the expression expr. When called with arguments y and prec (where the variable y is supposed to be already initialized), the function mpfr_const_something sets the precision of y to a suitable precision and stores in it an approximate value of c such that |y-c| <= |c|*2^(1-prec).
- When no filename filename is given or if default is given as filename, the source code produced by implementconstant is printed on standard output. Otherwise, when filename is given as a string of characters, the source code is output to a file named filename. If that file cannot be opened and/or written to, implementconstant fails and has no other effect.
- When functionname is given as an argument to implementconstant and functionname evaluates to a string of characters, the default name for the C function const_something is replaced by functionname. When default is given as functionname, the default name is used nevertheless, as if no functionname argument were given. When choosing a character sequence for functionname, the user should keep attention to the fact that functionname must be a valid C identifier in order to enable error-free compilation of the produced code.
- If expr refers to a constant defined with libraryconstant, the produced code uses the external code implementing this constant. The user should keep in mind that it is up to them to make sure the symbol for that external code can get resolved when the newly generated code is to be loaded.
- If a subexpression of expr evaluates to 0, implementconstant will most likely fail with an error message.
- implementconstant is unable to implement constant expressions expr that contain procedure-based functions, i.e. functions created from Sollya procedures using the function construct. If expr contains such a procedure-based function, implementconstant prints a warning and fails silently. The reason for this lack of functionality is that the produced C source code, which is supposed to be compiled, would have to call back to the Sollya interpreter in order to evaluate the procedure-based function.
- Similarily, implementconstant is currently unable to implement constant expressions expr that contain library-based functions, i.e. functions dynamically bound to Sollya using the library construct. If expr contains such a library-based function, implementconstant prints a warning and fails silently. Support for this feature is in principle feasible from a technical standpoint and might be added in a future release of Sollya.
- Currently, non-differentiable functions such as double, doubledouble, tripledouble, single, halfprecision, quad, doubleextended, floor, ceil, nearestint are not supported by implementconstant. If implementconstant encounters one of them, a warning message is displayed and no code is produced. However, if autosimplify equals on, it is possible that Sollya silently simplifies subexpressions of expr containing such functions and that implementconstant successfully produces code for evaluating expr.
- While it produces an MPFR-based C source code for expr, implementconstant takes architectural and system-dependent parameters into account. For example, it checks whether literal constants figuring in expr can be represented on a C long int type or if they must be stored in a different manner not to affect their accuracy. These tests, performed by Sollya during execution of implementconstant, depend themselves on the architecture Sollya is running on. Users should keep this matter in mind, especially when trying to compile source code on one machine whilst it has been produced on another.

#include <mpfr.h>

void

const_something (mpfr_ptr y, mp_prec_t prec)

{

/* Declarations */

mpfr_t tmp1;

mpfr_t tmp2;

mpfr_t tmp3;

mpfr_t tmp4;

mpfr_t tmp5;

mpfr_t tmp6;

mpfr_t tmp7;

/* Initializations */

mpfr_init2 (tmp2, prec+5);

mpfr_init2 (tmp1, prec+3);

mpfr_init2 (tmp4, prec+8);

mpfr_init2 (tmp3, prec+7);

mpfr_init2 (tmp6, prec+11);

mpfr_init2 (tmp7, prec+11);

mpfr_init2 (tmp5, prec+11);

/* Core */

mpfr_set_prec (tmp2, prec+4);

mpfr_set_ui (tmp2, 1, MPFR_RNDN);

mpfr_set_prec (tmp1, prec+3);

mpfr_exp (tmp1, tmp2, MPFR_RNDN);

mpfr_set_prec (tmp4, prec+8);

mpfr_set_ui (tmp4, 2, MPFR_RNDN);

mpfr_set_prec (tmp3, prec+7);

mpfr_log (tmp3, tmp4, MPFR_RNDN);

mpfr_set_prec (tmp6, prec+11);

mpfr_set_ui (tmp6, 1, MPFR_RNDN);

mpfr_set_prec (tmp7, prec+11);

mpfr_set_ui (tmp7, 10, MPFR_RNDN);

mpfr_set_prec (tmp5, prec+11);

mpfr_div (tmp5, tmp6, tmp7, MPFR_RNDN);

mpfr_set_prec (tmp4, prec+7);

mpfr_sqrt (tmp4, tmp5, MPFR_RNDN);

mpfr_set_prec (tmp2, prec+5);

mpfr_div (tmp2, tmp3, tmp4, MPFR_RNDN);

mpfr_set_prec (y, prec+3);

mpfr_add (y, tmp1, tmp2, MPFR_RNDN);

/* Cleaning stuff */

mpfr_clear(tmp1);

mpfr_clear(tmp2);

mpfr_clear(tmp3);

mpfr_clear(tmp4);

mpfr_clear(tmp5);

mpfr_clear(tmp6);

mpfr_clear(tmp7);

}

> readfile("sine_of_thirteen_seventeenth.c");

#include <mpfr.h>

void

const_something (mpfr_ptr y, mp_prec_t prec)

{

/* Declarations */

mpfr_t tmp1;

mpfr_t tmp2;

mpfr_t tmp3;

/* Initializations */

mpfr_init2 (tmp2, prec+6);

mpfr_init2 (tmp3, prec+6);

mpfr_init2 (tmp1, prec+6);

/* Core */

mpfr_set_prec (tmp2, prec+6);

mpfr_set_ui (tmp2, 13, MPFR_RNDN);

mpfr_set_prec (tmp3, prec+6);

mpfr_set_ui (tmp3, 17, MPFR_RNDN);

mpfr_set_prec (tmp1, prec+6);

mpfr_div (tmp1, tmp2, tmp3, MPFR_RNDN);

mpfr_set_prec (y, prec+2);

mpfr_sin (y, tmp1, MPFR_RNDN);

/* Cleaning stuff */

mpfr_clear(tmp1);

mpfr_clear(tmp2);

mpfr_clear(tmp3);

}

#include <mpfr.h>

void

arcsin_of_one_third_pi (mpfr_ptr y, mp_prec_t prec)

{

/* Declarations */

mpfr_t tmp1;

mpfr_t tmp2;

mpfr_t tmp3;

/* Initializations */

mpfr_init2 (tmp2, prec+8);

mpfr_init2 (tmp3, prec+8);

mpfr_init2 (tmp1, prec+8);

/* Core */

mpfr_set_prec (tmp2, prec+8);

mpfr_const_pi (tmp2, MPFR_RNDN);

mpfr_set_prec (tmp3, prec+8);

mpfr_set_ui (tmp3, 3, MPFR_RNDN);

mpfr_set_prec (tmp1, prec+8);

mpfr_div (tmp1, tmp2, tmp3, MPFR_RNDN);

mpfr_set_prec (y, prec+2);

mpfr_asin (y, tmp1, MPFR_RNDN);

/* Cleaning stuff */

mpfr_clear(tmp1);

mpfr_clear(tmp2);

mpfr_clear(tmp3);

}

> readfile("constant_code.c");

#include <mpfr.h>

void

magic_constant (mpfr_ptr y, mp_prec_t prec)

{

/* Initializations */

/* Core */

mpfr_set_prec (y, prec);

mpfr_set_ui (y, 3, MPFR_RNDN);

}

> bashexecute("gcc -shared -o libraryconstantexample libraryconstantexample.o -lgmp -lmpfr");

> euler_gamma = libraryconstant("./libraryconstantexample");

> implementconstant(euler_gamma^(1/3));

#include <mpfr.h>

void

const_something (mpfr_ptr y, mp_prec_t prec)

{

/* Declarations */

mpfr_t tmp1;

/* Initializations */

mpfr_init2 (tmp1, prec+1);

/* Core */

euler_gamma (tmp1, prec+1);

mpfr_set_prec (y, prec+2);

mpfr_root (y, tmp1, 3, MPFR_RNDN);

/* Cleaning stuff */

mpfr_clear(tmp1);

}

- The command implementpoly implements the polynomial polynomial in range
range as a function called functionname in C code
using double, double-double and triple-double arithmetic in a way that
the rounding error (estimated at its first order) is bounded by error bound.
The produced code is output in a file named filename. The
argument format indicates the double, double-double or triple-double
format of the variable in which the polynomial varies, influencing
also in the signature of the C function.

If a seventh or eighth argument proof filename is given and if this argument evaluates to a variable of type string, the command implementpoly will produce a Gappa proof that the rounding error is less than the given bound. This proof will be output in Gappa syntax in a file name proof filename.

The command implementpoly returns the polynomial that has been implemented. As the command implementpoly tries to adapt the precision needed in each evaluation step to its strict minimum and as it applies renormalization to double-double and triple-double precision coefficients to bring them to a round-to-nearest expansion form, the returned polynomial may differ from the polynomial polynomial. Nevertheless the difference will be small enough that the rounding error bound with regard to the polynomial polynomial (estimated at its first order) will be less than the given error bound.

If a seventh argument honor coefficient precisions is given and evaluates to a variable honorcoeffprec of type honorcoeffprec, implementpoly will honor the precision of the given polynomial polynomials. This means if a coefficient needs a double-double or a triple-double to be exactly stored, implementpoly will allocate appropriate space and use a double-double or triple-double operation even if the automatic (heuristic) determination implemented in command implementpoly indicates that the coefficient could be stored on less precision or, respectively, the operation could be performed with less precision. The use of honorcoeffprec has advantages and disadvantages. If the polynomial polynomial given has not been determined by a process considering directly polynomials with floating-point coefficients, honorcoeffprec should not be indicated. The implementpoly command can then determine the needed precision using the same error estimation as used for the determination of the precisions of the operations. Generally, the coefficients will get rounded to double, double-double and triple-double precision in a way that minimizes their number and respects the rounding error bound error bound. Indicating honorcoeffprec may in this case short-circuit most precision estimations leading to sub-optimal code. On the other hand, if the polynomial polynomial has been determined with floating-point precisions in mind, honorcoeffprec should be indicated because such polynomials often are very sensitive in terms of error propagation with regard to their coefficients' values. Indicating honorcoeffprec prevents the implementpoly command from rounding the coefficients and altering by many orders of magnitude the approximation error of the polynomial with regard to the function it approximates.

The implementer behind the implementpoly command makes some assumptions on its input and verifies them. If some assumption cannot be verified, the implementation will not succeed and implementpoly will evaluate to a variable error of type error. The same behaviour is observed if some file is not writable or some other side-effect fails, e.g. if the implementer runs out of memory.

As error estimation is performed only on the first order, the code produced by the implementpoly command should be considered valid iff a Gappa proof has been produced and successfully run in Gappa.

1 + x^2 * (-0.166666666666666657414808128123695496469736099243164 + x^2 * 8.3333333333333332176851016015461937058717012405395e-3)

> readfile("implementation.c");

#define p_coeff_0h 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000e+00

#define p_coeff_2h -1.66666666666666657414808128123695496469736099243164062500000000000000000000000000e-01

#define p_coeff_4h 8.33333333333333321768510160154619370587170124053955078125000000000000000000000000e-03

void p(double *p_resh, double x) {

double p_x_0_pow2h;

p_x_0_pow2h = x * x;

double p_t_1_0h;

double p_t_2_0h;

double p_t_3_0h;

double p_t_4_0h;

double p_t_5_0h;

p_t_1_0h = p_coeff_4h;

p_t_2_0h = p_t_1_0h * p_x_0_pow2h;

p_t_3_0h = p_coeff_2h + p_t_2_0h;

p_t_4_0h = p_t_3_0h * p_x_0_pow2h;

p_t_5_0h = p_coeff_0h + p_t_4_0h;

*p_resh = p_t_5_0h;

}

1 + x^2 * (-0.166666666666666657414808128123695496469736099243164 + x^2 * 8.3333333333333332176851016015461937058717012405395e-3)

> q = implementpoly(1 - simplify(TD(1/6)) * x^2,[-1b-10;1b-10],1b-60,DD,"p","implementation.c");

Warning: at least one of the coefficients of the given polynomial has been rounded in a way

that the target precision can be achieved at lower cost. Nevertheless, the implemented polynomial

is different from the given one.

> printexpansion(q);

0x3ff0000000000000 + x^2 * 0xbfc5555555555555

> r = implementpoly(1 - simplify(TD(1/6)) * x^2,[-1b-10;1b-10],1b-60,DD,"p","implementation.c",honorcoeffprec);

Warning: the infered precision of the 2th coefficient of the polynomial is greater than

the necessary precision computed for this step. This may make the automatic determination

of precisions useless.

> printexpansion(r);

0x3ff0000000000000 + x^2 * (0xbfc5555555555555 + 0xbc65555555555555 + 0xb905555555555555)

> q = implementpoly(p,[-1/2;1/2],1b-60,D,"p","implementation.c",honorcoeffprec,"implementation.gappa");

> if (q != p) then print("During implementation, rounding has happened.") else print("Polynomial implemented as given.");

Polynomial implemented as given.

See also: honorcoeffprec, roundcoefficients, double, doubledouble, tripledouble, readfile, printexpansion, error, remez, fpminimax, taylor

Go back to the list of commands- I is an interval.
- x is a real number.

- Returns the lower bound of the interval I. Each bound of an interval has its own precision, so this command is exact, even if the current precision is too small to represent the bound.
- When called on a real number x, inf considers it as an interval formed of a single point: [x, x]. In other words, inf behaves like the identity.

1

> inf(0);

0

> I=[0.111110000011111_2; 1];

> inf(I);

1.11110000011111_2 * 2^(-1)

> prec=12!;

> inf(I);

1.11110000011111_2 * 2^(-1)

- f is a function.
- I is an interval.
- filename (optional) is the name of the file into a proof will be saved.
- IList (optional) is a list of intervals to be excluded.

- infnorm(f,range) computes an interval bounding the infinity norm of the given function f on the interval I, e.g. computes an interval J such that max {|f(x)|, x in I} C J.
- If filename is given, a proof in English will be produced (and stored in file called filename) proving that max {|f(x)|, x in I} C J.
- If a list IList of intervals I1, ... ,In is given, the infinity norm will be computed on I \ (I1 U ... U I2).
- The function f is assumed to be at least twice continuous on I. More generally, if f is Ck, global variables hopitalrecursions and taylorrecursions must have values not greater than k.
- If the interval is reduced to a single point, the result of infnorm is an interval containing the exact absolute value of f at this point.
- If the interval is not bound, the result will be [0, +Inf] which is correct but perfectly useless. infnorm is not meant to be used with infinite intervals.
- The result of this command depends on the global variables prec, diam,
taylorrecursions and hopitalrecursions. The contribution of each variable is
not easy even to analyse.
- The algorithm uses interval arithmetic with precision prec. The precision should thus be set high enough to ensure that no critical cancellation will occur.
- When an evaluation is performed on an interval [a,b], if the result is considered being too large, the interval is split into [a,(a+b)/2] and [(a+b)/2,b] and so on recursively. This recursion step is not performed if the (b-a) < d*|I| where d is the value of variable diam. In other words, diam controls the minimum length of an interval during the algorithm.
- To perform the evaluation of a function on an interval, Taylor's rule is applied, e.g. f([a,b]) C f(m) + [a-m, b-m]*f'([a,b]) where m=(a+b)/2. This rule is recursively applied n times where n is the value of variable taylorrecursions. Roughly speaking, the evaluations will avoid decorrelation up to order n.
- When a function of the form g/h has to be evaluated on an interval [a,b] and when g and h vanish at a same point z of the interval, the ratio may be defined even if the expression g(z)/h(z)=0/0 does not make any sense. In this case, L'Hopital's rule may be used and (g/h)([a,b]) C (g'/h')([a,b]). Since the same can occur with the ratio g'/h', the rule is applied recursively. The variable hopitalrecursions controls the number of recursion steps.

- The algorithm used for this command is quite complex to be explained here.
Please find a complete description in the following article:

S. Chevillard and C. Lauter

A certified infinity norm for the implementation of elementary functions

LIP Research Report number RR2007-26

http://prunel.ccsd.cnrs.fr/ensl-00119810

- Users should be aware about the fact that the algorithm behind infnorm is inefficient in most cases and that other, better suited algorithms, such as supnorm, are available inside Sollya. As a matter of fact, while infnorm is maintained for compatibility reasons with legacy Sollya codes, users are advised to avoid using infnorm in new Sollya scripts and to replace it, where possible, by the supnorm command.

[2.00855369231876677409285296545817178969879078385537e1;2.00855369231876677409285296545817178969879078385544e1]

[2.00855369231876677409285296545817178969879078385537e1;2.00855369231876677409285296545817178969879078385544e1]

[2.00855369231876677409285296545817178969879078385537e1;2.00855369231876677409285296545817178969879078385544e1]

[2.00855369231876677409285296545817178969879078385537e1;2.00855369231876677409285296545817178969879078385544e1]

[2.71828182845904523536028747135266249775724709369989;2.71828182845904523536028747135266249775724709369998]

[0;@Inf@]

See also: prec, diam, hopitalrecursions, dirtyinfnorm, checkinfnorm, supnorm, findzeros, diff, taylorrecursions, autodiff, numberroots, taylorform

Go back to the list of commands- expr represents a constant expression
- range1 and range2 represent ranges (intervals)

- When its first operand is a constant expression expr, the operator in evaluates to true iff the constant value of the expression expr is contained in the interval range1.
- When both its operands are ranges (intervals), the operator in evaluates to true iff all values in range1 are contained in the interval range2.
- in is also used as a keyword for loops over the different elements of a list.

true

> 4 in [-1;1];

false

> 0 in sin([-17;17]);

true

true

> [2;3] in [4;5];

false

> [2;3] in [2.5;5];

false

1

2

3

4

5

- integer represents the machine integer type for declarations
of external procedures externalproc.

Remark that in contrast to other indicators, type indicators like integer cannot be handled outside the externalproc context. In particular, they cannot be assigned to variables.

- f is a function.
- I is an interval.

- integral(f,I) returns an interval J such that the exact value of the integral of f on I lies in J.
- This command is safe but very inefficient. Use dirtyintegral if you just want an approximate value.
- The result of this command depends on the global variable diam. The method used is the following: I is cut into intervals of length not greater then d*|I| where d is the value of global variable diam. On each small interval J, an evaluation of f by interval is performed. The result is multiplied by the length of J. Finally all values are summed.

-0.54402111088936981340474766185137728168364301291621

> integral(cos(x),[0;10]);

[-0.54710197983579690224097637163525943075698599257332;-0.54094015130013183848150540881373370744053741191728]

> diam=1e-5!;

> integral(cos(x),[0;10]);

[-0.54432915685955427101857780295936956775293876382777;-0.54371306401249969508039644221927489010425803173555]

- ident is a name.

- isbound(ident) returns a boolean value indicating whether the name ident is used or not to represent a variable. It returns true when ident is the name used to represent the global variable or if the name is currently used to refer to a (possibly local) variable.
- When a variable is defined in a block and has not been defined outside, isbound returns true when called inside the block, and false outside. Note that isbound returns true as soon as a variable has been declared with var, even if no value is actually stored in it.
- If ident1 is bound to a variable and if ident2 refers to the global variable, the command rename(ident2, ident1) hides the value of ident1 which becomes the global variable. However, if the global variable is again renamed, ident1 gets its value back. In this case, isbound(ident1) returns true. If ident1 was not bound before, isbound(ident1) returns false after that ident1 has been renamed.

false

> isbound(f);

false

> isbound(g);

false

> f=sin(x);

> isbound(x);

true

> isbound(f);

true

> isbound(g);

false

false

> { var a; isbound(a); };

true

> isbound(a);

false

> isbound(x);

true

> rename(x,y);

> isbound(x);

false

> f=sin(y);

> rename(y,x);

> f;

sin(x)

> x;

x

> isbound(x);

true

> rename(x,y);

> isbound(x);

true

> x;

1

- function represents a function
- constant represents a constant point

- isevaluable applied to function function and a constant constant returns a boolean indicating whether or not a subsequent call to evaluate on the same function function and constant constant will produce a numerical result or NaN. This means isevaluable returns false iff evaluate will return NaN.

true

> print(evaluate(sin(pi * 1/x), 0.75));

-0.86602540378443864676372317075293618347140262690518

true

> print(evaluate(sin(pi * 1/x), 0.5));

[-1.72986452514381269516508615031098129542836767991679e-12715;7.5941198201187963145069564314525661706039084390067e-12716]

false

> print(evaluate(sin(pi * 1/x), 0));

[@NaN@;@NaN@]

See also: evaluate

Go back to the list of commands- expr1 and expr2 represent constant expressions

- The operator <= evaluates to true iff its operands expr1 and expr2 evaluate to two floating-point numbers a1 respectively a2 with the global precision prec and a1 is less than or equal to a2. The user should be aware of the fact that because of floating-point evaluation, the operator <= is not exactly the same as the mathematical operation "less-than-or-equal-to".

false

> 5 <= 5;

true

> 5 <= 6;

true

> exp(2) <= exp(1);

false

> log(1) <= exp(2);

true

The precision has been set to 12 bits.

> 16385.1 <= 16384.1;

true

- L is a list.
- s is a string.

12

5

0

@Inf@

- The command libraryconstant lets you extend the set of mathematical constants known to Sollya. By default, the only mathematical constant known by Sollya is pi. For particular applications, one may want to manipulate other constants, such as Euler's gamma constant, for instance.
- libraryconstant makes it possible to let Sollya know about new constants. In order to let it know, you have to provide an implementation of the constant you are interested in. This implementation is a C file containing a function of the form: void my_ident(mpfr_t result, mp_prec_t prec) The semantic of this function is the following: it is an implementation of the constant in arbitrary precision. my_ident(result, prec) shall set the precision of the variable result to a suitable precision (the variable is assumed to be already initialized) and store in result an approximate value of the constant with a relative error not greater than 2^(1-prec). More precisely, if c is the exact value of the constant, the value stored in result should satisfy |result-c| <= 2^(1-prec)*|c|.
- You can include sollya.h in your implementation and use library functionnalities of Sollya for your implementation. However, this requires to have compiled Sollya with -fPIC in order to make the Sollya executable code position independent and to use a system on with programs, using dlopen to open dynamic routines can dynamically open themselves.
- To bind your constant into Sollya, you must use the same identifier as the function name used in your implementation file (my_ident in the previous example). Once the function code has been bound to an identifier, you can use a simple assignment to assign the bound identifier to yet another identifier. This way, you may use convenient names inside Sollya even if your implementation environment requires you to use a less convenient name.
- Once your constant is bound, it is considered by Sollya as an infinitely accurate constant (i.e. a 0-ary function, exactly like pi).

> bashexecute("gcc -shared -o libraryconstantexample libraryconstantexample.o -lgmp -lmpfr");

> euler_gamma = libraryconstant("./libraryconstantexample");

> prec = 20!;

> euler_gamma;

0.577215

> prec = 100!;

> euler_gamma;

0.577215664901532860606512090082

> midpointmode = on;

Midpoint mode has been activated.

> [euler_gamma];

0.577215664901532860606512090~0/1~

- The command library lets you extend the set of mathematical functions known to Sollya. By default, Sollya knows the most common mathematical functions such as exp, sin, erf, etc. Within Sollya, these functions may be composed. This way, Sollya should satisfy the needs of a lot of users. However, for particular applications, one may want to manipulate other functions such as Bessel functions, or functions defined by an integral or even a particular solution of an ODE.
- library makes it possible to let Sollya know about new functions. In order to let it know, you have to provide an implementation of the function you are interested in. This implementation is a C file containing a function of the form: int my_ident(sollya_mpfi_t result, sollya_mpfi_t op, int n) The semantic of this function is the following: it is an implementation of the function and its derivatives in interval arithmetic. my_ident(result, I, n) shall store in result an enclosure of the image set of the n-th derivative of the function f over I: f^(n)(I) C result.
- The integer value returned by the function implementation currently has no meaning.
- You do not need to provide a working implementation for any n. Most functions of Sollya requires a relevant implementation only for f, f' and f''. For higher derivatives, its is not so critical and the implementation may just store [-Inf, +Inf] in result whenever n>2.
- Note that you should respect somehow interval-arithmetic standards in your implementation: result has its own precision and you should perform the intermediate computations so that result is as tight as possible.
- You can include sollya.h in your implementation and use library functionnalities of Sollya for your implementation. However, this requires to have compiled Sollya with -fPIC in order to make the Sollya executable code position independent and to use a system on with programs, using dlopen to open dynamic routines can dynamically open themselves.
- To bind your function into Sollya, you must use the same identifier as the function name used in your implementation file (my_ident in the previous example). Once the function code has been bound to an identifier, you can use a simple assignment to assign the bound identifier to yet another identifier. This way, you may use convenient names inside Sollya even if your implementation environment requires you to use a less convenient name.

> bashexecute("gcc -shared -o libraryexample libraryexample.o -lgmp -lmpfr");

> myownlog = library("./libraryexample");

> evaluate(log(x), 2);

0.69314718055994530941723212145817656807550013436024

> evaluate(myownlog(x), 2);

0.69314718055994530941723212145817656807550013436024

- log10 is the decimal logarithm defined by: log10(x) = log(x)/log(10).
- It is defined only for x in [0; +Inf].

- log1p is the function defined by log1p(x) = log(1+x).
- It is defined only for x in [-1; +Inf].

See also: log

Go back to the list of commands- log2 is the binary logarithm defined by: log2(x) = log(x)/log(2).
- It is defined only for x in [0; +Inf].

- log is the natural logarithm defined as the inverse of the exponential function: log(y) is the unique real number x such that exp(x)=y.
- It is defined only for y in [0; +Inf].

- expr1 and expr2 represent constant expressions

- The operator < evaluates to true iff its operands expr1 and expr2 evaluate to two floating-point numbers a1 respectively a2 with the global precision prec and a1 is less than a2. The user should be aware of the fact that because of floating-point evaluation, the operator < is not exactly the same as the mathematical operation "less-than".

false

> 5 < 5;

false

> 5 < 6;

true

> exp(2) < exp(1);

false

> log(1) < exp(2);

true

The precision has been set to 12 bits.

> 16384.1 < 16385.1;

false

- x is a dyadic number.

> e=exponent(a);

> m=mantissa(a);

> m;

411775

> a-m*2^e;

0

- expr are constant expressions.
- l is a list of constant expressions.

- max determines which of a given set of constant expressions expr has maximum value. To do so, max tries to increase the precision used for evaluation until it can decide the ordering or some maximum precision is reached. In the latter case, a warning is printed indicating that there might actually be another expression that has a greater value.
- Even though max determines the maximum expression by evaluation, it returns the expression that is maximum as is, i.e. as an expression tree that might be evaluated to any accuracy afterwards.
- max can be given either an arbitrary number of constant expressions in argument or a list of constant expressions. The list however must not be end-elliptic.
- Users should be aware that the behavior of max follows the IEEE 754-2008 standard with respect to NaNs. In particular, a NaN given as the first argument will not be promoted as a result unless the other argument is a NaN. This means that NaNs may seem to disappear during computations.

1.48413159102576603421115580040552279623487667593878e2

> max(17);

17

> max(l);

1.48413159102576603421115580040552279623487667593878e2

exp(17)

> print(max(17 + log2(13)/log2(9),17 + log(13)/log(9)));

Warning: maximum computation relies on floating-point result that is faithfully evaluated and different faithful roundings toggle the result.

17 + log2(13) / log2(9)

- I is an interval.
- x is a real number.

- Returns the middle of the interval I. If the middle is not exactly representable at the current precision, the value is returned as an unevaluated expression.
- When called on a real number x, mid considers it as an interval formed of a single point: [x, x]. In other words, mid behaves like the identity.

2

> mid(17);

17

- activation value enables or disables the mode.

- midpointmode is a global variable. When its value is off, intervals are displayed as usual (in the form [a;b]). When its value is on, and if a and b have the same first significant digits, the interval in displayed in a way that lets one immediately see the common digits of the two bounds.
- This mode is supported only with display set to decimal. In other modes of display, midpointmode value is simply ignored.

> b = round(Pi,30,RU);

> d = [a,b];

> d;

[3.1415926516056060791015625;3.1415926553308963775634765625]

> midpointmode=on!;

> d;

0.314159265~1/6~e1

- expr are constant expressions.
- l is a list of constant expressions.

- min determines which of a given set of constant expressions expr has minimum value. To do so, min tries to increase the precision used for evaluation until it can decide the ordering or some maximum precision is reached. In the latter case, a warning is printed indicating that there might actually be another expression that has a lesser value.
- Even though min determines the minimum expression by evaluation, it returns the expression that is minimum as is, i.e. as an expression tree that might be evaluated to any accuracy afterwards.
- min can be given either an arbitrary number of constant expressions in argument or a list of constant expressions. The list however must not be end-elliptic.
- Users should be aware that the behavior of min follows the IEEE 754-2008 standard with respect to NaNs. In particular, a NaN given as the first argument will not be promoted as a result unless the other argument is a NaN. This means that NaNs may seem to disappear during computations.

-1.3862943611198906188344642429163531361510002687205

> min(17);

17

> min(l);

-1.3862943611198906188344642429163531361510002687205

sin(62)

> print(min(17 + log2(13)/log2(9),17 + log(13)/log(9)));

Warning: minimum computation relies on floating-point result that is faithfully evaluated and different faithful roundings toggle the result.

17 + log(13) / log(9)

- function1 and function2 represent functions
- interval1 and interval2 represent intervals (ranges)
- constant represents a constant or constant expression

- - represents the subtraction (function) on reals. The expression function1 - function2 stands for the function composed of the subtraction function and the two functions function1 and function2, where function1 is the subtrahend and function2 the subtractor.
- - can be used for interval arithmetic on intervals (ranges). - will evaluate to an interval that safely encompasses all images of the subtraction function with arguments varying in the given intervals. Any combination of intervals with intervals or constants (resp. constant expressions) is supported. However, it is not possible to represent families of functions using an interval as one argument and a function (varying in the free variable) as the other one.
- - stands also for the negation function.

3

-2 + x

0

cos(x) - exp(x)

[-3;-1]

> [1;2] - 17;

[-16;-15]

> 13 - [-4;17];

[-4;17]

-exp(x)

> -13;

-13

> -[13;17];

[-17;-13]

- function1 and function2 represent functions
- interval1 and interval2 represent intervals (ranges)
- constant represents a constant or constant expression

- * represents the multiplication (function) on reals. The expression function1 * function2 stands for the function composed of the multiplication function and the two functions function1 and function2.
- * can be used for interval arithmetic on intervals (ranges). * will evaluate to an interval that safely encompasses all images of the multiplication function with arguments varying in the given intervals. Any combination of intervals with intervals or constants (resp. constant expressions) is supported. However, it is not possible to represent families of functions using an interval as one argument and a function (varying in the free variable) as the other one.

10

x * 2

x^2

sin(x) * exp(x) + exp(x) * cos(x)

[3;8]

> [1;2] * 17;

[17;34]

> 13 * [-4;17];

[-52;221]

- nearestint is defined as usual: nearestint(x) is the integer nearest to x, with the special rule that the even integer is chosen if there exist two integers equally near to x.
- It is defined for every real number x.

- expr1 and expr2 represent expressions

- The operator != evaluates to true iff its operands expr1 and
expr2 are syntactically unequal and both different from error or
constant expressions that are not constants and that evaluate to two
different floating-point number with the global precision prec. The
user should be aware of the fact that because of floating-point
evaluation, the operator != is not exactly the same as the
negation of the mathematical equality.

Note that the expressions !(expr1 != expr2) and expr1 == expr2 do not evaluate to the same boolean value. See error for details.

false

> "Hello" != "Salut";

true

> "Hello" != 5;

true

> 5 + x != 5 + x;

false

false

> asin(1) * 2 != pi;

false

> exp(5) != log(4);

true

true

The precision has been set to 12 bits.

> 16384.1 != 16385.1;

false

false

- The command nop does nothing. This means it is an explicit parse element in the Sollya language that finally does not produce any result or side-effect.
- The command nop may take an optional positive integer argument n. The argument controls how much (useless) integer additions Sollya performs while doing nothing. With this behaviour, nop can be used for calibration of timing tests.
- The keyword nop is implicit in some procedure definitions. Procedures without imperative body get parsed as if they had an imperative body containing one nop statement.

> succ;

proc(n)

{

nop;

return (n) + (1);

}

> succ(5);

6

- expr represents a boolean expression

true

false

- p is a polynomial.
- I is an interval.

- numberroots rigorously computes the number of roots of polynomial the p in the interval I. The technique used is Sturm's algorithm. The value returned is not just a numerical estimation of the number of roots of p in I: it is the exact number of roots.
- The command findzeros computes safe enclosures of all the zeros of a function, without forgetting any, but it is not guaranteed to separate them all in distinct intervals. numberroots is more accurate since it guarantees the exact number of roots. However, it does not compute them. It may be used, for instance, to certify that findzeros did not put two distinct roots in the same interval.
- Multiple roots are counted only once.
- The interval I must be bounded. The algorithm cannot handle unbounded intervals. Moreover, the interval is considered as a closed interval: if one (or both) of the endpoints of I are roots of p, they are counted.
- The argument p can be any expression, but if Sollya fails to prove that it is a polynomial an error is produced. Also, please note that if the coefficients of p or the endpoints of I are not exactly representable, they are first numerically evaluated, before the algorithm is used. In that case, the counted number of roots corresponds to the rounded polynomial on the rounded interval *and not* to the exact parameters given by the user. A warning is displayed to inform the user.

1

> findzeros(1+x-x^2, [1,2]);

[|[1.617919921875;1.6180419921875]|]

2

> numberroots(x^2, [-1,1]);

1

> numberroots(x-pi, [0,4]);

Warning: the 0th coefficient of the polynomial is neither a floating point

constant nor can be evaluated without rounding to a floating point constant.

Will faithfully evaluate it with the current precision (165 bits)

1

> numberroots(1+x-x^2, [0, @Inf@]);

Warning: the given interval must have finite bounds.

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

error

> numberroots(exp(x), [0, 1]);

Warning: the given function must be a polynomial in this context.

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

error

- expr represents an expression

- If expr represents a fraction expr1/expr2, numerator(expr)
returns the numerator of this fraction, i.e. expr1.

If expr represents something else, numerator(expr) returns the expression itself, i.e. expr.

Note that for all expressions expr, numerator(expr) / denominator(expr) is equal to expr.

5

exp(x)

> b = numerator(a)/denominator(a);

> print(a);

5 / 3

> print(b);

5 / 3

> b = numerator(a)/denominator(a);

> print(a);

exp(x / 3)

> print(b);

exp(x / 3)

- off is a special value used to deactivate certain functionnalities of Sollya.
- As any value it can be affected to a variable and stored in lists.

Canonical automatic printing output has been activated.

> p=1+x+x^2;

> mode=off;

> p;

1 + x + x^2

> canonical=mode;

Canonical automatic printing output has been deactivated.

> p;

1 + x * (1 + x)

See also: on, autosimplify, canonical, timing, fullparentheses, midpointmode, rationalmode, roundingwarnings, timing, dieonerrormode

Go back to the list of commands- on is a special value used to activate certain functionnalities of Sollya.
- As any value it can be affected to a variable and stored in lists.

> mode=on;

> p;

1 + x * (1 + x)

> canonical=mode;

Canonical automatic printing output has been activated.

> p;

1 + x + x^2

See also: off, autosimplify, canonical, timing, fullparentheses, midpointmode, rationalmode, roundingwarnings, timing, dieonerrormode

Go back to the list of commands- expr1 and expr2 represent boolean expressions

false

true

- string represents a character sequence

- parse(string) parses the character sequence string containing
an expression built on constants and base functions.

If the character sequence does not contain a well-defined expression, a warning is displayed indicating a syntax error and parse returns a error of type error. - The character sequence to be parsed by parse may contain commands that return expressions, including parse itself. Those commands get executed after the string has been parsed. parse(string) will return the expression computed by the commands contained in the character sequence string.

exp(x)

> print("The string", text, "gives", parse(text));

The string remez(exp(x),5,[-1;1]) gives 1.00004475029055070643077052482053398765426158966754 + x * (1.00003834652983970735244541124504033817544233075356 + x * (0.49919698262882986492168824494240374771969012861297 + x * (0.166424656075155194415920597322727380932279602909199 + x * (4.37936963873280470271257566207183496659575464236489e-2 + x * 8.7381910388065551140158420278330960479960476713376e-3))))

> parse("5 + * 3");

Warning: syntax error, unexpected MULTOKEN. Will try to continue parsing (expecting ";"). May leak memory.

Warning: the string "5 + * 3" could not be parsed by the miniparser.

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

error

- The use of perturb in the command externalplot enables the addition
of some random noise around each sampling point in externalplot.

See externalplot for details.

> bashexecute("gcc -shared -o externalplotexample externalplotexample.o -lgmp -lmpfr");

> externalplot("./externalplotexample",relative,exp(x),[-1/2;1/2],12,perturb);

- pi is the constant Pi, defined as half the period of sine and cosine.
- In Sollya, pi is considered a 0-ary function. This way, the constant is not evaluated at the time of its definition but at the time of its use. For instance, when you define a constant or a function relating to Pi, the current precision at the time of the definition does not matter. What is important is the current precision when you evaluate the function or the constant value.
- Remark that when you define an interval, the bounds are first evaluated and then the interval is defined. In this case, pi will be evaluated as any other constant value at the definition time of the interval, thus using the current precision at this time.

> a = 2*pi;

> a;

Warning: rounding has happened. The value displayed is a faithful rounding of the true result.

6.283

> prec=20!;

> a;

Warning: rounding has happened. The value displayed is a faithful rounding of the true result.

6.283187

Display mode is binary numbers.

> prec=12!;

> d = [pi; 5];

> d;

[1.1001001_2 * 2^(1);1.01_2 * 2^(2)]

> prec=20!;

> d;

[1.1001001_2 * 2^(1);1.01_2 * 2^(2)]

- f1, ..., fn are functions to be plotted.
- L is a list of functions to be plotted.
- I is the interval where the functions have to be plotted.
- name is a string representing the name of a file.

- This command plots one or several functions f1, ... ,fn on an interval I. Functions can be either given as parameters of plot or as a list L which elements are functions. The functions are drawn on the same plot with different colors.
- If L contains an element that is not a function (or a constant), an error occurs.
- plot relies on the value of global variable points. Let n be the value of this variable. The algorithm is the following: each function is evaluated at n evenly distributed points in I. At each point, the computed value is a faithful rounding of the exact value with a sufficiently high precision. Each point is finally plotted. This should avoid numerical artefacts such as critical cancellations.
- You can save the function plot either as a data file or as a postscript file.
- If you use argument file with a string name, Sollya will save a data file called name.dat and a gnuplot directives file called name.p. Invoking gnuplot on name.p will plot the data stored in name.dat.
- If you use argument postscript with a string name, Sollya will save a postscript file called name.eps representing your plot.
- If you use argument postscriptfile with a string name, Sollya will produce the corresponding name.dat, name.p and name.eps.
- This command uses gnuplot to produce the final plot. If your terminal is not graphic (typically if you use Sollya through ssh without -X) gnuplot should be able to detect that and produce an ASCII-art version on the standard output. If it is not the case, you can either store the plot in a postscript file to view it locally, or use asciiplot command.
- If every function is constant, plot will not plot them but just display their value.
- If the interval is reduced to a single point, plot will just display the value of the functions at this point.

1

0.84147098480789650665250232163029899962256306079837

0.84147098480789650665250232163029899962256306079837

0.54030230586813971740093660744297660373231042061792

- function1 and function2 represent functions
- interval1 and interval2 represent intervals (ranges)
- constant represents a constant or constant expression

- + represents the addition (function) on reals. The expression function1 + function2 stands for the function composed of the addition function and the two functions function1 and function2.
- + can be used for interval arithmetic on intervals (ranges). + will evaluate to an interval that safely encompasses all images of the addition function with arguments varying in the given intervals. Any combination of intervals with intervals or constants (resp. constant expressions) is supported. However, it is not possible to represent families of functions using an interval as one argument and a function (varying in the free variable) as the other one.

3

2 + x

x * 2

cos(x) + exp(x)

[4;6]

> [1;2] + 17;

[18;19]

> 13 + [-4;17];

[9;30]

- n represents the number of points

- points is a global variable. Its value represents the number of points used in numerical algorithms of Sollya (namely dirtyinfnorm, dirtyintegral, dirtyfindzeros, plot).

> points=10;

The number of points has been set to 10.

> dirtyfindzeros(f, [0;1]);

[|0, 0.318309886183790671537767526745028724068919291480918|]

> points=100;

The number of points has been set to 100.

> dirtyfindzeros(f, [0;1]);

[|0, 2.4485375860291590118289809749617594159147637806224e-2, 3.97887357729738339422209408431285905086149114351147e-2, 4.54728408833986673625382181064326748669884702115589e-2, 5.3051647697298445256294587790838120678153215246819e-2, 6.3661977236758134307553505349005744813783858296183e-2, 7.9577471545947667884441881686257181017229822870229e-2, 0.106103295394596890512589175581676241356306430493638, 0.159154943091895335768883763372514362034459645740459, 0.318309886183790671537767526745028724068919291480918|]

- postscriptfile is a special value used in commands plot and externalplot to save the result of the command in a data file and a postscript file.
- As any value it can be affected to a variable and stored in lists.

> name="plotSinCos";

> plot(sin(x),0,cos(x),[-Pi,Pi],savemode, name);

- postscript is a special value used in commands plot and externalplot to save the result of the command in a postscript file.
- As any value it can be affected to a variable and stored in lists.

> name="plotSinCos";

> plot(sin(x),0,cos(x),[-Pi,Pi],savemode, name);

- function1 and function2 represent functions
- interval1 and interval2 represent intervals (ranges)
- constant represents a constant or constant expression

- ^ represents the power (function) on reals. The expression function1 ^ function2 stands for the function composed of the power function and the two functions function1 and function2, where function1 is the base and function2 the exponent. If function2 is a constant integer, ^ is defined on negative values of function1. Otherwise ^ is defined as exp(y * log(x)).
- Note that whenever several ^ are composed, the priority goes to the last ^. This corresponds to the natural way of thinking when a tower of powers is written on a paper. Thus, 2^3^5 is interpreted as 2^(3^5).
- ^ can be used for interval arithmetic on intervals (ranges). ^ will evaluate to an interval that safely encompasses all images of the power function with arguments varying in the given intervals. If the intervals given contain points where the power function is not defined, infinities and NaNs will be produced in the output interval. Any combination of intervals with intervals or constants (resp. constant expressions) is supported. However, it is not possible to represent families of functions using an interval as one argument and a function (varying in the free variable) as the other one.

25

x^2

4.1152263374485596707818930041152263374485596707818e-3

@NaN@

sin(x)^exp(x) * ((cos(x) * exp(x)) / sin(x) + exp(x) * log(sin(x)))

1.4134776518227074636666380005943348126619871175005e73

> (2^3)^5;

32768

> 2^(3^5);

1.4134776518227074636666380005943348126619871175005e73

[1;1.60000000000000000000000000000000000000000000000007e1]

> [1;2] ^ 17;

[1;131072]

> 13 ^ [-4;17];

[3.501277966457757081334687160813696999404782745702e-5;8650415919381337933]

- x is a dyadic number.

> precision(a);

19

> m=mantissa(a);

> ceil(log2(m));

19

> prec=50;

The precision has been set to 50 bits.

> dirtyinfnorm(exp(x),[1;2]);

1.110110001110011001001011100011010100110111011011_2 * 2^(2)

> prec=100;

The precision has been set to 100 bits.

> dirtyinfnorm(exp(x),[1;2]);

1.11011000111001100100101110001101010011011101101011011100110000110011101000111011101000100000011011_2 * 2^(2)

- x is an object of any type.
- L is a list (possibly empty).

- .: adds the element x at the beginning of the list L.
- Note that since x may be of any type, it can be in particular a list.

[|1, 2, 3, 4|]

[|[|1, 2, 3|], 4, 5, 6|]

[|1|]

- constant represents a constant

- Prints a constant value as a hexadecimal number on 16 hexadecimal
digits. The hexadecimal number represents the integer equivalent to
the 64 bit memory representation of the constant considered as a
double precision number.

If the constant value does not hold on a double precision number, it is first rounded to the nearest double precision number before displayed. A warning is displayed in this case.

0x4008000000000000

> verbosity = 1!;

> printdouble(exp(5));

Warning: the given expression is not a constant but an expression to evaluate. A faithful evaluation will be used.

Warning: rounding down occurred before printing a value as a double.

0x40628d389970338f

- polynomial represents the polynomial to be printed

- The command printexpansion prints the polynomial polynomial in Horner form
writing its coefficients as expansions of double precision
numbers. The double precision numbers themselves are displayed in
hexadecimal memory notation (see printdouble).

If some of the coefficients of the polynomial polynomial are not floating-point constants but constant expressions, they are evaluated to floating-point constants using the global precision prec. If a rounding occurs in this evaluation, a warning is displayed.

If the exponent range of double precision is not sufficient to display all the mantissa bits of a coefficient, the coefficient is displayed rounded and a warning is displayed.

If the argument polynomial does not a polynomial, nothing but a warning or a newline is displayed. Constants can be displayed using printexpansion since they are polynomials of degree 0.

0x3ff0000000000000 + x * (0x3ff0000000000000 + x * (0x3fe0000000000000 + x * ((0x3fc5555555555555 + 0x3c65555555555555) + x * ((0x3fa5555555555555 + 0x3c45555555555555) + x * (0x3f81111111111111 + 0x3c01111111111111)))))

(0x3ff0002eec908ce9 + 0xbc7df99eb225af5b + 0xb8d55834b08b1f18) + x * ((0x3ff0002835917719 + 0x3c6d82c073b25ebf + 0xb902cf062b54b7b6 + 0x35b0000000000000) + x * ((0x3fdff2d7e6a9c5e9 + 0xbc7b09a95b0d520f + 0xb915b639add55731 + 0x35a8000000000000) + x * ((0x3fc54d67338ba09f + 0x3c4867596d0631cf + 0xb8ef0756bdb4af6e) + x * ((0x3fa66c209b825167 + 0x3c45ec5b6655b076 + 0xb8d8c125286400bc) + x * (0x3f81e55425e72ab4 + 0x3c263b25a1bf597b + 0xb8c843e0401dadd0 + 0x3540000000000000)))))

> prec = 3500!;

> printexpansion(pi);

(0x400921fb54442d18 + 0x3ca1a62633145c07 + 0xb92f1976b7ed8fbc + 0x35c4cf98e804177d + 0x32631d89cd9128a5 + 0x2ec0f31c6809bbdf + 0x2b5519b3cd3a431b + 0x27e8158536f92f8a + 0x246ba7f09ab6b6a9 + 0xa0eedd0dbd2544cf + 0x1d779fb1bd1310ba + 0x1a1a637ed6b0bff6 + 0x96aa485fca40908e + 0x933e501295d98169 + 0x8fd160dbee83b4e0 + 0x8c59b6d799ae131c + 0x08f6cf70801f2e28 + 0x05963bf0598da483 + 0x023871574e69a459 + 0x8000000005702db3 + 0x8000000000000000)

Warning: the expansion is not complete because of the limited exponent range of double precision.

Warning: rounding occurred while printing.

See also: printdouble, horner, print, prec, remez, taylor, roundcoefficients, fpminimax, implementpoly

Go back to the list of commands- expr represents an expression
- filename represents a character sequence indicating a file name

- print(expr1,...,exprn) prints the expressions expr1 through
exprn separated by spaces and followed by a newline.

If a second argument filename is given after a single ">", the displaying is not output on the standard output of Sollya but if in the file filename that get newly created or overwritten. If a double ">>" is given, the output will be appended to the file filename.

The global variables display, midpointmode and fullparentheses have some influence on the formatting of the output (see display, midpointmode and fullparentheses).

Remark that if one of the expressions expri given in argument is of type string, the character sequence expri evaluates to is displayed. However, if expri is of type list and this list contains a variable of type string, the expression for the list is displayed, i.e. all character sequences get displayed surrounded by double quotes ("). Nevertheless, escape sequences used upon defining character sequences are interpreted immediately.

x + 2 + exp(sin(x))

> print("Hello","world");

Hello world

> print("Hello","you", 4 + 3, "other persons.");

Hello you 7 other persons.

Hello

> print([|"Hello"|]);

[|"Hello"|]

> s = "Hello";

> print(s,[|s|]);

Hello [|"Hello"|]

> t = "Hello\tyou";

> print(t,[|t|]);

Hello you [|"Hello\tyou"|]

> readfile("foo.sol");

x + 2 + exp(sin(x))

Display mode is decimal numbers.

> a = evaluate(sin(pi * x), 0.25);

> b = evaluate(sin(pi * x), [0.25; 0.25 + 1b-50]);

> print(a);

0.70710678118654752440084436210484903928483593768847

> display = binary;

Display mode is binary numbers.

> print(a);

1.01101010000010011110011001100111111100111011110011001001000010001011001011111011000100110110011011101010100101010111110100111110001110101101111011000001011101010001_2 * 2^(-1)

> display = hexadecimal;

Display mode is hexadecimal numbers.

> print(a);

0xb.504f333f9de6484597d89b3754abe9f1d6f60ba88p-4

> display = dyadic;

Display mode is dyadic numbers.

> print(a);

33070006991101558613323983488220944360067107133265b-165

> display = powers;

Display mode is dyadic numbers in integer-power-of-2 notation.

> print(a);

33070006991101558613323983488220944360067107133265 * 2^(-165)

> display = decimal;

Display mode is decimal numbers.

> midpointmode = off;

Midpoint mode has been deactivated.

> print(b);

[0.70710678118654752440084436210484903928483593768844;0.70710678118654949743721782517557347782646274417048]

> midpointmode = on;

Midpoint mode has been activated.

> print(b);

0.7071067811865~4/5~

> display = dyadic;

Display mode is dyadic numbers.

> print(b);

[2066875436943847413332748968013809022504194195829b-161;16535003495550825444196237019385936414432675156571b-164]

> display = decimal;

Display mode is decimal numbers.

> autosimplify = off;

Automatic pure tree simplification has been deactivated.

> fullparentheses = off;

Full parentheses mode has been deactivated.

> print(x + x * ((x + 1) + 1));

x + x * (x + 1 + 1)

> fullparentheses = on;

Full parentheses mode has been activated.

> print(x + x * ((x + 1) + 1));

x + (x * ((x + 1) + 1))

See also: write, printexpansion, printdouble, printsingle, printxml, readfile, autosimplify, display, midpointmode, fullparentheses, evaluate, rationalmode

Go back to the list of commands- constant represents a constant

- Prints a constant value as a hexadecimal number on 8 hexadecimal
digits. The hexadecimal number represents the integer equivalent to
the 32 bit memory representation of the constant considered as a
single precision number.

If the constant value does not hold on a single precision number, it is first rounded to the nearest single precision number before it is displayed. A warning is displayed in this case.

0x40400000

> verbosity = 1!;

> printsingle(exp(5));

Warning: the given expression is not a constant but an expression to evaluate. A faithful evaluation will be used.

Warning: rounding down occurred before printing a value as a simple.

0x431469c5

- expr represents a functional expression
- filename represents a character sequence indicating a file name

- printxml(expr) prints the functional expression expr as a tree of
MathML Content Definition Markups. This XML tree can be re-read in
external tools or by usage of the readxml command.

If a second argument filename is given after a single ">", the MathML tree is not output on the standard output of Sollya but if in the file filename that get newly created or overwritten. If a double ">" is given, the output will be appended to the file filename.

<?xml version="1.0" encoding="UTF-8"?>

<!-- generated by sollya: http://sollya.gforge.inria.fr/ -->

<!-- syntax: printxml(...); example: printxml(x^2-2*x+5); -->

<?xml-stylesheet type="text/xsl" href="http://sollya.gforge.inria.fr/mathmlc2p-web.xsl"?>

<?xml-stylesheet type="text/xsl" href="mathmlc2p-web.xsl"?>

<!-- This stylesheet allows direct web browsing of MathML-c XML files (http:// or file://) -->

<math xmlns="http://www.w3.org/1998/Math/MathML">

<semantics>

<annotation-xml encoding="MathML-Content">

<lambda>

<bvar><ci> x </ci></bvar>

<apply>

<apply>

<plus/>

<apply>

<plus/>

<ci> x </ci>

<cn type="integer" base="10"> 2 </cn>

</apply>

<apply>

<exp/>

<apply>

<sin/>

<ci> x </ci>

</apply>

</apply>

</apply>

</apply>

</lambda>

</annotation-xml>

<annotation encoding="sollya/text">(x + 1b1) + exp(sin(x))</annotation>

</semantics>

</math>

- identifier represents the name of the procedure to be defined and assigned
- formal parameter1, formal parameter2 through formal parameter n represent identifiers used as formal parameters
- formal list parameter represents an identifier used as a formal parameter for the list of an arbitrary number of parameters
- procedure body represents the imperative statements in the body of the procedure
- expression represents the expression procedure shall evaluate to

- The procedure keyword allows for defining and assigning procedures in the Sollya language. It is an abbreviation to a procedure definition using proc with the same formal parameters, procedure body and return-expression followed by an assignment of the procedure (object) to the identifier identifier. In particular, all rules concerning local variables declared using the var keyword apply for procedure.

> succ(5);

6

> 3 + succ(0);

4

> succ;

proc(n)

{

nop;

return (n) + (1);

}

> myprint("Lyon","Nancy","Beaverton","Coye-la-Foret","Amberg","Nizhny Novgorod","Cluj-Napoca");

Lyon

Nancy

Beaverton

Coye-la-Foret

Amberg

Nizhny Novgorod

Cluj-Napoca

- formal parameter1, formal parameter2 through formal parameter n represent identifiers used as formal parameters
- formal list parameter represents an identifier used as a formal parameter for the list of an arbitrary number of parameters
- procedure body represents the imperative statements in the body of the procedure
- expression represents the expression proc shall evaluate to

- The proc keyword allows for defining procedures in the Sollya language. These procedures are common Sollya objects that can be applied to actual parameters after definition. Upon such an application, the Sollya interpreter applies the actual parameters to the formal parameters formal parameter1 through formal parameter n (resp. builds up the list of arguments and applies it to the list formal list parameter) and executes the procedure body. The procedure applied to actual parameters evaluates then to the expression expression in the return statement after the procedure body or to void, if no return statement is given (i.e. a return void statement is implicitly given).
- Sollya procedures defined by proc have no name. They can be bound to an identifier by assigning the procedure object a proc expression produces to an identifier. However, it is possible to use procedures without giving them any name. For instance, Sollya procedures, i.e. procedure objects, can be elements of lists. They can even be given as an argument to other internal Sollya procedures. See also procedure on this subject.
- Upon definition of a Sollya procedure using proc, no type check
is performed. More precisely, the statements in procedure body are
merely parsed but not interpreted upon procedure definition with
proc. Type checks are performed once the procedure is applied to
actual parameters or to void. At this time, if the procedure was
defined using several different formal parameters formal parameter 1
through formal parameter n, it is checked whether the number of
actual parameters corresponds to the number of formal parameters. If
the procedure was defined using the syntax for a procedure with an
arbitrary number of parameters by giving a formal list parameter,
the number of actual arguments is not checked but only a list <formal
list parameter> of appropriate length is built up. Type checks are
further performed upon execution of each statement in procedure body
and upon evaluation of the expression expression to be returned.

Procedures defined by proc containing a quit or restart command cannot be executed (i.e. applied). Upon application of a procedure, the Sollya interpreter checks beforehand for such a statement. If one is found, the application of the procedure to its arguments evaluates to error. A warning is displayed. Remark that in contrast to other type or semantic correctness checks, this check is really performed before interpreting any other statement in the body of the procedure. - Through the var keyword it is possible to declare local variables and thus to have full support of recursive procedures. This means a procedure defined using proc may contain in its procedure body an application of itself to some actual parameters: it suffices to assign the procedure (object) to an identifier with an appropriate name.
- Sollya procedures defined using proc may return other procedures. Further procedure body may contain assignments of locally defined procedure objects to identifiers. See var for the particular behaviour of local and global variables.
- The expression expression returned by a procedure is evaluated with regard to Sollya commands, procedures and external procedures. Simplification may be performed. However, an application of a procedure defined by proc to actual parameters evaluates to the expression expression that may contain the free global variable or that may be composed.

> succ(5);

6

> 3 + succ(0);

4

> succ;

proc(n)

{

nop;

return (n) + (1);

}

> add(5,6);

11

> add;

proc(m, n)

{

var res;

res := (m) + (n);

return res;

}

> verbosity = 1!;

> add(3);

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

error

> add(true,false);

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

Warning: the given expression or command could not be handled.

error

> succ(5);

6

> succ(x);

1 + x

> hey();

Hello world.

> print(hey());

Hello world.

void

> hey;

proc()

{

print("Hello world.");

return void;

}

> fac(5);

120

> fac(11);

39916800

> fac;

proc(n)

{

var res;

if (n) == (0) then

res := 1

else

res := (n) * (fac((n) - (1)));

return res;

}

> (myprocs[0])(5,6);

11

> (myprocs[1])(5,6);

-1

> succ = proc(n) { return n + 1; };

> pred = proc(n) { return n - 1; };

> applier = proc(p,n) { return p(n); };

> applier(succ,5);

6

> applier(pred,5);

4

> myquit = proc(n) { print(n); quit; };

> myquit;

proc(n)

{

print(n);

quit;

return void;

}

> myquit(5);

Warning: a quit or restart command may not be part of a procedure body.

The procedure will not be executed.

Warning: an error occurred while executing a procedure.

Warning: the given expression or command could not be handled.

error

> printsucc(5);

Successor of 5 is 6

> makeadd(4);

n = 4

proc(m, n)

{

nop;

return (n) + (m);

}

> (makeadd(4))(5,6);

n = 4

11

> sumall;

proc(L = ...)

{

var acc, i;

acc = 0;

for i in L do

acc = (acc) + (i);

return acc;

}

> sumall();

0

> sumall(2);

2

> sumall(2,5);

7

> sumall(2,5,7,9,16);

39

> sumall @ [|1,...,10|];

55

- quad is both a function and a constant.
- As a function, it rounds its argument to the nearest IEEE 754 quad precision (i.e. IEEE754-2008 binary128) number. Subnormal numbers are supported as well as standard numbers: it is the real rounding described in the standard.
- As a constant, it symbolizes the quad precision format. It is used in contexts when a precision format is necessary, e.g. in the commands round and roundcoefficients. It is not supported for implementpoly. See the corresponding help pages for examples.

> QD(0.1);

1.100110011001100110011001100110011001100110011001100110011001100110011001100110011001100110011001100110011001101_2 * 2^(-4)

> QD(4.17);

1.000010101110000101000111101011100001010001111010111000010100011110101110000101000111101011100001010001111010111_2 * 2^(2)

> QD(1.011_2 * 2^(-16493));

1.1_2 * 2^(-16493)

See also: halfprecision, single, double, doubleextended, doubledouble, tripledouble, roundcoefficients, implementpoly, fpminimax, round, printsingle

Go back to the list of commands- The command quit, when executed, stops the execution of a Sollya
script and leaves the Sollya interpreter unless the quit command
is executed in a Sollya script read into a main Sollya script by
execute or #include.

Upon exiting the Sollya interpreter, all state is thrown away, all memory is deallocated, all bound libraries are unbound and the temporary files produced by plot and externalplot are deleted.

If the quit command does not lead to a halt of the Sollya interpreter, a warning is displayed.

- range represents the range type for declarations
of external procedures by means of externalproc.

Remark that in contrast to other indicators, type indicators like range cannot be handled outside the externalproc context. In particular, they cannot be assigned to variables.

- x is a number to approximate.
- n is a integer (representing a format).

- rationalapprox(x,n) returns a constant function of the form a/b where a and b are integers. The value a/b is an approximation of x. The quality of this approximation is determined by the parameter n that indicates the number of correct bits that a/b should have.
- The command is not safe in the sense that it is not ensured that the error between a/b and x is less than 2^(-n).
- The following algorithm is used: x is first rounded downwards and upwards to a format of n bits, thus obtaining an interval [xl,xu]. This interval is then developped into a continued fraction as far as the representation is the same for every elements of [xl,xu]. The corresponding fraction is returned.
- Since rational numbers are not a primitive object of Sollya, the fraction is returned as a constant function. This can be quite amazing, because Sollya immediately simplifies a constant function by evaluating it when the constant has to be displayed. To avoid this, you can use print (that displays the expression representing the constant and not the constant itself) or the commands numerator and denominator.

> pi50 = rationalapprox(Pi,50);

> pi100 = rationalapprox(Pi,100);

> print( pi10, ": ", simplify(floor(-log2(abs(pi10-Pi)/Pi))), "bits." );

22 / 7 : 11 bits.

> print( pi50, ": ", simplify(floor(-log2(abs(pi50-Pi)/Pi))), "bits." );

90982559 / 28960648 : 50 bits.

> print( pi100, ": ", simplify(floor(-log2(abs(pi100-Pi)/Pi))), "bits." );

4850225745369133 / 1543874804974140 : 101 bits.

> b=rationalapprox(a,4);

> numerator(b); denominator(b);

1

10

> print(simplify(floor(-log2(abs((b-a)/a)))), "bits.");

166 bits.

- activation value controls if rational arithmetic should be used or not

- rationalmode is a global variable. When its value is off, which is the default, Sollya will not use rational arithmetic to simplify expressions. All computations, including the evaluation of constant expressions given on the Sollya prompt, will be performed using floating-point and interval arithmetic. Constant expressions will be approximated by floating-point numbers, which are in most cases faithful roundings of the expressions, when shown at the prompt.
- When the value of the global variable rationalmode is on, Sollya will use rational arithmetic when simplifying expressions. Constant expressions, given at the Sollya prompt, will be simplified to rational numbers and displayed as such when they are in the set of the rational numbers. Otherwise, flaoting-point and interval arithmetic will be used to compute a floating-point approximation, which is in most cases a faithful rounding of the constant expression.

> 19/17 + 3/94;

1.1495619524405506883604505632040050062578222778473

> rationalmode=on!;

> 19/17 + 3/94;

1837 / 1598

> exp(19/17 + 3/94);

3.15680977395514136754709208944824276340328162814418

> rationalmode=on!;

> exp(19/17 + 3/94);

3.15680977395514136754709208944824276340328162814418

> round(Pi,20,RD);

1.1001001000011111101_2 * 2^(1)

- filename represents a character sequence indicating a file name

> t = readfile("myfile.txt");

> t;

Hello world

> readfile("afile.txt");

Warning: the file "afile.txt" could not be opened for reading.

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

error

- filename represents a character sequence indicating a file name

- readxml(filename) reads the first occurrence of a lambda
application with one bounded variable on applications of the supported
basic functions in file filename and returns it as a Sollya
functional expression.

If the file filename does not contain a valid MathML-Content tree, readxml tries to find an "annotation encoding" markup of type "sollya/text". If this annotation contains a character sequence that can be parsed by parse, readxml returns that expression. Otherwise readxml displays a warning and returns an error variable of type error.

2 + x + exp(sin(x))

- The use of relative in the command externalplot indicates that during
plotting in externalplot a relative error is to be considered.

See externalplot for details. - Used with fpminimax, relative indicates that fpminimax must try to minimize
the relative error.

See fpminimax for details. - When given in argument to supnorm, relative indicates that a relative error
is to be considered for supremum norm computation.

See supnorm for details.

> bashexecute("gcc -shared -o externalplotexample externalplotexample.o -lgmp -lmpfr");

> externalplot("./externalplotexample",absolute,exp(x),[-1/2;1/2],12,perturb);

- f is the function to be approximated
- n is the degree of the polynomial that must approximate f
- L is a list of monomials that can be used to represent the polynomial that must approximate f
- range is the interval where the function must be approximated
- w (optional) is a weight function. Default is 1.
- quality (optional) is a parameter that controls the quality of the returned polynomial p, with respect to the exact minimax p*. Default is 1e-5.

- remez computes an approximation of the function f with respect to the weight function w on the interval range. More precisely, it searches a polynomial p such that ||p*w-f|| is (almost) minimal among all polynomials p of a certain form. The norm is the infinity norm, e.g. ||g|| = max {|g(x)|, x in range}.
- If w=1 (the default case), it consists in searching the best polynomial approximation of f with respect to the absolute error. If f=1 and w is of the form 1/g, it consists in searching the best polynomial approximation of g with respect to the relative error.
- If n is given, the polynomial p is searched among the polynomials with degree not greater than n. If L is given, the polynomial p is searched as a linear combination of monomials X^k where k belongs to L. L may contain ellipses but cannot be end-elliptic.
- The polynomial is obtained by a convergent iteration called Remez' algorithm. The algorithm computes a sequence p1, ..., pk, ... such that ek = ||pk*w-f|| converges towards the optimal value e. The algorithm is stopped when the relative error between ek and e is less than quality.
- Note: the algorithm may not converge in certain cases. Moreover, it may converge towards a polynomial that is not optimal. These cases correspond to the cases when the Haar condition is not fulfilled. See [Cheney - Approximation theory] for details.

> degree(p);

5

> dirtyinfnorm(p-exp(x),[0;1]);

1.12956984638214536849843017679626063762687501534126e-6

> canonical=on!;

> p;

0.99999999994393749280444571988532724907643631727381 + -0.49999999571557467737204931630836834563663039748203 * x^2 + 4.16666132335010905188253972212748718651775241902969e-2 * x^4 + -1.38865291475286141707180658383176799662601691348739e-3 * x^6 + 2.437267919111162694221738667927916761689966804242e-5 * x^8

> p2 = remez(exp(x),5,[0;1],default,1e-10);

> p3 = remez(exp(x),5,[0;1],default,1e-15);

> dirtyinfnorm(p1-exp(x),[0;1]);

1.12956984638214536849843017679626063762687501534126e-6

> dirtyinfnorm(p2-exp(x),[0;1]);

1.12956980227478687332174207517728389861926659249056e-6

> dirtyinfnorm(p3-exp(x),[0;1]);

1.12956980227478687332174207517728389861926659249056e-6

- ident1 is the current name of the free variable.
- ident2 is a fresh name.

- rename permits a change of the name of the free variable. Sollya can handle only one free variable at a time. The first time in a session that an unbound name is used in a context where it can be interpreted as a free variable, the name is used to represent the free variable of Sollya. In the following, this name can be changed using rename.
- Be careful: if ident2 has been set before, its value will be lost. Use the command isbound to know if ident2 is already used or not.
- If ident1 is not the current name of the free variable, an error occurs.
- If rename is used at a time when the name of the free variable has not been defined, ident1 is just ignored and the name of the free variable is set to ident2.

> f;

sin(x)

> rename(x,y);

> f;

sin(y)

> f=sin(x);

> rename(x,a);

> a;

a

> f;

sin(a)

> f=sin(x);

> rename(y,z);

Warning: the current free variable is named "x" and not "y". Can only rename the free variable.

The last command will have no effect.

> isbound(x);

false

> isbound(y);

true

- The command restart brings Sollya back to its initial state. All
current state is abandoned, all libraries unbound and all memory freed.

The restart command has no effect when executed inside a Sollya script read into a main Sollya script using execute. It is executed in a Sollya script included by a #include macro.

Using the restart command in nested elements of imperative programming like for or while loops is possible. Since in most cases abandoning the current state of Sollya means altering a loop invariant, warnings for the impossibility of continuing a loop may follow unless the state is rebuilt.

exp(x)

> a = 3;

> restart;

The tool has been restarted.

> print(x);

x

> a;

Warning: the identifier "a" is neither assigned to, nor bound to a library function nor external procedure, nor equal to the current free variable.

Will interpret "a" as "x".

x

exp(x)

> for i from 1 to 10 do {

print(i);

if (i == 5) then restart;

};

1

2

3

4

5

The tool has been restarted.

Warning: the tool has been restarted inside a for loop.

The for loop will no longer be executed.

exp(x)

> a = 3;

> for i from 1 to 10 do {

print(i);

if (i == 5) then {

restart;

i = 7;

};

};

1

2

3

4

5

The tool has been restarted.

8

9

10

> print(x);

x

> a;

Warning: the identifier "a" is neither assigned to, nor bound to a library function nor external procedure, nor equal to the current free variable.

Will interpret "a" as "x".

x

- expression represents the expression to be returned

- The keyword return allows for returning the (evaluated) expression
expression at the end of a begin-end-block ({}-block) used as a
Sollya procedure body. See proc for further details concerning
Sollya procedure definitions.

Statements for returning expressions using return are only possible at the end of a begin-end-block used as a Sollya procedure body. Only one return statement can be given per begin-end-block. - If at the end of a procedure definition using proc no return statement is given, a return void statement is implicitely added. Procedures, i.e. procedure objects, when printed out in Sollya defined with an implicit return void statement are displayed with this statement explicitly given.

> succ(5);

6

> succ;

proc(n)

{

var res;

res := (n) + (1);

return res;

}

> hey("world");

Hello world

> hey;

proc(s)

{

print("Hello", s);

return void;

}

- L is a list.

[| |]

[|4, 1, 2, 5, 3, 2|]

> round(Pi,20,RN);

1.100100100001111111_2 * 2^(1)

- p is a function. Usually a polynomial.
- L is a list of formats.

- If p is a polynomial and L a list of floating-point formats, roundcoefficients(p,L) rounds each coefficient of p to the corresponding format in L.
- If p is not a polynomial, roundcoefficients does not do anything.
- If L contains other elements than HP, halfprecision, SG, single, D, double, DE, doubleextended, DD, doubledouble, QD, quad, TD and tripledouble, an error occurs.
- The coefficients in p corresponding to X^i is rounded to the format L[i]. If L does not contain enough elements (e.g. if length(L) < degree(p)+1), a warning is displayed. However, the coefficients corresponding to an element of L are rounded. The trailing coefficients (that do not have a corresponding element in L) are kept with their own precision. If L contains too much elements, the trailing useless elements are ignored. In particular L may be end-elliptic in which case roundcoefficients has the natural behavior.

> display=binary!;

> roundcoefficients(p,[|DD,D,D|]);

1.010110111111000010101000101100010100010101110110100101010011010101011111101110001010110001000000010011101_2 * 2^(1) + x * (1.110110001110011001001011100011010100110111011010111_2 * 2^(2) + x * (1.010000010101111001011011111101101111101100010000011_2 * 2^(4)))

> roundcoefficients(p,[|DD,D...|]);

1.010110111111000010101000101100010100010101110110100101010011010101011111101110001010110001000000010011101_2 * 2^(1) + x * (1.110110001110011001001011100011010100110111011010111_2 * 2^(2) + x * (1.010000010101111001011011111101101111101100010000011_2 * 2^(4)))

> display=binary!;

> f;

sin(x * (1.0101101111110000101010001011000101000101011101101001010100110101010111111011100010101100010000000100111001111010011110011110001110110001011100111000101100000111101_2 * 2^(1)))

> roundcoefficients(f,[|D...|]);

sin(x * (1.0101101111110000101010001011000101000101011101101001010100110101010111111011100010101100010000000100111001111010011110011110001110110001011100111000101100000111101_2 * 2^(1)))

> verbosity=1!;

> display=binary!;

> roundcoefficients(p,[|DD,D|]);

Warning: the number of the given formats does not correspond to the degree of the given polynomial.

Warning: the 0th coefficient of the given polynomial does not evaluate to a floating-point constant without any rounding.

Will evaluate the coefficient in the current precision in floating-point before rounding to the target format.

Warning: the 1th coefficient of the given polynomial does not evaluate to a floating-point constant without any rounding.

Will evaluate the coefficient in the current precision in floating-point before rounding to the target format.

Warning: rounding may have happened.

1.010110111111000010101000101100010100010101110110100101010011010101011111101110001010110001000000010011101_2 * 2^(1) + x * (1.110110001110011001001011100011010100110111011010111_2 * 2^(2) + x * (1.01000001010111100101101111110110111110110001000001011111001011010100101111011111110001010011011101000100110000111010001110010000010110000101100000111001011100101001_2 * 2^(4)))

See also: halfprecision, single, double, doubleextended, doubledouble, quad, tripledouble, fpminimax, remez, implementpoly, subpoly

Go back to the list of commands- range represents a range in which an exact value lies

- Let range be a range of values, determined by some approximation
process, safely bounding an unknown value v. The command
roundcorrectly(range) determines a precision such that for this precision,
rounding to the nearest any value in range yields to the same
result, i.e. to the correct rounding of v.

If no such precision exists, a warning is displayed and roundcorrectly evaluates to NaN.

1.01_2

> printbinary(roundcorrectly([1.00001_2; 1.001_2]));

1_2

@NaN@

- activation value controls if warnings should be shown or not

- roundingwarnings is a global variable. When its value is on, warnings are emitted in appropriate verbosity modes (see verbosity) when roundings occur. When its value is off, these warnings are suppressed.
- This mode depends on a verbosity of at least 1. See verbosity for more details.
- Default is on when the standard input is a terminal and off when Sollya input is read from a file.

> roundingwarnings = on;

Rounding warning mode has been activated.

> exp(0.1);

Warning: Rounding occurred when converting the constant "0.1" to floating-point with 165 bits.

If safe computation is needed, try to increase the precision.

Warning: rounding has happened. The value displayed is a faithful rounding of the true result.

1.1051709180756476248117078264902466682245471947375

> roundingwarnings = off;

Rounding warning mode has been deactivated.

> exp(0.1);

1.1051709180756476248117078264902466682245471947375

- x is a constant to be rounded.
- n is the precision of the target format.
- format is the name of a supported floating-point format.
- mode is the desired rounding mode.

- If used with an integer parameter n, round(x,n,mode) rounds x to a floating-point number with precision n, according to rounding-mode mode.
- If used with a format parameter format, round(x,format,mode) rounds x to a floating-point number in the floating-point format format, according to rounding-mode mode.
- Subnormal numbers are handled for the case when format is one of halfprecision, single, double, doubleextended, doubledouble, quad or tripledouble. Otherwise, when format is an integer, round does not take any exponent range into consideration, i.e. typically uses the full exponent range of the underlying MPFR library.

> round(Pi,20,RN);

1.100100100001111111_2 * 2^(1)

0x417709348c0ea4f9

> printdouble(D(exp(17)));

0x417709348c0ea4f9

> a=2^(-1100);

> round(a,53,RN);

1_2 * 2^(-1100)

> round(a,D,RN);

0

> double(a);

0

See also: RN, RD, RU, RZ, halfprecision, single, double, doubleextended, doubledouble, quad, tripledouble, roundcoefficients, roundcorrectly, printdouble, printsingle, ceil, floor, nearestint

Go back to the list of commands> round(Pi,20,RU);

1.100100100001111111_2 * 2^(1)

> round(Pi,20,RZ);

1.1001001000011111101_2 * 2^(1)

- function represents the function to be considered
- start represents a value around which the search is to be performed
- preimage precision represents the precision (discretization) for the eligible preimage values
- steps represents the binary logarithm (log2) of the number of search steps to be performed
- format represents the format the image of the function is to be rounded to
- error bound represents a upper bound on the relative rounding error when rounding the image
- list of functions represents the functions to be considered
- list of formats represents the respective formats the images of the functions are to be rounded to
- list of error bounds represents a upper bound on the relative rounding error when rounding the image

- The command searchgal searches for a preimage z of function
function or a list of functions list of functions such that
z is a floating-point number with preimage precision
significant mantissa bits and the image y of the function,
respectively each image yi of the functions, rounds to
format format respectively to the corresponding format in list of format
with a relative rounding error less than error bound
respectively the corresponding value in list of error bounds. During
this search, at most 2^steps attempts are made. The search
starts with a preimage value equal to start. This value is then
increased and decreased by 1 ulp in precision preimage precision
until a value is found or the step limit is reached.

If the search finds an appropriate preimage z, searchgal evaluates to a list containing this value. Otherwise, searchgal evaluates to an empty list.

[| |]

> searchgal(log(x),2,53,18,DD,1b-112);

[|2.0000000000384972054234822280704975128173828125|]

> s = searchgal(f,2,53,18,DD,1b-112);

> if (s != [||]) then {

v = s[0];

print("The rounding error is 2^(",evaluate(log2(abs(DD(f)/f - 1)),v),")");

} else print("No value found");

The rounding error is 2^( -1.12106878438809380148206984258358542322113874177832e2 )

[|1.00000000000159494639717649988597258925437927246094|]

- function represents the expression to be simplified

- The command simplifysafe simplifies the expression given in argument
representing the function function. The command simplifysafe does not
endanger the safety of computations even in Sollya's floating-point
environment: the function returned is mathematically equal to the
function function.

Remark that the simplification provided by simplifysafe is not perfect: they may exist simpler equivalent expressions for expressions returned by simplifysafe.

8 + 6 * x

(pi) / 2

(pi) / 2 - (pi) / 4 * 2

- function represents the expression to be simplified

sin(3.14159265358979323846264338327950288419716939937508 * x)

> print(simplify(erf(exp(3) + x * log(4))));

erf(2.00855369231876677409285296545817178969879078385544e1 + x * 1.3862943611198906188344642429163531361510002687205)

> t = erf(0.5);

> s = simplify(erf(0.5));

> prec = 200!;

> t;

0.5204998778130465376827466538919645287364515757579637000588058

> s;

0.52050018310546875

- single is both a function and a constant.
- As a function, it rounds its argument to the nearest IEEE 754 single precision (i.e. IEEE754-2008 binary32) number. Subnormal numbers are supported as well as standard numbers: it is the real rounding described in the standard.
- As a constant, it symbolizes the single precision format. It is used in contexts when a precision format is necessary, e.g. in the commands round and roundcoefficients. In is not supported for implementpoly. See the corresponding help pages for examples.

> SG(0.1);

1.10011001100110011001101_2 * 2^(-4)

> SG(4.17);

1.000010101110000101001_2 * 2^(2)

> SG(1.011_2 * 2^(-1073));

0

See also: halfprecision, double, doubleextended, doubledouble, quad, tripledouble, roundcoefficients, implementpoly, round, printsingle

Go back to the list of commands- sinh is the usual hyperbolic sine function: sinh(x) = (exp(x)-exp(-x))/2.
- It is defined for every real number x.

- sin is the usual sine function.
- It is defined for every real number x.

- L is a list.

[| |]

> sort([|2,3,5,2,1,4|]);

[|1, 2, 2, 3, 4, 5|]

- sqrt is the square root, e.g. the inverse of the function square: sqrt(y) is the unique positive x such that x^2=y.
- It is defined only for x in [0; +Inf].

- string represents the string type for declarations
of external procedures by means of externalproc.

Remark that in contrast to other indicators, type indicators like string cannot be handled outside the externalproc context. In particular, they cannot be assigned to variables.

- polynomial represents the polynomial the coefficients are taken from
- list represents the list of monomials to be taken

- subpoly extracts the coefficients of a polynomial polynomial and builds up a
new polynomial out of those coefficients associated to monomial degrees figuring in
the list list.

If polynomial represents a function that is not a polynomial, subpoly returns 0.

If list is a list that is end-elliptic, let be j the last value explicitly specified in the list. All coefficients of the polynomial associated to monomials greater or equal to j are taken.

> s = subpoly(p,[|1,3,5|]);

> print(p);

1 + x * (1 + x * (0.5 + x * (1 / 6 + x * (1 / 24 + x / 120))))

> print(s);

x * (1 + x^2 * (1 / 6 + x^2 / 120))

> subpoly(p,[|1,3,5...|]);

x * (0.99986632946591986997581285958052433296267358727229 + x^2 * (-0.330304785504861260596093435534236137298206064685038 + x^2 * (0.180159294636523467997437751178959039617773054107393 + x * (-1.21704858321866028906175835649390114311877360260197e-14 + x * (-8.5156350833702702996505336803770858918120961566741e-2 + x * (1.39681284176342339364451388757935358048374878126733e-14 + x * (2.0845114175434561643018447784809880955983412532269e-2 + x * (-5.6810131012579436265697622426011349460288598691964e-15))))))))

0

- f is a function.
- g is a function.
- t is a real number.

- substitute(f, g) produces the function (f o g) : x -> f(g(x)).
- substitute(f, t) is the constant f(t). Note that the constant is represented by its expression until it has been evaluated (exactly the same way as if you type the expression f replacing instances of the free variable by t).
- If f is stored in a variable F, the effect of the commands substitute(F,g) or substitute(F,t) is absolutely equivalent to writing F(g) resp. F(t).

> g=cos(x);

> substitute(f,g);

sin(cos(x))

> f(g);

sin(cos(x))

> f=sin(x);

> substitute(f,a);

0.84147098480789650665250232163029899962256306079837

> f(a);

0.84147098480789650665250232163029899962256306079837

See also: evaluate

Go back to the list of commands- supnorm(p, f, I, errorType, accuracy) tries to compute an interval bound r=[l,u] for the supremum norm of the error function epsilon_absolute=p-f (when errorType evaluates to absolute) or epsilon_relative=p/f-1 (when errorType evaluates to relative), over the interval I, such that sup{|epsilon(x)|, x in I} C r and 0<|u/l-1|< accuracy. If supnorm succeeds in computing a suitable interval r, which it returns, that interval is guaranteed to contain the supremum norm value and to satisfy the required quality. Otherwise, supnorm evaluates to error, displaying a corresponding error message. These failure cases are rare and basically happen only for functions which are too complicated.
- Roughly speaking, supnorm is based on taylorform to obtain a higher degree polynomial approximation for f. This process is coupled with an a posteriori validation of a potential supremum norm upper bound. The validation is based on showing a certain polynomial the problem gets reduced to does not vanish. In cases when this process alone does not succeed, for instance because taylorform is unable to compute a sufficiently good approximation to f, supnorm falls back to bisecting the working interval until safe supremum norm bounds can be computed with the required accuracy or until the width of the subintervals becomes less than diam times the original interval I, in which case supnorm fails.
- The algorithm used for supnorm is quite complex, which makes it impossible
to explain it here in further detail.
Please find a complete description in the following article:

Sylvain Chevillard, John Harrison, Mioara Joldes, Christoph Lauter

Efficient and accurate computation of upper bounds of approximation errors

Journal of Theoretical Computer Science (TCS), 2010

LIP Research Report number RR LIP2010-2

http://prunel.ccsd.cnrs.fr/ensl-00445343/fr/

- In practical cases, supnorm should be able to automatically handle removable discontinuities that relative errors might have. This means that usually, if f vanishes at a point x0 in the interval considered, the approximation polynomial p is designed such that it also vanishes at the same point with a multiplicity large enough. Hence, although f vanishes, epsilon_relative=p/f-1 may be defined by continuous extension at such points x0, called removable discontinuities (see Example 3).

> midpointmode=on!;

> supnorm(p, exp(x), [-1;1], absolute, 2^(-40));

0.452055210438~2/7~e-4

> midpointmode=on!;

> d = [1;2];

> f = exp(cos(x)^2 + 1);

> p = remez(1,15,d,1/f,1e-40);

> theta=1b-60;

> prec=default!;

> mode=relative;

> supnorm(p,f,d,mode,theta);

0.30893006200251428~5/6~e-13

> mode=relative;

> theta=1b-135;

> d = [-1b-2;1b-2];

> f = expm1(x);

> p = x * (1 + x * ( 2097145 * 2^(-22) + x * ( 349527 * 2^(-21) + x * (87609 * 2^(-21) + x * 4369 * 2^(-19)))));

> theta=1b-40;

> supnorm(p,f,d,mode,theta);

0.98349131972~2/3~e-7

See also: dirtyinfnorm, infnorm, checkinfnorm, absolute, relative, taylorform, autodiff, numberroots, diam

Go back to the list of commands- I is an interval.
- x is a real number.

- Returns the upper bound of the interval I. Each bound of an interval has its own precision, so this command is exact, even if the current precision is too small to represent the bound.
- When called on a real number x, sup considers it as an interval formed of a single point: [x, x]. In other words, sup behaves like the identity.

3

> sup(5);

5

> I=[0; 0.111110000011111_2];

> sup(I);

1.11110000011111_2 * 2^(-1)

> prec=12!;

> sup(I);

1.11110000011111_2 * 2^(-1)

- L is a list.

[|2, 3|]

> tail([|1,2...|]);

[|2...|]

See also: head

Go back to the list of commands- tanh is the hyperbolic tangent function, defined by tanh(x) = sinh(x)/cosh(x).
- It is defined for every real number x.

- tan is the tangent function, defined by tan(x) = sin(x)/cos(x).
- It is defined for every real number x that is not of the form n*Pi+Pi/2 where n is an integer.

- f is the function to be approximated.
- n is the degree of the polynomial that must approximate f.
- x0 is the point (it can be a real number or an interval) where the Taylor exansion of the function is to be considered.
- I is the interval over which the function is to be approximated. If this parameter is omitted, the behavior is changed (see detailed description below).
- errorType (optional) is the type of error to be considered. See the detailed description below. Default is absolute.

- taylorform computes an approximation polynomial and an interval error bound
for function f. More precisely, it returns a list L=[p, coeffErrors, Delta]
where:
- p is an approximation polynomial of degree n which is roughly speaking a numerical Taylor expansion of f at the point x0.
- coeffsErrors is a list of n+1 intervals. Each interval coeffsErrors[i] contains an enclosure of all the errors accumulated when computing the i-th coefficient of p.
- Delta is an interval that provides a bound for the approximation error between p and f. Its significance depends on the errorType considered.

- The polynomial p and the bound Delta are obtained using Taylor Models principles.
- Please note that x0 can be an interval. In general, it is meant to be a small interval approximating a non representable value. If x0 is given as a constant expression, it is first numerically evaluated (leading to a faithful rounding xapprox0 at precision prec), and it is then replaced by the (exactly representable) point-interval [xapprox0, xapprox0]. In particular, it is not the same to call taylorform with x0 = pi and with x0 = [pi], for instance. In general, if the point around which one desires to compute the polynomial is not exactly representable, one should preferably use a small interval for x0.
- More formally, the mathematical property ensured by the algorithm may be
stated as follows. For all xi0 in x0, there exist (small) values
eps[i] in coeffsErrors[i] such that:

If errorType is absolute, for all x in I, there exists delta in Delta such that f(x) - p(x-xi0) = sum{i=0...n} eps[i]*(x-xi0)^i + delta.

If errorType is relative, for all x in I, there exists delta in Delta such that f(x) - p(x-xi0) = sum{i=0...n} eps[i]*(x-xi0)^i + delta*(x-xi0)^(n+1). - It is also possible to use a large interval for x0, though it is not obvious to give an intuitive sense to the result of taylorform in that case. A particular case that might be interesting is when x0=I in relative mode. In that case, denoting by p_i the coefficient of p of order i, the interval p_i + coeffsError[i] gives an enclosure of f^(i)(I)/i!. However, the command autodiff is more convenient for computing such enclosures.
- When the interval I is not given, the approximated Taylor polynomial is computed but no remainder is produced. In that case the returned list is L=[p, coeffErrors].
- The relative case is especially useful when functions with removable singularities are considered. In such a case, this routine is able to compute a finite remainder bound, provided that the expansion point given is the problematic removable singularity point.
- The algorithm does not guarantee that by increasing the degree of the approximation, the remainder bound will become smaller. Moreover, it may even become larger due to the dependecy phenomenon present with interval arithmetic. In order to reduce this phenomenon, a possible solution is to split the definition domain I into several smaller intervals.
- The command taylor also computes a Taylor polynomial of a function. However it does not provide a bound on the remainder. Besides, taylor is a somehow symbolic command: each coefficient of the Taylor polynomial is computed exactly and returned as an expression tree exactly equal to theoretical value. It is henceforth much more inefficient than taylorform and taylorform should be prefered if only numercial (yet safe) computations are required. The same difference exists between commands diff and autodiff.

> p=TL[0];

> Delta=TL[2];

> errors=TL[1];

> for epsi in errors do epsi;

[0;0]

[0;0]

[0;5.3455294201843912922810729343029637576303937602101e-51]

[0;0]

[-3.3409558876152445576756705839393523485189961001313e-52;3.3409558876152445576756705839393523485189961001313e-52]

[0;0]

[-1.04404871487976392427364705748104760891218628129103e-53;1.04404871487976392427364705748104760891218628129103e-53]

[0;0]

[-1.63132611699963113167757352731413688892529106451724e-55;1.63132611699963113167757352731413688892529106451724e-55]

[0;0]

[-1.91171029335894273243465647732125416670932546623114e-57;1.91171029335894273243465647732125416670932546623114e-57]

> p; Delta;

1 + x^2 * (-0.16666666666666666666666666666666666666666666666667 + x^2 * (8.3333333333333333333333333333333333333333333333333e-3 + x^2 * (-1.98412698412698412698412698412698412698412698412698e-4 + x^2 * (2.75573192239858906525573192239858906525573192239859e-6 + x^2 * (-2.50521083854417187750521083854417187750521083854419e-8)))))

[-1.6135797443886066084999806203254010793747502812764e-10;1.6135797443886066084999806203254010793747502812764e-10]

> p=TL[0];

> Delta=TL[2];

> p; Delta;

1 + x * (1 + x * (0.5 + x * (0.16666666666666666666666666666666666666666666666667 + x * (4.1666666666666666666666666666666666666666666666667e-2 + x * (8.3333333333333333333333333333333333333333333333333e-3 + x * (1.38888888888888888888888888888888888888888888888889e-3 + x * (1.98412698412698412698412698412698412698412698412698e-4 + x * (2.4801587301587301587301587301587301587301587301587e-5 + x * (2.75573192239858906525573192239858906525573192239859e-6 + x * 2.7557319223985890652557319223985890652557319223986e-7)))))))))

[-2.31142719641187619441242534182684745832539555102969e-8;2.7312660755642474420206278018039434042553645532164e-8]

> TL2 = taylorform(exp(x), 10, [log2(10)], [-1,1], absolute);

> TL1==TL2;

false

> TL2 = taylorform(exp(x), 3, 0, relative);

> TL1[0]==TL2[0];

true

> TL1[1]==TL2[1];

true

> length(TL1);

3

> length(TL2);

2

> TL = taylorform(f, 3, x0);

> T1 = TL[0];

> T2 = taylor(f, 3, x0);

> print(coeff(T1, 2));

-1.35914091422952261768014373567633124887862354684999

> print(coeff(T2, 2));

-(0.5 * exp(1))

- n represents the number of recursions

- taylorrecursions is a global variable. Its value represents the number of steps of recursion that are used when applying Taylor's rule. This rule is applied by the interval evaluator present in the core of Sollya (and particularly visible in commands like infnorm).
- To improve the quality of an interval evaluation of a function f, in particular when there are problems of decorrelation), the evaluator of Sollya uses Taylor's rule: f([a,b]) C f(m) + [a-m, b-m]*f'([a,b]) where m=(a+b)/2. This rule can be applied recursively. The number of step in this recursion process is controlled by taylorrecursions.
- Setting taylorrecursions to 0 makes Sollya use this rule only once; setting it to 1 makes Sollya use the rule twice, and so on. In particular: the rule is always applied at least once.

> p=remez(f,3,[0;1]);

> taylorrecursions=0;

The number of recursions for Taylor evaluation has been set to 0.

> evaluate(f-p, [0;1]);

[-0.46839364816303627522963565754743169862357620487739;0.46947781754667086491682464997088054443583003517779]

> taylorrecursions=1;

The number of recursions for Taylor evaluation has been set to 1.

> evaluate(f-p, [0;1]);

[-0.13813111495387910066337940912697015317218647208804;0.13921528433751369035056840155041899898444030238844]

- function represents the function to be expanded
- degree represents the degree of the expansion to be delivered
- point represents the point in which the function is to be developped

- The command taylor returns an expression that is a Taylor expansion
of function function in point point having the degree degree.

Let f be the function function, t be the point point and n be the degree degree. Then, taylor(function,degree,point) evaluates to an expression mathematically equal to f(t) + f'(t) * x + ... + 1/(n!) * f[n](t) * x^n. In other words, if p(x) denotes the polynomial returned by taylor, p(x-t) is the Taylor polynomial of degree n of f developped at point t.

Remark that taylor evaluates to 0 if the degree degree is negative.

exp(1) + x * (exp(1) + x * (0.5 * exp(1) + x * exp(1) / 6))

x * (1 + x^2 * (1 / 6 + x^2 * (9 / 120 + x^2 * 225 / 5040)))

x * (1 / sqrt((pi) / 4) + x^2 * ((sqrt((pi) / 4) * 4 / (pi) * (-2)) / 6 + x^2 * (sqrt((pi) / 4) * 4 / (pi) * 12) / 120))

- code is the code to be timed.

- time permits timing a Sollya instruction, resp. a begin-end block of Sollya instructions. The timing value, measured in seconds, is returned as a Sollya constant (and not merely displayed as for timing). This permits performing computations of the timing measurement value inside Sollya.
- The extended nop command permits executing a defined number of useless instructions. Taking the ratio of the time needed to execute a certain Sollya instruction and the time for executing a nop therefore gives a way to abstract from the speed of a particular machine when evaluating an algorithm's performance.

> write(t,"s were spent computing p = ",p,"\n");

0.143369000000000000019824419883462951474939472973347s were spent computing p = -3.3426550293345171908513995127407122194691200059639e-17 + x * (0.99999999973628359955372011464713121003442988167693 + x * (7.8802751877302786684499343799047732495568873819693e-16 + x * (-0.166666661386013037032912982196741385680498698107285 + x * (-5.3734444911159112186289355138557504839692987221233e-15 + x * (8.3333037186548537651002133031675072810009327877148e-3 + x * (1.33797221389218815884112341005509831429347230871284e-14 + x * (-1.98344863018277416493268155154158924422004290239026e-4 + x * (-1.3789116451286674170531616441916183417598709732816e-14 + x * (2.6876259495430304684251822024896210963401672262005e-6 + x * 5.0282378350010211058128384123578805586173782863605e-15)))))))))

The error is 2^(log2(2.39602467695631727848641768186659313738474584992648e-11))

0.27866300000000000000165839564303382758225779980421 s were spent

> write(~(t-10),"s of execution overhead.\n");

1.79400000000000045541348470123921288177371025085449e-3s of execution overhead.

> write("This ratio = ", ratio, " should somehow be independent of the type of machine.\n");

This ratio = 3.56569315332201930141279648066430730792081206162027 should somehow be independent of the type of machine.

- activation value controls if timing should be performed or not

> timing=on;

Timing has been activated.

> p=remez(sin(x),10,[-1;1]);

Information: Remez: computing the matrix spent 1 ms

Information: Remez: computing the quality of approximation spent 7 ms

Information: Remez: computing the matrix spent 1 ms

Information: Remez: computing the quality of approximation spent 7 ms

Information: Remez: computing the matrix spent 1 ms

Information: Remez: computing the quality of approximation spent 7 ms

Information: computing a minimax approximation spent 142 ms

Information: assignment spent 142 ms

Information: full execution of the last parse chunk spent 142 ms

- tripledouble is both a function and a constant.
- As a function, it rounds its argument to the nearest number that can be written as the sum of three double precision numbers.
- The algorithm used to compute tripledouble(x) is the following: let xh = double(x), let xm = double(x - xh) and let xl = double(x - xh - xm). Return the number xh + xm + xl. Note that if the current precision is not sufficient to represent exactly xh + xm + xl, a rounding will occur and the result of tripledouble(x) will be useless.
- As a constant, it symbolizes the triple-double precision format. It is used in contexts when a precision format is necessary, e.g. in the commands roundcoefficients and implementpoly. See the corresponding help pages for examples.

> a = 1+ 2^(-55)+2^(-115);

> TD(a);

1.00000000000000002775557561562891353466491600711096

> prec=110!;

> TD(a);

Warning: double rounding occurred on invoking the triple-double rounding operator.

Try to increase the working precision.

1.000000000000000027755575615628913

See also: halfprecision, single, double, doubleextended, doubledouble, quad, roundcoefficients, implementpoly, fpminimax, printexpansion

Go back to the list of commands- true is the usual boolean value.

false

> 2>1;

true

- identifier1, identifier2,... , identifiern represent variable identifiers

- The keyword var allows for the declaration of local variables
identifier1 through identifiern in a begin-end-block ({}-block).
Once declared as a local variable, an identifier will shadow
identifiers declared in higher scopes and undeclared identifiers
available at top-level.

Variable declarations using var are only possible in the beginning of a begin-end-block. Several var statements can be given. Once another statement is given in a begin-end-block, no more var statements can be given.

Variables declared by var statements are dereferenced as error until they are assigned a value.

exp(x)

> a = 3;

> {var a, b; a=5; b=3; {var a; var b; b = true; a = 1; a; b;}; a; b; };

1

true

5

3

> a;

3

- n controls the amount of information to be displayed

- verbosity accepts any integer value. At level 0, commands do not display anything on standard output. Note that very critical information may however be displayed on standard error.
- Default level is 1. It displays important information such as warnings when roundings happen.
- For higher levels more information is displayed depending on the command.

> 1.2+"toto";

error

> verbosity=1!;

> 1.2+"toto";

Warning: Rounding occurred when converting the constant "1.2" to floating-point with 165 bits.

If safe computation is needed, try to increase the precision.

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

error

> verbosity=2!;

> 1.2+"toto";

Warning: Rounding occurred when converting the constant "1.2" to floating-point with 165 bits.

If safe computation is needed, try to increase the precision.

Warning: at least one of the given expressions or a subexpression is not correctly typed

or its evaluation has failed because of some error on a side-effect.

Information: the expression or a partial evaluation of it has been the following:

(1.19999999999999999999999999999999999999999999999999) + ("toto")

error

- The variable void represents the functional result of a
side-effect or an empty argument. It is used only in combination with
the applications of procedures or identifiers bound through
externalproc to external procedures.

The void result produced by a procedure or an external procedure is not printed at the prompt. However, it is possible to print it out in a print statement or in complex data types such as lists.

The void argument is implicit when giving no argument to a procedure or an external procedure when applied. It can nevertheless be given explicitly. For example, suppose that foo is a procedure or an external procedure with a void argument. Then foo() and foo(void) are correct calls to foo. Here, a distinction must be made for procedures having an arbitrary number of arguments. In this case, an implicit void as the only parameter to a call of such a procedure gets converted into an empty list of arguments, an explicit void gets passed as-is in the formal list of parameters the procedure receives. - void is used also as a type identifier for externalproc. Typically, an external procedure taking void as an argument or returning void is bound with a signature void -> some type or some type -> void. See externalproc for more details.

void

> void;

> hey;

proc()

{

print("Hello world.");

return void;

}

> hey();

Hello world.

> hey(void);

Hello world.

> print(hey());

Hello world.

void

> bashexecute("gcc -fPIC -shared -o externalprocvoidexample externalprocvoidexample.o");

> externalproc(foo, "./externalprocvoidexample", void -> void);

> foo;

foo(void) -> void

> foo();

Hello from the external world.

> foo(void);

Hello from the external world.

> print(foo());

Hello from the external world.

void

> blub(1);

Argument list: [|1|]

> blub();

Argument list: [| |]

> blub(void);

Argument list: [|void|]

- function represents the function to be considered
- preimage precision represents the precision of the preimages
- preimage exponent range represents the exponents in the preimage format
- image precision represents the precision of the format the images are to be rounded to
- error bound represents the upper bound for the search w.r.t. the relative rounding error
- filename represents a character sequence containing a filename

prec = 165

x = 1.99999988079071044921875 f(x) = 7.3890552520751953125 eps = 4.5998601423446695596184695493764120138001954979037e-9 = 2^(-27.695763)

x = 2 f(x) = 7.38905620574951171875 eps = 1.44563608749673018122228379395533417878125150587072e-8 = 2^(-26.043720)

- expr represents an expression
- filename represents a character sequence indicating a file name

- write(expr1,...,exprn) prints the expressions expr1 through
exprn. The character sequences corresponding to the expressions are
concatenated without any separator. No newline is displayed at the
end. In contrast to print, write expects the user to give all
separators and newlines explicitly.

If a second argument filename is given after a single ">", the displaying is not output on the standard output of Sollya but if in the file filename that get newly created or overwritten. If a double ">>" is given, the output will be appended to the file filename.

The global variables display, midpointmode and fullparentheses have some influence on the formatting of the output (see display, midpointmode and fullparentheses).

Remark that if one of the expressions expri given in argument is of type string, the character sequence expri evaluates to is displayed. However, if expri is of type list and this list contains a variable of type string, the expression for the list is displayed, i.e. all character sequences get displayed surrounded by quotes ("). Nevertheless, escape sequences used upon defining character sequences are interpreted immediately.

> write("Hello\n");

x + 2 + exp(sin(x))Hello

> write("Hello","world\n");

Helloworld

> write("Hello","you", 4 + 3, "other persons.\n");

Helloyou7other persons.

Hello

> write([|"Hello"|],"\n");

[|"Hello"|]

> s = "Hello";

> write(s,[|s|],"\n");

Hello[|"Hello"|]

> t = "Hello\tyou";

> write(t,[|t|],"\n");

Hello you[|"Hello\tyou"|]

> readfile("foo.sol");

x + 2 + exp(sin(x))

See also: print, printexpansion, printdouble, printsingle, printxml, readfile, autosimplify, display, midpointmode, fullparentheses, evaluate, roundingwarnings, autosimplify

Go back to the list of commands