Révision 176 pobysoPythonSage/src/sageSLZ/sageSLZ.sage
| sageSLZ.sage (revision 176) | ||
|---|---|---|
| 112 | 112 |
if numberAsRR == 0: |
| 113 | 113 |
return RR(-infinity) |
| 114 | 114 |
else: |
| 115 |
return floor(numberAsRR.abs().log2()) |
|
| 115 |
realField = numberAsRR.parent() |
|
| 116 |
numberLog2 = numberAsRR.abs().log2() |
|
| 117 |
floorNumberLog2 = floor(numberLog2) |
|
| 118 |
## Do not get caught by rounding of log2() both ways. |
|
| 119 |
## When numberLog2 is an integer, compare numberAsRR |
|
| 120 |
# with 2^numberLog2. |
|
| 121 |
if floorNumberLog2 == numberLog2: |
|
| 122 |
if numberAsRR.abs() < realField(2^floorNumberLog2): |
|
| 123 |
return floorNumberLog2 - 1 |
|
| 124 |
else: |
|
| 125 |
return floorNumberLog2 |
|
| 126 |
else: |
|
| 127 |
return floorNumberLog2 |
|
| 116 | 128 |
# End slz_compute_binade |
| 117 | 129 |
|
| 118 | 130 |
# |
| ... | ... | |
| 480 | 492 |
if lowerBoundSa > upperBoundSa: |
| 481 | 493 |
return None |
| 482 | 494 |
RRR = lowerBoundSa.parent() |
| 483 |
intervalShrinkConstFactorSa = RRR('0.8')
|
|
| 495 |
intervalShrinkConstFactorSa = RRR('0.9')
|
|
| 484 | 496 |
absoluteErrorTypeSo = pobyso_absolute_so_so() |
| 485 | 497 |
currentRangeSo = pobyso_bounds_to_range_sa_so(lowerBoundSa, upperBoundSa) |
| 486 | 498 |
currentUpperBoundSa = upperBoundSa |
| ... | ... | |
| 498 | 510 |
sollya_lib_clear_obj(polySo) |
| 499 | 511 |
sollya_lib_clear_obj(intervalCenterSo) |
| 500 | 512 |
sollya_lib_clear_obj(maxErrorSo) |
| 501 |
shrinkFactorSa = RRR('5')/(maxErrorSa/approxPrecSa).log2().abs()
|
|
| 502 |
#shrinkFactorSa = 1.5/(maxErrorSa/approxPrecSa) |
|
| 513 |
# Very empirical schrinking factor. |
|
| 514 |
shrinkFactorSa = 1 / (maxErrorSa/approxPrecSa).log2().abs() |
|
| 515 |
print "Shrink factor:", shrinkFactorSa, intervalShrinkConstFactorSa |
|
| 503 | 516 |
#errorRatioSa = approxPrecSa/maxErrorSa |
| 504 | 517 |
#print "Error ratio: ", errorRatioSa |
| 505 | 518 |
if shrinkFactorSa > intervalShrinkConstFactorSa: |
| ... | ... | |
| 540 | 553 |
#print "Sollya prec:", |
| 541 | 554 |
#pobyso_autoprint(sollya_lib_get_prec(None)) |
| 542 | 555 |
sollya_lib_clear_obj(absoluteErrorTypeSo) |
| 543 |
return((polySo, currentRangeSo, intervalCenterSo, maxErrorSo))
|
|
| 556 |
return (polySo, currentRangeSo, intervalCenterSo, maxErrorSo)
|
|
| 544 | 557 |
# End slz_compute_polynomial_and_interval |
| 545 | 558 |
|
| 546 | 559 |
def slz_compute_reduced_polynomial(matrixRow, |
| ... | ... | |
| 795 | 808 |
degreeSo = pobyso_constant_from_int_sa_so(degreeSa) |
| 796 | 809 |
scaledBoundsSo = pobyso_bounds_to_range_sa_so(scaledLowerBoundSa, |
| 797 | 810 |
scaledUpperBoundSa) |
| 798 |
# Compute the first Taylor expansion. |
|
| 799 |
print "Computing first Taylor expansion..." |
|
| 800 |
(polySo, boundsSo, intervalCenterSo, maxErrorSo) = \ |
|
| 801 |
slz_compute_polynomial_and_interval(functionSo, degreeSo, |
|
| 802 |
scaledLowerBoundSa, scaledUpperBoundSa, |
|
| 803 |
approxPrecSa, internalSollyaPrecSa) |
|
| 804 |
print "...done." |
|
| 805 |
if polySo is None: |
|
| 806 |
print "slz_get_intervals_and_polynomials: Aborting and returning None!" |
|
| 807 |
if precChangedSa: |
|
| 808 |
pobyso_set_prec_so_so(currentSollyaPrecSo) |
|
| 809 |
sollya_lib_clear_obj(currentSollyaPrecSo) |
|
| 810 |
sollya_lib_clear_obj(functionSo) |
|
| 811 |
sollya_lib_clear_obj(degreeSo) |
|
| 812 |
sollya_lib_clear_obj(scaledBoundsSo) |
|
| 813 |
return None |
|
| 811 |
|
|
| 814 | 812 |
realIntervalField = RealIntervalField(max(lowerBoundSa.parent().precision(), |
| 815 | 813 |
upperBoundSa.parent().precision())) |
| 816 |
boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField) |
|
| 817 |
errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo) |
|
| 818 |
print "First approximation error:", errorSa.n(digits=10) |
|
| 819 |
# If the error and interval are OK a the first try, just return. |
|
| 820 |
if boundsSa.endpoints()[1] >= scaledUpperBoundSa: |
|
| 821 |
# Change variable stuff in Sollya x -> x0-x. |
|
| 822 |
changeVarExpressionSo = sollya_lib_build_function_sub( \ |
|
| 823 |
sollya_lib_build_function_free_variable(), \ |
|
| 814 |
currentScaledLowerBoundSa = scaledLowerBoundSa |
|
| 815 |
currentScaledUpperBoundSa = scaledUpperBoundSa |
|
| 816 |
while True: |
|
| 817 |
## Compute the first Taylor expansion. |
|
| 818 |
print "Computing a Taylor expansion..." |
|
| 819 |
(polySo, boundsSo, intervalCenterSo, maxErrorSo) = \ |
|
| 820 |
slz_compute_polynomial_and_interval(functionSo, degreeSo, |
|
| 821 |
currentScaledLowerBoundSa, |
|
| 822 |
currentScaledUpperBoundSa, |
|
| 823 |
approxPrecSa, internalSollyaPrecSa) |
|
| 824 |
print "...done." |
|
| 825 |
## If slz_compute_polynomial_and_interval fails, it returns None. |
|
| 826 |
# This value goes to the first variable: polySo. Other variables are |
|
| 827 |
# not assigned and should not be tested. |
|
| 828 |
if polySo is None: |
|
| 829 |
print "slz_get_intervals_and_polynomials: Aborting and returning None!" |
|
| 830 |
if precChangedSa: |
|
| 831 |
pobyso_set_prec_so_so(currentSollyaPrecSo) |
|
| 832 |
sollya_lib_clear_obj(currentSollyaPrecSo) |
|
| 833 |
sollya_lib_clear_obj(functionSo) |
|
| 834 |
sollya_lib_clear_obj(degreeSo) |
|
| 835 |
sollya_lib_clear_obj(scaledBoundsSo) |
|
| 836 |
return None |
|
| 837 |
## Add to the result array. |
|
| 838 |
### Change variable stuff in Sollya x -> x0-x. |
|
| 839 |
changeVarExpressionSo = \ |
|
| 840 |
sollya_lib_build_function_sub( \ |
|
| 841 |
sollya_lib_build_function_free_variable(), |
|
| 824 | 842 |
sollya_lib_copy_obj(intervalCenterSo)) |
| 825 |
polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo) |
|
| 843 |
polyVarChangedSo = \ |
|
| 844 |
sollya_lib_evaluate(polySo, changeVarExpressionSo) |
|
| 845 |
#### Get rid of the variable change Sollya stuff. |
|
| 826 | 846 |
sollya_lib_clear_obj(changeVarExpressionSo) |
| 827 |
resultArray.append((polySo, polyVarChangedSo, boundsSo, \
|
|
| 847 |
resultArray.append((polySo, polyVarChangedSo, boundsSo, |
|
| 828 | 848 |
intervalCenterSo, maxErrorSo)) |
| 829 |
if internalSollyaPrecSa != currentSollyaPrecSa: |
|
| 830 |
pobyso_set_prec_sa_so(currentSollyaPrecSa) |
|
| 831 |
sollya_lib_clear_obj(currentSollyaPrecSo) |
|
| 832 |
sollya_lib_clear_obj(functionSo) |
|
| 833 |
sollya_lib_clear_obj(degreeSo) |
|
| 834 |
sollya_lib_clear_obj(scaledBoundsSo) |
|
| 835 |
#print "Approximation error:", errorSa |
|
| 836 |
return resultArray |
|
| 837 |
# The returned interval upper bound does not reach the requested upper |
|
| 838 |
# upper bound: compute the next upper bound. |
|
| 839 |
# The following ratio is always >= 1 |
|
| 840 |
currentErrorRatio = approxPrecSa / errorSa |
|
| 841 |
# Starting point for the next upper bound. |
|
| 842 |
currentScaledUpperBoundSa = boundsSa.endpoints()[1] |
|
| 843 |
boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0] |
|
| 844 |
# Compute the increment. |
|
| 845 |
if currentErrorRatio > RR('1000'): # ]1.5, infinity[
|
|
| 846 |
currentScaledUpperBoundSa += \ |
|
| 847 |
currentErrorRatio * boundsWidthSa * 2 |
|
| 848 |
else: # [1, 1.5] |
|
| 849 |
currentScaledUpperBoundSa += \ |
|
| 850 |
(RR('1.0') + currentErrorRatio.log() / 500) * boundsWidthSa
|
|
| 851 |
# Take into account the original interval upper bound. |
|
| 852 |
if currentScaledUpperBoundSa > scaledUpperBoundSa: |
|
| 853 |
currentScaledUpperBoundSa = scaledUpperBoundSa |
|
| 854 |
# Compute the other expansions. |
|
| 855 |
while boundsSa.endpoints()[1] < scaledUpperBoundSa: |
|
| 856 |
currentScaledLowerBoundSa = boundsSa.endpoints()[1] |
|
| 857 |
(polySo, boundsSo, intervalCenterSo, maxErrorSo) = \ |
|
| 858 |
slz_compute_polynomial_and_interval(functionSo, degreeSo, |
|
| 859 |
currentScaledLowerBoundSa, |
|
| 860 |
currentScaledUpperBoundSa, |
|
| 861 |
approxPrecSa, |
|
| 862 |
internalSollyaPrecSa) |
|
| 849 |
boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField) |
|
| 863 | 850 |
errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo) |
| 864 |
if errorSa < approxPrecSa: |
|
| 865 |
# Change variable stuff |
|
| 851 |
print "Computed approximation error:", errorSa.n(digits=10) |
|
| 852 |
# If the error and interval are OK a the first try, just return. |
|
| 853 |
if (boundsSa.endpoints()[1] >= scaledUpperBoundSa) and \ |
|
| 854 |
(errorSa <= approxPrecSa): |
|
| 855 |
if precChangedSa: |
|
| 856 |
pobyso_set_prec_sa_so(currentSollyaPrecSa) |
|
| 857 |
sollya_lib_clear_obj(currentSollyaPrecSo) |
|
| 858 |
sollya_lib_clear_obj(functionSo) |
|
| 859 |
sollya_lib_clear_obj(degreeSo) |
|
| 860 |
sollya_lib_clear_obj(scaledBoundsSo) |
|
| 866 | 861 |
#print "Approximation error:", errorSa |
| 867 |
changeVarExpressionSo = sollya_lib_build_function_sub( |
|
| 868 |
sollya_lib_build_function_free_variable(), |
|
| 869 |
sollya_lib_copy_obj(intervalCenterSo)) |
|
| 870 |
polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo) |
|
| 871 |
sollya_lib_clear_obj(changeVarExpressionSo) |
|
| 872 |
resultArray.append((polySo, polyVarChangedSo, boundsSo, \ |
|
| 873 |
intervalCenterSo, maxErrorSo)) |
|
| 874 |
boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField) |
|
| 875 |
# Compute the next upper bound. |
|
| 876 |
# The following ratio is always >= 1 |
|
| 862 |
return resultArray |
|
| 863 |
## The returned interval upper bound does not reach the requested upper |
|
| 864 |
# upper bound: compute the next upper bound. |
|
| 865 |
## The following ratio is always >= 1. If errorSa, we may want to |
|
| 866 |
# enlarge the interval |
|
| 877 | 867 |
currentErrorRatio = approxPrecSa / errorSa |
| 878 |
# Starting point for the next upper bound. |
|
| 879 |
currentScaledUpperBoundSa = boundsSa.endpoints()[1] |
|
| 868 |
## --|--------------------------------------------------------------|-- |
|
| 869 |
# --|--------------------|-------------------------------------------- |
|
| 870 |
# --|----------------------------|------------------------------------ |
|
| 871 |
## Starting point for the next upper bound. |
|
| 880 | 872 |
boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0] |
| 881 | 873 |
# Compute the increment. |
| 882 |
if currentErrorRatio > RR('1000'): # ]1.5, infinity[
|
|
| 883 |
currentScaledUpperBoundSa += \ |
|
| 884 |
currentErrorRatio * boundsWidthSa * 2 |
|
| 885 |
else: # [1, 1.5] |
|
| 886 |
currentScaledUpperBoundSa += \ |
|
| 887 |
(RR('1.0') + currentErrorRatio.log()/500) * boundsWidthSa
|
|
| 888 |
#print "currentErrorRatio:", currentErrorRatio |
|
| 889 |
#print "currentScaledUpperBoundSa", currentScaledUpperBoundSa |
|
| 890 |
# Test for insufficient precision. |
|
| 891 |
if currentScaledUpperBoundSa == scaledLowerBoundSa: |
|
| 874 |
newBoundsWidthSa = \ |
|
| 875 |
((currentErrorRatio.log() / 10) + 1) * boundsWidthSa |
|
| 876 |
currentScaledLowerBoundSa = boundsSa.endpoints()[1] |
|
| 877 |
currentScaledUpperBoundSa = boundsSa.endpoints()[1] + newBoundsWidthSa |
|
| 878 |
# Take into account the original interval upper bound. |
|
| 879 |
if currentScaledUpperBoundSa > scaledUpperBoundSa: |
|
| 880 |
currentScaledUpperBoundSa = scaledUpperBoundSa |
|
| 881 |
if currentScaledUpperBoundSa == currentScaledLowerBoundSa: |
|
| 892 | 882 |
print "Can't shrink the interval anymore!" |
| 893 | 883 |
print "You should consider increasing the Sollya internal precision" |
| 894 | 884 |
print "or the polynomial degree." |
| 895 | 885 |
print "Giving up!" |
| 896 |
if internalSollyaPrecSa != currentSollyaPrecSa:
|
|
| 886 |
if precChangedSa:
|
|
| 897 | 887 |
pobyso_set_prec_sa_so(currentSollyaPrecSa) |
| 898 | 888 |
sollya_lib_clear_obj(currentSollyaPrecSo) |
| 899 | 889 |
sollya_lib_clear_obj(functionSo) |
| 900 | 890 |
sollya_lib_clear_obj(degreeSo) |
| 901 | 891 |
sollya_lib_clear_obj(scaledBoundsSo) |
| 902 | 892 |
return None |
| 903 |
if currentScaledUpperBoundSa > scaledUpperBoundSa: |
|
| 904 |
currentScaledUpperBoundSa = scaledUpperBoundSa |
|
| 905 |
if internalSollyaPrecSa > currentSollyaPrecSa: |
|
| 906 |
pobyso_set_prec_so_so(currentSollyaPrecSo) |
|
| 907 |
sollya_lib_clear_obj(currentSollyaPrecSo) |
|
| 908 |
sollya_lib_clear_obj(functionSo) |
|
| 909 |
sollya_lib_clear_obj(degreeSo) |
|
| 910 |
sollya_lib_clear_obj(scaledBoundsSo) |
|
| 911 |
return(resultArray) |
|
| 893 |
# Compute the other expansions. |
|
| 894 |
# Test for insufficient precision. |
|
| 912 | 895 |
# End slz_get_intervals_and_polynomials |
| 913 | 896 |
|
| 914 |
|
|
| 915 | 897 |
def slz_interval_scaling_expression(boundsInterval, expVar): |
| 916 | 898 |
""" |
| 917 | 899 |
Compute the scaling expression to map an interval that spans at most |
Formats disponibles : Unified diff