Skip to content

Commit d939cb5

Browse files
Tom St Denissjaeckel
Tom St Denis
authored and
sjaeckel
committedOct 19, 2010
added libtomfloat-0.02
·
0.02
1 parent 6c3b48a commit d939cb5

17 files changed

+299
-35
lines changed
 

‎changes.txt

+16-2
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,17 @@
1-
May 5th, 2004 v0.01 -- wrote base of LTF
1+
June 21st, 2004
2+
v0.02 -- Added missing objects to the makefile [oops]
3+
-- fixed up mpf_add and mpf_sub to be more reliable about the precision
4+
-- added required limited ln function mpf_ln_l (domain 0 < x <= 2)
5+
It's still incomplete as it converges slowly (and therefore yields incorrect results)
6+
-- Added mpf_ln and mpf_atan
7+
-- Added short-circuits to sin, cos, invsqrt, atan, ln and sqrt [huge speedup]
8+
-- Optimized mpf_sqrt and mpf_invsqrt by using quick estimates. Fixed
9+
circular dependency as well (mpf_sqrt requires mpf_invsqrt but not vice versa)
10+
++ Note: No further releases are planned for a while. I encourage interested
11+
parties to fork this code base and extend it!
12+
13+
May 5th, 2004
14+
v0.01 -- wrote base of LTF
215

3-
Anytime v0.00 -- no LTF existed.
16+
Anytime
17+
v0.00 -- no LTF existed.

‎demos/ex1.c

+25-1
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,9 @@ void draw(mp_float *a)
1010
int main(void)
1111
{
1212
mp_float a, b, c, d, e;
13-
mpf_init_multi(96, &a, &b, &c, &d, &e, NULL);
13+
int err;
14+
15+
mpf_init_multi(100, &a, &b, &c, &d, &e, NULL);
1416

1517
mpf_const_d(&a, 1); draw(&a);
1618
mpf_const_d(&b, 2); draw(&b);
@@ -21,6 +23,12 @@ int main(void)
2123
mpf_sub(&b, &c, &e); printf("2 - 3 =="); draw(&e);
2224
mpf_mul(&b, &c, &e); printf("2 * 3 == "); draw(&e);
2325
mpf_div(&b, &c, &e); printf("2 / 3 == "); draw(&e);
26+
mpf_add_d(&b, 3, &e); printf("2 + 3 == "); draw(&e);
27+
mpf_sub_d(&b, 3, &e); printf("2 - 3 =="); draw(&e);
28+
mpf_mul_d(&b, 3, &e); printf("2 * 3 == "); draw(&e);
29+
mpf_div_d (&b, 3, &e); printf("2 / 3 == "); draw(&e);
30+
mpf_const_d(&e, 0); mpf_add_d(&e, 1, &e); printf("0 + 1 == "); draw(&e);
31+
mpf_const_d(&e, 0); mpf_sub_d(&e, 1, &e); printf("0 - 1 == "); draw(&e);
2432
printf("\n");
2533
mpf_invsqrt(&d, &e); printf("1/sqrt(4) == 1/2 == "); draw(&e);
2634
mpf_invsqrt(&c, &e); printf("1/sqrt(3) == "); draw(&e);
@@ -29,6 +37,8 @@ int main(void)
2937
mpf_inv(&c, &e); printf("1/3 == "); draw(&e);
3038
mpf_inv(&d, &e); printf("1/4 == "); draw(&e);
3139
printf("\n");
40+
mpf_const_pi(&e); printf("Pi == "); draw(&e);
41+
printf("\n");
3242
mpf_const_e(&e); printf("e == "); draw(&e);
3343
mpf_exp(&c, &e); printf("e^3 == "); draw(&e);
3444
mpf_sqrt(&e, &e); printf("sqrt(e^3) == "); draw(&e);
@@ -46,7 +56,21 @@ int main(void)
4656
mpf_tan(&b, &e); printf("tan(2) == "); draw(&e);
4757
mpf_tan(&c, &e); printf("tan(3) == "); draw(&e);
4858
mpf_tan(&d, &e); printf("tan(4) == "); draw(&e);
59+
mpf_inv(&a, &e); mpf_atan(&e, &e); printf("atan(1/1) == "); draw(&e);
60+
mpf_inv(&b, &e); mpf_atan(&e, &e); printf("atan(1/2) == "); draw(&e);
61+
mpf_inv(&c, &e); mpf_atan(&e, &e); printf("atan(1/3) == "); draw(&e);
62+
mpf_inv(&d, &e); mpf_atan(&e, &e); printf("atan(1/4) == "); draw(&e);
4963
printf("\n");
64+
#define lntest(x) if ((err = mpf_const_ln_d(&e, x)) != MP_OKAY) { printf("Failed ln(%3d), %d\n", x, err); } else { printf("ln(%3d) == ", x); draw(&e); };
65+
lntest(0);
66+
lntest(1);
67+
lntest(2);
68+
lntest(4);
69+
lntest(8);
70+
lntest(17);
71+
lntest(1000);
72+
lntest(100000);
73+
lntest(250000);
5074
return 0;
5175
}
5276

