Skip to content

Commit ffd0519

Browse files
guihuazCopilot
andauthored
Fix #222: Wrong result for e.g. √(2.25) - 1.5 #222 (#2424)
* add snap zero in add/sub/log * add a comment * Update src/CalcManager/Ratpack/exp.cpp Co-authored-by: Copilot <[email protected]> * Update src/CalcManager/Ratpack/rat.cpp Co-authored-by: Copilot <[email protected]> * Update src/CalcManager/Ratpack/rat.cpp Co-authored-by: Copilot <[email protected]> --------- Co-authored-by: Copilot <[email protected]>
1 parent aab8915 commit ffd0519

File tree

11 files changed

+221
-84
lines changed

11 files changed

+221
-84
lines changed

src/CalcManager/Ratpack/exp.cpp

Lines changed: 31 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,7 @@ void exprat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
7878
const int32_t intpwr = rattoi32(pint, radix, precision);
7979
ratpowi32(&pwr, intpwr, precision);
8080

81-
subrat(px, pint, precision);
81+
_subrat(px, pint, precision);
8282

8383
// It just so happens to be an integral power of e.
8484
if (rat_gt(*px, rat_negsmallest, precision) && rat_lt(*px, rat_smallest, precision))
@@ -97,7 +97,7 @@ void exprat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
9797

9898
//-----------------------------------------------------------------------------
9999
//
100-
// FUNCTION: lograt, _lograt
100+
// FUNCTION: lograt, _lograt, __lograt
101101
//
102102
// ARGUMENTS: x PRAT representation of number to logarithim
103103
//
@@ -119,10 +119,13 @@ void exprat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
119119
// Number is scaled between one and e_to_one_half prior to taking the
120120
// log. This is to keep execution time from exploding.
121121
//
122+
// lograt tries to snap to zero. Use _lograt inside ratpak by default.
123+
// __lograt is part of _lograt private implementation and should not be used.
124+
//
122125
//
123126
//-----------------------------------------------------------------------------
124127

125-
void _lograt(PRAT* px, int32_t precision)
128+
void __lograt(PRAT* px, int32_t precision)
126129

127130
{
128131
CREATETAYLOR();
@@ -149,7 +152,7 @@ void _lograt(PRAT* px, int32_t precision)
149152
DESTROYTAYLOR();
150153
}
151154

152-
void lograt(_Inout_ PRAT* px, int32_t precision)
155+
void _lograt(_Inout_ PRAT* px, int32_t precision)
153156

154157
{
155158
PRAT pwr = nullptr; // pwr is the large scaling factor.
@@ -192,18 +195,18 @@ void lograt(_Inout_ PRAT* px, int32_t precision)
192195
while (rat_gt(*px, e_to_one_half, precision))
193196
{
194197
divrat(px, e_to_one_half, precision);
195-
addrat(&offset, rat_one, precision);
198+
_addrat(&offset, rat_one, precision);
196199
}
197200

198-
_lograt(px, precision);
201+
__lograt(px, precision);
199202

200203
// Add the large and small scaling factors, take into account
201204
// small scaling was done in e_to_one_half chunks.
202205
divrat(&offset, rat_two, precision);
203-
addrat(&pwr, offset, precision);
206+
_addrat(&pwr, offset, precision);
204207

205208
// And add the resulting scaling factor to the answer.
206-
addrat(px, pwr, precision);
209+
_addrat(px, pwr, precision);
207210

208211
trimit(px, precision);
209212

@@ -217,6 +220,18 @@ void lograt(_Inout_ PRAT* px, int32_t precision)
217220
destroyrat(pwr);
218221
}
219222

223+
void lograt(_Inout_ PRAT* px, int32_t precision)
224+
225+
{
226+
PRAT a = nullptr;
227+
DUPRAT(a, *px);
228+
229+
_lograt(px, precision);
230+
231+
_snaprat(px, a, nullptr, precision);
232+
destroyrat(a);
233+
}
234+
220235
void log10rat(_Inout_ PRAT* px, int32_t precision)
221236

222237
{
@@ -235,8 +250,8 @@ bool IsEven(PRAT x, uint32_t radix, int32_t precision)
235250
DUPRAT(tmp, x);
236251
divrat(&tmp, rat_two, precision);
237252
fracrat(&tmp, radix, precision);
238-
addrat(&tmp, tmp, precision);
239-
subrat(&tmp, rat_one, precision);
253+
_addrat(&tmp, tmp, precision);
254+
_subrat(&tmp, rat_one, precision);
240255
if (rat_lt(tmp, rat_zero, precision))
241256
{
242257
bRet = true;
@@ -339,11 +354,11 @@ void powratNumeratorDenominator(_Inout_ PRAT* px, _In_ PRAT y, uint32_t radix, i
339354
DUPRAT(roundedResult, originalResult);
340355
if (roundedResult->pp->sign == -1)
341356
{
342-
subrat(&roundedResult, rat_half, precision);
357+
_subrat(&roundedResult, rat_half, precision);
343358
}
344359
else
345360
{
346-
addrat(&roundedResult, rat_half, precision);
361+
_addrat(&roundedResult, rat_half, precision);
347362
}
348363
intrat(&roundedResult, radix, precision);
349364

@@ -423,7 +438,7 @@ void powratcomp(_Inout_ PRAT* px, _In_ PRAT y, uint32_t radix, int32_t precision
423438
{
424439
PRAT pxint = nullptr;
425440
DUPRAT(pxint, *px);
426-
subrat(&pxint, rat_one, precision);
441+
_subrat(&pxint, rat_one, precision);
427442
if (rat_gt(pxint, rat_negsmallest, precision) && rat_lt(pxint, rat_smallest, precision) && (sign == 1))
428443
{
429444
// *px is one, special case a 1 return.
@@ -442,12 +457,12 @@ void powratcomp(_Inout_ PRAT* px, _In_ PRAT y, uint32_t radix, int32_t precision
442457
// If power is an integer let ratpowi32 deal with it.
443458
PRAT iy = nullptr;
444459
DUPRAT(iy, y);
445-
subrat(&iy, podd, precision);
460+
_subrat(&iy, podd, precision);
446461
int32_t inty = rattoi32(iy, radix, precision);
447462

448463
PRAT plnx = nullptr;
449464
DUPRAT(plnx, *px);
450-
lograt(&plnx, precision);
465+
_lograt(&plnx, precision);
451466
mulrat(&plnx, iy, precision);
452467
if (rat_gt(plnx, rat_max_exp, precision) || rat_lt(plnx, rat_min_exp, precision))
453468
{
@@ -515,7 +530,7 @@ void powratcomp(_Inout_ PRAT* px, _In_ PRAT y, uint32_t radix, int32_t precision
515530
sign = 1;
516531
}
517532

518-
lograt(px, precision);
533+
_lograt(px, precision);
519534
mulrat(px, y, precision);
520535
exprat(px, radix, precision);
521536
}

src/CalcManager/Ratpack/fact.cpp

Lines changed: 19 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,6 @@
1515
//-----------------------------------------------------------------------------
1616
#include "ratpak.h"
1717

18-
#define ABSRAT(x) (((x)->pp->sign = 1), ((x)->pq->sign = 1))
1918
#define NEGATE(x) ((x)->pp->sign *= -1)
2019

2120
//-----------------------------------------------------------------------------
@@ -73,35 +72,35 @@ void _gamma(PRAT* pn, uint32_t radix, int32_t precision)
7372

7473
// Find the best 'A' for convergence to the required precision.
7574
a = i32torat(radix);
76-
lograt(&a, precision);
75+
_lograt(&a, precision);
7776
mulrat(&a, ratprec, precision);
7877

7978
// Really is -ln(n)+1, but -ln(n) will be < 1
8079
// if we scale n between 0.5 and 1.5
81-
addrat(&a, rat_two, precision);
80+
_addrat(&a, rat_two, precision);
8281
DUPRAT(tmp, a);
83-
lograt(&tmp, precision);
82+
_lograt(&tmp, precision);
8483
mulrat(&tmp, *pn, precision);
85-
addrat(&a, tmp, precision);
86-
addrat(&a, rat_one, precision);
84+
_addrat(&a, tmp, precision);
85+
_addrat(&a, rat_one, precision);
8786

8887
// Calculate the necessary bump in precision and up the precision.
8988
// The following code is equivalent to
9089
// precision += ln(exp(a)*pow(a,n+1.5))-ln(radix));
9190
DUPRAT(tmp, *pn);
9291
one_pt_five = i32torat(3L);
9392
divrat(&one_pt_five, rat_two, precision);
94-
addrat(&tmp, one_pt_five, precision);
93+
_addrat(&tmp, one_pt_five, precision);
9594
DUPRAT(term, a);
9695
powratcomp(&term, tmp, radix, precision);
9796
DUPRAT(tmp, a);
9897
exprat(&tmp, radix, precision);
9998
mulrat(&term, tmp, precision);
100-
lograt(&term, precision);
99+
_lograt(&term, precision);
101100
const auto ratRadix = i32torat(radix);
102101
DUPRAT(tmp, ratRadix);
103-
lograt(&tmp, precision);
104-
subrat(&term, tmp, precision);
102+
_lograt(&tmp, precision);
103+
_subrat(&term, tmp, precision);
105104
precision += rattoi32(term, radix, precision);
106105

107106
// Set up initial terms for series, refer to series in above comment block.
@@ -118,10 +117,10 @@ void _gamma(PRAT* pn, uint32_t radix, int32_t precision)
118117
DUPRAT(sum, rat_one);
119118
divrat(&sum, *pn, precision);
120119
DUPRAT(tmp, *pn);
121-
addrat(&tmp, rat_one, precision);
120+
_addrat(&tmp, rat_one, precision);
122121
DUPRAT(term, a);
123122
divrat(&term, tmp, precision);
124-
subrat(&sum, term, precision);
123+
_subrat(&sum, term, precision);
125124

126125
DUPRAT(err, ratRadix);
127126
NEGATE(ratprec);
@@ -134,7 +133,7 @@ void _gamma(PRAT* pn, uint32_t radix, int32_t precision)
134133
// Loop until precision is reached, or asked to halt.
135134
while (!zerrat(term) && rat_gt(term, err, precision))
136135
{
137-
addrat(pn, rat_two, precision);
136+
_addrat(pn, rat_two, precision);
138137

139138
// WARNING: mixing numbers and rationals here.
140139
// for speed and efficiency.
@@ -146,22 +145,22 @@ void _gamma(PRAT* pn, uint32_t radix, int32_t precision)
146145
divrat(&factorial, a2, precision);
147146

148147
DUPRAT(tmp, *pn);
149-
addrat(&tmp, rat_one, precision);
148+
_addrat(&tmp, rat_one, precision);
150149
destroyrat(term);
151150
createrat(term);
152151
DUPNUM(term->pp, count);
153152
DUPNUM(term->pq, num_one);
154-
addrat(&term, rat_one, precision);
153+
_addrat(&term, rat_one, precision);
155154
mulrat(&term, tmp, precision);
156155
DUPRAT(tmp, a);
157156
divrat(&tmp, term, precision);
158157

159158
DUPRAT(term, rat_one);
160159
divrat(&term, *pn, precision);
161-
subrat(&term, tmp, precision);
160+
_subrat(&term, tmp, precision);
162161

163162
divrat(&term, factorial, precision);
164-
addrat(&sum, term, precision);
163+
_addrat(&sum, term, precision);
165164
ABSRAT(term);
166165
}
167166

@@ -214,7 +213,7 @@ void factrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
214213
while (rat_gt(*px, rat_zero, precision) && (LOGRATRADIX(*px) > -precision))
215214
{
216215
mulrat(&fact, *px, precision);
217-
subrat(px, rat_one, precision);
216+
_subrat(px, rat_one, precision);
218217
}
219218

220219
// Added to make numbers 'close enough' to integers use integer factorial.
@@ -226,13 +225,13 @@ void factrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
226225

227226
while (rat_lt(*px, neg_rat_one, precision))
228227
{
229-
addrat(px, rat_one, precision);
228+
_addrat(px, rat_one, precision);
230229
divrat(&fact, *px, precision);
231230
}
232231

233232
if (rat_neq(*px, rat_zero, precision))
234233
{
235-
addrat(px, rat_one, precision);
234+
_addrat(px, rat_one, precision);
236235
_gamma(px, radix, precision);
237236
mulrat(px, fact, precision);
238237
}

src/CalcManager/Ratpack/itrans.cpp

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -95,7 +95,7 @@ void asinrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
9595

9696
// Avoid the really bad part of the asin curve near +/-1.
9797
DUPRAT(phack, *px);
98-
subrat(&phack, rat_one, precision);
98+
_subrat(&phack, rat_one, precision);
9999
// Since *px might be epsilon near zero we must set it to zero.
100100
if (rat_le(phack, rat_smallest, precision) && rat_ge(phack, rat_negsmallest, precision))
101101
{
@@ -109,7 +109,7 @@ void asinrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
109109
{
110110
if (rat_gt(*px, rat_one, precision))
111111
{
112-
subrat(px, rat_one, precision);
112+
_subrat(px, rat_one, precision);
113113
if (rat_gt(*px, rat_smallest, precision))
114114
{
115115
throw(CALC_E_DOMAIN);
@@ -122,11 +122,11 @@ void asinrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
122122
DUPRAT(pret, *px);
123123
mulrat(px, pret, precision);
124124
(*px)->pp->sign *= -1;
125-
addrat(px, rat_one, precision);
125+
_addrat(px, rat_one, precision);
126126
rootrat(px, rat_two, radix, precision);
127127
_asinrat(px, precision);
128128
(*px)->pp->sign *= -1;
129-
addrat(px, pi_over_two, precision);
129+
_addrat(px, pi_over_two, precision);
130130
destroyrat(pret);
131131
}
132132
else
@@ -214,7 +214,7 @@ void acosrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
214214
(*px)->pp->sign = sgn;
215215
asinrat(px, radix, precision);
216216
(*px)->pp->sign *= -1;
217-
addrat(px, pi_over_two, precision);
217+
_addrat(px, pi_over_two, precision);
218218
}
219219
}
220220

@@ -297,15 +297,15 @@ void atanrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
297297
tmpx->pp->sign = sgn;
298298
tmpx->pq->sign = 1;
299299
DUPRAT(*px, pi_over_two);
300-
subrat(px, tmpx, precision);
300+
_subrat(px, tmpx, precision);
301301
destroyrat(tmpx);
302302
}
303303
else
304304
{
305305
(*px)->pp->sign = sgn;
306306
DUPRAT(tmpx, *px);
307307
mulrat(&tmpx, *px, precision);
308-
addrat(&tmpx, rat_one, precision);
308+
_addrat(&tmpx, rat_one, precision);
309309
rootrat(&tmpx, rat_two, radix, precision);
310310
divrat(px, tmpx, precision);
311311
destroyrat(tmpx);
@@ -322,6 +322,6 @@ void atanrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
322322
}
323323
if (rat_gt(*px, pi_over_two, precision))
324324
{
325-
subrat(px, pi, precision);
325+
_subrat(px, pi, precision);
326326
}
327327
}

src/CalcManager/Ratpack/itransh.cpp

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -59,10 +59,10 @@ void asinhrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
5959
PRAT ptmp = nullptr;
6060
DUPRAT(ptmp, (*px));
6161
mulrat(&ptmp, *px, precision);
62-
addrat(&ptmp, rat_one, precision);
62+
_addrat(&ptmp, rat_one, precision);
6363
rootrat(&ptmp, rat_two, radix, precision);
64-
addrat(px, ptmp, precision);
65-
lograt(px, precision);
64+
_addrat(px, ptmp, precision);
65+
_lograt(px, precision);
6666
destroyrat(ptmp);
6767
}
6868
else
@@ -113,10 +113,10 @@ void acoshrat(_Inout_ PRAT* px, uint32_t radix, int32_t precision)
113113
PRAT ptmp = nullptr;
114114
DUPRAT(ptmp, (*px));
115115
mulrat(&ptmp, *px, precision);
116-
subrat(&ptmp, rat_one, precision);
116+
_subrat(&ptmp, rat_one, precision);
117117
rootrat(&ptmp, rat_two, radix, precision);
118-
addrat(px, ptmp, precision);
119-
lograt(px, precision);
118+
_addrat(px, ptmp, precision);
119+
_lograt(px, precision);
120120
destroyrat(ptmp);
121121
}
122122
}
@@ -143,11 +143,11 @@ void atanhrat(_Inout_ PRAT* px, int32_t precision)
143143
{
144144
PRAT ptmp = nullptr;
145145
DUPRAT(ptmp, (*px));
146-
subrat(&ptmp, rat_one, precision);
147-
addrat(px, rat_one, precision);
146+
_subrat(&ptmp, rat_one, precision);
147+
_addrat(px, rat_one, precision);
148148
divrat(px, ptmp, precision);
149149
(*px)->pp->sign *= -1;
150-
lograt(px, precision);
150+
_lograt(px, precision);
151151
divrat(px, rat_two, precision);
152152
destroyrat(ptmp);
153153
}

src/CalcManager/Ratpack/logic.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -246,7 +246,7 @@ void modrat(_Inout_ PRAT* pa, _In_ PRAT b)
246246

247247
if (needAdjust && !zerrat(*pa))
248248
{
249-
addrat(pa, b, BASEX);
249+
_addrat(pa, b, BASEX);
250250
}
251251

252252
// Get *pa back in the integer over integer form.

0 commit comments

Comments
 (0)