Allegro.cc - Online Community

Allegro.cc Forums » Programming Questions » Looking for math library

This thread is locked; no one can reply to it. rss feed Print
 1   2   3   4 
Looking for math library
Edgar Reynaldo
Major Reynaldo
May 2007
avatar

I'm looking for a decent math library. It has to satisfy the following conditions :

1) When scanned, "#.#" has to equal #.# - ie... "5.1" has to equal 5.1 when scanned from a string. Because currently if I use sscanf on "5.1" to get a double and then printf the double back out with high precision, it is 5.0999999999999996 instead of 5.1, and this occurs far too often with other numbers as well, even if it is only off by 4e-16 or so. I want to be able to move from string to double and back without losing any precision.

2) It should support at least as many math functions as the C standard library.

3) C or C++ only.

4) Non GPL if possible, but maybe you can convince me otherwise if
a) My code doesn't have to be GPL too.
b) I'm not forced to release all of my source code, even if I might anyway.

5) Optional - extra features?

I looked at The wikipedia list article of numerical libraries, but about the only thing that looked good at first glance was the GNU MultiPrecision Library. The GNU Scientific Library doesn't seem to offer precision greater than doubles do, only a truckload of math functions.

The problem with GMP is that it uses the LGPL, and I still haven't wrapped my head around what you can and can't do with the (L)GPL. I wish there was a (L)GPL for dummies guide around somewhere. :P

So, any recommendations / experience to share with me on this matter?

SiegeLord
Member #7,827
October 2006
avatar

I want to be able to move from string to double and back without losing any precision.

This is impossible. You need to have decimal fraction storage to do that, and double does not provide that.

I only know of GMP and GSL too, but generally I've never needed the additional precision of the first.