‎float.pdf

-360 Bytes
Binary file not shown.

‎float.tex

+1-1
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@
4949
\begin{document}
5050
\frontmatter
5151
\pagestyle{empty}
52-
\title{LibTomFloat User Manual \\ v0.01}
52+
\title{LibTomFloat User Manual \\ v0.02}
5353
\author{Tom St Denis \\ tomstdenis@iahu.ca}
5454
\maketitle
5555
This text and the library are hereby placed in the public domain. This book has been formatted for B5 [176x250] paper using the \LaTeX{} {\em book}

‎makefile

+3-3
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ default: libtomfloat.a
66

77
CFLAGS += -Os -Wall -W -I./
88

9-
VERSION=0.01
9+
VERSION=0.02
1010

1111
#default files to install
1212
LIBNAME=libtomfloat.a
@@ -24,7 +24,7 @@ DATAPATH=/usr/share/doc/libtomfloat/pdf
2424
OBJECTS = \
2525
mpf_init.o mpf_clear.o mpf_init_multi.o mpf_clear_multi.o mpf_init_copy.o \
2626
\
27-
mpf_copy.o mpf_exch.o \
27+
mpf_copy.o mpf_exch.o mpf_abs.o mpf_neg.o \
2828
\
2929
mpf_cmp.o mpf_cmp_d.o \
3030
\
@@ -70,7 +70,7 @@ install: libtomfloat.a
7070

