The Inbox: Kernel-nice.903.mcz
Locked
 
|
Nicolas Cellier uploaded a new version of Kernel to project The Inbox:
http://source.squeak.org/inbox/Kernel-nice.903.mcz ==================== Summary ==================== Name: Kernel-nice.903 Author: nice Time: 14 February 2015, 2:19:44.184 am UUID: d4f06790-d47b-4cb0-b9ff-9e0ff7dea9fc Ancestors: Kernel-nice.902 Introduce two new alternatives to integer division: #ratio: and #residue: are like #quo: and #rem: except that they round the quotient to nearest integer (tie to even) instead of truncating (note that // and \\ floor the quotient...). The second thing that they do differently is that they compute the exact ratio and exact residue when given a pair of Float. The third thing that they do differently is that they first coerce a pair of numbers to the highest generality before attempting any operation. =============== Diff against Kernel-nice.902 =============== Item was added: + ----- Method: Float>>fratio: (in category 'mathematical functions') ----- + fratio: aFloat + "Compute the exact ratio of self divided by aFloat (rounded to nearest integer, tie to even). + The result is either a Float in case of null (to wear correct sign) or exceptional value, + or an Integer (required in order to avoid overflow and inexact rounding)." + + | ex ey sx sy xa ya n result sign | + (self isFinite and: [self ~= 0.0 and: [aFloat isFinite and: [aFloat ~= 0.0]]]) + ifFalse: + ["purge cases of exceptional and trivial values" + ^self/aFloat]. + + xa := self abs. + ya := aFloat abs. + ex := xa exponent. + ey := ya exponent. + sx := xa significand. + sy := ya significand. + n := ex - ey. + sign := self significand*aFloat significand. + n < -1 ifTrue: [^sign copySignTo: 0.0]. + result := 0. + "subtract as many y significand as we can from x significand" + [(n >= 0 and: [sx >= sy]) + ifTrue: + [sx := sx - sy. + result := 1<<n+result]. + sx = 0.0 or: [n < 0]] + whileFalse: + [n := n - 1. + sx := sx * 2.0]. + (sx > sy or: [result odd and: [sx = sy]]) + ifTrue: + ["Round to nearest integer, tie to even" + result := result + 1]. + ^sign copySignTo: (result = 0 ifTrue: [0.0] ifFalse: [result])! Item was added: + ----- Method: Float>>fresidue: (in category 'mathematical functions') ----- + fresidue: aFloat + "Compute the exact residue of self divided by aFloat (rounded to nearest integer, tie to even). + This operation emulates libm remainder(self,aFloat) and is always exact. + The residue is between aFloat negated / 2 and aFloat / 2." + + | ex ey sx sy xa ya n sign odd roundToZero | + (self isFinite and: [self ~= 0.0 and: [aFloat isFinite and: [aFloat ~= 0.0]]]) + ifFalse: + ["purge cases of exceptional values" + ^self*aFloat / (self * aFloat)]. + + xa := self abs. + ya := aFloat abs. + ex := xa exponent. + ey := ya exponent. + sx := xa significand. + sy := ya significand. + n := ex - ey. + n < -1 ifTrue: [^self]. + sign := self. + roundToZero := true. + odd := false. + "subtract as many y significand as we can from x significand" + [(n >= 0 and: [sx >= sy]) + ifTrue: + [sx := sx - sy. + n = 0 ifTrue: [odd := true]. + sx = 0.0 ifTrue: [^sign copySignTo: 0.0]]. + n < 0] + whileFalse: + [n := n - 1. + sx := sx * 2.0]. + (sx > sy or: [odd and: [sx = sy]]) + ifTrue: + ["Round to nearest integer, tie to even" + sx := sx * 0.5 - sy. + sign := sign negated] + ifFalse: + [sx := sx * 0.5]. + ^(sign copySignTo: sx) timesTwoPower: ey! Item was added: + ----- Method: Float>>ratio: (in category 'arithmetic') ----- + ratio: aNumber + aNumber isFloat ifTrue: [^self fratio: aNumber]. + ^aNumber adaptToFloat: self andSend: #ratio:! Item was added: + ----- Method: Float>>residue: (in category 'arithmetic') ----- + residue: aNumber + aNumber isFloat ifTrue: [^self fresidue: aNumber]. + ^aNumber adaptToFloat: self andSend: #residue:! Item was added: + ----- Method: Fraction>>ratio: (in category 'arithmetic') ----- + ratio: aNumber + aNumber isFraction ifTrue: [^super ratio: aNumber]. + ^aNumber adaptToFraction: self andSend: #ratio:! Item was added: + ----- Method: Fraction>>residue: (in category 'arithmetic') ----- + residue: aNumber + aNumber isFraction ifTrue: [^super residue: aNumber]. + ^aNumber adaptToFraction: self andSend: #residue:! Item was added: + ----- Method: Fraction>>roundedTieToEven (in category 'truncation and round off') ----- + roundedTieToEven + "Faster than super." + | result | + denominator = 2 ifFalse: [^self rounded]. + result := numerator abs bitShift: -1. + result := result + (result bitAt: 1). + ^numerator positive + ifTrue: [result] + ifFalse: [result negated]! Item was added: + ----- Method: Integer>>ratio: (in category 'arithmetic') ----- + ratio: aNumber + aNumber isFraction ifTrue: [^super ratio: aNumber]. + ^aNumber adaptToInteger: self andSend: #ratio:! Item was added: + ----- Method: Integer>>residue: (in category 'arithmetic') ----- + residue: aNumber + aNumber isFraction ifTrue: [^super residue: aNumber]. + ^aNumber adaptToInteger: self andSend: #residue:! Item was added: + ----- Method: Integer>>roundedTieToEven (in category 'truncation and round off') ----- + roundedTieToEven + ^self! Item was added: + ----- Method: Number>>ratio: (in category 'arithmetic') ----- + ratio: aNumber + "Integer quotient defined by division with rounding, tie to nearest even. + residue: answers the remainder from this division." + + ^(self / aNumber) roundedTieToEven! Item was added: + ----- Method: Number>>residue: (in category 'arithmetic') ----- + residue: aNumber + "Remainder defined in terms of ratio:. + Answer a Number between aNumber/-2 and aNumber/2" + + ^self - ((self ratio: aNumber) * aNumber)! Item was added: + ----- Method: Number>>roundedTieToEven (in category 'truncation and round off') ----- + roundedTieToEven + "Answer the integer nearest the receiver. + In case of perfect tie, round to nearest even. + Unlike rounded (round away from zero), this mode of rounding is unbiased." + + | sign truncated fractionPart | + sign := self sign. + truncated := self truncated. + fractionPart := self - self truncated. + ^(fractionPart abs > (1/2) or: [fractionPart = (sign/2) and: [truncated odd]]) + ifTrue: [truncated + sign] + ifFalse: [truncated]! Item was added: + ----- Method: ScaledDecimal>>ratio: (in category 'arithmetic') ----- + ratio: aNumber + "Answer the integer quotient after dividing the receiver by aNumber + with rounding, tie to even." + (aNumber isKindOf: ScaledDecimal) ifTrue: [^super ratio: aNumber]. + ^aNumber adaptToScaledDecimal: self andSend: #ratio:! Item was added: + ----- Method: ScaledDecimal>>residue: (in category 'arithmetic') ----- + residue: aNumber + (aNumber isKindOf: ScaledDecimal) ifTrue: [^super residue: aNumber]. + ^aNumber adaptToScaledDecimal: self andSend: #residue:! |
Re: The Inbox: Kernel-nice.903.mcz
|
|
2015-02-14 2:20 GMT+01:00 <[hidden email]>: Nicolas Cellier uploaded a new version of Kernel to project The Inbox: Hi folks, I've already bent Squeak/Pharo handling of Float a bit, and I'd like to increase further Smalltalk dialects compliance to language independent arithmetic standard (ISO/IEC 10967). The subject of the day is state of the art quotient and remainder operations for Floats. I
noticed that the libm function fmod which computes the remainder of
division of two floats (quotient being the truncated fraction) is and
should be exact. Our equivalent Float version of #rem: is not, and I had the feeling that it ain't good. They are not only exact, they also support extreme Floats without overflow, like (Float fmax quo: Float fmin*5). (Float fmax rem: Float fmin*5). Once plugged, some results might unfortunately look surprising like: (1 to: 100) count: [:i | ((i/100.0) fquo: 0.01) ~= i]. -> 58 (1 to: 100) count: [:i | ((i/100.0) fmod: 0.01) ~= 0]. -> 93 The seven exact multiple are of course (0 to: 6) collect:[:i | 0.01 timesTwoPower: i] The others just tell the awful truth about floating point...
While with the old implementation was a little less surprising:(1 to: 100) count: [:i | ((i/100.0) / 0.01) truncated ~= i]. -> 6 (1 to: 100) count: [:i | (((i/100.0) / 0.01) truncated * 0.01 - (i/100.0)) ~= 0]. -> 13 Then I realized that language independent arithmetic standard and also older floating point standard (IEEE 754), do rather define quotient and remainder (ratio and residue) in term of rounded rather than truncated quotient. So I implemented fratio: and fresidue: as prescribed by the standard (with an unbiased roundedTieToEven). And these two operations better fit our expectations, at least for the ratio:(1 to: 100) count: [:i | ((i/100.0) fratio: 0.01) ~= i]. -> 0 The residues are exact allways small: (1 to: 100) count: [:i | | a b fa fb ea eb | a := 1/100. fa := a asFloat. ea := a - fa asFraction. "exact, float, and error" b := i*a. fb := b asFloat.. eb := b - fb asFraction. [eb <= (0.5*fb ulp) and: [ea <= (0.5*fa ulp)]] assert. (fb fresidue: fa) abs > (ea*i-eb) abs]. -> 0 The naive implementation works well too for the ratio: (1 to: 100) count: [:i | ((i/100.0) / 0.01) rounded ~= i]. -> 0 but of course with inexact residue: (1 to: 100) count: [:i | | a b fa fb ea eb | a := 1/100. fa := a asFloat. ea := a - fa asFraction. "exact, float, and error" b := i*a. fb := b asFloat.. eb := b - fb asFraction. [eb <= (0.5*fb ulp) and: [ea <= (0.5*fa ulp)]] assert. (fb - ((fb/fa) rounded*fb)) abs > (ea*i-eb) abs]. -> 99 Rounded quotient is exactly what we want for Float. The rounding problems still exist, but they are rejected near exact tie. But that means implementing yet another quo/rem message (#ratio: #residue: using terms of the standard) while we already have two pairs... Do you think that such extension ratio/residue is valuable? Do you think that modification of quo:/rem: is bearable? The alternative would be to provide an external package (eventually with overrides) to turn Squeak/Pharo behavior into standard, but I don't like the idea that the kernel does not respect the standards by default... ------------------------------------------------ About // and \\: The problem with \\ is that it can't always be exact when expressed as Float. (b \\ a) can be near a when b abs < a abs, a*b < 0. So when b is of very small magnitude relative to a, remainder requires many bits. For example Float fmin negated asFraction \\ 1.0 asFraction would exactly be (1.0 asFraction - Float fmin asFraction) and thus require a significand of 1074 bits. We don't have ArbitraryPrecisionFloat in trunk, and Fraction is overkill... We can't either let the exact remainder be rounded to nearest Float because it will be equal to the divider (it should be less than in magnitude). So, there is no easy solution for this case in Float... But note that current naive implementation suffers exactly from same problem and would answer 1.0 (the divider). ------------------------------------------------ I've also asked peers for a good reason to implement an exact #rem: but knowledgeable peers are unreachable yet ;) http://cs.stackexchange.com/questions/24362/why-does-floating-point-modulus-exactness-matters So maybe it's time to push this stuff |
«
Return to Squeak - Dev
|
1 view|%1 views
| Free forum by Nabble | Edit this page |