"For in much wisdom is much grief: and he that increases knowledge increases sorrow."-Ecclesiastes 1:18
[SiegeLord's Abode][Codes]:[DAllegro5]:[RustAllegro]

Matthew Leverton
Supreme Loser
January 1999
avatar

You won't have a problem with LGPL libraries corrupting your source code as long as you dynamically link against them.

Arthur Kalliokoski
Second in Command
February 2005
avatar

Because currently if I use sscanf on "5.1" to get a double and then printf the double back out with high precision, it is 5.0999999999999996 instead of 5.1,

#include <stdio.h>

double f;
int i;

int main(void)
{
  for(i=0;i<100;i++)
  {
    printf("%0.1lf %lf\n",f,f);  //print the single decimal place version followed by the ordinary version.
    f += 0.11111;
  }
  return 0;
}

They all watch too much MSNBC... they get ideas.

decepto
Member #7,102
April 2006
avatar

I love this library, even when I'm not programming OpenGL. It's just great all around.

--------------------------------------------------
Boom!

Edgar Reynaldo
Major Reynaldo
May 2007
avatar

SiegeLord said:

This is impossible. You need to have decimal fraction storage to do that, and double does not provide that.

I know double doesn't do that, and that was kind of the point.

You won't have a problem with LGPL libraries corrupting your source code as long as you dynamically link against them.

Dynamic linking is fine with me. Am I supposed to mention the library and provide a copy of the GPL and LGPL along with it? Can I just provide a link to the GPL, LGPL, and the url of the library used?

@Arthur Kalliokoski
I don't really know what you are trying to say with that code. Using a precision specifier with printf only allows you to specify the number of digits after the decimal point, not the actual number of significant digits. Ie 12.3 has 3 significant digits, not 1. It also does rounding, which I want to minimize as much as possible.

decepto said:

I love this library [glm.g-truc.net], even when I'm not programming OpenGL. It's just great all around.

It looks decent to start off, but there's no guarantee of precision with their two high precision data types :

typedef highp_float_t highp_float

High precision floating-point numbers.

There is no guarantee on the actual precision. From GLSL 1.30.8 specification

Definition at line 54 of file type_float.hpp.
typedef detail::highp_int_t highp_int

High precision signed integer.

There is no guarantee on the actual precision. From GLSL 1.30.8 specification.

Definition at line 62 of file type_int.hpp.

What I want most is to be able to have a number string convert to an exactly corresponding decimal data type and back again. I can limit the precision later if necessary.

The decNumber library is looking good too, except for the GPL bit. If I use a GPL library in my program, is my program then automatically bound by the GPL as well?

Arthur Kalliokoski
Second in Command
February 2005
avatar

Using a precision specifier with printf only allows you to specify the number of digits after the decimal point, not the actual number of significant digits. Ie 12.3 has 3 significant digits, not 1. It also does rounding, which I want to minimize as much as possible.

You want to truncate digits? Why? If you want to eliminate rounding, multiply by the power of 10 you want to have the precision up to, floor() it, then divide again.
Or you could simply sprintf this number to a buffer, use strchr() to find the decimal point, use the index to the sought number for strcpy, etc. Sure, it'd be slower than pure math operations, but a library would slow things down as well.

They all watch too much MSNBC... they get ideas.

tobing
Member #5,213
November 2004
avatar

If I use a GPL library in my program, is my program then automatically bound by the GPL as well?

Yes, you have to publish your source code as well under GPL. GPL is infectious.

When I found this: http://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic I was surprised how many such libraries there are... I have no experiences with any of them, however.

Arthur Kalliokoski
Second in Command
February 2005
avatar

I have a 32 bit assembler arbitrary precision library you can have if you want, but I haven't been motivated enough to finish the 64 bit version. To convert to assembly code for Windows, you'd have to add an underscore prefix to the global namespaces. I never did quite finish the transcendental functions either, but it does the fourbanger calculator stuff just fine, as well as print to ascii buffers and scan from ascii buffers.

They all watch too much MSNBC... they get ideas.

james_lohr
Member #1,947
February 2002

You need to tell us what you want this for, because most likely you don't actually know what you really want; it is very unlikely that you need a fully blown Decimal library for whatever it is you are doing. Chances are that all you really need is to write your own double.ToString() method.

Edgar Reynaldo
Major Reynaldo
May 2007
avatar

You want to truncate digits? Why? If you want to eliminate rounding, multiply by the power of 10 you want to have the precision up to, floor() it, then divide again.
Or you could simply sprintf this number to a buffer, use strchr() to find the decimal point, use the index to the sought number for strcpy, etc. Sure, it'd be slower than pure math operations, but a library would slow things down as well.

No, I said I want to minimize rounding - ie maintain precision up to an arbitrary point, maybe around 30 significant digits or so. Multiplying and dividing by 10 won't help fix the failure in precision of statements like 0.1 or 1/10.

Also, super speed is not essential in this application. It won't be computationally expensive over long periods of time.

Tobing said:

Yes, you have to publish your source code as well under GPL. GPL is infectious.

<runs for antibiotics>

I have a 32 bit assembler arbitrary precision library you can have if you want, but I haven't been motivated enough to finish the 64 bit version. To convert to assembly code for Windows, you'd have to add an underscore prefix to the global namespaces. I never did quite finish the transcendental functions either, but it does the fourbanger calculator stuff just fine, as well as print to ascii buffers and scan from ascii buffers.

I appreciate the offer, but I'm going to need trig functions.

You need to tell us what you want this for, because most likely you don't actually know what you really want; it is very unlikely that you need a fully blown Decimal library for whatever it is you are doing. Chances are that all you really need is to write your own double.ToString() method.

I've already written a FormatDouble function, and it works fine. It retains the maximum amount of precision available from a double and strips leading and trailing zeros and optimizes the exponent for string length consideration.

What I'm using it for is a console/plotter program that lets you enter equations as strings and then evaluate / simplify / derive / integrate? / plot them.

The LUA interpreter has better precision than a C++ double - it doesn't have any problem with 5.1 or 0.1, or most other common (relatively small) numbers and has precision of somewhere around 14 digits.

BTW, keep comments like 'I don't actually know what I really want' to yourself. I've already stated several times that I want more precision than doubles afford, and I want perfect string to decimal data type and back to string conversions.

Audric
Member #907
January 2001

The restrictions of LGPL are super-light. If you don't need modify the "official" DLL, you only need to provide a readme like this one:

doc/README-SDL.txt said:

The following file:

SDL.dll

is the runtime environment for the SDL library.

The Simple DirectMedia Layer (SDL for short) is a cross-platfrom library
designed to make it easy to write multi-media software, such as games and
emulators.

The Simple DirectMedia Layer library source code is available from:
http://www.libsdl.org/

This library is distributed under the terms of the GNU LGPL license:
http://www.gnu.org/copyleft/lesser.html

A little more complex example where I had to compile the DLL myself: since I didn't have to modifiy the sources, I simply point to the specific versions used :

doc/README-SDL_image.txt said:

The following file

SDL_image.dll

is the runtime environment for the library SDL_image 1.2.10 by Sam Lantiga and Mattias Engdegård.

This library is distributed under the terms of the GNU LGPL license:
http://www.gnu.org/copyleft/lesser.html

The source is available from the libraries page at the SDL website:
http://www.libsdl.org/

=========================

Compiled by yrizoud on 16/06/2010 with static builds of :
libjpeg (JPEG library http://www.ijg.org/ release 6b of 27-Mar-1998 : jpegsr6.zip)
libtiff (TIFF library ftp://ftp.sgi.com/graphics/tiff/ : tiff-v3.4-tar.gz)
And dynamic version of :
libpng (http://www.libpng.org/pub/png/libpng.html) v1.4.2 --> libpng14-14.dll requirement

Edgar Reynaldo
Major Reynaldo
May 2007
avatar

james_lohr
Member #1,947
February 2002

It makes it all too apparent what a dated language C++ is; Java has had BigDecimal for over a decade, and "decimal" is a basic type in C#.

Kitty Cat
Member #2,815
October 2002
avatar

It makes it all too apparent what a dated language C++ is; Java has had BigDecimal for over a decade, and "decimal" is a basic type in C#.

Except it isn't a basic type. It may be built-in to the language, but the CPU/FPU doesn't handle it directly. Software has to break it down, do the work, then reconstruct it.

C/C++'s basic types are what the processor deals with. You don't even get float and double if your target machine doesn't have an FPU (unless your compiler offers FPU emulation, but that's not required as part of the language).

At most, it could be part of libc/libstdc++, but there are already libraries which implement it just fine.

--
"Do not meddle in the affairs of cats, for they are subtle and will pee on your computer." -- Bruce Graham

Arthur Kalliokoski
Second in Command
February 2005
avatar

IIRC, the lcc-win32 compiler has arbitrary precision arithmetic, but lcc-win32 is slow and you can't use it for commercial programs.

They all watch too much MSNBC... they get ideas.

orz
Member #565
August 2000

1) When scanned, "#.#" has to equal #.# - ie... "5.1" has to equal 5.1 when scanned from a string. Because currently if I use sscanf on "5.1" to get a double and then printf the double back out with high precision, it is 5.0999999999999996 instead of 5.1, and this occurs far too often with other numbers as well, even if it is only off by 4e-16 or so. I want to be able to move from string to double and back without losing any precision.

I don't think you should ask for something like that without specifying what your precision requirements are.

If you only care about destringifaction+stringification not producing the original value, simply limiting every decimal stringification to 14 digits meets your requirements well over a fairly wide range of numbers, though that's actually a decrease in precision.

On the other hand if you need purely decimal math, that's a very different thing (though it would solve the above problem as long as all stringifications were to decimal). Likewise if you need infinite precision.

Evert
Member #794
November 2000
avatar

Java has had BigDecimal for over a decade, and "decimal" is a basic type in C#.

Since computers use binary rather than decimal, I would consider that a rather useless feature.
I mean, sure, you can't express decimal farctions exactly in a binary representation, but they're both equally arbitrary anyway. You can't express a fraction of a third exactly in either, for instance.

Kitty Cat
Member #2,815
October 2002
avatar

orz said:

If you only care about destringifaction+stringification not producing the original value, simply limiting every decimal stringification to 14 digits meets your requirements well over a fairly wide range of numbers, though that's actually a decrease in precision.

Or you can use %a (C99), which uses a hexadecimal notation for floating-point values. The default precision allows you to distinguish between two double values.

a, A   (C99; not in SUSv2) For a conversion, the double argument  is  converted  to
       hexadecimal notation (using the letters abcdef) in the style [-]0xh.hhhhp±d;
       for A conversion the prefix 0X, the letters ABCDEF, and the exponent separa‐
       tor P is used.  There is one hexadecimal digit before the decimal point, and
       the number of digits after it is equal to the precision.  The default preci‐
       sion suffices for an exact representation of the value if an exact represen‐
       tation in base 2 exists and otherwise is sufficiently large  to  distinguish
       values  of  type  double.  The digit before the decimal point is unspecified
       for nonnormalized numbers, and nonzero but otherwise unspecified for normal‐
       ized numbers.

--
"Do not meddle in the affairs of cats, for they are subtle and will pee on your computer." -- Bruce Graham

William Labbett
Member #4,486
March 2004
avatar

Evert said:

I mean, sure, you can't express decimal farctions exactly in a binary representation, but they're both equally arbitrary anyway. You can't express a fraction of a third exactly in either, for instan

...let alone irrational numbers.

What about the gnu big number library ? Not sure if that would help.

Audric
Member #907
January 2001

The decimal type in C# is not for mathematicians, but accountants.

Arthur Kalliokoski
Second in Command
February 2005
avatar

Audric said:

The decimal type in C# is not for mathematicians, but accountants.

Isn't the '#.#' specifier from COBOL? It's meant for accountants and $PHB's.

They all watch too much MSNBC... they get ideas.

Edgar Reynaldo
Major Reynaldo
May 2007
avatar

orz said:

I don't think you should ask for something like that without specifying what your precision requirements are.

My precision requirements are that string numbers have an exact representation when converted to data. With my current implementation using doubles, my command line calculator / function parser is too stupid to return 51 for 5.1/0.1. That's just not acceptable to me, nor should it be acceptable to anyone using a calculator program.

Kitty Cat said:

Or you can use %a (C99), which uses a hexadecimal notation for floating-point values. The default precision allows you to distinguish between two double values.

Kitty Cat said:

The default precision suffices for an exact representation of the value if an exact representation in base 2 exists and otherwise is sufficiently large to distinguish values of type double.

This is the basis of the problem - most base 2 floating point representations are just incapable of representing decimals correctly. If they used base 2 to represent the numerals and used a separate number to represent the decimal point, then there shouldn't ever be a problem with it, but you can't represent every number with powers of 1/2.

What about the gnu big number library ? Not sure if that would help.

That's the same thing as the GNU Multi Precision library. I'm looking into it.

verthex
Member #11,340
September 2009
avatar

Theres LIDIA but I'm not sure what it could do for you.

Arthur Kalliokoski
Second in Command
February 2005
avatar

With my current implementation using doubles, my command line calculator / function parser is too stupid to return 51 for 5.1/0.1.

#SelectExpand
1#include <stdio.h> 2#include <stdlib.h> 3#include <errno.h> 4 5int main(int argc, char **argv) 6{ 7 double parm1; 8 double parm2; 9 double result; 10 int intresult; 11 int operator; 12 char *end1; 13 14 if(argc != 2) 15 { 16 fprintf(stderr,"I want one string in quotes to parse for four banger math, e.g. \"5/4\"\n"); 17 return 0; 18 } 19 20 errno = 0; 21 parm1=strtod(argv[1],&end1); 22 if(errno == ERANGE) 23 { 24 fprintf(stderr,"First number out of range\n"); 25 return 1; 26 } 27 28 operator = (char)*end1; 29 30 end1++; 31 32 parm2=strtod(end1,0); 33 34 switch(operator) 35 { 36 case '+': 37 result = parm1 + parm2; 38 break; 39 case '-': 40 result = parm1 - parm2; 41 break; 42 case '*': 43 result = parm1 * parm2; 44 break; 45 case '/': 46 result = parm1 / parm2; 47 break; 48 default: 49 fprintf(stderr,"Operator '%c' isn't allowed\n",operator); 50 return 1; 51 } 52 53 intresult = result; 54 if(intresult == result) 55 printf("answer is %d\n",intresult); 56 else 57 printf("answer is %lg\n",result); 58 59 return 0; 60}

[EDIT]

Forgot the range check on second number, but you know what I mean...

[EDIT2]

The arbitrary precision lib I mentioned earlier does have transcendental functions, but some of them could be optimized with precalculated numbers in a file for identities, which I haven't bothered with.

I tried it just to see how fast it was on this 3.1Ghz Athlon II, but somehow argc & argv were getting corrupted. After hours of fiddling with gdb, valgrind and Electric Fence, I tried the Watcom compiler, which worked fine (although Watcom crashes on an fflush() :P).

Printing gives out results like:

count = 7021

Input number    0.7021
Sin    0.6458224341509764981664646025062778090792175100125560862802365350598448
Cos    0.7634876446592358739582226598267228849622487768704818372243868793930424
Tan    0.8458845911504221153467969183513439373684976836083291157721866966077248

count = 7022

Input number    0.7022
Sin    0.6458987796862030057491918266634106923511167321066618731945697836968761
Cos    0.7634230585984901932657298394753035699212712280249010525166211455076285
Tan    0.846056157737701829266851733611519926546263348925174306300998247937501

count = 7023

Input number    0.7023
Sin    0.6459751187624417218523788222528204136464410918821905557397273206694932
Cos    0.7633584649035139329501939060000511971689937336425396753570132724535036
Tan    0.8462277533584315588906546332567283389122943855797825858950950746135725

and here's how long it takes to print out 10000 loops with 80 bytes in the integral part and 80 bytes in the fractional part. The first one is redirecting to results.txt, where the output above was copied from.

pepsi@fractal_comet:/home/prog/bignums 02:49 AM $ time t > results.txt

real    0m59.132s
user    0m59.021s
sys     0m0.087s
pepsi@fractal_comet:/home/prog/bignums 02:50 AM $ time t noprint

real    0m0.038s
user    0m0.038s
sys     0m0.000s
pepsi@fractal_comet:/home/prog/bignums 02:50 AM $ 

The second timing result is how long it takes to calculate the trig funcs (1/30 of a second for all 30000 results) but not printing to decimal (quite a chore in itself)

[EDIT5]

I forgot to show the input number for the trig stuff... edited all the above

[EDIT6]

Recompiling with 800 bytes of integral and decimal bytes with 700 decimal positions printed gives this for 100 loops

pepsi@fractal_comet:/home/prog/bignums 03:02 AM $ time t > results.txt

real    4m9.194s
user    4m8.927s
sys     0m0.093s
pepsi@fractal_comet:/home/prog/bignums 03:07 AM $ time t noprint

real    0m0.004s
user    0m0.003s
sys     0m0.000s

noprint is faster than before because only 100 loops

count = 42

Input number    0.42
Sin    0.4077604530595701859727871580863423156477998312272199837618432020183270288765336837439844031983341253779690576147116094979039277815379680051349489526859390640113444468728511707118838605394853620887365370981415293797665495368793143468270256405581832984876247762177846669463896709533912328054270139386495849713716053377167824292341132590924517670116080193774565000394404325176018879137001116501917905458262931179583116640073280342696907129516730267817169709076756630570369478068356496088806733562791025935518005717823975399647433911242565646284290408108965129580005045630325009645600351052252234931946795935673759622200090939484176528145403518844626429027118291618481207436517990068199382319726055853613
Cos    0.9130889403123082724360887896656672808656036779177695185367711006995955759007027812268532869219485750392032411157021985701743449718712839893630769060956987500869050561130530606277722475164033362715938822340958362676305487363307447311046606036732548594394880660779345469357363927362599236287787967113779712125417724313266162079002839798171344659068700365238628373106988573793108841436286988550471435158583090282785551299521559463196485941733155093326870424140702825354668014736012094116310173282063294354797176618106265223818473429154560119875265514193935088602078450686698516654784487568476528092182943197051363283336770971514484898843637969744795863457522041097293912185666895010296292417843629169566
Tan    0.4465725462845951080354340492493297863664702530228539431793643752790652121798457119308281093711337117351824302405127265043941078743201413802742986904431097907418452657673429572010439544368058258668460621337750507856743137674138128805462521543172690127434221187279537588377878373065777689804867923360071570888272946017486343012618629153112816524990960304473577139058272000642627832629543035542998512407586811092722803660439984439880202739952675782181893386116210218071303871409986764914583800844104115292506997715868582578179441467373632677384179542325076687281260631102061490708090583107700757824466976889133117431057537299832543055583282997302563566765119476476011271217072270990986942660758601073145

They all watch too much MSNBC... they get ideas.

 1   2   3   4 


Go to: