Révision 218
| pobysoPythonSage/src/pobyso.py (revision 218) | ||
|---|---|---|
| 90 | 90 |
return(sollya_lib_absolute(None)) |
| 91 | 91 |
|
| 92 | 92 |
def pobyso_autoprint(arg): |
| 93 |
sollya_lib_autoprint(arg,None) |
|
| 93 |
sollya_lib_autoprint(arg, None)
|
|
| 94 | 94 |
|
| 95 | 95 |
def pobyso_autoprint_so_so(arg): |
| 96 | 96 |
sollya_lib_autoprint(arg,None) |
| ... | ... | |
| 369 | 369 |
## Precision. If a precision is given, convert the polynomial |
| 370 | 370 |
# into the right polynomial field. If not convert it straight |
| 371 | 371 |
# to Sollya. |
| 372 |
if not precSa is None: |
|
| 373 |
RRR = RealField(precSa) |
|
| 374 |
## Create a Sage polynomial in the "right" precision. |
|
| 375 |
P_RRR = RRR[polySa.variables()[0]] |
|
| 376 |
polyFloatSa = P_RRR(polySa) |
|
| 377 |
else: |
|
| 378 |
polyFloatSa = polySa |
|
| 372 |
sollyaPrecChanged = False |
|
| 373 |
(curSollyaPrecSo, curSollyaPrecSa) = pobyso_get_prec_so_so_sa() |
|
| 374 |
if precSa is None: |
|
| 379 | 375 |
precSa = polySa.parent().base_ring().precision() |
| 376 |
if (precSa != curSollyaPrecSa): |
|
| 377 |
precSo = pobyso_constant_from_int(precSa) |
|
| 378 |
pobyso_set_prec_so_so(precSo) |
|
| 379 |
sollya_lib_clear_obj(precSo) |
|
| 380 |
sollyaPrecChanged = True |
|
| 380 | 381 |
## Get exponents and coefficients. |
| 381 |
exponentsSa = polyFloatSa.exponents()
|
|
| 382 |
coefficientsSa = polyFloatSa.coefficients()
|
|
| 382 |
exponentsSa = polySa.exponents()
|
|
| 383 |
coefficientsSa = polySa.coefficients()
|
|
| 383 | 384 |
## Build the polynomial. |
| 384 | 385 |
polySo = None |
| 385 | 386 |
for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa): |
| ... | ... | |
| 393 | 394 |
monomialSo = sollya_lib_build_function_pow( |
| 394 | 395 |
sollya_lib_build_function_free_variable(), |
| 395 | 396 |
exponentSo) |
| 397 |
polyTermSo = sollya_lib_build_function_mul(coefficientSo, |
|
| 398 |
monomialSo) |
|
| 396 | 399 |
if polySo is None: |
| 397 |
polySo = sollya_lib_build_function_mul(coefficientSo, |
|
| 398 |
monomialSo) |
|
| 400 |
polySo = polyTermSo |
|
| 399 | 401 |
else: |
| 400 |
polyTermSo = sollya_lib_build_function_mul(coefficientSo, |
|
| 401 |
monomialSo) |
|
| 402 | 402 |
polySo = sollya_lib_build_function_add(polySo, polyTermSo) |
| 403 |
if sollyaPrecChanged: |
|
| 404 |
pobyso_set_prec_so_so(curSollyaPrecSo) |
|
| 405 |
sollya_lib_clear_obj(curSollyaPrecSo) |
|
| 403 | 406 |
return polySo |
| 404 | 407 |
# End pobyso_float_poly_sa_so |
| 405 | 408 |
|
| 406 | 409 |
def pobyso_float_poly_so_sa(polySo, realFieldSa=None): |
| 407 | 410 |
""" |
| 408 | 411 |
Convert a Sollya polynomial into a Sage floating-point polynomial. |
| 409 |
We assume that the polynomial is in canonical form. |
|
| 410 | 412 |
If no realField is given, a RealField corresponding to the maximum |
| 411 | 413 |
precision of the coefficients is internally computed. |
| 412 | 414 |
The real field is not returned but can be easily retrieved from |
| ... | ... | |
| 415 | 417 |
- (optional) compute the RealField of the coefficients; |
| 416 | 418 |
- convert the Sollya expression into a Sage expression; |
| 417 | 419 |
- convert the Sage expression into a Sage polynomial |
| 418 |
TODO: the canonical thing for the polynomial. |
|
| 419 | 420 |
""" |
| 420 | 421 |
if realFieldSa is None: |
| 421 | 422 |
expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo) |
| 422 |
print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa |
|
| 423 |
#print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
|
|
| 423 | 424 |
realFieldSa = RealField(expressionPrecSa) |
| 424 | 425 |
#print "Sollya expression before...", |
| 425 | 426 |
#pobyso_autoprint(polySo) |
| 426 | 427 |
|
| 427 | 428 |
expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, |
| 428 | 429 |
realFieldSa) |
| 429 |
#print "...Sollya expression after.",
|
|
| 430 |
#print "...Sollya expression after." |
|
| 430 | 431 |
#pobyso_autoprint(polySo) |
| 431 | 432 |
polyVariableSa = expressionSa.variables()[0] |
| 432 | 433 |
polyRingSa = realFieldSa[str(polyVariableSa)] |
| ... | ... | |
| 575 | 576 |
The precision of the floating-point number returned is that of the Sollya |
| 576 | 577 |
constant. |
| 577 | 578 |
""" |
| 579 |
#print "Before computing precision of variable..." |
|
| 580 |
#pobyso_autoprint(constExpSo) |
|
| 578 | 581 |
precisionSa = pobyso_get_prec_of_constant_so_sa(constExpSo) |
| 582 |
#print "precisionSa:", precisionSa |
|
| 579 | 583 |
## If the expression can not be exactly converted, None is returned. |
| 580 | 584 |
# In this case opt for the Sollya current expression. |
| 581 | 585 |
if precisionSa is None: |
| ... | ... | |
| 641 | 645 |
|
| 642 | 646 |
def pobyso_get_head_function_so_sa(expressionSo): |
| 643 | 647 |
functionType = c_int(0) |
| 644 |
sollya_lib_get_head_function(byref(functionType), expressionSo, None)
|
|
| 648 |
sollya_lib_get_head_function(byref(functionType), expressionSo) |
|
| 645 | 649 |
return int(functionType.value) |
| 646 | 650 |
|
| 647 | 651 |
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ): |
| ... | ... | |
| 813 | 817 |
Tries to find a precision to represent ctExpSo without rounding. |
| 814 | 818 |
If not possible, returns None. |
| 815 | 819 |
""" |
| 820 |
#print "Entering pobyso_get_prec_of_constant_so_sa..." |
|
| 816 | 821 |
prec = c_int(0) |
| 817 | 822 |
retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None) |
| 818 | 823 |
if retc == 0: |
| 824 |
#print "pobyso_get_prec_of_constant_so_sa failed." |
|
| 819 | 825 |
return None |
| 826 |
#print "...exiting pobyso_get_prec_of_constant_so_sa." |
|
| 820 | 827 |
return int(prec.value) |
| 821 | 828 |
|
| 822 | 829 |
def pobyso_get_prec_of_range_so_sa(rangeSo): |
| ... | ... | |
| 876 | 883 |
#print "This is a constant" |
| 877 | 884 |
return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa) |
| 878 | 885 |
elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE: |
| 879 |
#print "This is free variable" |
|
| 886 |
#print "This is the free variable"
|
|
| 880 | 887 |
return eval(sollyaLibFreeVariableName) |
| 881 | 888 |
else: |
| 882 | 889 |
print "Unexpected" |
| ... | ... | |
| 1123 | 1130 |
maxPrecSo, |
| 1124 | 1131 |
maxErrSo): |
| 1125 | 1132 |
print "Input arguments:" |
| 1126 |
pobyso_autoprint(polySo) |
|
| 1127 |
pobyso_autoprint(funcSo) |
|
| 1128 |
pobyso_autoprint(icSo) |
|
| 1129 |
pobyso_autoprint(intervalSo) |
|
| 1130 |
pobyso_autoprint(itpSo) |
|
| 1131 |
pobyso_autoprint(ftpSo) |
|
| 1132 |
pobyso_autoprint(maxPrecSo) |
|
| 1133 |
pobyso_autoprint(maxErrSo) |
|
| 1134 |
print "________________" |
|
| 1133 |
#pobyso_autoprint(polySo)
|
|
| 1134 |
#pobyso_autoprint(funcSo)
|
|
| 1135 |
#pobyso_autoprint(icSo)
|
|
| 1136 |
#pobyso_autoprint(intervalSo)
|
|
| 1137 |
#pobyso_autoprint(itpSo)
|
|
| 1138 |
#pobyso_autoprint(ftpSo)
|
|
| 1139 |
#pobyso_autoprint(maxPrecSo)
|
|
| 1140 |
#pobyso_autoprint(maxErrSo)
|
|
| 1141 |
#print "________________"
|
|
| 1135 | 1142 |
|
| 1136 | 1143 |
## Higher order function see: |
| 1137 | 1144 |
# http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm |
| ... | ... | |
| 1148 | 1155 |
|
| 1149 | 1156 |
# |
| 1150 | 1157 |
degreeSa = pobyso_polynomial_degree_so_sa(polySo) |
| 1158 |
print "degreeSa:", degreeSa |
|
| 1151 | 1159 |
ratio = precision_decay_ratio_function(degreeSa) |
| 1160 |
print "ratio:", ratio |
|
| 1152 | 1161 |
itpSa = pobyso_constant_from_int_so_sa(itpSo) |
| 1162 |
print "itpsSa:", itpSa |
|
| 1153 | 1163 |
ftpSa = pobyso_constant_from_int_so_sa(ftpSo) |
| 1164 |
print "ftpSa:", ftpSa |
|
| 1154 | 1165 |
maxPrecSa = pobyso_constant_from_int_so_sa(maxPrecSo) |
| 1166 |
print "maxPrecSa:", maxPrecSa |
|
| 1155 | 1167 |
maxErrSa = pobyso_get_constant_as_rn_so_sa(maxErrSo) |
| 1156 |
resPolySo = pobyso_constant_0_sa_so()
|
|
| 1168 |
print "maxErrSa:", maxErrSa
|
|
| 1157 | 1169 |
lastResPolySo = None |
| 1170 |
lastInfNormSo = None |
|
| 1171 |
print "About to enter the while loop..." |
|
| 1158 | 1172 |
while True: |
| 1173 |
resPolySo = pobyso_constant_0_sa_so() |
|
| 1159 | 1174 |
pDeltaSa = ftpSa - itpSa |
| 1160 | 1175 |
for indexSa in reversed(xrange(0,degreeSa+1)): |
| 1161 |
print "Index:", indexSa |
|
| 1176 |
#print "Index:", indexSa
|
|
| 1162 | 1177 |
indexSo = pobyso_constant_from_int_sa_so(indexSa) |
| 1163 | 1178 |
#pobyso_autoprint(indexSo) |
| 1164 | 1179 |
#print ratio(indexSa) |
| ... | ... | |
| 1177 | 1192 |
#pobyso_autoprint(cmonrSo) |
| 1178 | 1193 |
resPolySo = sollya_lib_build_function_add(resPolySo, |
| 1179 | 1194 |
cmonrSo) |
| 1180 |
pobyso_autoprint(resPolySo) |
|
| 1195 |
#pobyso_autoprint(resPolySo)
|
|
| 1181 | 1196 |
# End for index |
| 1182 | 1197 |
freeVarSo = sollya_lib_build_function_free_variable() |
| 1183 | 1198 |
changeVarSo = sollya_lib_sub(freeVarSo, icSo) |
| 1184 | 1199 |
resPolyCvSo = sollya_lib_evaluate(resPolySo, changeVarSo) |
| 1185 |
errFuncSo = sollya_lib_build_function_sub(funcSo, resPolyCvSo)
|
|
| 1186 |
infNormSo = sollya_lib_dirtyinform(errFuncSo, intervalSo)
|
|
| 1187 |
print "Infnorm (Sollya):", ; pobyso_autoprint(infNormSo)
|
|
| 1200 |
errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
|
|
| 1201 |
resPolyCvSo)
|
|
| 1202 |
infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
|
|
| 1188 | 1203 |
cerrSa = pobyso_get_constant_as_rn_so_sa(infNormSo) |
| 1189 |
sollya_lib_clear_obj(resPolyCvSo)
|
|
| 1190 |
sollya_lib_clear_obj(infNormSo)
|
|
| 1191 |
print "Infnorm (Sage):", cerrSa |
|
| 1204 |
print "Infnorm (Sollya):", pobyso_autoprint(infNormSo)
|
|
| 1205 |
sollya_lib_clear_obj(errFuncSo)
|
|
| 1206 |
#print "Infnorm (Sage):", cerrSa
|
|
| 1192 | 1207 |
if (cerrSa > maxErrSa): |
| 1193 | 1208 |
print "Error is too large." |
| 1194 |
if not lastResPolySo is None:
|
|
| 1209 |
if lastResPolySo is None: |
|
| 1195 | 1210 |
print "Enlarging prec." |
| 1196 | 1211 |
ntpSa = floor(ftpSa + ftpSa/50) |
| 1197 | 1212 |
## Can't enlarge (numerical) |
| 1198 | 1213 |
if ntpSa == ftpSa: |
| 1199 |
sollya_lib_clear_obj(lastResPolySo) |
|
| 1200 | 1214 |
sollya_lib_clear_obj(resPolySo) |
| 1201 | 1215 |
return None |
| 1202 | 1216 |
## Can't enlarge (not enough precision left) |
| 1203 | 1217 |
if ntpSa > maxPrecSa: |
| 1204 |
sollya_lib_clear_obj(lastResPolySo) |
|
| 1205 | 1218 |
sollya_lib_clear_obj(resPolySo) |
| 1206 | 1219 |
return None |
| 1207 | 1220 |
ftpSa = ntpSa |
| 1208 |
lastResPolySo = resPolySo |
|
| 1209 | 1221 |
continue |
| 1210 | 1222 |
## One enlargement took place. |
| 1211 | 1223 |
else: |
| 1224 |
print "Exit with the last before last polynomial." |
|
| 1212 | 1225 |
sollya_lib_clear_obj(resPolySo) |
| 1213 |
return lastResPolySo |
|
| 1214 |
# cerrSa <= maxErrSa: scrap more bits. |
|
| 1226 |
sollya_lib_clear_obj(infNormSo) |
|
| 1227 |
return (lastResPolySo, lastInfNormSo) |
|
| 1228 |
# cerrSa <= maxErrSa: scrap more bits, possibly. |
|
| 1215 | 1229 |
else: |
| 1216 |
print "Shrinking prec." |
|
| 1217 |
ntpSa = floor(ftpSa - ftpSa/50) |
|
| 1218 |
## Can't shrink (numerical) |
|
| 1219 |
if ntpSa == ftpSa: |
|
| 1230 |
print "Error is too small" |
|
| 1231 |
if cerrSa <= (maxErrSa/2): |
|
| 1232 |
print "Shrinking prec." |
|
| 1233 |
ntpSa = floor(ftpSa - ftpSa/50) |
|
| 1234 |
## Can't shrink (numerical) |
|
| 1235 |
if ntpSa == ftpSa: |
|
| 1236 |
if not lastResPolySo is None: |
|
| 1237 |
sollya_lib_clear_obj(lastResPolySo) |
|
| 1238 |
if not lastInfNormSo is None: |
|
| 1239 |
sollya_lib_clear_obj(lastInfNormSo) |
|
| 1240 |
return (resPolySo, infNormSo) |
|
| 1241 |
## Can't shrink (not enough precision left) |
|
| 1242 |
if ntpSa <= itpSa: |
|
| 1243 |
if not lastResPolySo is None: |
|
| 1244 |
sollya_lib_clear_obj(lastResPolySo) |
|
| 1245 |
if not lastInfNormSo is None: |
|
| 1246 |
sollya_lib_clear_obj(lastInfNormSo) |
|
| 1247 |
return (resPolySo, infNormSo) |
|
| 1248 |
ftpSa = ntpSa |
|
| 1220 | 1249 |
if not lastResPolySo is None: |
| 1221 | 1250 |
sollya_lib_clear_obj(lastResPolySo) |
| 1222 |
return resPolySo |
|
| 1223 |
## Can't enlarge (not enough precision left) |
|
| 1224 |
if ntpSa <= itpSa: |
|
| 1251 |
if not lastInfNormSo is None: |
|
| 1252 |
sollya_lib_clear_obj(lastInfNormSo) |
|
| 1253 |
lastResPolySo = resPolySo |
|
| 1254 |
lastInfNormSo = infNormSo |
|
| 1255 |
continue |
|
| 1256 |
else: # Error is not that small, just return |
|
| 1225 | 1257 |
if not lastResPolySo is None: |
| 1226 | 1258 |
sollya_lib_clear_obj(lastResPolySo) |
| 1227 |
return resPolySo
|
|
| 1228 |
ftpSa = ntpSa
|
|
| 1229 |
continue
|
|
| 1259 |
if not lastInfNormSo is None:
|
|
| 1260 |
sollya_lib_clear_obj(lastInfNormSo)
|
|
| 1261 |
return (resPolySo, infNormSo)
|
|
| 1230 | 1262 |
# End wile True |
| 1231 | 1263 |
return None |
| 1232 | 1264 |
# End pobyso_polynomial_coefficients_progressive_truncate_so_so. |
| pobysoPythonSage/src/sollya_lib.sage (revision 218) | ||
|---|---|---|
| 42 | 42 |
sollya_lib_copy_obj = sollya.sollya_lib_copy_obj |
| 43 | 43 |
sollya_lib_cos = sollya.sollya_lib_cos |
| 44 | 44 |
sollya_lib_degree = sollya.sollya_lib_degree |
| 45 |
sollya_lib_dirtyinform = sollya.sollya_lib_dirtyinfnorm |
|
| 45 |
sollya_lib_dirtyinfnorm = sollya.sollya_lib_dirtyinfnorm
|
|
| 46 | 46 |
sollya_lib_end_elliptic_list= \ |
| 47 | 47 |
sollya.sollya_lib_end_elliptic_list |
| 48 | 48 |
sollya_lib_error=sollya.sollya_lib_error |
| ... | ... | |
| 129 | 129 |
try: |
| 130 | 130 |
sollya_lib_add.argtypes = [c_ulong, c_ulong] |
| 131 | 131 |
sollya_lib_append.argstypes = [c_ulong, c_ulong] |
| 132 |
sollya_lib_autoprint.argtypes = [c_ulong] |
|
| 132 |
#sollya_lib_autoprint.argtypes = [c_ulong]
|
|
| 133 | 133 |
sollya_lib_build_function_add.argtypes = [c_int, c_int] |
| 134 | 134 |
sollya_lib_constant.argtypes = [c_ulong] |
| 135 | 135 |
sollya_lib_constant_from_int64.argtypes = [c_long] |
| 136 | 136 |
sollya_lib_constant_from_mpq.argtypes = [c_ulong] |
| 137 | 137 |
sollya_lib_constant.restype = c_ulong |
| 138 |
sollya_lib_dirtyinfnorm.argstypes = [c_ulong, c_ulong] |
|
| 138 | 139 |
sollya_lib_end_elliptic_list.argtypes = [c_ulong, c_int] |
| 139 | 140 |
sollya_lib_get_constant.argtypes = [c_ulong, c_ulong] |
| 140 | 141 |
sollya_lib_get_constant_as_int.argtypes = [POINTER(c_int), c_int] |
| 141 | 142 |
sollya_lib_get_constant_as_int64.argtypes= [POINTER(c_long), c_int] |
| 142 | 143 |
sollya_lib_get_function_arity.argtypes = [POINTER(c_int), c_int] |
| 143 |
sollya_lib_get_head_function.argtypes = [POINTER(c_int), c_int]
|
|
| 144 |
sollya_lib_get_head_function.argtypes = [POINTER(c_int), c_ulong]
|
|
| 144 | 145 |
sollya_lib_get_interval_from_range.argtypes = [c_ulong, c_ulong] |
| 145 | 146 |
sollya_lib_get_list_elements.argtypes = [POINTER(POINTER(c_longlong)), \ |
| 146 | 147 |
POINTER(c_int),\ |
| 147 | 148 |
POINTER(c_int), c_int] |
| 148 |
sollya_lib_get_prec_of_constant.argtypes = [POINTER(c_int), c_int]
|
|
| 149 |
sollya_lib_get_prec_of_constant.argtypes = [POINTER(c_int), c_ulong]
|
|
| 149 | 150 |
sollya_lib_get_prec_of_range.argtypes = [POINTER(c_int), c_int] |
| 150 | 151 |
sollya_lib_get_subfunctions.argtypes = [c_int, POINTER(c_int), \ |
| 151 | 152 |
POINTER(c_int), POINTER(c_int), POINTER(c_int), POINTER(c_int), \ |
| pobysoPythonSage/src/sageSLZ/sageSLZ.sage (revision 218) | ||
|---|---|---|
| 494 | 494 |
# End slz_compute_integer_polynomial_modular_roots_2. |
| 495 | 495 |
# |
| 496 | 496 |
def slz_compute_polynomial_and_interval(functionSo, degreeSo, lowerBoundSa, |
| 497 |
upperBoundSa, approxPrecSa,
|
|
| 498 |
sollyaPrecSa=None):
|
|
| 497 |
upperBoundSa, approxAccurSa,
|
|
| 498 |
precSa=None):
|
|
| 499 | 499 |
""" |
| 500 | 500 |
Under the assumptions listed for slz_get_intervals_and_polynomials, compute |
| 501 | 501 |
a polynomial that approximates the function on a an interval starting |
| ... | ... | |
| 511 | 511 |
- an range (an interval, not in the sense of number given as an interval); |
| 512 | 512 |
- the center of the interval; |
| 513 | 513 |
- the maximum error in the approximation of the input functionSo by the |
| 514 |
output polynomial ; this error <= approxPrecSaS.
|
|
| 514 |
output polynomial ; this error <= approxAccurSaS.
|
|
| 515 | 515 |
|
| 516 | 516 |
""" |
| 517 |
#print"In slz_compute_polynomial_and_interval..." |
|
| 517 | 518 |
## Superficial argument check. |
| 518 | 519 |
if lowerBoundSa > upperBoundSa: |
| 519 | 520 |
return None |
| 521 |
## Change Sollya precision, if requested. |
|
| 522 |
sollyaPrecChanged = False |
|
| 523 |
(sollyaPrecSo, sollyaPrecSa) = pobyso_get_prec_so_so_sa() |
|
| 524 |
if precSa is None: |
|
| 525 |
precSa = ceil((RR('1.5') * abs(RR(approxAccurSa).log2())) / 64) * 64
|
|
| 526 |
#print "Computed internal precision:", precSa |
|
| 527 |
if precSa < 192: |
|
| 528 |
precSa = 192 |
|
| 529 |
if precSa != sollyaPrecSa: |
|
| 530 |
precSo = pobyso_constant_from_int_sa_so(precSa) |
|
| 531 |
pobyso_set_prec_so_so(precSo) |
|
| 532 |
sollya_lib_clear_obj(precSo) |
|
| 533 |
sollyaPrecChanged = True |
|
| 520 | 534 |
RRR = lowerBoundSa.parent() |
| 521 | 535 |
intervalShrinkConstFactorSa = RRR('0.9')
|
| 522 | 536 |
absoluteErrorTypeSo = pobyso_absolute_so_so() |
| ... | ... | |
| 531 | 545 |
currentRangeSo, |
| 532 | 546 |
absoluteErrorTypeSo) |
| 533 | 547 |
maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo) |
| 534 |
while maxErrorSa > approxPrecSa:
|
|
| 548 |
while maxErrorSa > approxAccurSa:
|
|
| 535 | 549 |
print "++Approximation error:", maxErrorSa.n() |
| 536 | 550 |
sollya_lib_clear_obj(polySo) |
| 537 | 551 |
sollya_lib_clear_obj(intervalCenterSo) |
| 538 | 552 |
sollya_lib_clear_obj(maxErrorSo) |
| 539 | 553 |
# Very empirical shrinking factor. |
| 540 |
shrinkFactorSa = 1 / (maxErrorSa/approxPrecSa).log2().abs()
|
|
| 554 |
shrinkFactorSa = 1 / (maxErrorSa/approxAccurSa).log2().abs()
|
|
| 541 | 555 |
print "Shrink factor:", \ |
| 542 | 556 |
shrinkFactorSa.n(), \ |
| 543 | 557 |
intervalShrinkConstFactorSa |
| 544 | 558 |
|
| 545 |
#errorRatioSa = approxPrecSa/maxErrorSa
|
|
| 559 |
#errorRatioSa = approxAccurSa/maxErrorSa
|
|
| 546 | 560 |
#print "Error ratio: ", errorRatioSa |
| 547 | 561 |
# Make sure interval shrinks. |
| 548 | 562 |
if shrinkFactorSa > intervalShrinkConstFactorSa: |
| ... | ... | |
| 566 | 580 |
print "Use either or both a higher polynomial degree or a higher", |
| 567 | 581 |
print "internal precision." |
| 568 | 582 |
print "Aborting!" |
| 569 |
return (None, None, None, None) |
|
| 583 |
if sollyaPrecChanged: |
|
| 584 |
pobyso_set_prec_so_so(sollyaPrecSo) |
|
| 585 |
sollya_lib_clear_obj(sollyaPrecSo) |
|
| 586 |
return None |
|
| 570 | 587 |
currentRangeSo = pobyso_bounds_to_range_sa_so(currentLowerBoundSa, |
| 571 | 588 |
currentUpperBoundSa) |
| 572 | 589 |
# print "New interval:", |
| ... | ... | |
| 583 | 600 |
#print "Max errorSa:", maxErrorSa |
| 584 | 601 |
#print "Sollya prec:", |
| 585 | 602 |
#pobyso_autoprint(sollya_lib_get_prec(None)) |
| 603 |
# End while |
|
| 586 | 604 |
sollya_lib_clear_obj(absoluteErrorTypeSo) |
| 605 |
if sollyaPrecChanged: |
|
| 606 |
pobyso_set_prec_so_so(sollyaPrecSo) |
|
| 607 |
sollya_lib_clear_obj(sollyaPrecSo) |
|
| 587 | 608 |
return (polySo, currentRangeSo, intervalCenterSo, maxErrorSo) |
| 588 | 609 |
# End slz_compute_polynomial_and_interval |
| 610 |
|
|
| 611 |
def slz_compute_polynomial_and_interval_01(functionSo, degreeSo, lowerBoundSa, |
|
| 612 |
upperBoundSa, approxAccurSa, |
|
| 613 |
sollyaPrecSa=None): |
|
| 614 |
""" |
|
| 615 |
Under the assumptions listed for slz_get_intervals_and_polynomials, compute |
|
| 616 |
a polynomial that approximates the function on a an interval starting |
|
| 617 |
at lowerBoundSa and finishing at a value that guarantees that the polynomial |
|
| 618 |
approximates with the expected precision. |
|
| 619 |
The interval upper bound is lowered until the expected approximation |
|
| 620 |
precision is reached. |
|
| 621 |
The polynomial, the bounds, the center of the interval and the error |
|
| 622 |
are returned. |
|
| 623 |
OUTPUT: |
|
| 624 |
A tuple made of 4 Sollya objects: |
|
| 625 |
- a polynomial; |
|
| 626 |
- an range (an interval, not in the sense of number given as an interval); |
|
| 627 |
- the center of the interval; |
|
| 628 |
- the maximum error in the approximation of the input functionSo by the |
|
| 629 |
output polynomial ; this error <= approxAccurSaS. |
|
| 630 |
|
|
| 631 |
""" |
|
| 632 |
print"In slz_compute_polynomial_and_interval..." |
|
| 633 |
## Superficial argument check. |
|
| 634 |
if lowerBoundSa > upperBoundSa: |
|
| 635 |
return None |
|
| 636 |
## Change Sollya precision, if requested. |
|
| 637 |
if not sollyaPrecSa is None: |
|
| 638 |
sollyaPrecSo = pobyso_get_prec_so() |
|
| 639 |
pobyso_set_prec_sa_so(sollyaPrecSa) |
|
| 640 |
else: |
|
| 641 |
sollyaPrecSa = pobyso_get_prec_so_sa() |
|
| 642 |
sollyaPrecSo = None |
|
| 643 |
RRR = lowerBoundSa.parent() |
|
| 644 |
intervalShrinkConstFactorSa = RRR('0.9')
|
|
| 645 |
absoluteErrorTypeSo = pobyso_absolute_so_so() |
|
| 646 |
currentRangeSo = pobyso_bounds_to_range_sa_so(lowerBoundSa, upperBoundSa) |
|
| 647 |
currentUpperBoundSa = upperBoundSa |
|
| 648 |
currentLowerBoundSa = lowerBoundSa |
|
| 649 |
# What we want here is the polynomial without the variable change, |
|
| 650 |
# since our actual variable will be x-intervalCenter defined over the |
|
| 651 |
# domain [lowerBound-intervalCenter , upperBound-intervalCenter]. |
|
| 652 |
(polySo, intervalCenterSo, maxErrorSo) = \ |
|
| 653 |
pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, |
|
| 654 |
currentRangeSo, |
|
| 655 |
absoluteErrorTypeSo) |
|
| 656 |
maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo) |
|
| 657 |
while maxErrorSa > approxAccurSa: |
|
| 658 |
print "++Approximation error:", maxErrorSa.n() |
|
| 659 |
sollya_lib_clear_obj(polySo) |
|
| 660 |
sollya_lib_clear_obj(intervalCenterSo) |
|
| 661 |
sollya_lib_clear_obj(maxErrorSo) |
|
| 662 |
# Very empirical shrinking factor. |
|
| 663 |
shrinkFactorSa = 1 / (maxErrorSa/approxAccurSa).log2().abs() |
|
| 664 |
print "Shrink factor:", \ |
|
| 665 |
shrinkFactorSa.n(), \ |
|
| 666 |
intervalShrinkConstFactorSa |
|
| 667 |
|
|
| 668 |
#errorRatioSa = approxAccurSa/maxErrorSa |
|
| 669 |
#print "Error ratio: ", errorRatioSa |
|
| 670 |
# Make sure interval shrinks. |
|
| 671 |
if shrinkFactorSa > intervalShrinkConstFactorSa: |
|
| 672 |
actualShrinkFactorSa = intervalShrinkConstFactorSa |
|
| 673 |
#print "Fixed" |
|
| 674 |
else: |
|
| 675 |
actualShrinkFactorSa = shrinkFactorSa |
|
| 676 |
#print "Computed",shrinkFactorSa,maxErrorSa |
|
| 677 |
#print shrinkFactorSa, maxErrorSa |
|
| 678 |
#print "Shrink factor", actualShrinkFactorSa |
|
| 679 |
currentUpperBoundSa = currentLowerBoundSa + \ |
|
| 680 |
(currentUpperBoundSa - currentLowerBoundSa) * \ |
|
| 681 |
actualShrinkFactorSa |
|
| 682 |
#print "Current upper bound:", currentUpperBoundSa |
|
| 683 |
sollya_lib_clear_obj(currentRangeSo) |
|
| 684 |
# Check what is left with the bounds. |
|
| 685 |
if currentUpperBoundSa <= currentLowerBoundSa or \ |
|
| 686 |
currentUpperBoundSa == currentLowerBoundSa.nextabove(): |
|
| 687 |
sollya_lib_clear_obj(absoluteErrorTypeSo) |
|
| 688 |
print "Can't find an interval." |
|
| 689 |
print "Use either or both a higher polynomial degree or a higher", |
|
| 690 |
print "internal precision." |
|
| 691 |
print "Aborting!" |
|
| 692 |
if not sollyaPrecSo is None: |
|
| 693 |
pobyso_set_prec_so_so(sollyaPrecSo) |
|
| 694 |
sollya_lib_clear_obj(sollyaPrecSo) |
|
| 695 |
return None |
|
| 696 |
currentRangeSo = pobyso_bounds_to_range_sa_so(currentLowerBoundSa, |
|
| 697 |
currentUpperBoundSa) |
|
| 698 |
# print "New interval:", |
|
| 699 |
# pobyso_autoprint(currentRangeSo) |
|
| 700 |
#print "Second Taylor expansion call." |
|
| 701 |
(polySo, intervalCenterSo, maxErrorSo) = \ |
|
| 702 |
pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, |
|
| 703 |
currentRangeSo, |
|
| 704 |
absoluteErrorTypeSo) |
|
| 705 |
#maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo, RRR) |
|
| 706 |
#print "Max errorSo:", |
|
| 707 |
#pobyso_autoprint(maxErrorSo) |
|
| 708 |
maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo) |
|
| 709 |
#print "Max errorSa:", maxErrorSa |
|
| 710 |
#print "Sollya prec:", |
|
| 711 |
#pobyso_autoprint(sollya_lib_get_prec(None)) |
|
| 712 |
# End while |
|
| 713 |
sollya_lib_clear_obj(absoluteErrorTypeSo) |
|
| 714 |
itpSo = pobyso_constant_from_int_sa_so(floor(sollyaPrecSa/3)) |
|
| 715 |
ftpSo = pobyso_constant_from_int_sa_so(floor(2*sollyaPrecSa/3)) |
|
| 716 |
maxPrecSo = pobyso_constant_from_int_sa_so(sollyaPrecSa) |
|
| 717 |
approxAccurSo = pobyso_constant_sa_so(RR(approxAccurSa)) |
|
| 718 |
print "About to call polynomial rounding with:" |
|
| 719 |
print "polySo: ", ; pobyso_autoprint(polySo) |
|
| 720 |
print "functionSo: ", ; pobyso_autoprint(functionSo) |
|
| 721 |
print "intervalCenterSo: ", ; pobyso_autoprint(intervalCenterSo) |
|
| 722 |
print "currentRangeSo: ", ; pobyso_autoprint(currentRangeSo) |
|
| 723 |
print "itpSo: ", ; pobyso_autoprint(itpSo) |
|
| 724 |
print "ftpSo: ", ; pobyso_autoprint(ftpSo) |
|
| 725 |
print "maxPrecSo: ", ; pobyso_autoprint(maxPrecSo) |
|
| 726 |
print "approxAccurSo: ", ; pobyso_autoprint(approxAccurSo) |
|
| 727 |
(roundedPolySo, roundedPolyMaxErrSo) = \ |
|
| 728 |
pobyso_polynomial_coefficients_progressive_round_so_so(polySo, |
|
| 729 |
functionSo, |
|
| 730 |
intervalCenterSo, |
|
| 731 |
currentRangeSo, |
|
| 732 |
itpSo, |
|
| 733 |
ftpSo, |
|
| 734 |
maxPrecSo, |
|
| 735 |
approxAccurSo) |
|
| 736 |
|
|
| 737 |
sollya_lib_clear_obj(polySo) |
|
| 738 |
sollya_lib_clear_obj(maxErrorSo) |
|
| 739 |
sollya_lib_clear_obj(itpSo) |
|
| 740 |
sollya_lib_clear_obj(ftpSo) |
|
| 741 |
sollya_lib_clear_obj(approxAccurSo) |
|
| 742 |
if not sollyaPrecSo is None: |
|
| 743 |
pobyso_set_prec_so_so(sollyaPrecSo) |
|
| 744 |
sollya_lib_clear_obj(sollyaPrecSo) |
|
| 745 |
print "1: ", ; pobyso_autoprint(roundedPolySo) |
|
| 746 |
print "2: ", ; pobyso_autoprint(currentRangeSo) |
|
| 747 |
print "3: ", ; pobyso_autoprint(intervalCenterSo) |
|
| 748 |
print "4: ", ; pobyso_autoprint(roundedPolyMaxErrSo) |
|
| 749 |
return (roundedPolySo, currentRangeSo, intervalCenterSo, roundedPolyMaxErrSo) |
|
| 750 |
# End slz_compute_polynomial_and_interval_01 |
|
| 751 |
|
|
| 589 | 752 |
def slz_compute_reduced_polynomial(matrixRow, |
| 590 | 753 |
knownMonomials, |
| 591 | 754 |
var1, |
| ... | ... | |
| 996 | 1159 |
scaledLowerBoundSa, scaledUpperBoundSa = \ |
| 997 | 1160 |
scaledUpperBoundSa, scaledLowerBoundSa |
| 998 | 1161 |
#print "Switching!" |
| 999 |
print "Approximation precision: ", RR(approxPrecSa)
|
|
| 1162 |
print "Approximation accuracy: ", RR(approxAccurSa)
|
|
| 1000 | 1163 |
# Prepare the arguments for the Taylor expansion computation with Sollya. |
| 1001 | 1164 |
functionSo = \ |
| 1002 | 1165 |
pobyso_parse_string_sa_so(fff._assume_str().replace('_SAGE_VAR_', ''))
|
| ... | ... | |
| 1015 | 1178 |
slz_compute_polynomial_and_interval(functionSo, degreeSo, |
| 1016 | 1179 |
currentScaledLowerBoundSa, |
| 1017 | 1180 |
currentScaledUpperBoundSa, |
| 1018 |
approxPrecSa, internalSollyaPrecSa)
|
|
| 1181 |
approxAccurSa, internalSollyaPrecSa)
|
|
| 1019 | 1182 |
print "...done." |
| 1020 | 1183 |
## If slz_compute_polynomial_and_interval fails, it returns None. |
| 1021 | 1184 |
# This value goes to the first variable: polySo. Other variables are |
| ... | ... | |
| 1046 | 1209 |
print "Computed approximation error:", errorSa.n(digits=10) |
| 1047 | 1210 |
# If the error and interval are OK a the first try, just return. |
| 1048 | 1211 |
if (boundsSa.endpoints()[1] >= scaledUpperBoundSa) and \ |
| 1049 |
(errorSa <= approxPrecSa):
|
|
| 1212 |
(errorSa <= approxAccurSa):
|
|
| 1050 | 1213 |
if precChangedSa: |
| 1051 | 1214 |
pobyso_set_prec_sa_so(currentSollyaPrecSa) |
| 1052 | 1215 |
sollya_lib_clear_obj(currentSollyaPrecSo) |
| ... | ... | |
| 1059 | 1222 |
# upper bound: compute the next upper bound. |
| 1060 | 1223 |
## The following ratio is always >= 1. If errorSa, we may want to |
| 1061 | 1224 |
# enlarge the interval |
| 1062 |
currentErrorRatio = approxPrecSa / errorSa
|
|
| 1225 |
currentErrorRatio = approxAccurSa / errorSa
|
|
| 1063 | 1226 |
## --|--------------------------------------------------------------|-- |
| 1064 | 1227 |
# --|--------------------|-------------------------------------------- |
| 1065 | 1228 |
# --|----------------------------|------------------------------------ |
| ... | ... | |
| 1223 | 1386 |
- [3]: the interval center; |
| 1224 | 1387 |
- [4]: the approximation error. |
| 1225 | 1388 |
|
| 1226 |
The function return a 5 elements tuple: formed with all the |
|
| 1389 |
The function returns a 5 elements tuple: formed with all the
|
|
| 1227 | 1390 |
input elements converted into their Sollya counterpart. |
| 1228 | 1391 |
""" |
| 1229 |
polynomialSa = pobyso_get_poly_so_sa(polyRangeCenterErrorSo[0]) |
|
| 1230 |
polynomialChangedVarSa = pobyso_get_poly_so_sa(polyRangeCenterErrorSo[1]) |
|
| 1392 |
polynomialSa = pobyso_float_poly_so_sa(polyRangeCenterErrorSo[0]) |
|
| 1393 |
#print "Polynomial after first conversion: ", pobyso_autoprint(polyRangeCenterErrorSo[1]) |
|
| 1394 |
polynomialChangedVarSa = pobyso_float_poly_so_sa(polyRangeCenterErrorSo[1]) |
|
| 1231 | 1395 |
intervalSa = \ |
| 1232 | 1396 |
pobyso_get_interval_from_range_so_sa(polyRangeCenterErrorSo[2]) |
| 1233 | 1397 |
centerSa = \ |
Formats disponibles : Unified diff