Unums or BigDecimals+Prolog?

[burs...@gmail.com]

Hi,

Discovered a couple of gems. One gem are the BigDecimals.
It seems that they not only offer unlimited precision
calculations, as can be seen here:

?- X is 0d1.23 * 0d4.56.
X = 0d5.6088

But all the arithmetic Java methods such as multiple()
etc.. have a variant where we can pass a math context
as an additional argument which defines some rounding.

These additional functions can be made accessible via
an evaluable functions vpa/2. This evaluable function is
special, since it does't directly evaluate its first

argument. Instead it crawls through the first argument,
and replaces all unlimited precision calculations by
some limited precision calculations:

?- X is vpa(0d1.23 * 0d4.56, 3).
X = 0d5.61

What can be done except +, -, *, / with these variant
methods? Well we embarked on implementing sin, cos, atan,
exp, log and are almost done, here are some examples:

?- X is vpa(exp(100000), 32).
X = 0d2.8066633604261231793183858327690E+43429

?- X is vpa(log(exp(100000)), 30).
X = 0d100000.000000000000000000000001

?- X is vpa(log(exp(100000)), 31).
X = 0d99999.99999999999999999999999995

?- X is vpa(log(exp(100000)), 32).
X = 0d100000.00000000000000000000000000

The BigDecimals have an amazing breadth when the mantissas
is resricted by a given precision. The exponent can be easily
bigger than what the new Windows 10 calculator allows.

Have a Nice day!

Short description: Variable Precission Functions in Prolog
https://plus.google.com/u/0/b/103259555581227445618/+JekejekeCh/posts/dLpzvJ2BrRj

Sourc code: Gist Source Code of sin, cos, atan, etc..
https://gist.github.com/jburse/a209b1fab191964adbd3eb15539b188c#file-trans-p

[j4n bur53]

Judging from some results published elsewhere(*),
we are kind of 3x times slower than MATLAB:

% VPA (MATLAB R2016b x64):
>> digits(25);
>> A = vpa(200*(rand(2000)-0.5));

>> tic; log(A); toc;
Elapsed time is 1161.7028 seconds.

I get the following figures, where we recompute constants
log(2) and log(10) internally each time (sic!):

Jekejeke Prolog 2, Runtime Library 1.1.8
(c) 1985-2017, XLOG Technologies GmbH, Switzerland

?- time((between(1, 2000, _), X is 100*random,
Y is vpa(log(X), 34), fail; true)).
% Uptime 1,889 ms, GC Time 20 ms, Thread Cpu Time 1,859 ms
Yes

?- X is 1.889*2000.
X = 3778.0

I guess with some constant caching or precomputation
we could improve the runtime drastically.

(*)
http://www.advanpix.com/documentation/users-manual/#Symbolic_Math_Toolbox_VPA_vs_Multiprecision_Computing_Toolbox

[burs...@gmail.com]

With cached constants even competitive with
MATLAB. I now get the following:

?- time((between(1, 2000, _), X is 100*random,
Y is vpa(log(X), 34), fail; true)).

% Uptime 554 ms, GC Time 8 ms, Thread Cpu Time 562 ms
Yes

?- X is 0.554*2000.
X = 1108.0

?- listing(_:vpa_cache/_).
vpa_cache(log5div4, 34, 0d0.2231435513142097557662950903098345).
vpa_cache(log32div25, 34, 0d0.2468600779315257978846419408385076).

But there are of course more factors, the
tests were not run on the same machine etc..

[burs...@gmail.com]

Here some results for same machine runnings,
I get a factor 3x slower of the Prolog implementation
compared with the MATLAB VPA implementation:

>> version
ans = 9.1.0.441655 (R2016b)

>> digits(25);
>> A = vpa(100*(rand(1,2000)));
>> tic; log(A); toc;
Elapsed time is 0.148637 seconds.
>> 0.148637 * 2000
ans = 297.274

And the Advanpix library has definitively the faster
logarithm functions, I get the following results so far:

>> mp.Digits(34);
>> A = mp(100*rand(1,2000));
>> tic; log(A); toc;
Elapsed time is 0.000820 seconds.
>> 0.000820 * 2000
ans = 1.64

A check shows that computation is indeed 34 digits:

>> A(1,1)
ans = 95.09760483334387970444367965683341
>> B = log(A);
>> B(1,1)
ans = 4.554903783464923127304549242118646

Note our own computation is not yet optimized precision
wise, i.e no extra digits yet internally, so we have:

Jekejeke Prolog 2, Runtime Library 1.1.8
(c) 1985-2017, XLOG Technologies GmbH, Switzerland

?- X is vpa(log(95.09760483334387970444367965683341), 34).
X = 0d4.554903783464923127304549242118647

[burs...@gmail.com]

In the C programming language, one would probably go with
the types _Decimal32, _Decimal64 and _Decimal128:

6.13 Decimal Floating Types
https://gcc.gnu.org/onlinedocs/gcc/Decimal-Float.html

But the Java class is more generic, its just a math context
that says, what decimal you want, but its a single class BigDecimal

There is also a ISO draft which includes sin, cos, etc.. and
exotic things like sinpi, cospi, etc..:

C11 extension ISO/IEC TS 18661-4:2015
http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1950.pdf

Am Freitag, 17. Februar 2017 14:50:39 UTC+1 schrieb burs...@gmail.com:

[burs...@gmail.com]

Practical usage of n = 2, i.e. a circle is two half-turns,
is seen in the system-R package. The advantage is that the

quadrant decision, which is usually needed in computing them,
can be made more precisie, the description is:

cospi(x), sinpi(x), and tanpi(x), compute cos(pi*x),
sin(pi*x), and tan(pi*x).

cospi(x), sinpi(x), and tanpi(x) are accurate
for x values which are multiples of a half.

https://stat.ethz.ch/R-manual/R-devel/library/base/html/Trig.html

[burs...@gmail.com]

It seems that the argumemt reduction for sin(),
cos() etc.. isn't trivial. Here is a test case:

Welcome to SWI-Prolog (threaded, 64 bits, version 7.5.1)
SWI-Prolog comes with ABSOLUTELY NO WARRANTY. This is free software.

?- X is sin(10^22).
X = 0.46261304076460175.

?- X is cos(10^22).
X = -0.8865603050636369.

And here is another result relying on Java Math:

Jekejeke Prolog 2, Runtime Library 1.1.7
(c) 1985-2016, XLOG Technologies GmbH, Switzerland

?- X is sin(10^22).
X = -0.8522008497671888

?- X is cos(10^22).
X = 0.523214785395139

The general wishdom is that the correct result is

X = -0.852200849...

X = 0.523214785...

See also Table 1: sin(x) and cos(x) for x = 1022 radians here:

ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit
K. C. Ng and the members of the FP group of SunPro
Works in progress, March 24, 1992
https://www.csee.umbc.edu/~phatak/645/supl/Ng-ArgReduction.pdf