Many hyperlinks are disabled.
Use anonymous login
to enable hyperlinks.
Overview
| Comment: | Only use special mp_sqrt() code when double format/tommath format are exactly what's expected. Otherwise, use original always-working tommath code. Simplify overflow check in bignum expononent code, not using bignums where it's not necessary. Don't overallocate bignums when using wideint's only. |
|---|---|
| Downloads: | Tarball | ZIP archive | SQL archive |
| Timelines: | family | ancestors | descendants | both | core-8-6-branch |
| Files: | files | file ages | folders |
| SHA3-256: | 3f35b523558b49db296fa34756dad9bd |
| User & Date: | jan.nijtmans 2019-04-11 20:09:27 |
Context
|
2019-04-17
| ||
| 14:28 | Isolate tests of [info frame] results from testing environment. check-in: 572f113bbb user: dgp tags: core-8-6-branch | |
|
2019-04-11
| ||
| 20:37 | Merge 8.6 check-in: 2a6c012bff user: jan.nijtmans tags: core-8-branch | |
| 20:09 | Only use special mp_sqrt() code when double format/tommath format are exactly what's expected. Other... check-in: 3f35b52355 user: jan.nijtmans tags: core-8-6-branch | |
|
2019-04-09
| ||
| 19:39 | merge 8.5 check-in: 4cb9044dfa user: sebres tags: core-8-6-branch | |
Changes
Changes to generic/tclExecute.c.
6028
6029
6030
6031
6032
6033
6034
6035
6036
6037
6038
6039
6040
6041
6042
6043
6044
6045
6046
6047
6048
6049
6050
6051
6052
6053
6054
6055
6056
6057
6058
6059
....
8980
8981
8982
8983
8984
8985
8986
8987
8988
8989
8990
8991
8992
8993
8994
8995
8996
8997
8998
8999
9000
9001
9002
9003
9004
9005
9006
9007
9008
9009
|
if (GetNumberFromObj(NULL, OBJ_AT_TOS, &ptr1, &type1) != TCL_OK) {
type1 = 0;
} else if (type1 == TCL_NUMBER_LONG) {
/* value is between LONG_MIN and LONG_MAX */
/* [string is integer] is -UINT_MAX to UINT_MAX range */
int i;
if (Tcl_GetIntFromObj(NULL, OBJ_AT_TOS, &i) != TCL_OK) {
type1 = TCL_NUMBER_WIDE;
}
#ifndef TCL_WIDE_INT_IS_LONG
} else if (type1 == TCL_NUMBER_WIDE) {
/* value is between WIDE_MIN and WIDE_MAX */
/* [string is wideinteger] is -UWIDE_MAX to UWIDE_MAX range */
int i;
if (Tcl_GetIntFromObj(NULL, OBJ_AT_TOS, &i) == TCL_OK) {
type1 = TCL_NUMBER_LONG;
}
#endif
} else if (type1 == TCL_NUMBER_BIG) {
/* value is an integer outside the WIDE_MIN to WIDE_MAX range */
/* [string is wideinteger] is -UWIDE_MAX to UWIDE_MAX range */
Tcl_WideInt w;
if (Tcl_GetWideIntFromObj(NULL, OBJ_AT_TOS, &w) == TCL_OK) {
type1 = TCL_NUMBER_WIDE;
}
}
TclNewIntObj(objResultPtr, type1);
TRACE(("\"%.20s\" => %d\n", O2S(OBJ_AT_TOS), type1));
NEXT_INST_F(1, 1, 1);
................................................................................
wResult = oddExponent ? -Exp64Value[base] : Exp64Value[base];
WIDE_RESULT(wResult);
}
}
#endif
overflowExpon:
Tcl_TakeBignumFromObj(NULL, value2Ptr, &big2);
if ((big2.used > 1)
#if DIGIT_BIT > 28
|| ((big2.used == 1) && (big2.dp[0] >= (1<<28)))
#endif
) {
mp_clear(&big2);
Tcl_SetObjResult(interp, Tcl_NewStringObj(
"exponent too large", -1));
return GENERAL_ARITHMETIC_ERROR;
}
Tcl_TakeBignumFromObj(NULL, valuePtr, &big1);
mp_init(&bigResult);
mp_expt_d(&big1, big2.dp[0], &bigResult);
mp_clear(&big1);
mp_clear(&big2);
BIG_RESULT(&bigResult);
}
case INST_ADD:
case INST_SUB:
case INST_MULT:
case INST_DIV:
|
|
|
|
>
|
<
|
|
<
<
<
|
<
|
6028
6029
6030
6031
6032
6033
6034
6035
6036
6037
6038
6039
6040
6041
6042
6043
6044
6045
6046
6047
6048
6049
6050
6051
6052
6053
6054
6055
6056
6057
6058
6059
....
8980
8981
8982
8983
8984
8985
8986
8987
8988
8989
8990
8991
8992
8993
8994
8995
8996
8997
8998
8999
9000
9001
9002
9003
9004
9005
|
if (GetNumberFromObj(NULL, OBJ_AT_TOS, &ptr1, &type1) != TCL_OK) {
type1 = 0;
} else if (type1 == TCL_NUMBER_LONG) {
/* value is between LONG_MIN and LONG_MAX */
/* [string is integer] is -UINT_MAX to UINT_MAX range */
int i;
if (TclGetIntFromObj(NULL, OBJ_AT_TOS, &i) != TCL_OK) {
type1 = TCL_NUMBER_WIDE;
}
#ifndef TCL_WIDE_INT_IS_LONG
} else if (type1 == TCL_NUMBER_WIDE) {
/* value is between WIDE_MIN and WIDE_MAX */
/* [string is wideinteger] is -UWIDE_MAX to UWIDE_MAX range */
int i;
if (TclGetIntFromObj(NULL, OBJ_AT_TOS, &i) == TCL_OK) {
type1 = TCL_NUMBER_LONG;
}
#endif
} else if (type1 == TCL_NUMBER_BIG) {
/* value is an integer outside the WIDE_MIN to WIDE_MAX range */
/* [string is wideinteger] is -UWIDE_MAX to UWIDE_MAX range */
Tcl_WideInt w;
if (TclGetWideIntFromObj(NULL, OBJ_AT_TOS, &w) == TCL_OK) {
type1 = TCL_NUMBER_WIDE;
}
}
TclNewIntObj(objResultPtr, type1);
TRACE(("\"%.20s\" => %d\n", O2S(OBJ_AT_TOS), type1));
NEXT_INST_F(1, 1, 1);
................................................................................
wResult = oddExponent ? -Exp64Value[base] : Exp64Value[base];
WIDE_RESULT(wResult);
}
}
#endif
overflowExpon:
if ((TclGetWideIntFromObj(NULL, value2Ptr, &w2) != TCL_OK)
|| (value2Ptr->typePtr != &tclIntType)
|| (Tcl_WideUInt)w2 >= (1<<28)) {
Tcl_SetObjResult(interp, Tcl_NewStringObj(
"exponent too large", -1));
return GENERAL_ARITHMETIC_ERROR;
}
Tcl_TakeBignumFromObj(NULL, valuePtr, &big1);
mp_init(&bigResult);
mp_expt_d(&big1, w2, &bigResult);
mp_clear(&big1);
BIG_RESULT(&bigResult);
}
case INST_ADD:
case INST_SUB:
case INST_MULT:
case INST_DIV:
|
Changes to generic/tclTomMathInterface.c.
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
...
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
|
unsigned long v;
mp_digit *p;
/*
* Allocate enough memory to hold the largest possible long
*/
status = mp_init_size(a,
(CHAR_BIT * sizeof(long) + DIGIT_BIT - 1) / DIGIT_BIT);
if (status != MP_OKAY) {
Tcl_Panic("initialization failure in TclBNInitBignumFromLong");
}
/*
* Convert arg to sign and magnitude.
*/
................................................................................
int status;
mp_digit *p;
/*
* Allocate enough memory to hold the largest possible Tcl_WideUInt.
*/
status = mp_init_size(a,
(CHAR_BIT * sizeof(Tcl_WideUInt) + DIGIT_BIT - 1) / DIGIT_BIT);
if (status != MP_OKAY) {
Tcl_Panic("initialization failure in TclBNInitBignumFromWideUInt");
}
a->sign = 0;
/*
|
|
<
|
<
|
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
...
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
|
unsigned long v;
mp_digit *p;
/*
* Allocate enough memory to hold the largest possible long
*/
status = mp_init(a);
if (status != MP_OKAY) {
Tcl_Panic("initialization failure in TclBNInitBignumFromLong");
}
/*
* Convert arg to sign and magnitude.
*/
................................................................................
int status;
mp_digit *p;
/*
* Allocate enough memory to hold the largest possible Tcl_WideUInt.
*/
status = mp_init(a);
if (status != MP_OKAY) {
Tcl_Panic("initialization failure in TclBNInitBignumFromWideUInt");
}
a->sign = 0;
/*
|
Changes to libtommath/bn_mp_set_double.c.
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 |
*
* SPDX-License-Identifier: Unlicense
*/
#if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559)
int mp_set_double(mp_int *a, double b)
{
uint64_t frac;
int exp, res;
union {
double dbl;
uint64_t bits;
} cast;
cast.dbl = b;
exp = (int)((unsigned)(cast.bits >> 52) & 0x7FFU);
frac = (cast.bits & ((1ULL << 52) - 1ULL)) | (1ULL << 52);
if (exp == 0x7FF) { /* +-inf, NaN */
|
| | |
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 |
*
* SPDX-License-Identifier: Unlicense
*/
#if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559)
int mp_set_double(mp_int *a, double b)
{
unsigned long long frac;
int exp, res;
union {
double dbl;
unsigned long long bits;
} cast;
cast.dbl = b;
exp = (int)((unsigned)(cast.bits >> 52) & 0x7FFU);
frac = (cast.bits & ((1ULL << 52) - 1ULL)) | (1ULL << 52);
if (exp == 0x7FF) { /* +-inf, NaN */
|
Changes to libtommath/bn_mp_sqrt.c.
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 .. 33 34 35 36 37 38 39 40 41 42 43 44 45 46 .. 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 .. 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 |
* additional optimizations in place.
*
* SPDX-License-Identifier: Unlicense
*/
#ifndef NO_FLOATING_POINT
#include <math.h>
#endif
/* this function is less generic than mp_n_root, simpler and faster */
int mp_sqrt(const mp_int *arg, mp_int *ret)
{
int res;
mp_int t1, t2;
int i, j, k;
#ifndef NO_FLOATING_POINT
volatile double d;
mp_digit dig;
#endif
/* must be positive */
if (arg->sign == MP_NEG) {
return MP_VAL;
................................................................................
}
/* easy out */
if (mp_iszero(arg) == MP_YES) {
mp_zero(ret);
return MP_OKAY;
}
i = (arg->used / 2) - 1;
j = 2 * i;
if ((res = mp_init_size(&t1, i+2)) != MP_OKAY) {
return res;
}
................................................................................
goto E2;
}
for (k = 0; k < i; ++k) {
t1.dp[k] = (mp_digit) 0;
}
#ifndef NO_FLOATING_POINT
/* Estimate the square root using the hardware floating point unit. */
d = 0.0;
for (k = arg->used-1; k >= j; --k) {
d = ldexp(d, DIGIT_BIT) + (double)(arg->dp[k]);
}
................................................................................
} else {
t1.used = i+1;
t1.dp[i] = ((mp_digit) d) - 1;
}
#else
/* Estimate the square root as having 1 in the most significant place. */
t1.used = i + 2;
t1.dp[i+1] = (mp_digit) 1;
t1.dp[i] = (mp_digit) 0;
#endif
/* t1 > 0 */
if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
goto E1;
}
|
> > > < > > > < < | > | < < > | > > > > > |
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 .. 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 .. 53 54 55 56 57 58 59 60 61 62 63 64 65 66 .. 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 |
* additional optimizations in place. * * SPDX-License-Identifier: Unlicense */ #ifndef NO_FLOATING_POINT #include <math.h> #if (DIGIT_BIT != 28) || (FLT_RADIX != 2) || (DBL_MANT_DIG != 53) || (DBL_MAX_EXP != 1024) #define NO_FLOATING_POINT #endif #endif /* this function is less generic than mp_n_root, simpler and faster */ int mp_sqrt(const mp_int *arg, mp_int *ret) { int res; mp_int t1, t2; #ifndef NO_FLOATING_POINT int i, j, k; volatile double d; mp_digit dig; #endif /* must be positive */ if (arg->sign == MP_NEG) { return MP_VAL; ................................................................................ } /* easy out */ if (mp_iszero(arg) == MP_YES) { mp_zero(ret); return MP_OKAY; } #ifndef NO_FLOATING_POINT i = (arg->used / 2) - 1; j = 2 * i; if ((res = mp_init_size(&t1, i+2)) != MP_OKAY) { return res; } ................................................................................ goto E2; } for (k = 0; k < i; ++k) { t1.dp[k] = (mp_digit) 0; } /* Estimate the square root using the hardware floating point unit. */ d = 0.0; for (k = arg->used-1; k >= j; --k) { d = ldexp(d, DIGIT_BIT) + (double)(arg->dp[k]); } ................................................................................ } else { t1.used = i+1; t1.dp[i] = ((mp_digit) d) - 1; } #else if ((res = mp_init_copy(&t1, arg)) != MP_OKAY) { return res; } if ((res = mp_init(&t2)) != MP_OKAY) { goto E2; } /* First approx. (not very bad for large arg) */ mp_rshd(&t1, t1.used/2); #endif /* t1 > 0 */ if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) { goto E1; } |