Comments (20)
moritz
commented on August 16, 2024
Here is a first shot at such a function. I'm neither a good C hacker nor very well versed with libtommath, so beware of bugs -- the few cases I tested seemed to work though.
#include "tommath.h"
#include <math.h>
double mp_get_double(mp_int *a) {
double d = 0.0;
double sign = SIGN(a) == MP_NEG ? -1.0 : 1.0;
if (USED(a) == 0)
return d;
if (USED(a) == 1)
return sign * (double) mp_get_int(a);
int i;
for (i = USED(a) - 1; DIGIT(a, i) == 0 && i > 0; i--) {
/* do nothing */
}
d = (double) DIGIT(a, i);
i--;
if (i == -1) {
return sign * d;
}
d *= pow(2.0, DIGIT_BIT);
d += (double) DIGIT(a, i);
d *= pow(2.0, DIGIT_BIT * i);
return sign * d;
}
from libtommath.
gerdr
commented on August 16, 2024
Depending on the value of DIGIT_BIT, moritz' code loses precision. Barring off-by-one errors, the following code should work in general:
#include <tommath.h>
#include <math.h>
#define PRECISION 53
double mp_get_double(mp_int *a)
{
static const int NEED_DIGITS = (PRECISION + 2 * DIGIT_BIT - 2) / DIGIT_BIT;
static const double DIGIT_MULTI = (mp_digit)1 << DIGIT_BIT;
int i, limit;
double d = 0.0;
mp_clamp(a);
i = USED(a);
limit = i <= NEED_DIGITS ? 0 : i - NEED_DIGITS;
while (i-- > limit) {
d += DIGIT(a, i);
d *= DIGIT_MULTI;
}
if(SIGN(a) == MP_NEG)
d *= -1.0;
d *= pow(2.0, i * DIGIT_BIT);
return d;
}
from libtommath.
hikari-no-yume
commented on August 16, 2024
I'd find this useful for my patch to add native bigints to PHP's Zend Engine.
from libtommath.
hikari-no-yume
commented on August 16, 2024
The reverse, exporting an mp_int to a double would also be handy, to save me from implementing it myself. Especially since I don't want a dependency on LibTomMath's internals... abstraction is my friend.
from libtommath.
czurnieden
commented on August 16, 2024
The reverse [...] would also be handy,
That would be mp_set_double?
#include <math.h>
#include <tommath.h>
int mp_set_double(mp_int * c, double d){
int exp, res, sign;
double frac;
sign = (d < 0)?MP_NEG:MP_ZPOS;
frac = frexp(abs(d), &exp);
if(exp <= 0){
return MP_OKAY;
}
if(exp == 1 && frac < .5){
return MP_OKAY;
}
mp_zero(c);
if(frac == 0){
c->sign = sign;
return MP_OKAY;
}
while(exp-- >= 0){
frac *= 2.0;
if(frac >= 1.0){
if ((res = mp_add_d(c, 1, c)) != MP_OKAY) {
return res;
}
frac -= 1.0;
}
if(exp > 0){
if ((res = mp_mul_2d(c, 1, c)) != MP_OKAY) {
return res;
}
}
}
c->sign = sign;
return MP_OKAY;
}Try it out with e.g.:
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char **argv){
char *eptr, retstr;
double v;
mp_int r;
v = strtod(argv[1],&eptr);
mp_init(&r);
mp_set_double(&r,v);
mp_toradix(&r, &retstr, 10);
printf("%f = %s\n",v,&retstr);
exit(EXIT_SUCCESS);
}The problem might be the rounding mode — it doesn't respect it (doesn't even read it ;-) ) it always rounds to zero (truncates). That's the reason the code is listed here instead of being committed regularly: any code handling floating point values must respect rounding modes.
PS: did the FFT work?
from libtommath.
czurnieden
commented on August 16, 2024
It seems as if
nearbyint()would do here, so
No, it doesn't. Silly me! Sorry.
from libtommath.
czurnieden
commented on August 16, 2024
Nothing in TV, so here wethey are ;-)
https://github.com/czurnieden/libtommath/blob/master/bn_mp_get_double.c
https://github.com/czurnieden/libtommath/blob/master/bn_mp_set_double.c
Please check.
from libtommath.
hikari-no-yume
commented on August 16, 2024
@czurnieden You do a == HUGE_VAL check. On systems where it's available (C99), perhaps better to do == INFINITY?
from libtommath.
czurnieden
commented on August 16, 2024
You do a
== HUGE_VALcheck. On systems where it's available (C99), perhaps better to do== INFINITY?
I just bracketed all c99 stuff out before I came back here ;-)
But one more or less ...
BTW: does the FFT work? I cannot repair something when I do not know that it is broken (a variation of "works for me") save how it is broken.
from libtommath.
sjaeckel
commented on August 16, 2024
I just had an idea how we could include that functionality in ltm without adding a dependency to libm ...
We provide a macro with the implementation, so someone who wants to use it puts that macro in his code and then he can use the implementation...
I only proposed this so someone else can propose something better ;-)
I revived this because I think that it's useful and we should just do it somehow (but without adding a dependency to libm)
from libtommath.
czurnieden
commented on August 16, 2024
mp_get_double (convert mp_int to double) does not need the libm but mp_set_double (double to mp_int) does use frexp (isnan, isinf, and nearbyint if available). It is not that complicated to write a frexp if you know the endinaness in advance, otherwise it is quite tedious (vid. e.g.: the implementation in the coreutils in coreutils-8.24/lib/frexp.c ), so it's a rewrite without frexp et al.?
from libtommath.
czurnieden
commented on August 16, 2024
Done (for big and little endian only), to be found under the same adress. Pleeeeease check before use.
from libtommath.
hikari-no-yume
commented on August 16, 2024
What's wrong with depending on libm? math.h is part of the C standard library.
from libtommath.
czurnieden
commented on August 16, 2024
@TazeTSchnitzel LibTomMath gets used quite widely, especially on embedded stuff. Those PICs do not necessarily have a FPU--some don't even support floats, they have to do it all with integers. So it might be already slow to use floating point math. On top of that: they don't have a lot of memory available, loading a library for one or two of its functions is a large overhead that might not be acceptable or it is even impossible to load in the first place because of lack of memory.
And while I'm writing that I get increasingly uncomfortable with the bit-pushing in my implementation of frexp(3). I think I will do it the long way instead. But not today ;-)
from libtommath.
hikari-no-yume
commented on August 16, 2024
I can understand the embedded use-case. Hmm.
Perhaps we could have an option to not build the floating-point parts of the library. Or maybe we could use libm if available... oh I don't know.
from libtommath.
czurnieden
commented on August 16, 2024
I did a frexp implementation the long way and without libm. Please run to test. Needs libm but only for the tests because I use libm's frexp to check against. The problem is the check for NaN and Infinity without libm. Some compilers might protest against the way I did it.
Code attached.
Or not?
Well ... sigh ... sorry for that mess.
double frexp_intern(double x, int *eptr)
{
int sign, exponent;
int i;
/*
* The exponent of an IEEE-754 double (binary64) is an 11-bit large integer
*/
double ap_2[11] = {
2.0000000000000000000000000000000000000,
4.0000000000000000000000000000000000000,
16.000000000000000000000000000000000000,
256.00000000000000000000000000000000000,
65536.000000000000000000000000000000000,
4294967296.0000000000000000000000000000,
18446744073709551616.000000000000000000,
3.4028236692093846346337460743176821146e38,
1.1579208923731619542357098500868790785e77,
1.3407807929942597099574024998205846128e154,
1.7976931348623157e308 // DBL_MAX
};
double ap_half[11] = {
0.50000000000000000000000000000000000000,
0.25000000000000000000000000000000000000,
0.062500000000000000000000000000000000000,
0.0039062500000000000000000000000000000000,
1.5258789062500000000000000000000000000e-5,
2.3283064365386962890625000000000000000e-10,
5.4210108624275221700372640043497085571e-20,
2.9387358770557187699218413430556141946e-39,
8.6361685550944446253863518628003995711e-78,
7.4583407312002067432909653154629338374e-155,
5.5626846462680034577255817933310101606e-309 // < DBL_MIN
};
// TODO: not every compiler might eat this check for Inf and NaN
// GCC-4.8.4 does
// TCC 0.9.25 does
// clang 3.4-1ubuntu3 (based on LLVM 3.4) does
if (x == 1.0 / 0.0) {
*eptr = 0;
return x;
}
if (x == 0.0 / 0.0) {
*eptr = 0;
return x;
}
if (x == 0.0) {
*eptr = 0;
return x;
}
exponent = 0.0;
/*
* Easier to work with positive values
*/
if (x < 0) {
x = -x;
sign = 1;
} else {
sign = 0;
}
if (x >= 1.0) {
/*
* Big steps
*/
for (i = 0; x >= ap_2[i]; i++) {
exponent += (1 << i);
x *= ap_half[i];
}
/*
* Small steps
*/
if (x < 0.5) {
while (x < 0.5) {
x *= 2.0;
exponent--;
}
} else {
while (x > 1.0) {
x /= 2.0;
exponent++;
}
}
} else {
/*
* Same as above, but in the opposite direction
*/
for (i = 0; x < ap_half[i]; i++) {
exponent -= (1 << i);
x *= ap_2[i];
}
if (x < 0.5) {
while (x < 0.5) {
x *= 2.0;
exponent--;
}
} else {
while (x > 1.0) {
x /= 2.0;
exponent++;
}
}
}
if (sign) {
x = -x;
}
*eptr = exponent;
return x;
}
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
// for the internal frexp() only
#include <math.h>
int main()
{
int i, e1, e2;
double d1, d2, t, f1, f2;
srand((unsigned int) time(NULL));
t = 0.0;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("0 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("0 MANT NAN:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
// NaN
// Will fail MANT-test because NaN != NaN
t = 0.0 / 0.0;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("1 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("1 MANT NAN:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
//exit(EXIT_FAILURE);
}
// Infinity
t = 1e308;
t *= t;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("2 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("2 MANT FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
// -Infinity
t = -1e308;
t *= t;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("2 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("2 MANT FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
// subnormal (denormal) numbers
t = 4.e-324;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("2 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("2 MANT FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
for (i = 0; i < 1000000; i++) {
d1 = (double) rand();
d2 = (double) rand();
t = d1;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("3 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("3 MANT FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
t = 1 / d1;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("4 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("4 MANT FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
t = d2 / d1;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("5 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("5 MANT FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
t = d1 / d2;
f1 = frexp_intern(t, &e1);
f2 = frexp(t, &e2);
if (e1 != e2) {
printf("6 EXPO FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
if (f1 != f2) {
printf("6 MANT FAIL:t = %g, e1 = %d, e2 = %d, f1 = %g, f2 = %g\n", t, e1,
e2, f1, f2);
exit(EXIT_FAILURE);
}
}
exit(EXIT_SUCCESS);
}from libtommath.
tomstdenis
commented on August 16, 2024
This is not really a 0.43 candidate so let's not merge it into develop at all for now. As for the general comment ... It'd be nice to not have any float code "standard" in LTM. I'd be ok with bundling it with the package but not with making it a core API component.
To be fair though embedded should be using TFM unless they have some specific requirement from LTM.
from libtommath.
sjaeckel
commented on August 16, 2024
I also don't mind delaying that one
from libtommath.
minad
commented on August 16, 2024
@sjaeckel I guess this can be closed.
from libtommath.
sjaeckel
commented on August 16, 2024
closed via #123
from libtommath.
Related Issues (20)
- Build fails on Linux/Sparc with 32 bit userspace HOT 3
- CryptAcquireContextW and CryptGenRandom are deprecated HOT 3
- Potentially lossy conversion in s_read_wincsp HOT 1
- Dependency on dead code elimination HOT 4
- s_read_wincsp leaks handle HOT 2
- Is the memory representation stable? HOT 6
- mp_is_square says 0 is not a square HOT 3
- Library can not be used in android (arm64) HOT 24
- Next release with SPM support
- tommath_c89.h versus other solutions for backward compatibility with <C99 HOT 2
- should i>=x.used better? HOT 14
- does it support cross-build? HOT 1
- typos HOT 6
- s_mp_invmod_odd returns wrong result for negative numbers HOT 2
- speed please!? HOT 1
- Compiler & base OS versions
- mp_fwrite should not write \0 character HOT 3
- cmake HOT 2
- mp_exptmod incorrect result since version 0.32 HOT 3
- No more tommath.pdf ?
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from libtommath.