On Wed, Aug 14, 2024 at 04:34:23PM +0200, Waldek Hebisch wrote:
> I am looking more into output convertion for floats and current
> implementation looks somewhat irregular.
I have a little correction to what I wrote above: default output
precision is determined by bit length of mantissa. This means
that there will be some discrepancy converting this to number
of digits. As implemented, number of digits may be smaller.
Usiually this is OK, because value is rounded which leads to
simpler output. But numbers with exponent 0 are in fact
exact integers, so we want to print them in full accuracy.
So special case for exponent 0 is justified.
> I think about following changes:
> - make slightly different rule in "general" convertion: use fixed
> form when there are at most 5 leading zeros before most significant
> digit (this is current rule) or there are at most 2 trailing
> zeros before the dot (this is change compared to current rule),
> - with specified precision print given number of digits in
> "general" convertion (this is change, current code may skip trailing
> zeros in such case)
> - with default precision (that is taken from floating point precision)
> trim trailng zeros (this is mostly current behaviour, but I would
> like to eliminate exceptions).
Attached is my current code, it implements changed rule for general
output. Changes in test result are quite small. As you can see
code is quite different than the current one. To use 'fixed_mr'
in FloatingPointSystem I had to export it. I did not want to
export it from Float (which is exposed), so it had to go to
a separate new package. 'fixed_mr' needs several support routines,
so they had to go to the new package. To avoid fiddling with
global variables all needed information is passed via parameters.
In particular there is a new package exporting operations which
beside normal arguments also take requested precision as an
extra argument. Once I was doing that I also implemented some
my ideas for improvements to elementary functions, in particular
logarithm is computed solving equation 'x = exp(y)' (it is currently
called 'log_newton' and 'log_newton2', I probably need a better
name as this is different than Newton method).
Anyway, there is new file containg new routines. Currently changes
in behaviour are small, but it is easy to add new fields to mode
record to allow more variation. AFAICS possibly controlable
behaviour includes:
- suppressing (or not) trailing digits
- putting spaces (or not) around "E" in exponential notation
- we could easily make separator configurable
We could probably add option for different "floating" notation,
so that dot goes after first significant digit. In principle
we could also put dot in other positions, but I am not sure
if that would be good idea (I suspect that exponential notation
with dot in different place than befor or after first significant
digit are less readable).
--
Waldek Hebisch
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/fricas-devel/Zr-JAnt_JjU8_QXY%40fricas.org.
diff -Naur '--exclude=.git' trunk/src/algebra/float.spad trunk.p1/src/algebra/float.spad
--- trunk/src/algebra/float.spad 2024-04-28 15:52:10.132572747 +0200
+++ trunk.p1/src/algebra/float.spad 2024-08-16 17:50:28.253002665 +0200
@@ -166,8 +166,6 @@
atanSeries : % -> % -- atan(x) by taylor series |x| < 1/2
atanInverse : I -> % -- atan(1/n) for n an integer > 1
expInverse : I -> % -- exp(1/n) for n an integer
- expSeries : % -> % -- exp(x) by taylor series |x| < 1/2
- logSeries : % -> % -- log(x) by taylor series 1/2 < x < 2
sinSeries : % -> % -- sin(x) by taylor series |x| < 1/2
cosSeries : % -> % -- cos(x) by taylor series |x| < 1/2
piRamanujan : () -> % -- pi using Ramanujans series
@@ -394,105 +392,16 @@
inc(5-p); r := log((x+1)/(1-x))/2; dec(5-p)
normalize r
- log x ==
- negative? x => error1("log: negative argument", x)
- zero? x => error "log 0 generated"
- p := bits(); inc 5
- -- apply log(x) = n log 2 + log(x/2^n) so that 1/2 < x < 2
- if (n := order x) < 0 then n := n+1
- l := if n = 0 then 0 else (x := shift(x, -n); n * log2)
- -- speed the series convergence by finding m and k such that
- -- | exp(m/2^k) x - 1 | < 1 / 2 ^ O(sqrt p)
- -- write log(exp(m/2^k) x) as m/2^k + log x
- k := ISQRT (p-100)::I quo 3
- if k > 1 then
- k := max(1, k+order(x-1))
- inc k
- ek := expInverse (2^k::N)
- dec(p quo 2); m := order square(x, k); inc(p quo 2)
- m := (6847196937 * m) quo 9878417065 -- m := m log 2
- x := x * ek ^ (-m)
- l := l + [m, -k]
- l := l + logSeries x
- bits p
- normalize l
+ log(x) == normalize(log(x::Rep, bits())$FloatElementaryFunctions)
- logSeries x ==
- -- log(x) = 2 y (1 + y^2/3 + y^4/5 ...) for y = (x-1) / (x+1)
- -- given 1/2 < x < 2 on input we have -1/3 < y < 1/3
- p := bits() + (g := LENGTH bits() + 3)
- inc g; y := (x-1)/(x+1); dec g
- s : I := d : I := shift(1, p)
- z := times(y, y)
- t := m := shift2(z.mantissa, z.exponent+p)
- for i in 3.. by 2 while t ~= 0 repeat
- s := s + t quo i
- t := m * t quo d
- y * [s, 1-p]
+ log2() == normalize(log2(bits())$FloatElementaryFunctions)
- L2 : StoredConstant := [1, 1]
- log2() ==
- -- log x = 2 * sum( ((x-1)/(x+1))^(2*k+1)/(2*k+1), k=1.. )
- -- log 2 = 2 * sum( 1/9^k / (2*k+1), k=0..n ) / 3
- n := bits() :: N
- n <= L2.precision => normalize L2.value
- n := n + LENGTH n + 3 -- guard bits
- s : I := shift(1, n+1) quo 3
- t : I := s quo 9
- for k in 3.. by 2 while t ~= 0 repeat
- s := s + t quo k
- t := t quo 9
- L2 := [bits(), [s, -n]]
- normalize L2.value
-
- L10 : StoredConstant := [1, [1, 1]]
- log10() ==
- -- log x = 2 * sum( ((x-1)/(x+1))^(2*k+1)/(2*k+1), k=0.. )
- -- log 5/4 = 2 * sum( 1/81^k / (2*k+1), k=0.. ) / 9
- n := bits() :: N
- n <= L10.precision => normalize L10.value
- n := n + LENGTH n + 5 -- guard bits
- s : I := shift(1, n+1) quo 9
- t : I := s quo 81
- for k in 3.. by 2 while t ~= 0 repeat
- s := s + t quo k
- t := t quo 81
- -- We have log 10 = log 5 + log 2 and log 5/4 = log 5 - 2 log 2
- inc 2; L10 := [bits(), [s, -n] + 3*log2]; dec 2
- normalize L10.value
+ log10() == normalize(log10(bits())$FloatElementaryFunctions)
log2(x) == (inc 2; r := log(x)/log2; dec 2; normalize r)
log10(x) == (inc 2; r := log(x)/log10; dec 2; normalize r)
- exp(x) ==
- -- exp(n+x) = exp(1)^n exp(x) for n such that |x| < 1
- p := bits(); inc 5; e1 : % := 1
- if (n := wholePart x) ~= 0 then
- inc LENGTH n; e1 := exp1 ^ n; dec LENGTH n
- x := fractionPart x
- if zero? x then (bits p; return normalize e1)
- -- make |x| < O( 2^(-sqrt p) ) < 1/2 to speed series convergence
- -- by repeated use of the formula exp(2*x/2) = exp(x/2)^2
- -- results in an overall running time of O( sqrt p M(p) )
- k := ISQRT (p-100)::I quo 3
- k := max(0, 2 + k + order x)
- if k > 0 then (inc k; x := shift(x, -k))
- e := expSeries x
- if k > 0 then e := square(e, k)
- bits p
- e * e1
-
- expSeries x ==
- -- exp(x) = 1 + x + x^2/2 + ... + x^i/i! valid for all x
- p := bits() + LENGTH bits() + 1
- s : I := d : I := shift(1, p)
- t : I := n : I := shift2(x.mantissa, x.exponent+p)
- for i in 2.. while t ~= 0 repeat
- s := s + t
- t := (n * t) quo i
- t := t quo d
- normalize [s, -p]
-
+ exp(x) == normalize(exp(x::Rep, bits())$FloatElementaryFunctions)
Ress ==> Record(highn : I, matr : Matrix(I))
@@ -717,78 +626,6 @@
bits p
normalize y
- -- Utility routines for conversion to decimal
- ceilLength10 : I -> I
- chop10 : (%, I) -> %
- convert10 : (%, I) -> %
- floorLength10 : I -> I
- length10 : I -> I
- normalize10 : (%, I) -> %
- quotient10 : (%, %, I) -> %
- power10 : (%, I, I) -> %
- times10 : (%, %, I) -> %
-
- convert10(x, d) ==
- m := x.mantissa; e := x.exponent
- --!! deal with bits here
- b := bits(); (q, r) := divide(abs e, b)
- b := 2^b::N; r := 2^r::N
- -- compute 2^e = b^q * r
- h := power10([b, 0], q, d+5)
- h := chop10([r*h.mantissa, h.exponent], d+5)
- if e < 0 then h := quotient10([m, 0], h, d)
- else times10([m, 0], h, d)
-
- ceilLength10 n == 146::I * length(n) quo 485 + 1
- floorLength10 n == 643::I * length(n) quo 2136
- length10 n ==
- ln : Integer := LENGTH(n := abs n)
- upper := 76573 * ln quo 254370
- lower := 21306 * (ln-1) quo 70777
- upper = lower => upper + 1
- n := n quo (10^lower::N)
- while n >= 10 repeat
- n := n quo 10
- lower := lower + 1
- lower + 1
-
- chop10(x, p) ==
- e : I := floorLength10 x.mantissa - p
- if e > 0 then x := [x.mantissa quo 10^e::N, x.exponent+e]
- x
- normalize10(x, p) ==
- ma := x.mantissa
- ex := x.exponent
- e : I := length10 ma - p
- if e > 0 then
- ma := ma quo 10^(e-1)::N
- ex := ex + e
- (ma, r) := divide(ma, 10)
- if r > 4 then
- ma := ma + 1
- if ma = 10^p::N then (ma := 1; ex := ex + p)
- [ma, ex]
- times10(x, y, p) == normalize10(times(x, y), p)
- quotient10(x, y, p) ==
- ew := floorLength10 y.mantissa - ceilLength10 x.mantissa + p + 2
- if ew < 0 then ew := 0
- mw := (x.mantissa * 10^ew::N) quo y.mantissa
- ew := x.exponent - y.exponent - ew
- normalize10([mw, ew], p)
- power10(x, n, d) ==
- x = 0 => 0
- n = 0 => 1
- n = 1 => x
- x = 1 => 1
- p : I := d + LENGTH n + 1
- e : I := n
- y : % := 1
- z : % := x
- repeat
- if odd? e then y := chop10(times(y, z), p)
- if (e := e quo 2) = 0 then return y
- z := chop10(times(z, z), p)
-
--------------------------------
-- Output routines for Floats --
--------------------------------
@@ -803,112 +640,7 @@
floating : % -> S
general : % -> S
- padFromLeft(s : S) : S ==
- zero? SPACING => s
- n : I := #s - 1
- t := new( (n + 1 + n quo SPACING) :: N , separator )
- for i in 0..n for j in minIndex t .. repeat
- t.j := s.(i + minIndex s)
- if (i+1) rem SPACING = 0 then j := j+1
- t
- padFromRight(s : S) : S ==
- SPACING = 0 => s
- n : I := #s - 1
- t := new( (n + 1 + n quo SPACING) :: N , separator )
- for i in n..0 by -1 for j in maxIndex t .. by -1 repeat
- t.j := s.(i + minIndex s)
- if (n-i+1) rem SPACING = 0 then j := j-1
- t
-
- fixed f ==
- d := if OUTPREC = -1 then digits()::I else OUTPREC
- zero? f =>
- OUTPREC = -1 => "0.0"
- concat("0", concat(".", padFromLeft new(d::N, zero)))
- zero? exponent f =>
- concat(padFromRight convert(mantissa f)@S,
- concat(".", padFromLeft new(d::N, zero)))
- negative? f => concat("-", fixed abs f)
- bl := LENGTH(f.mantissa) + f.exponent
- dd :=
- OUTPREC = -1 => d
- bl > 0 => (146*bl) quo 485 + 1 + d
- d
- g := convert10(abs f, dd)
- m := g.mantissa; e := g.exponent
- if OUTPREC ~= -1 then
- -- round g to OUTPREC digits after the decimal point
- l := length10 m
- if -e > OUTPREC and -e < 2*digits()::I then
- g := normalize10(g, l+e+OUTPREC)
- m := g.mantissa; e := g.exponent
- s := convert(m)@S; n := #s; o := e+n
- p := if OUTPREC = -1 then n::I else OUTPREC
- t : S
- if e >= 0 then
- s := concat(s, new(e::N, zero))
- t := ""
- else if o <= 0 then
- t := concat(new((-o)::N, zero), s)
- s := "0"
- else
- t := s(o + minIndex s .. n + minIndex s - 1)
- s := s(minIndex s .. o + minIndex s - 1)
- n := #t
- if OUTPREC = -1 then
- t := rightTrim(t, zero)
- if t = "" then t := "0"
- else if n > p then t := t(minIndex t .. p + minIndex t- 1)
- else t := concat(t, new((p-n)::N, zero))
- concat(padFromRight s, concat(".", padFromLeft t))
-
- floating f ==
- zero? f => "0.0"
- negative? f => concat("-", floating abs f)
- t:S := if zero? SPACING then "E" else " E "
- zero? exponent f =>
- s := convert(mantissa f)@S
- concat ["0.", padFromLeft s, t, convert(#s)@S]
- -- base conversion to decimal rounded to the requested precision
- d := if OUTPREC = -1 then digits()::I else OUTPREC
- g := convert10(f, d); m := g.mantissa; e := g.exponent
- -- I'm assuming that length10 m = # s given n > 0
- s := convert(m)@S; n := #s; o := e+n
- s := padFromLeft s
- concat ["0.", s, t, convert(o)@S]
-
- general(f) ==
- zero? f => "0.0"
- negative? f => concat("-", general abs f)
- d := if OUTPREC = -1 then digits()::I else OUTPREC
- zero? exponent f =>
- d := d + 1
- s := convert(mantissa f)@S
- OUTPREC ~= -1 and (e := #s) > d =>
- t:S := if zero? SPACING then "E" else " E "
- concat ["0.", padFromLeft s, t, convert(e)@S]
- concat(padFromRight(s), ".0")
- -- base conversion to decimal rounded to the requested precision
- g := convert10(f, d); m := g.mantissa; e := g.exponent
- -- I'm assuming that length10 m = # s given n > 0
- s := convert(m)@S; n := #s; o := n + e
- -- Note: at least one digit is displayed after the decimal point
- -- and trailing zeroes after the decimal point are dropped
- if o > 0 and o <= max(n, d) then
- -- fixed format: add trailing zeroes before the decimal point
- if o > n then s := concat(s, new((o-n)::N, zero))
- t := rightTrim(s(o + minIndex s .. n + minIndex s - 1), zero)
- if t = "" then t := "0" else t := padFromLeft t
- s := padFromRight s(minIndex s .. o + minIndex s - 1)
- concat(s, concat(".", t))
- else if o <= 0 and o >= -5 then
- -- fixed format: up to 5 leading zeroes after the decimal point
- concat("0.",padFromLeft concat(new((-o)::N,zero),rightTrim(s,zero)))
- else
- -- print using E format written 0.mantissa E exponent
- t := padFromLeft rightTrim(s, zero)
- s := if zero? SPACING then "E" else " E "
- concat ["0.", t, s, convert(e+n)@S]
+ mode_Rec ==> Record(default? : B, out_mode : S, spacing : I)
outputSpacing n ==
old_val := SPACING
@@ -936,18 +668,12 @@
OUTPREC := prec
convert(f) : S ==
- b : Integer :=
- OUTPREC = -1 and not zero? f =>
- bits(length(abs mantissa f)::PositiveInteger)
- 0
- s :=
- OUTMODE = "fixed" => fixed f
- OUTMODE = "floating" => floating f
- OUTMODE = "general" => general f
- empty()$String
- if b > 0 then bits(b::PositiveInteger)
- s = empty()$String => error "bad output mode"
- s
+ d : Integer :=
+ OUTPREC = -1 and not zero? f =>
+ max(1, 4004*(length(abs(mantissa(f))) - 1) quo 13301)
+ OUTPREC
+ mr := [OUTPREC = -1, OUTMODE, SPACING]$mode_Rec
+ convert_mr(f::Rep, d, mr)$FloatingPointConvertion
coerce(f) : OutputForm ==
f >= 0 => message(convert(f)@S)
diff -Naur '--exclude=.git' trunk/src/algebra/flopak.spad trunk.p1/src/algebra/flopak.spad
--- trunk/src/algebra/flopak.spad 1970-01-01 01:00:00.000000000 +0100
+++ trunk.p1/src/algebra/flopak.spad 2024-08-16 18:42:33.600558944 +0200
@@ -0,0 +1,437 @@
+)abbrev package FLOOPS FloatingPointOperations
+FloatingPointOperations : Exports == Implementation where
+ I ==> Integer
+ F ==> Record(mantissa : I, exponent : I)
+ Exports ==> with
+ plus : (F, F) -> F
+ ++ plus(x, y) computes \spad{x + y} with no rounding.
+ minus : (F, F) -> F
+ ++ minus(x, y) computes \spad{x - y} with no rounding.
+ times : (F, F) -> F
+ ++ times(x, y, b) computes \spad{x*y} with no rounding.
+
+ Implementation ==> add
+
+ shift(x, y) ==> ASH(x, y)$Lisp
+
+ plus(x, y) ==
+ (xm, xe) := x
+ (ym, ye) := y
+ xm = 0 => y
+ ym = 0 => x
+ xe = ye => [xm + ym, ye]
+ de := xe - ye
+ de > 0 => [shift(xm, de) + ym, ye]
+ [xm + shift(ym, -de), xe]
+
+ minus(x, y) ==
+ (ym, ye) := y
+ plus(x, [-ym, ye])
+
+ times(x, y) ==
+ (xm, xe) := x
+ (ym, ye) := y
+ [xm*ym, xe + ye]
+
+)abbrev package FELEM FloatElementaryFunctions
+FloatElementaryFunctions : Exports == Implementation where
+ I ==> Integer
+ PI ==> PositiveInteger
+ N ==> NonNegativeInteger
+ F ==> Record(mantissa : I, exponent : I)
+ DF ==> DoubleFloat
+ Exports ==> with
+ log2 : PI -> F
+ ++ log2(b) computes approximation to log(2) accurate to at least
+ ++ b bits.
+ log10 : PI -> F
+ ++ log10(b) computes approximation to log(10) accurate to at least
+ ++ b bits.
+ exp : (F, PI) -> F
+ ++ exp(x, b) computes approximation to exp(x) accurate to at least
+ ++ b bits.
+ log : (F, PI) -> F
+ ++ log(x, b) computes approximation to exp(x) accurate to at least
+ ++ b bits.
+ plus : (F, F, PI) -> F
+ ++ plus(x, y, b) computes approximation to \spad{x + y} accurate
+ ++ to at least b bits.
+ minus : (F, F, PI) -> F
+ ++ minus(x, y, b) computes approximation to \spad{x - y} accurate
+ ++ to at least b bits.
+ times : (F, F, PI) -> F
+ ++ times(x, y, b) computes approximation to \spad{x*y} accurate
+ ++ to at least b bits.
+ quotient : (F, F, PI) -> F
+ ++ quotient(x, y, b) computes approximation to \spad{x/y} accurate
+ ++ to at least b bits.
+ exp_series : (F, PI) -> F
+ log_series : (F, PI) -> F
+ Implementation ==> add
+
+ shift(x, y) ==> ASH(x, y)$Lisp
+
+ chop(x : F, b : PI) : F ==
+ (xm, xe) := x
+ l1 := length(xm)
+ l1 < b + 4 => x
+ db := l1 - b - 3
+ [shift(xm, -db), xe + db]
+
+ plus(x : F, y : F, b : PI) : F ==
+ (xm, xe) := x
+ (ym, ye) := y
+ de := xe - ye
+ dl := length(xm) - length(ym)
+ de + dl > b + 4 => x
+ de + dl < -b - 4 => y
+ chop(plus(x, y)$FloatingPointOperations, b)
+
+ minus(x : F, y : F, b : PI) : F ==
+ plus(x, [-y.mantissa, y.exponent], b)
+
+ times(x : F, y : F, b : PI) : F ==
+ chop(times(x, y)$FloatingPointOperations, b)
+
+ quotient(x : F, y : F, b : PI) : F ==
+ error "unimplemented"
+
+
+ -- compute exp(f) using Taylor series
+ exp_series(f : F, b : PI) : F ==
+ (m, e) := f
+ b1 := b + length(b) + 3
+ m := shift(m, e + b1)
+ -- start with f, we will add 1 at the end
+ s := tk := m
+ k := 2
+ repeat
+ tk := (m*tk) quo k
+ tk := shift(tk, -b1)
+ tk = 0 or tk = -1 => break
+ s := s + tk
+ k := k + 1
+ [shift(s + shift(1, b1), b - b1 + 3), -b - 3]
+
+ pow2k(x : F, k : I, b : PI) : F ==
+ for i in 1..k repeat
+ x := times(x, x, b)
+ x
+
+ exp(x : F, b : PI) : F ==
+ (m, e) := x
+ lm := length(m)
+ de := lm + e
+ de < -3 => exp_series(x, b)
+ k := de + 4
+ e1 := exp_series([m, e - k], b + k::PI)
+ chop(pow2k(e1, k, b + length(b)::PI + 4), b)
+
+ val_Rec ==> Record(val : F, prec : PI)
+
+ log2_val : val_Rec := [[235865763225513294137944142764154484399, -128],
+ 128]
+
+ log10_val : val_Rec := [[195882276370220766701334620675861842472, -126],
+ 125]
+
+ log_series(x : F, b : PI) : F ==
+ (m, e) := x
+ m := shift(m, b + e)
+ s := tk := m
+ k := 2
+ repeat
+ tk := (m*tk)
+ tk := shift(-tk, -b)
+ tk = 0 or tk = -1 => return [s, -b]
+ s := s + (tk quo k)
+ k := k + 1
+
+ log_newton2(y : F, x0 : F, b0 : PI, b : PI) : F ==
+ repeat
+ b1 : PI :=
+ b < 5*b0 => b
+ 3*b0
+ (x0m, x0e) := x0
+ y1 := exp([-x0m, x0e], b1 + 10)
+ yd := minus(times(y, y1, b1 + 10), [1, 0], b1 + 10)
+ xd := log_series(yd, b1 + 10)
+ x0 := plus(x0, xd, b1 + 2)
+ not(b1 < b) => return x0
+ b0 := b1
+
+ log2(b : PI) : F ==
+ p1 := log2_val.prec
+ v0 := log2_val.val
+ p1 = b => v0
+ b < p1 =>
+ p1 < b + 10 => v0
+ (m, e) := v0
+ de := p1 - b - 1
+ [shift(m, -de), e + de]
+ v0 := log_newton2([2, 0], v0, p1, b)
+ log2_val.val := v0
+ log2_val.prec := b
+ v0
+
+ log10(b : PI) : F ==
+ p1 := log10_val.prec
+ v0 := log10_val.val
+ p1 = b => v0
+ b < p1 =>
+ p1 < b + 10 => v0
+ (m, e) := v0
+ de := p1 - b - 1
+ [shift(m, -de), e + de]
+ v0 := log_newton2([10, 0], v0, p1, b)
+ log10_val.val := v0
+ log10_val.prec := b
+ v0
+
+ log_newton(x : F, b : PI) : F ==
+ (m, e) := x
+ m1 := m
+ e1 := e
+ if 50 < -e then
+ e1 := -50
+ m1 := shift(m, e - e1)
+ l0 := log(float(m1, e1, 2)$DF)$DF
+ m0 := mantissa(l0)
+ e0 := exponent(l0)
+ log_newton2(x, [m0, e0], 40, b)
+
+ log(x : F, b : PI) : F ==
+ (m, e) := x
+ not(0 < m) => error "argument to log must be positive"
+ lm := length(m)
+ m1 : I := shift(m, 6 - lm)
+ if m1 < 42 then lm := lm - 1
+ e1 := lm + e
+ e1 = 0 => log_newton(x, b)
+ l1 := log_newton([m, -lm], b + 3)
+ plus(l1, times([e1, 0], log2(b + 3), b + 3), b + 2)
+
+)abbrev package FOUT FloatingPointConvertion
+FloatingPointConvertion : Exports == Implementation where
+ B ==> Boolean
+ I ==> Integer
+ S ==> String
+ PI ==> PositiveInteger
+ RN ==> Fraction Integer
+ SF ==> DoubleFloat
+ N ==> NonNegativeInteger
+ F ==> Record(mantissa : I, exponent : I)
+ mode_Rec ==> Record(default? : B, out_mode : S, spacing : I)
+ Exports ==> with
+ fixed_mr : (F, I, mode_Rec) -> S
+ floating_mr : (F, I, mode_Rec) -> S
+ general_mr : (F, I, mode_Rec) -> S
+
+ convert_mr : (F, I, mode_Rec) -> S
+
+ convert_to_decimal :(F, I) -> F
+
+ Implementation ==> add
+
+ convert_mr(f : F, d : I, mr : mode_Rec) : S ==
+ mr.out_mode = "fixed" => fixed_mr(f, d, mr)
+ mr.out_mode = "floating" => floating_mr(f, d, mr)
+ mr.out_mode = "general" => general_mr(f, d, mr)
+ error "convert_mr: bad output mode"
+
+ convert_to_decimal(f : F, d : I) : F ==
+ ba0 := d*70777 quo 21306 + 1
+ m := f.mantissa
+ e := f.exponent
+ ea := abs(length(m) + e)
+ ba1 := (ba0 + length(ea) + 25)::PI
+ if length(m) > ba1 then
+ ed := length(m) - ba1
+ m := m quo 2^(ed::N)
+ e := e + ed
+ lm := log([m, 0], ba1)$FloatElementaryFunctions
+ le := times([e, 0], log2(ba1)$FloatElementaryFunctions, ba1
+ )$FloatElementaryFunctions
+ lf := plus(lm, le, ba1)$FloatElementaryFunctions
+ l10 := log10(ba1)$FloatElementaryFunctions
+ ed := lf.exponent - l10.exponent
+ n10 := lf.mantissa
+ d10 := l10.mantissa
+ if ed < 0 then
+ ed := -ed
+ d10 := d10*2^(ed::N)
+ else
+ n10 := n10*2^(ed::N)
+ e10 := n10 quo d10 - d
+ le10 := times([e10, 0], l10, ba1)$FloatElementaryFunctions
+ ldm := minus(lf, le10, ba1)$FloatElementaryFunctions
+ dm := exp(ldm, ba1)$FloatElementaryFunctions
+ dm.exponent >= 0 =>
+ error "impossible"
+ (q, r) := divide(dm.mantissa, 2^((-dm.exponent - 1)::N))
+ dfm1 := q quo 2
+ if odd?(q) and r > 0 then dfm1 := dfm1 + 1
+ if dfm1 >= 10^(d::N) then
+ e10 := e10 + 1
+ (q1, r1) := divide(q, 10)
+ dfm1 := q1 quo 2
+ if odd?(q1) and (r1 > 0 or r > 0) then dfm1 := dfm1 + 1
+ [dfm1, e10]
+
+ --------------------------------
+ -- Output routines for Floats --
+ --------------------------------
+ zero ==> char("0")
+ separator ==> underscore()$Character
+
+ insert_separators_from_left(s : S, mr : mode_Rec) : S ==
+ sp := mr.spacing
+ zero?(sp) => s
+ n : I := #s
+ t := new((n + (n - 1) quo sp)::N, separator)
+ j : I := 1
+ for i in 1..n repeat
+ t(j) := s(i)
+ if i rem sp = 0 then j := j + 1
+ j := j + 1
+ t
+
+ insert_separators_from_right(s : S, mr : mode_Rec) : S ==
+ sp := mr.spacing
+ zero?(sp) => s
+ n : I := #s
+ t := new((n + (n - 1) quo sp)::N, separator)
+ j : I := #t
+ for i in n..1 by -1 repeat
+ t(j) := s(i)
+ if (n - i + 1) rem sp = 0 then j := j - 1
+ j := j - 1
+ t
+
+ split_digits(s : S, top_digit : I) : List(S) ==
+ s2 := s
+ n1 := #s
+ s1 :=
+ top_digit > 0 =>
+ not(top_digit < n1) =>
+ s2 := "0"
+ concat(s, new((top_digit - n1)::N, zero))
+ s2 := s(top_digit + 1..n1)
+ s(1..top_digit)
+ "0"
+ if top_digit < 0 then
+ s2 := concat(new((-top_digit)::N, zero), s2)
+ [s1, s2]
+
+ fixed_default(f : F, d : I, mr : mode_Rec) : S ==
+ df := convert_to_decimal(f, d)
+ (dm, de) := df
+ s := convert(dm)@S
+ n0 := #s
+ top_digit := de + n0
+ s := rightTrim(s, zero)
+ (s1, s2) := split_digits(s, top_digit)
+ s1 := concat(insert_separators_from_right(s1, mr), ".")
+ s2 := insert_separators_from_left(s2, mr)
+ concat(s1, s2)
+
+ fixed_mr(f : F, d : I, mr : mode_Rec) : S ==
+ m := f.mantissa
+ e := f.exponent
+ m = 0 =>
+ mr.default? => "0.0"
+ concat("0", concat(".",
+ insert_separators_from_left(new(d::N, zero), mr)))
+ m < 0 => concat("-", fixed_mr([-m, e], d, mr))
+ even?(m) and e < 0 =>
+ lb := length(m)
+ if -e < lb then
+ lb := -e
+ i : I := 0
+ while i <= lb and not(bit?(m, i)) repeat
+ i := i + 1
+ m := shift(m, -i)
+ fixed_mr([m, e + i], d, mr)
+ not(e < 0) =>
+ nm :=
+ e = 0 => m
+ shift(m::N, e)
+ s1 := insert_separators_from_right(convert(nm)@S, mr)
+ s2 : S :=
+ mr.default? => ".0"
+ concat(".", insert_separators_from_left(new(d::N, zero), mr))
+ concat(s1, s2)
+ mr.default? => fixed_default(f, d, mr)
+ e := -e
+ top_bit := length(m) - e
+ (top_bit + 2)*12655 < -d*42039 =>
+ fixed_mr([0, 0], d, mr)
+ m := m*5^(e::N)
+ if e > d then
+ d10 := 10^((e - d)::N)
+ (q, r) := divide(m, d10)
+ d10h := d10 quo 2
+ m :=
+ r < d10h => q
+ r > d10h => q + 1
+ even?(q) => q
+ q + 1
+ e := d
+ s := convert(m)@S
+ n0 := #s
+ top_digit := n0 - e
+ (s1, s2) := split_digits(s, top_digit)
+ if d = 0 then
+ s2 := ""
+ s1 := concat(insert_separators_from_right(s1, mr), ".")
+ if e < d then
+ s2 := concat(s2, new((d - e)::N, zero))
+ s2 := insert_separators_from_left(s2, mr)
+ concat(s1, s2)
+
+ floating_mr(f : F, d : I, mr : mode_Rec) : S ==
+ f.mantissa = 0 => "0.0"
+ f.mantissa < 0 => concat("-", floating_mr([-f.mantissa,
+ f.exponent], d, mr))
+ es : S := if zero?(mr.spacing) then "E" else " E "
+ zero? f.exponent =>
+ s := convert(f.mantissa)@S
+ concat(["0.", insert_separators_from_left(s, mr), es,
+ convert(#s)@S])
+ (m, e) := convert_to_decimal(f, d)
+ s := convert(m)@S
+ n0 := #s
+ if mr.default? then
+ s := rightTrim(s, zero)
+ s := concat("0.", insert_separators_from_left(s, mr))
+ concat(s, concat(es, convert(e + n0)@S))
+
+ general_mr(f : F, d : I, mr : mode_Rec) : S ==
+ (m, e) := f
+ m = 0 => "0.0"
+ m < 0 =>
+ concat("-", general_mr([-m, e], d, mr))
+ e = 0 =>
+ s := convert(m)@S
+ top_digit := #s
+ mr.default? or top_digit < d + 3 =>
+ concat(insert_separators_from_right(s, mr), ".0")
+ s1 := concat("0.", insert_separators_from_left(s(1..d), mr))
+ es : S := if zero?(mr.spacing) then "E" else " E "
+ concat(s1, concat(es, convert(top_digit - d)@S))
+ (m, e) := convert_to_decimal(f, d)
+ s := convert(m)@S
+ n0 := #s
+ top_digit := e + n0
+ if mr.default? then
+ s := rightTrim(s, zero)
+ if -6 < top_digit and top_digit < d + 2 then
+ (s1, s2) := split_digits(s, top_digit)
+ s1 := concat(insert_separators_from_right(s1, mr), ".")
+ s2 := insert_separators_from_left(s2, mr)
+ concat(s1, s2)
+ else
+ es : S := if zero?(mr.spacing) then "E" else " E "
+ s1 := concat("0.", insert_separators_from_left(s, mr))
+ s2 := concat(es, convert(top_digit)@S)
+ concat(s1, s2)
diff -Naur '--exclude=.git' trunk/src/algebra/flopak.spad.bb trunk.p1/src/algebra/flopak.spad.bb
diff -Naur '--exclude=.git' trunk/src/algebra/Makefile.in trunk.p1/src/algebra/Makefile.in
--- trunk/src/algebra/Makefile.in 2024-07-01 02:50:56.455724693 +0200
+++ trunk.p1/src/algebra/Makefile.in 2024-08-15 23:26:18.590651315 +0200
@@ -42,7 +42,7 @@
eigen elemntry elfuts equation1 error \
evalut expexpan export3D expps expr2ups expr extred \
facutil fdalg ffact ffcat ffdoms fffg \
- files float fmod \
+ files float flopak fmod \
fmt fmt1d fmt2d fmtjfricas fmtlatex fmtmathjax \
fname fnla \
fortcat fortmac fortout fortpak fortran forttyp \
@@ -147,7 +147,7 @@
EXPR2UPS EXPRODE EXPR EXPRTUBE EXPUPXS EXTRED FACTCAT FACTFUNC \
FACUTIL FAKEPOL FAMR FARRAY \
FAXF FBICPO FCDCPO FCOMP FCPAK1 FCTOOL FC FDALG FDIV2 FDIVCAT \
- FDIV FDCPO FELFUN FEVALAB FEXPR FFCAT2 FFCAT FFCGP \
+ FDIV FDCPO FELEM FELFUN FEVALAB FEXPR FFCAT2 FFCAT FFCGP \
FFCG FFCGX FFF FFHOM FFIELDC FFINTBAS \
FFNBP FFNB FFNBX FFPOLY2 FFPOLY FFP \
FFSLPE FF FFX FGLMICPK FGROUP FGRPH FIELD \
@@ -155,11 +155,11 @@
SPADLIST3=\
FLAGG FLALG FLASORT FLINEXP FLIOUFUN FLOATCP FLOATRP \
- FLOAT FMAGMA FMCAT FMC FMFUN FMOEBF FMONOID \
+ FLOAT FLOOPS FMAGMA FMCAT FMC FMFUN FMOEBF FMONOID \
FM FM2 FMCF2 FMTC \
FMT1D FMT2D FMTCAT FMTLATEX FMTMJAX FMTOUT \
FNAME FNCAT FNGRPH FNLA FOP \
- FORDER FORMAT FORMCAT FORTCAT FORTFN FORTFORM FORTRAN \
+ FORDER FORMAT FORMCAT FORTCAT FORTFN FORTFORM FORTRAN FOUT \
FPARFRAC FPATMAB FPC FPOSET FPS FR2 \
FRAC2 FRAC FRAMALG FRETRCT FRIDEAL2 FRIDEAL \
FRIMOD FRMOD FRNAAF2 FRNAALG FR FRUTIL FS2EXPXP \
diff -Naur '--exclude=.git' trunk/src/algebra/sf.spad trunk.p1/src/algebra/sf.spad
--- trunk/src/algebra/sf.spad 2024-07-23 01:48:02.241534107 +0200
+++ trunk.p1/src/algebra/sf.spad 2024-08-16 18:30:02.517450645 +0200
@@ -160,13 +160,11 @@
digits() == max(1, 4004 * (bits()-1) quo 13301)::PositiveInteger
toString x == convert(x)@String
+
toString(x : %, n : NonNegativeInteger) : String ==
- x1 := round x
- x0 := abs(x-x1)
- tenn := float(1, n::Integer, 10)
- x2 := round(tenn*x0)
- res:List String := [convert(retract(x1))@String,".",convert(retract(x2))@String]
- concat res
+ fixed_mr([mantissa(x), exponent(x)], n,
+ [false, "fixed", 0])$FloatingPointConvertion
+
)if false
\section{domain DFLOAT DoubleFloat}