7171
clean:
7272
rm -f $(OBJECTS) libtomfloat.a *~ demos/*.o demos/*~ ex1
73-
rm -f float.aux float.dvi float.log float.idx float.lof float.out float.toc
73+
rm -f float.aux float.dvi float.log float.idx float.lof float.out float.toc float.ilg float.ind float.pdf
7474

7575
zipup: clean manual
7676
cd .. ; rm -rf ltf* libtomfloat-$(VERSION) ; mkdir libtomfloat-$(VERSION) ; \

‎mpf_add.c

+14
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,20 @@ int mpf_add(mp_float *a, mp_float *b, mp_float *c)
1818
mp_float tmp, *other;
1919
long diff;
2020

21+
if (mpf_iszero(a)) {
22+
diff = c->radix;
23+
if ((err = mpf_copy(b, c)) != MP_OKAY) {
24+
return err;
25+
}
26+
return mpf_normalize_to(c, diff);
27+
} else if (mpf_iszero(b)) {
28+
diff = c->radix;
29+
if ((err = mpf_copy(a, c)) != MP_OKAY) {
30+
return err;
31+
}
32+
return mpf_normalize_to(c, diff);
33+
}
34+
2135
if (a->exp < b->exp) {
2236
/* tmp == a normalize to b's exp */
2337
if ((err = mpf_init_copy(a, &tmp)) != MP_OKAY) {

‎mpf_atan.c

+79
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,85 @@
1212
*/
1313
#include <tomfloat.h>
1414

15+
/* y = y - (tan(y) - x)/(tan(y+0.1)-tan(x)) */
16+
1517
int mpf_atan(mp_float *a, mp_float *b)
1618
{
19+
mp_float oldval, tmp, tmpx, res, sqr;
20+
int oddeven, ires, err, itts;
21+
long n;
22+
23+
/* ensure -1 <= a <= 1 */
24+
if ((err = mpf_cmp_d(a, -1, &ires)) != MP_OKAY) {
25+
return err;
26+
}
27+
if (ires == MP_LT) {
28+
return MP_VAL;
29+
}
30+
31+
if ((err = mpf_cmp_d(a, 1, &ires)) != MP_OKAY) {
32+
return err;
33+
}
34+
if (ires == MP_GT) {
35+
return MP_VAL;
36+
}
37+
38+
/* easy out if a == 0 */
39+
if (mpf_iszero(a) == MP_YES) {
40+
return mpf_const_d(b, 1);
41+
}
42+
43+
/* now a != 0 */
44+
45+
/* initialize temps */
46+
if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
47+
return err;
48+
}
49+
50+
/* initlialize temps */
51+
/* res = 0 */
52+
/* tmpx = 1/a */
53+
if ((err = mpf_inv(a, &tmpx)) != MP_OKAY) { goto __ERR; }
54+
55+
/* sqr = a^2 */
56+
if ((err = mpf_sqr(a, &sqr)) != MP_OKAY) { goto __ERR; }
57+
58+
/* this is the denom counter. Goes up by two per pass */
59+
n = 1;
60+
61+
/* we alternate between adding and subtracting */
62+
oddeven = 0;
63+
64+
/* get number of iterations */
65+
itts = mpf_iterations(b);
66+
67+
while (itts-- > 0) {
68+
if ((err = mpf_copy(&res, &oldval)) != MP_OKAY) { goto __ERR; }
69+
70+
/* compute 1/(2n-1) */
71+
if ((err = mpf_const_d(&tmp, (2*n++ - 1))) != MP_OKAY) { goto __ERR; }
72+
if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY) { goto __ERR; }
73+
74+
/* now multiply a into tmpx twice */
75+
if ((err = mpf_mul(&tmpx, &sqr, &tmpx)) != MP_OKAY) { goto __ERR; }
76+
77+
/* now multiply the two */
78+
if ((err = mpf_mul(&tmpx, &tmp, &tmp)) != MP_OKAY) { goto __ERR; }
79+
80+
/* now depending on if this is even or odd we add/sub */
81+
oddeven ^= 1;
82+
if (oddeven == 1) {
83+
if ((err = mpf_add(&res, &tmp, &res)) != MP_OKAY) { goto __ERR; }
84+
} else {
85+
if ((err = mpf_sub(&res, &tmp, &res)) != MP_OKAY) { goto __ERR; }
86+
}
87+
88+
if (mpf_cmp(&oldval, &res) == MP_EQ) {
89+
break;
90+
}
91+
}
92+
mpf_exch(&res, b);
93+
__ERR: mpf_clear_multi(&oldval, &tmpx, &tmp, &res, &sqr, NULL);
94+
return err;
95+
1796
}

‎mpf_const_ln_d.c

+14
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,20 @@
1515
int mpf_const_ln_d(mp_float *a, long b)
1616
{
1717
int err;
18+
19+
/* test input */
20+
if (b < 0) {
21+
return MP_VAL;
22+
}
23+
24+
if (b == 0) {
25+
return mpf_const_d(a, 1);
26+
}
27+
28+
if (b == 1) {
29+
return mpf_const_d(a, 0);
30+
}
31+
1832
if ((err = mpf_const_d(a, b)) != MP_OKAY) {
1933
return err;
2034
}

‎mpf_cos.c

+8-3
Original file line numberDiff line numberDiff line change
@@ -15,11 +15,11 @@
1515
/* using cos x == \sum_{n=0}^{\infty} ((-1)^n/(2n)!) * x^2n */
1616
int mpf_cos(mp_float *a, mp_float *b)
1717
{
18-
mp_float tmpovern, tmp, tmpx, res, sqr;
18+
mp_float oldval, tmpovern, tmp, tmpx, res, sqr;
1919
int oddeven, err, itts;
2020
long n;
2121
/* initialize temps */
22-
if ((err = mpf_init_multi(b->radix, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
22+
if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
2323
return err;
2424
}
2525

@@ -40,6 +40,7 @@ int mpf_cos(mp_float *a, mp_float *b)
4040
itts = mpf_iterations(b);
4141

4242
while (itts-- > 0) {
43+
if ((err = mpf_copy(&res, &oldval)) != MP_OKAY) { goto __ERR; }
4344
/* compute 1/(2n)! from 1/(2(n-1))! by multiplying by (1/n)(1/(n+1)) */
4445
if ((err = mpf_const_d(&tmp, ++n)) != MP_OKAY) { goto __ERR; }
4546
if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY) { goto __ERR; }
@@ -62,8 +63,12 @@ int mpf_cos(mp_float *a, mp_float *b)
6263
} else {
6364
if ((err = mpf_sub(&res, &tmp, &res)) != MP_OKAY) { goto __ERR; }
6465
}
66+
67+
if (mpf_cmp(&res, &oldval) == MP_EQ) {
68+
break;
69+
}
6570
}
6671
mpf_exch(&res, b);
67-
__ERR: mpf_clear_multi(&tmpx, &tmpovern, &tmp, &res, &sqr, NULL);
72+
__ERR: mpf_clear_multi(&oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL);
6873
return err;
6974
}

‎mpf_exp.c

+9-4
Original file line numberDiff line numberDiff line change
@@ -16,12 +16,12 @@
1616
/* compute b = e^a using e^x == \sum_{n=0}^{\infty} {1 \over n!}x^n */
1717
int mpf_exp(mp_float *a, mp_float *b)
1818
{
19-
mp_float tmpx, tmpovern, tmp, res;
19+
mp_float oldval, tmpx, tmpovern, tmp, res;
2020
int err, itts;
2121
long n;
2222

2323
/* initialize temps */
24-
if ((err = mpf_init_multi(b->radix, &tmpx, &tmpovern, &tmp, &res, NULL)) != MP_OKAY) {
24+
if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmpovern, &tmp, &res, NULL)) != MP_OKAY) {
2525
return err;
2626
}
2727

