--* From postmaster%watson.vnet.ibm.com@yktvmv.watson.ibm.com Fri Jul 15 09:54:43 1994 --* Received: from yktvmv-ob.watson.ibm.com by asharp.watson.ibm.com (AIX 3.2/UCB 5.64/930311) --* id AA21494; Fri, 15 Jul 1994 09:54:43 -0400 --* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3) --* with BSMTP id 4845; Fri, 15 Jul 94 09:54:45 EDT --* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0" --* id --* for asbugs@watson; Fri, 15 Jul 94 09:54:45 -0400 --* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0" --* id 2921; Fri, 15 Jul 1994 09:54:44 EDT --* Received: from spadserv.watson.ibm.com by yktvmv.watson.ibm.com --* (IBM VM SMTP V2R3) with TCP; Fri, 15 Jul 94 09:54:43 EDT --* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524) --* id AA12439; Fri, 15 Jul 1994 09:54:53 -0400 --* Date: Fri, 15 Jul 1994 09:54:53 -0400 --* From: rmc@spadserv.watson.ibm.com --* X-External-Networks: yes --* Message-Id: <9407151354.AA12439@spadserv.watson.ibm.com> --* To: asbugs@watson.ibm.com --* Subject: [1] seg fault -gloop [rratio.as][36.0] --@ Fixed by: PI Mon Aug 1 18:40:55 EDT 1994 --@ Tested by: none --@ Summary: General fixes interpreter & libraries. --+ asharp -gloop --+ %1 >> #include "aslib" --+ %2 >> #library rr "rratio.aso" --+ %3 >> SI ==> SingleInteger; --+ %4 >> R ==> RoundedRatio(SI, 1000); --+ %5 >> import from SI, rr, R; --+ %6 >> x:R := 1@R; --+ #1 (Error) Program fault (segmentation violation). -- -- First new type: Rounded Rationals. RMC 25/5/1994 -- -- Basic idea: use continued fractions to keep rationals short. -- -- Lessons. -- 6) Compile with asharp -Fo RoundedRatio.as for later use. #include "aslib.as" -- The parameter MAXDENOM ensures that denominators -- are not too big. Currently numerators are allowed to -- be MUCH larger---probably this is a mistake. RoundedRatio( I: IntegerNumberSystem, MAXDENOM: I): Join(Field, OrderedRing) with { -- the following extra operations: +: (I, %) -> %; -- These eight operations allow easy use of -: (I, %) -> %; -- embedded integers. *: (I, %) -> %; /: (I, %) -> %; +: (%, I) -> %; -: (%, I) -> %; *: (%, I) -> %; /: (%, I) -> %; /: (%, %) -> %; \: (%, %) -> %; /: (I, I) -> %; inv: % -> %; numer: % -> I; denom: % -> I; -- coerce: I -> %; quo: (%, %) -> %; -- needed if we want to rem: (%, %) -> %; -- pretend we're a EuclideanDomain gcd: (%, %) -> %; divide: (%, %) -> (%, %); sqrt: (%) -> %; root3: () -> %; export from I; } == add { Rep ==> Record(numer: I, denom: I); import from Rep; default a, b: %; ratio(n: I, d: I): % == per [n, d]; -- Use continued fractions to keep denominators small. round( x: Rep) : % == { local a0,a1,q,r,p0,p1,p2,q0,q1,q2,s : I; a0 := abs x.numer; a1 := abs x.denom; s := sign(x.numer)*sign(x.denom); (q,r) := divide(a0,a1); a0 := a1; a1 := r; p0 := 1; p1 := q; q0 := 0; q1 := 1; while r ~= 0 repeat { (q,r) := divide(a0,a1); a0 := a1; a1 := r; p2 := q*p1 + p0; q2 := q*q1 + q0; q2 > MAXDENOM => break; p0 := p1; p1 := p2; q0 := q1; q1 := q2; } per [s* abs p1, abs q1] } -- Public Part -- -- 1) Needed for all types (from BasicType) (stream: TextWriter) << (z: %) : TextWriter == stream << numer z << "/" << denom z ; -- 2) Operations necessary to qualify as an OrderedRing -- (though we don't guarantee the ring axioms hold) (a: %) = (b: %): Boolean == numer a = numer b and denom a = denom b; (a: %) ~= (b: %): Boolean == ~(a = b); (a: %) < (b: %): Boolean == denom b * numer a < denom a * numer b; (a: %) <= (b: %): Boolean == denom b * numer a <= denom a * numer b; (a: %) > (b: %): Boolean == denom b * numer a > denom a * numer b; (a: %) >= (b: %): Boolean == denom b * numer a >= denom a * numer b; max(a: %, b: %): % == if a > b then a else b; min(a: %, b: %): % == if a < b then a else b; sign (a: %): % == ratio(sign numer a, 1); abs (a: %): % == if negative? a then -a else a; zero? (a: %): Boolean == zero? numer a; negative? (a: %): Boolean == negative? numer a; positive? (a: %): Boolean == positive? numer a; -- 3) Operations necessary to qualify as an ArithmeticSystem -- Default coercion used for coerce(n: Integer) : % -- unless I=Integer. Decide that at compile time. -- coerce(n: I) :% == ratio(n, 1); coerce(n: Integer) :% == ratio(n::I, 1); -- craps out if I=SingleInteger coerce(n: SingleInteger) :% == ratio(n::I, 1); 0: % == ratio(0,1); 1: % == ratio(1,1); +(a: %): % == a; -(a: %): % == ratio(-numer a, denom a); (a: %) + (b: %): % == { z := [numer a * denom b + numer b * denom a, denom a * denom b ]; round z; } (a: %) - (b: %): % == a + (-b); (a: %) * (b: %): % == round( [ numer a * numer b, denom a * denom b] ); -- Use binary powering, with rounding at each step. (alpha: %) ^ (n: Integer): % == { s: % := alpha; m: Integer := abs n; len: SingleInteger := length m; c: % := 1; { m=0 => c := 1; m=1 => c := s; -- otherwise (loop invariant `answer' = c*s^v, v=m initially) { for i:SingleInteger in 1..(len-1) repeat { if bit(m,i-1) then c := c*s; -- v <== v - 1 s := s*s; -- v <== v/2 } c := c*s; -- bit(m, length(m) - 1) = 1 always. } } if n >= 0 then c else round [denom c, numer c] } -- 4) Extra operations. numer(a: %): I == rep(a).numer; denom(a: %): I == rep(a).denom; (a: %) quo (b: %): % == a/b; (a: %) rem (b: %): % == 0; gcd(a:%, b:%): % == 1; divide(a:%, b:%): (%, %) == (a/b, 0); -- We experiment here with allowing denormalized numbers to appear -- temporarily as a result of "inv". inv(a: %): % == per( [denom a, numer a] ); -- The results of other operations, however, are rounded. (n: I) / (d: I): % == round( [n, d] ) ; (a: %) / (b: %): % == round( [ denom b * numer a, numer b * denom a] ); (a: %) \ (b: %): % == round( [ denom a * numer b, denom b * numer a] ); (n: I) + (b: %): % == ratio(n, 1) + b; (a: %) + (n: I): % == a + ratio(n, 1); (n: I) - (b: %): % == ratio(n, 1) - b; (a: %) - (n: I): % == a - ratio(n, 1); (n: I) * (b: %): % == ratio(n, 1) * b; (a: %) * (n: I): % == a * ratio(n, 1); (n: I) / (b: %): % == ratio(n, 1) / b; (a: %) / (n: I): % == a / ratio(n, 1); root3() : % == { -- sqrt(3) = 1 + [1, 2, 1, 2, 1, 2, ... ] local p0,p1,p2,q0,q1,q2: I; p0 := 1; p1 := 1; q0 := 0; q1 := 1; while true repeat { p2 := p1 + p0; q2 := q1 + q0; q2 > MAXDENOM => break; p0 := p1; p1 := p2; q0 := q1; q1 := q2; p2 := 2*p1 + p0; q2 := 2*q1 + q0; q2 > MAXDENOM => break; p0 := p1; p1 := p2; q0 := q1; q1 := q2; } per [p1, q1] } sqrt(x:%) :% == { x < 0 => error "negative square root"; guess: % := 1; residual: % := guess * guess - x; while abs residual > 100/(MAXDENOM * MAXDENOM) repeat { guess := guess - residual/(2 * guess); residual := guess * guess - x; } guess } }