@@ -36,8 +36,9 @@ int mpf_exp(mp_float *a, mp_float *b)
3636
itts = mpf_iterations(b);
3737

3838
while (itts-- > 0) {
39+
if ((err = mpf_copy(&res, &oldval)) != MP_OKAY) { goto __ERR; }
40+
3941
/* compute 1/n! as 1/(n-1)! * 1/n */
40-
// hack: this won't be portable for n>127
4142
if ((err = mpf_const_d(&tmp, n++)) != MP_OKAY) { goto __ERR; }
4243
if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY) { goto __ERR; }
4344
if ((err = mpf_mul(&tmp, &tmpovern, &tmpovern)) != MP_OKAY) { goto __ERR; }
@@ -48,10 +49,14 @@ int mpf_exp(mp_float *a, mp_float *b)
4849
/* multiply and sum them */
4950
if ((err = mpf_mul(&tmpovern, &tmpx, &tmp)) != MP_OKAY) { goto __ERR; }
5051
if ((err = mpf_add(&tmp, &res, &res)) != MP_OKAY) { goto __ERR; }
52+
53+
if (mpf_cmp(&oldval, &res) == MP_EQ) {
54+
break;
55+
}
5156
}
5257

5358
mpf_exch(&res, b);
54-
__ERR: mpf_clear_multi(&tmpx, &tmpovern, &tmp, &res, NULL);
59+
__ERR: mpf_clear_multi(&oldval, &tmpx, &tmpovern, &tmp, &res, NULL);
5560
return err;
5661
}
5762

‎mpf_invsqrt.c

+13-8
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
/* using newtons method we have 1/sqrt(x) = Y_{n+1} = y_n * ((3 - xy^2_n)/2) */
1616
int mpf_invsqrt(mp_float *a, mp_float *b)
1717
{
18-
mp_float tmp1, tmp2, const_3;
18+
mp_float oldval, tmp1, tmp2, const_3;
1919
int err, itts;
2020

2121
/* ensure a is not zero or negative */
@@ -27,7 +27,7 @@ int mpf_invsqrt(mp_float *a, mp_float *b)
2727
itts = mpf_iterations(b);
2828

2929
/* init temps */
30-
if ((err = mpf_init_multi(b->radix, &tmp1, &tmp2, &const_3, NULL)) != MP_OKAY) {
30+
if ((err = mpf_init_multi(b->radix, &oldval, &tmp1, &tmp2, &const_3, NULL)) != MP_OKAY) {
3131
return err;
3232
}
3333

@@ -36,13 +36,13 @@ int mpf_invsqrt(mp_float *a, mp_float *b)
3636

3737
/* tmp1 == reasonable guess at sqrt */
3838
if ((err = mpf_copy(a, &tmp1)) != MP_OKAY) { goto __ERR; }
39-
40-
/* negate exponent and halve */
41-
tmp1.radix = b->radix;
42-
if ((err = mp_sqrt(&(tmp1.mantissa), &(tmp1.mantissa))) != MP_OKAY) { goto __ERR; }
43-
if ((err = mpf_normalize(&tmp1)) != MP_OKAY) { goto __ERR; }
39+
mp_rshd(&(tmp1.mantissa), tmp1.mantissa.used>>1);
40+
if ((err = mpf_normalize_to(&tmp1, b->radix)) != MP_OKAY) { goto __ERR; }
4441

4542
while (itts-- > 0) {
43+
/* grap copy of tmp1 for early out */
44+
if ((err = mpf_copy(&tmp1, &oldval)) != MP_OKAY) { goto __ERR; }
45+
4646
/* first tmp2 = y^2 == tmp1^2 */
4747
if ((err = mpf_sqr(&tmp1, &tmp2)) != MP_OKAY) { goto __ERR; }
4848
/* multiply by x, tmp1 * a */
@@ -53,10 +53,15 @@ int mpf_invsqrt(mp_float *a, mp_float *b)
5353
if ((err = mpf_div_2(&tmp2, &tmp2)) != MP_OKAY) { goto __ERR; }
5454
/* multiply by y_n and feedback */
5555
if ((err = mpf_mul(&tmp1, &tmp2, &tmp1)) != MP_OKAY) { goto __ERR; }
56+
57+
/* early out if stable */
58+
if (mpf_cmp(&oldval, &tmp1) == MP_EQ) {
59+
break;
60+
}
5661
}
5762

5863
mpf_exch(&tmp1, b);
59-
__ERR: mpf_clear_multi(&tmp1, &tmp2, &const_3, NULL);
64+
__ERR: mpf_clear_multi(&oldval, &tmp1, &tmp2, &const_3, NULL);
6065
return err;
6166
}
6267

‎mpf_ln.c

+66
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,72 @@
1212
*/
1313
#include <tomfloat.h>
1414

15+
/*
16+
17+
Using the newton approximation y1 = y - (e^y - x)/e^y
18+
19+
Which converges quickly and we can reuse e^y so we only calc it once per loop
20+
21+
*/
1522
int mpf_ln(mp_float *a, mp_float *b)
1623
{
24+
mp_float oldval, tmpey, tmpy, val;
25+
int itts, err;
26+
long k;
27+
28+
/* ensure positive */
29+
if (a->mantissa.sign == MP_NEG) {
30+
return MP_VAL;
31+
}
32+
33+
/* easy out for 0 */
34+
if (mpf_iszero(a) == MP_YES) {
35+
return mpf_const_d(b, 1);
36+
}
37+
38+
/* initialize temps */
39+
if ((err = mpf_init_multi(b->radix, &oldval, &tmpey, &tmpy, &val, NULL)) != MP_OKAY) {
40+
return err;
41+
}
42+
43+
/* initial guess */
44+
if ((err = mpf_const_e(&val)) != MP_OKAY) { goto __ERR; }
45+
if ((err = mpf_sqr(&val, &val)) != MP_OKAY) { goto __ERR; }
46+
if ((err = mpf_inv(&val, &tmpey)) != MP_OKAY) { goto __ERR; }
47+
if ((err = mpf_copy(a, &tmpy)) != MP_OKAY) { goto __ERR; }
48+
if ((err = mpf_normalize_to(&tmpy, b->radix)) != MP_OKAY) { goto __ERR; }
49+
50+
/* divide out e's */
51+
k = 0;
52+
while (mpf_cmp(&tmpy, &val) == MP_GT) {
53+
++k;
54+
if ((err = mpf_mul(&tmpy, &tmpey, &tmpy)) != MP_OKAY) { goto __ERR; }
55+
}
56+
if ((err = mpf_const_d(&tmpy, k*2)) != MP_OKAY) { goto __ERR; }
57+
58+
/* number of iterations */
59+
itts = mpf_iterations(b);
60+
61+
while (itts--) {
62+
if ((err = mpf_copy(&tmpy, &oldval)) != MP_OKAY) { goto __ERR; }
63+
64+
/* get e^y and save it */
65+
if ((err = mpf_exp(&tmpy, &tmpey)) != MP_OKAY) { goto __ERR; }
66+
67+
/* now compute e^y - x */
68+
if ((err = mpf_sub(&tmpey, a, &val)) != MP_OKAY) { goto __ERR; }
69+
70+
/* now compute (e^y - x) / e^y */
71+
if ((err = mpf_div(&val, &tmpey, &val)) != MP_OKAY) { goto __ERR; }
72+
73+
/* y = y - (e^y - x)/e^y */
74+
if ((err = mpf_sub(&tmpy, &val, &tmpy)) != MP_OKAY) { goto __ERR; }
75+
76+
if (mpf_cmp(&tmpy, &oldval) == MP_EQ) {
77+
break;
78+
}
79+
}
80+
mpf_exch(&tmpy, b);
81+
__ERR: mpf_clear_multi(&oldval, &tmpey, &tmpy, &val, NULL);
82+
return err;
1783
}

‎mpf_normalize.c

+15-2
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,9 @@
1414

1515
int mpf_normalize(mp_float *a)
1616
{
17-
long cb, diff;
17+
long cb, diff;
18+
int err;
19+
mp_digit c;
1820

1921
/* sanity */
2022
if (a->radix < 2) {
@@ -25,7 +27,18 @@ int mpf_normalize(mp_float *a)
2527
if (cb > a->radix) {
2628
diff = cb - a->radix;
2729
a->exp += diff;
28-
return mp_div_2d(&(a->mantissa), diff, &(a->mantissa), NULL);
30+
31+
/* round it, add 1 after shift if diff-1'th bit is 1 */
32+
c = a->mantissa.dp[diff/DIGIT_BIT] & (1U<<(diff%DIGIT_BIT));
33+
if ((err = mp_div_2d(&(a->mantissa), diff, &(a->mantissa), NULL)) != MP_OKAY) {
34+
return err;
35+
}
36+
37+
if (c != 0) {
38+
return mp_add_d(&(a->mantissa), 1, &(a->mantissa));
39+
} else {
40+
return MP_OKAY;
41+
}
2942
} else if (cb < a->radix) {
3043
if (mp_iszero(&(a->mantissa)) == MP_YES) {
3144
return mpf_const_0(a);

‎mpf_sin.c

+9-3
Original file line numberDiff line numberDiff line change
@@ -15,12 +15,12 @@
1515
/* using sin x == \sum_{n=0}^{\infty} ((-1)^n/(2n+1)!) * x^(2n+1) */
1616
int mpf_sin(mp_float *a, mp_float *b)
1717
{
18-
mp_float tmpovern, tmp, tmpx, res, sqr;
18+
mp_float oldval, tmpovern, tmp, tmpx, res, sqr;
1919
int oddeven, err, itts;
2020
long n;
2121

2222
/* initialize temps */
23-
if ((err = mpf_init_multi(b->radix, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
23+
if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
2424
return err;
2525
}
2626

@@ -49,6 +49,8 @@ int mpf_sin(mp_float *a, mp_float *b)
4949
itts = mpf_iterations(b);
5050

5151
while (itts-- > 0) {
52+
if ((err = mpf_copy(&res, &oldval)) != MP_OKAY) { goto __ERR; }
53+
5254
/* compute 1/(2n)! from 1/(2(n-1))! by multiplying by (1/n)(1/(n+1)) */
5355
if ((err = mpf_const_d(&tmp, ++n)) != MP_OKAY) { goto __ERR; }
5456
if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY) { goto __ERR; }
@@ -71,8 +73,12 @@ int mpf_sin(mp_float *a, mp_float *b)
7173
} else {
7274
if ((err = mpf_sub(&res, &tmp, &res)) != MP_OKAY) { goto __ERR; }
7375
}
76+
77+
if (mpf_cmp(&res, &oldval) == MP_EQ) {
78+
break;
79+
}
7480
}
7581
mpf_exch(&res, b);
76-
__ERR: mpf_clear_multi(&tmpx, &tmpovern, &tmp, &res, &sqr, NULL);
82+
__ERR: mpf_clear_multi(&oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL);
7783
return err;
7884
}

‎mpf_sqrt.c

+11-6
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
/* compute using sqrt(x) = y_{n+1} = 0.5 * (y_{n} + x/y_{n}) */
1616
int mpf_sqrt(mp_float *a, mp_float *b)
1717
{
18-
mp_float tmp, res;
18+
mp_float oldval, tmp, res;
1919
int err, itts;
2020

2121
/* make sure it's positive */
@@ -24,31 +24,36 @@ int mpf_sqrt(mp_float *a, mp_float *b)
2424
}
2525

2626
/* let's roll */
27-
if ((err = mpf_init_multi(b->radix, &tmp, &res, NULL)) != MP_OKAY) {
27+
if ((err = mpf_init_multi(b->radix, &oldval, &tmp, &res, NULL)) != MP_OKAY) {
2828
return err;
2929
}
3030

3131
/* get copy of a and make reasonable guestimte towards the sqrt */
3232
if ((err = mpf_copy(a, &res)) != MP_OKAY) { goto __ERR; }
33-
if ((err = mp_sqrt(&res.mantissa, &res.mantissa)) != MP_OKAY) { goto __ERR; }
33+
mp_rshd(&(res.mantissa), res.mantissa.used >> 1);
3434
res.exp = res.exp / 2;
35-
res.radix = b->radix;
36-
if ((err = mpf_normalize(&res)) != MP_OKAY) { goto __ERR; }
35+
if ((err = mpf_normalize_to(&res, b->radix)) != MP_OKAY) { goto __ERR; }
3736

3837
/* number of iterations */
3938
itts = mpf_iterations(b);
4039

4140
while (itts--) {
41+
if ((err = mpf_copy(&res, &oldval)) != MP_OKAY) { goto __ERR; }
42+
4243
/* compute x/res */
4344
if ((err = mpf_div(a, &res, &tmp)) != MP_OKAY) { goto __ERR; }
4445
/* res + x/res */
4546
if ((err = mpf_add(&res, &tmp, &res)) != MP_OKAY) { goto __ERR; }
4647
/* 0.5 * (res + x/res) */
4748
if ((err = mpf_div_2(&res, &res)) != MP_OKAY) { goto __ERR; }
49+
50+
if (mpf_cmp(&res, &oldval) == MP_EQ) {
51+
break;
52+
}
4853
}
4954

5055
mpf_exch(&res, b);
51-
__ERR: mpf_clear_multi(&tmp, &res, NULL);
56+
__ERR: mpf_clear_multi(&oldval, &tmp, &res, NULL);
5257
return err;
5358
}
5459

‎mpf_sub.c

+15-1
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,20 @@ int mpf_sub(mp_float *a, mp_float *b, mp_float *c)
1818
mp_float tmp;
1919
long diff;
2020

21+
if (mpf_iszero(a)) {
22+
diff = c->radix;
23+
if ((err = mpf_neg(b, c)) != MP_OKAY) {
24+
return err;
25+
}
26+
return mpf_normalize_to(c, diff);
27+
} else if (mpf_iszero(b)) {
28+
diff = c->radix;
29+
if ((err = mpf_copy(a, c)) != MP_OKAY) {
30+
return err;
31+
}
32+
return mpf_normalize_to(c, diff);
33+
}
34+
2135
if (a->exp < b->exp) {
2236
/* tmp == a normalize to b's exp */
2337
if ((err = mpf_init_copy(a, &tmp)) != MP_OKAY) {
@@ -40,7 +54,7 @@ int mpf_sub(mp_float *a, mp_float *b, mp_float *c)
4054
diff = a->exp - tmp.exp;
4155
tmp.exp = a->exp;
4256
if ((err = mp_div_2d(&(tmp.mantissa), diff, (&tmp.mantissa), NULL)) != MP_OKAY) { goto __TMP; }
43-
if ((err = mp_sub(&(a->mantissa), &(tmp.mantissa), &(c->mantissa))) != MP_OKAY) { goto __TMP; }
57+
if ((err = mp_sub(&(a->mantissa), &(tmp.mantissa), &(c->mantissa))) != MP_OKAY) { goto __TMP; }
4458
c->exp = a->exp;
4559
}
4660

‎tomfloat.h

+1-1
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ int mpf_init_copy(mp_float *a, mp_float *b);
3434
int mpf_copy(mp_float *src, mp_float *dest);
3535
void mpf_exch(mp_float *a, mp_float *b);
3636

37-
/* maintainers */
37+
/* maintainers/misc */
3838
int mpf_normalize(mp_float *a);
3939
int mpf_normalize_to(mp_float *a, long radix);
4040
int mpf_iterations(mp_float *a);

0 commit comments

Comments
 (0)
Please sign in to comment.