diff --git a/lang/m2/libm2/.distr b/lang/m2/libm2/.distr
index a91543de3..5c9d5afa5 100644
--- a/lang/m2/libm2/.distr
+++ b/lang/m2/libm2/.distr
@@ -6,6 +6,7 @@ Conversion.def
 FIFFEF.def
 InOut.def
 Makefile
+Mathlib.def
 MathLib0.def
 Processes.def
 RealInOut.def
diff --git a/lang/m2/libm2/FIFFEF.def b/lang/m2/libm2/FIFFEF.def
index 922332d12..ea49df35a 100644
--- a/lang/m2/libm2/FIFFEF.def
+++ b/lang/m2/libm2/FIFFEF.def
@@ -1,11 +1,12 @@
+(*$Foreign*)
 DEFINITION MODULE FIFFEF;
 
-	PROCEDURE FIF(arg1, arg2: REAL; VAR intres: REAL) : REAL;
+	PROCEDURE FIF(arg1, arg2: LONGREAL; VAR intres: LONGREAL) : LONGREAL;
 	(* multiplies arg1 and arg2, and returns the integer part of the
 	   result in "intres" and the fraction part as the function result.
 	*)
 
-	PROCEDURE FEF(arg: REAL; VAR exp: INTEGER) : REAL;
+	PROCEDURE FEF(arg: LONGREAL; VAR exp: INTEGER) : LONGREAL;
 	(* splits "arg" in mantissa and a base-2 exponent.
 	   The mantissa is returned, and the exponent is left in "exp".
 	*)
diff --git a/lang/m2/libm2/FIFFEF.e b/lang/m2/libm2/FIFFEF.e
index b7588e6a9..849cfc2b7 100644
--- a/lang/m2/libm2/FIFFEF.e
+++ b/lang/m2/libm2/FIFFEF.e
@@ -2,50 +2,45 @@
  mes 2,EM_WSIZE,EM_PSIZE
 
 #define ARG1    0
-#define ARG2    EM_FSIZE
-#define IRES    2*EM_FSIZE
+#define ARG2    EM_DSIZE
+#define IRES    2*EM_DSIZE
 
-; FIFFEF_FIF is called with three parameters:
+; FIF is called with three parameters:
 ;       - address of integer part result (IRES)
 ;       - float two (ARG2)
 ;       - float one (ARG1)
-; and returns an EM_FSIZE-byte floating point number
+; and returns an EM_DSIZE-byte floating point number
 ; Definition:
-;	PROCEDURE FIF(ARG1, ARG2: REAL; VAR IRES: REAL) : REAL;
+;	PROCEDURE FIF(ARG1, ARG2: LONGREAL; VAR IRES: LONGREAL) : LONGREAL;
 
- exp $FIFFEF_FIF
- pro $FIFFEF_FIF,0
+ exp $FIF
+ pro $FIF,0
  lal 0
- loi 2*EM_FSIZE
- fif EM_FSIZE
+ loi 2*EM_DSIZE
+ fif EM_DSIZE
  lal IRES
  loi EM_PSIZE
- sti EM_FSIZE
- ret EM_FSIZE
+ sti EM_DSIZE
+ ret EM_DSIZE
  end ?
 
 #define FARG    0
-#define ERES    EM_FSIZE
+#define ERES    EM_DSIZE
 
-; FIFFEF_FEF is called with two parameters:
+; FEF is called with two parameters:
 ;       - address of base 2 exponent result (ERES)
 ;       - floating point number to be split (FARG)
-; and returns an EM_FSIZE-byte floating point number (the mantissa)
+; and returns an EM_DSIZE-byte floating point number (the mantissa)
 ; Definition:
-;	PROCEDURE FEF(FARG: REAL; VAR ERES: integer): REAL;
+;	PROCEDURE FEF(FARG: LONGREAL; VAR ERES: integer): LONGREAL;
 
- exp $FIFFEF_FEF
- pro $FIFFEF_FEF,0
+ exp $FEF
+ pro $FEF,0
  lal FARG
- loi EM_FSIZE
- fef EM_FSIZE
+ loi EM_DSIZE
+ fef EM_DSIZE
  lal ERES
  loi EM_PSIZE
  sti EM_WSIZE
- ret EM_FSIZE
- end ?
-
- exp $FIFFEF
- pro $FIFFEF,0
- ret 0
+ ret EM_DSIZE
  end ?
diff --git a/lang/m2/libm2/LIST b/lang/m2/libm2/LIST
index 05d106755..9e2d8a06c 100644
--- a/lang/m2/libm2/LIST
+++ b/lang/m2/libm2/LIST
@@ -1,18 +1,19 @@
 tail_m2.a
+RealInOut.mod
 InOut.mod
 Terminal.mod
 TTY.mod
 ASCII.mod
-FIFFEF.e
 MathLib0.mod
+Mathlib.mod
 Processes.mod
 RealConver.mod
-RealInOut.mod
 Storage.mod
 Conversion.mod
 Semaphores.mod
 random.mod
 Strings.mod
+FIFFEF.e
 Arguments.c
 catch.c
 hol0.e
diff --git a/lang/m2/libm2/Makefile b/lang/m2/libm2/Makefile
index 1b2dbf2f1..77b7bc94b 100644
--- a/lang/m2/libm2/Makefile
+++ b/lang/m2/libm2/Makefile
@@ -4,7 +4,8 @@ DEFDIR = $(HOME)/lib/m2
 SOURCES =	ASCII.def FIFFEF.def MathLib0.def Processes.def \
 		RealInOut.def Storage.def Arguments.def Conversion.def \
 		random.def Semaphores.def Unix.def RealConver.def \
-		Strings.def InOut.def Terminal.def TTY.def
+		Strings.def InOut.def Terminal.def TTY.def \
+		Mathlib.def
 
 all:
 
diff --git a/lang/m2/libm2/MathLib0.def b/lang/m2/libm2/MathLib0.def
index cbae3a382..985a6da44 100644
--- a/lang/m2/libm2/MathLib0.def
+++ b/lang/m2/libm2/MathLib0.def
@@ -1,4 +1,8 @@
 DEFINITION MODULE MathLib0;
+(*
+	Exists for compatibility.
+	A more elaborate math lib can be found in Mathlib.def
+*)
 
 	PROCEDURE sqrt(x : REAL) : REAL;
 
diff --git a/lang/m2/libm2/MathLib0.mod b/lang/m2/libm2/MathLib0.mod
index 93c8a7864..d48833c47 100644
--- a/lang/m2/libm2/MathLib0.mod
+++ b/lang/m2/libm2/MathLib0.mod
@@ -1,326 +1,35 @@
 IMPLEMENTATION MODULE MathLib0;
-(*	Rewritten in Modula-2.
-	The originals came from the Pascal runtime library.
-*)
 
-FROM FIFFEF	IMPORT	FIF, FEF;
-
-CONST
-	HUGE =	1.701411733192644270E38;
-
-PROCEDURE sinus(arg: REAL; quad: INTEGER): REAL;
-
-(*	Coefficients for sin/cos are #3370 from Hart & Cheney (18.80D).
-*)
-CONST
-	twoopi	= 0.63661977236758134308;
-	p0	= 0.1357884097877375669092680E8;
-	p1	= -0.4942908100902844161158627E7;
-	p2	= 0.4401030535375266501944918E6;
-	p3	= -0.1384727249982452873054457E5;
-	p4	= 0.1459688406665768722226959E3;
-	q0	= 0.8644558652922534429915149E7;
-	q1	= 0.4081792252343299749395779E6;
-	q2	= 0.9463096101538208180571257E4;
-	q3	= 0.1326534908786136358911494E3;
-VAR
-	e, f: REAL;
-	ysq: REAL;
-	x,y: REAL;
-	k: INTEGER;
-	temp1, temp2: REAL;
-BEGIN
-	x := arg;
-	IF x < 0.0 THEN
-		x := -x;
-		quad := quad + 2;
-	END;
-	x := x*twoopi;	(*underflow?*)
-	IF x>32764.0 THEN
-		y := FIF(x, 10.0, e);
-		e := e + FLOAT(quad);
-		temp1 := FIF(0.25, e, f);
-		quad := TRUNC(e - 4.0*f);
-	ELSE
-		k := TRUNC(x);
-		y := x - FLOAT(k);
-		quad := (quad + k) MOD 4;
-	END;
-	IF ODD(quad) THEN
-		y := 1.0-y;
-	END;
-	IF quad > 1 THEN
-		y := -y;
-	END;
-
-	ysq := y*y;
-	temp1 := ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
-	temp2 := ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
-	RETURN temp1/temp2;
-END sinus;
+  IMPORT Mathlib;
 
 PROCEDURE cos(arg: REAL): REAL;
 BEGIN
-	IF arg < 0.0 THEN
-		arg := -arg;
-	END;
-	RETURN sinus(arg, 1);
+	RETURN Mathlib.cos(arg);
 END cos;
 
 PROCEDURE sin(arg: REAL): REAL;
 BEGIN
-	RETURN sinus(arg, 0);
+	RETURN Mathlib.sin(arg);
 END sin;
 
-(*
-	floating-point arctangent
-
-	arctan returns the value of the arctangent of its
-	argument in the range [-pi/2,pi/2].
-
-	coefficients are #5077 from Hart & Cheney. (19.56D)
-*)
-
-CONST
-	sq2p1	= 2.414213562373095048802E0;
-	sq2m1	= 0.414213562373095048802E0;
-	pio2	= 1.570796326794896619231E0;
-	pio4	= 0.785398163397448309615E0;
-	p4	= 0.161536412982230228262E2;
-	p3	= 0.26842548195503973794141E3;
-	p2	= 0.11530293515404850115428136E4;
-	p1	= 0.178040631643319697105464587E4;
-	p0	= 0.89678597403663861959987488E3;
-	q4	= 0.5895697050844462222791E2;
-	q3	= 0.536265374031215315104235E3;
-	q2	= 0.16667838148816337184521798E4;
-	q1	= 0.207933497444540981287275926E4;
-	q0	= 0.89678597403663861962481162E3;
-
-(*
-	xatan evaluates a series valid in the
-	range [-0.414...,+0.414...].
-*)
-
-PROCEDURE xatan(arg: REAL) : REAL;
-VAR
-	argsq, value: REAL;
-BEGIN
-	argsq := arg*arg;
-	value := ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0);
-	value := value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0);
-	RETURN value*arg;
-END xatan;
-
-PROCEDURE satan(arg: REAL): REAL;
-BEGIN
-	IF arg < sq2m1 THEN
-		RETURN xatan(arg);
-	ELSIF arg > sq2p1 THEN
-		RETURN pio2 - xatan(1.0/arg);
-	ELSE
-		RETURN pio4 + xatan((arg-1.0)/(arg+1.0));
-	END;
-END satan;
-
-(*
-	atan makes its argument positive and
-	calls the inner routine satan.
-*)
-
 PROCEDURE arctan(arg: REAL): REAL;
 BEGIN
-	IF arg>0.0 THEN
-		RETURN satan(arg);
-	ELSE
-		RETURN -satan(-arg);
-	END;
+	RETURN Mathlib.arctan(arg);
 END arctan;
 
-(*
-	sqrt returns the square root of its floating
-	point argument. Newton's method.
-*)
-
 PROCEDURE sqrt(arg: REAL): REAL;
-VAR
-	x, temp: REAL;
-	exp, i: INTEGER;
 BEGIN
-	IF arg <= 0.0 THEN
-		IF arg < 0.0 THEN
-			(* ??? *)
-			;
-		END;
-		RETURN 0.0;
-	END;
-	x := FEF(arg,exp);
-	(*
-	 * NOTE
-	 * this wont work on 1's comp
-	 *)
-	IF ODD(exp) THEN
-		x := 2.0 * x;
-		DEC(exp);
-	END;
-	temp := 0.5*(1.0 + x);
-
-	WHILE exp > 28 DO
-		temp := temp * 16384.0;
-		exp := exp - 28;
-	END;
-	WHILE exp < -28 DO
-		temp := temp / 16384.0;
-		exp := exp + 28;
-	END;
-	WHILE exp >= 2 DO
-		temp := temp * 2.0;
-		exp := exp - 2;
-	END;
-	WHILE exp <= -2 DO
-		temp := temp / 2.0;
-		exp := exp + 2;
-	END;
-	FOR i := 0 TO 4 DO
-		temp := 0.5*(temp + arg/temp);
-	END;
-	RETURN temp;
+	RETURN Mathlib.sqrt(arg);
 END sqrt;
 
-(*
-	ln returns the natural logarithm of its floating
-	point argument.
-
-	The coefficients are #2705 from Hart & Cheney. (19.38D)
-*)
 PROCEDURE ln(arg: REAL): REAL;
-CONST
-	log2	= 0.693147180559945309E0;
-	sqrto2	= 0.707106781186547524E0;
-	p0	= -0.240139179559210510E2;
-	p1	= 0.309572928215376501E2;
-	p2	= -0.963769093368686593E1;
-	p3	= 0.421087371217979714E0;
-	q0	= -0.120069589779605255E2;
-	q1	= 0.194809660700889731E2;
-	q2	= -0.891110902798312337E1;
-VAR
-	x,z, zsq, temp: REAL;
-	exp: INTEGER;
 BEGIN
-	IF arg <= 0.0 THEN
-		(* ??? *)
-		RETURN -HUGE;
-	END;
-	x := FEF(arg,exp);
-	IF x<sqrto2 THEN
-		x := x + x;
-		DEC(exp);
-	END;
-
-	z := (x-1.0)/(x+1.0);
-	zsq := z*z;
-
-	temp := ((p3*zsq + p2)*zsq + p1)*zsq + p0;
-	temp := temp/(((zsq + q2)*zsq + q1)*zsq + q0);
-	temp := temp*z + FLOAT(exp)*log2;
-	RETURN temp;
+	RETURN Mathlib.ln(arg);
 END ln;
 
-(*
-	exp returns the exponential function of its
-	floating-point argument.
-
-	The coefficients are #1069 from Hart and Cheney. (22.35D)
-*)
-
-PROCEDURE floor(d: REAL): REAL;
-BEGIN
-	IF d < 0.0 THEN
-		d := -d;
-		IF FIF(d, 1.0, d) # 0.0 THEN
-			d := d + 1.0;
-		END;
-		d := -d;
-	ELSE
-		IF FIF(d, 1.0, d) # 0.0 THEN
-			(* just ignore result of FIF *)
-			;
-		END;
-	END;
-	RETURN d;
-END floor;
-
-PROCEDURE ldexp(fr: REAL; exp: INTEGER): REAL;
-VAR
-	neg,i: INTEGER;
-BEGIN
-	neg := 1;
-	IF fr < 0.0 THEN
-		fr := -fr;
-		neg := -1;
-	END;
-	fr := FEF(fr, i);
-	exp := exp + i;
-	IF exp > 127 THEN
-		(* Too large. ??? *)
-		RETURN FLOAT(neg) * HUGE;
-	END;
-	IF exp < -127 THEN
-		RETURN 0.0;
-	END;
-	WHILE exp > 14 DO
-		fr := fr * 16384.0;
-		exp := exp - 14;
-	END;
-	WHILE exp < -14 DO
-		fr := fr / 16384.0;
-		exp := exp + 14;
-	END;
-	WHILE exp > 0 DO
-		fr := fr + fr;
-		DEC(exp);
-	END;
-	WHILE exp < 0 DO
-		fr := fr / 2.0;
-		INC(exp);
-	END;
-	RETURN FLOAT(neg) * fr;
-END ldexp;
-
 PROCEDURE exp(arg: REAL): REAL;
-CONST
-	p0	= 0.2080384346694663001443843411E7;
-	p1	= 0.3028697169744036299076048876E5;
-	p2	= 0.6061485330061080841615584556E2;
-	q0	= 0.6002720360238832528230907598E7;
-	q1	= 0.3277251518082914423057964422E6;
-	q2	= 0.1749287689093076403844945335E4;
-	log2e	= 1.4426950408889634073599247;
-	sqrt2	= 1.4142135623730950488016887;
-	maxf	= 10000.0;
-VAR
-	fract: REAL;
-	temp1, temp2, xsq: REAL;
-	ent: INTEGER;
 BEGIN
-	IF arg = 0.0 THEN
-		RETURN 1.0;
-	END;
-	IF arg < -maxf THEN
-		RETURN 0.0;
-	END;
-	IF arg > maxf THEN
-		(* result too large ??? *)
-		RETURN HUGE;
-	END;
-	arg := arg * log2e;
-	ent := TRUNC(floor(arg));
-	fract := (arg-FLOAT(ent)) - 0.5;
-	xsq := fract*fract;
-	temp1 := ((p2*xsq+p1)*xsq+p0)*fract;
-	temp2 := ((xsq+q2)*xsq+q1)*xsq + q0;
-	RETURN ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent);
+	RETURN Mathlib.exp(arg);
 END exp;
 
 PROCEDURE entier(x: REAL): INTEGER;
diff --git a/lang/m2/libm2/Mathlib.def b/lang/m2/libm2/Mathlib.def
new file mode 100644
index 000000000..979d1b696
--- /dev/null
+++ b/lang/m2/libm2/Mathlib.def
@@ -0,0 +1,66 @@
+DEFINITION MODULE Mathlib;
+
+  (* Some mathematical constants: *)
+
+  CONST
+	(* From:	Handbook of Mathematical Functions
+			Edited by M. Abramowitz and I.A. Stegun
+			National Bureau of Standards
+			Applied Mathematics Series 55
+	*)
+
+	pi	= 3.141592653589793238462643;
+	twicepi	= 6.283185307179586476925286;
+	halfpi	= 1.570796326794896619231322;
+	quartpi	= 0.785398163397448309615661;
+	e	= 2.718281828459045235360287;
+	ln2	= 0.693147180559945309417232;
+	ln10	= 2.302585092994045684017992;
+
+  (* basic functions *)
+
+  PROCEDURE pow(x: REAL; i: INTEGER): REAL;
+
+  PROCEDURE sqrt(x: REAL): REAL;
+
+  PROCEDURE exp(x: REAL): REAL;
+
+  PROCEDURE ln(x: REAL): REAL;	(* natural log *)
+
+  PROCEDURE log(x: REAL): REAL;	(* log with base 10 *)
+
+  (* trigonometric functions; arguments in radians *)
+
+  PROCEDURE sin(x: REAL): REAL;
+
+  PROCEDURE cos(x: REAL): REAL;
+
+  PROCEDURE tan(x: REAL): REAL;
+
+  PROCEDURE arcsin(x: REAL): REAL;
+
+  PROCEDURE arccos(x: REAL): REAL;
+
+  PROCEDURE arctan(x: REAL): REAL;
+
+  (* hyperbolic functions *)
+
+  PROCEDURE sinh(x: REAL): REAL;
+
+  PROCEDURE cosh(x: REAL): REAL;
+
+  PROCEDURE tanh(x: REAL): REAL;
+
+  PROCEDURE arcsinh(x: REAL): REAL;
+
+  PROCEDURE arccosh(x: REAL): REAL;
+
+  PROCEDURE arctanh(x: REAL): REAL;
+
+  (* conversions *)
+
+  PROCEDURE RadianToDegree(x: REAL): REAL;
+
+  PROCEDURE DegreeToRadian(x: REAL): REAL;
+
+END Mathlib.
diff --git a/lang/m2/libm2/Mathlib.mod b/lang/m2/libm2/Mathlib.mod
new file mode 100644
index 000000000..e1268bb32
--- /dev/null
+++ b/lang/m2/libm2/Mathlib.mod
@@ -0,0 +1,443 @@
+IMPLEMENTATION MODULE Mathlib;
+
+  FROM FIFFEF IMPORT FIF, FEF;
+
+	(* From:	Handbook of Mathematical Functions
+			Edited by M. Abramowitz and I.A. Stegun
+			National Bureau of Standards
+			Applied Mathematics Series 55
+	*)
+
+  CONST
+	OneRadianInDegrees	= 57.295779513082320876798155D;
+	OneDegreeInRadians	=  0.017453292519943295769237D;
+
+  (* basic functions *)
+
+  PROCEDURE pow(x: REAL; i: INTEGER): REAL;
+  BEGIN
+	RETURN SHORT(longpow(LONG(x), i));
+  END pow;
+
+  PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL;
+    VAR
+	val: LONGREAL;
+	ri: LONGREAL;
+  BEGIN
+	ri := FLOATD(i);
+	IF x < 0.0D THEN
+		val := longexp(longln(-x) * ri);
+		IF ODD(i) THEN RETURN -val;
+		ELSE RETURN val;
+		END;
+	ELSIF x = 0.0D THEN
+		RETURN 0.0D;
+	ELSE
+		RETURN longexp(longln(x) * ri);
+	END;
+  END longpow;
+
+  PROCEDURE sqrt(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longsqrt(LONG(x)));
+  END sqrt;
+
+  PROCEDURE longsqrt(x: LONGREAL): LONGREAL;
+    VAR
+	temp: LONGREAL;
+	exp, i: INTEGER;
+  BEGIN
+	IF x <= 0.0D THEN
+		IF x < 0.0D THEN
+			(* ??? *)
+			;
+		END;
+		RETURN 0.0D;
+	END;
+	temp := FEF(x,exp);
+	(*
+	 * NOTE
+	 * this wont work on 1's comp
+	 *)
+	IF ODD(exp) THEN
+		temp := 2.0D * temp;
+		DEC(exp);
+	END;
+	temp := 0.5D*(1.0D + temp);
+
+	WHILE exp > 28 DO
+		temp := temp * 16384.0D;
+		exp := exp - 28;
+	END;
+	WHILE exp < -28 DO
+		temp := temp / 16384.0D;
+		exp := exp + 28;
+	END;
+	WHILE exp >= 2 DO
+		temp := temp * 2.0D;
+		exp := exp - 2;
+	END;
+	WHILE exp <= -2 DO
+		temp := temp / 2.0D;
+		exp := exp + 2;
+	END;
+	FOR i := 0 TO 4 DO
+		temp := 0.5D*(temp + x/temp);
+	END;
+	RETURN temp;
+  END longsqrt;
+
+  PROCEDURE exp(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longexp(LONG(x)));
+  END exp;
+
+  PROCEDURE longexp(x: LONGREAL): LONGREAL;
+  (*
+   * n = floor(x / ln2), d = x / ln2 - n
+   * exp(x) = exp((x / ln2) * ln2) = exp((n + d) * ln2) =
+   * exp(n * ln2) * exp(d * ln2) = 2 ** n * exp(d * ln2)
+   *)
+    CONST
+	a1 = -0.9999999995D;
+	a2 =  0.4999999206D;
+	a3 = -0.1666653019D;
+	a4 =  0.0416573475D;
+	a5 = -0.0083013598D;
+	a6 =  0.0013298820D;
+	a7 = -0.0001413161D;
+    VAR
+	neg: BOOLEAN;
+	polval: LONGREAL;
+	n: LONGREAL;
+	n1 : INTEGER;
+  BEGIN
+	neg := x < 0.0D;
+	IF neg THEN
+		x := -x;
+	END;
+	x := FIF(x, 1.0D/LONG(ln2), n) * LONG(ln2);
+	polval := 1.0D /(1.0D + x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7)))))));
+	n1 := TRUNCD(n + 0.5D);
+	WHILE n1 >= 16 DO
+		polval := polval * 65536.0D;
+		n1 := n1 - 16;
+	END;
+	WHILE n1 > 0 DO
+		polval := polval * 2.0D;
+		DEC(n1);
+	END;
+	IF neg THEN RETURN 1.0D/polval; END;
+	RETURN polval;
+  END longexp;
+
+  PROCEDURE ln(x: REAL): REAL;	(* natural log *)
+  BEGIN
+	RETURN SHORT(longln(LONG(x)));
+  END ln;
+
+  PROCEDURE longln(x: LONGREAL): LONGREAL;	(* natural log *)
+    CONST
+	a1 =  0.9999964239D;
+	a2 = -0.4998741238D;
+	a3 =  0.3317990258D;
+	a4 = -0.2407338084D;
+	a5 =  0.1676540711D;
+	a6 = -0.0953293897D;
+	a7 =  0.0360884937D;
+	a8 = -0.0064535442D;
+    VAR
+	exp: INTEGER;
+	polval: LONGREAL;
+
+  BEGIN
+	IF x <= 0.0D THEN
+		(* ??? *)
+		RETURN 0.0D;
+	END;
+	x := FEF(x, exp);
+	WHILE x < 1.0D DO
+		x := x + x;
+		DEC(exp);
+	END;
+	x := x - 1.0D;
+	polval := x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*(a7+a8*x)))))));
+	RETURN polval + FLOATD(exp) * LONG(ln2);
+  END longln;
+
+  PROCEDURE log(x: REAL): REAL;	(* log with base 10 *)
+  BEGIN
+	RETURN SHORT(longlog(LONG(x)));
+  END log;
+
+  PROCEDURE longlog(x: LONGREAL): LONGREAL;	(* log with base 10 *)
+  BEGIN
+	RETURN longln(x)/LONG(ln10);
+  END longlog;
+
+  (* trigonometric functions; arguments in radians *)
+
+  PROCEDURE sin(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longsin(LONG(x)));
+  END sin;
+
+  PROCEDURE longsin(x: LONGREAL): LONGREAL;
+    CONST
+	a2  = -0.1666666664D;
+	a4  =  0.0083333315D;
+	a6  = -0.0001984090D;
+	a8  =  0.0000027526D;
+	a10 = -0.0000000239D;
+    VAR
+	xsqr: LONGREAL;
+	neg: BOOLEAN;
+  BEGIN
+	neg := FALSE;
+	IF x < 0.0D THEN
+		neg := TRUE;
+		x := -x;
+	END;
+	x := FIF(x, 1.0D / LONG(twicepi), (* dummy *) xsqr) * LONG(twicepi);
+	IF x >= LONG(pi) THEN
+		neg := NOT neg;
+		x := x - LONG(pi);
+	END;
+	IF x > LONG(halfpi) THEN
+		x := LONG(pi) - x;
+	END;
+	xsqr := x * x;
+	x := x * (1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*a10)))));
+	IF neg THEN RETURN -x; END;
+	RETURN x;
+  END longsin;
+
+  PROCEDURE cos(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longcos(LONG(x)));
+  END cos;
+
+  PROCEDURE longcos(x: LONGREAL): LONGREAL;
+    CONST
+	a2  = -0.4999999963D;
+	a4  =  0.0416666418D;
+	a6  = -0.0013888397D;
+	a8  =  0.0000247609D;
+	a10 = -0.0000002605D;
+    VAR
+	xsqr: LONGREAL;
+	neg: BOOLEAN;
+  BEGIN
+	neg := FALSE;
+	IF x < 0.0D THEN x := -x; END;
+	x := FIF(x, 1.0D / LONG(twicepi), (* dummy *) xsqr) * LONG(twicepi);
+	IF x >= LONG(pi) THEN
+		x := LONG(twicepi) - x;
+	END;
+	IF x > LONG(halfpi) THEN
+		neg := NOT neg;
+		x := LONG(pi) - x;
+	END;
+	xsqr := x * x;
+	x := 1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*a10))));
+	IF neg THEN RETURN -x; END;
+	RETURN x;
+  END longcos;
+
+  PROCEDURE tan(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longtan(LONG(x)));
+  END tan;
+
+  PROCEDURE longtan(x: LONGREAL): LONGREAL;
+    VAR cosinus: LONGREAL;
+  BEGIN
+	cosinus := longcos(x);
+	IF cosinus = 0.0D THEN
+		(* ??? *)
+		RETURN 0.0D;
+	END;
+	RETURN longsin(x)/cosinus;
+  END longtan;
+
+  PROCEDURE arcsin(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longarcsin(LONG(x)));
+  END arcsin;
+
+  PROCEDURE longarcsin(x: LONGREAL): LONGREAL;
+    CONST
+	a0 =  1.5707963050D;
+	a1 = -0.2145988016D;
+	a2 =  0.0889789874D;
+	a3 = -0.0501743046D;
+	a4 =  0.0308918810D;
+	a5 = -0.0170881256D;
+	a6 =  0.0066700901D;
+	a7 = -0.0012624911D;
+  BEGIN
+	IF x < 0.0D THEN x := -x; END;
+	IF x > 1.0D THEN
+		(* ??? *)
+		RETURN 0.0D;
+	END;
+	RETURN LONG(halfpi) -
+		longsqrt(1.0D - x)*(a0+x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7)))))));
+  END longarcsin;
+
+  PROCEDURE arccos(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longarccos(LONG(x)));
+  END arccos;
+
+  PROCEDURE longarccos(x: LONGREAL): LONGREAL;
+  BEGIN
+	RETURN LONG(halfpi) - longarcsin(x);
+  END longarccos;
+
+  PROCEDURE arctan(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longarctan(LONG(x)));
+  END arctan;
+
+  PROCEDURE longarctan(x: LONGREAL): LONGREAL;
+    CONST
+	a2  = -0.3333314528D;
+	a4  =  0.1999355085D;
+	a6  = -0.1420889944D;
+	a8  =  0.1065626393D;
+	a10 = -0.0752896400D;
+	a12 =  0.0429096318D;
+	a14 = -0.0161657367D;
+	a16 =  0.0028662257D;
+    VAR
+	xsqr: LONGREAL;
+	rev: BOOLEAN;
+	neg: BOOLEAN;
+  BEGIN
+	rev := FALSE;
+	neg := FALSE;
+	IF x < 0.0D THEN
+		neg := TRUE;
+		x := -x;
+	END;
+	IF x > 1.0D THEN
+		rev := TRUE;
+		x := 1.0D / x;
+	END;
+	xsqr := x * x;
+	x := x * (1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*(a10+xsqr*(a12+xsqr*(a14+xsqr*a16))))))));
+	IF rev THEN
+		x := LONG(quartpi) - x;
+	END;
+	IF neg THEN RETURN -x; END;
+	RETURN x;
+  END longarctan;
+
+  (* hyperbolic functions *)
+
+  PROCEDURE sinh(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longsinh(LONG(x)));
+  END sinh;
+
+  PROCEDURE longsinh(x: LONGREAL): LONGREAL;
+    VAR expx: LONGREAL;
+  BEGIN
+	expx := longexp(x);
+	RETURN (expx - 1.0D/expx)/2.0D;
+  END longsinh;
+
+  PROCEDURE cosh(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longcosh(LONG(x)));
+  END cosh;
+
+  PROCEDURE longcosh(x: LONGREAL): LONGREAL;
+    VAR expx: LONGREAL;
+  BEGIN
+	expx := longexp(x);
+	RETURN (expx + 1.0D/expx)/2.0D;
+  END longcosh;
+
+  PROCEDURE tanh(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longtanh(LONG(x)));
+  END tanh;
+
+  PROCEDURE longtanh(x: LONGREAL): LONGREAL;
+    VAR expx: LONGREAL;
+  BEGIN
+	expx := longexp(x);
+	RETURN (expx - 1.0D/expx) / (expx + 1.0D/expx);
+  END longtanh;
+
+  PROCEDURE arcsinh(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longarcsinh(LONG(x)));
+  END arcsinh;
+
+  PROCEDURE longarcsinh(x: LONGREAL): LONGREAL;
+    VAR neg: BOOLEAN;
+  BEGIN
+	neg := FALSE;
+	IF x < 0.0D THEN
+		neg := TRUE;
+		x := -x;
+	END;
+	x := longln(x + longsqrt(x*x+1.0D));
+	IF neg THEN RETURN -x; END;
+	RETURN x;
+  END longarcsinh;
+
+  PROCEDURE arccosh(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longarccosh(LONG(x)));
+  END arccosh;
+
+  PROCEDURE longarccosh(x: LONGREAL): LONGREAL;
+  BEGIN
+	IF x < 1.0D THEN
+		(* ??? *)
+		RETURN 0.0D;
+	END;
+	RETURN longln(x + longsqrt(x*x - 1.0D));
+  END longarccosh;
+
+  PROCEDURE arctanh(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longarctanh(LONG(x)));
+  END arctanh;
+
+  PROCEDURE longarctanh(x: LONGREAL): LONGREAL;
+  BEGIN
+	IF (x <= -1.0D) OR (x >= 1.0D) THEN
+		(* ??? *)
+		RETURN 0.0D;
+	END;
+	RETURN longln((1.0D + x)/(1.0D - x)) / 2.0D;
+  END longarctanh;
+
+  (* conversions *)
+
+  PROCEDURE RadianToDegree(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longRadianToDegree(LONG(x)));
+  END RadianToDegree;
+
+  PROCEDURE longRadianToDegree(x: LONGREAL): LONGREAL;
+  BEGIN
+	RETURN x * OneRadianInDegrees;
+  END longRadianToDegree;
+
+  PROCEDURE DegreeToRadian(x: REAL): REAL;
+  BEGIN
+	RETURN SHORT(longDegreeToRadian(LONG(x)));
+  END DegreeToRadian;
+
+  PROCEDURE longDegreeToRadian(x: LONGREAL): LONGREAL;
+  BEGIN
+	RETURN x * OneDegreeInRadians;
+  END longDegreeToRadian;
+
+END Mathlib.
diff --git a/lang/m2/libm2/RealConver.mod b/lang/m2/libm2/RealConver.mod
index 8d510e175..9af9b9205 100644
--- a/lang/m2/libm2/RealConver.mod
+++ b/lang/m2/libm2/RealConver.mod
@@ -2,21 +2,30 @@ IMPLEMENTATION MODULE RealConversions;
 
   FROM FIFFEF IMPORT FIF, FEF;
 
-  PROCEDURE RealToString(r: REAL;
+  PROCEDURE RealToString(arg: REAL;
+		width, digits: INTEGER;
+		VAR str: ARRAY OF CHAR;
+		VAR ok: BOOLEAN);
+  BEGIN
+	LongRealToString(LONG(arg), width, digits, str, ok);
+  END RealToString;
+
+  PROCEDURE LongRealToString(arg: LONGREAL;
 		width, digits: INTEGER;
 		VAR str: ARRAY OF CHAR;
 		VAR ok: BOOLEAN);
     VAR	pointpos: INTEGER;
 	i: CARDINAL;
 	ecvtflag: BOOLEAN;
-	intpart, fractpart: REAL;
+	r, intpart, fractpart: LONGREAL;
 	ind1, ind2 : CARDINAL;
 	sign: BOOLEAN;
 	tmp : CHAR;
 	ndigits: CARDINAL;
-	dummy, dig: REAL;
+	dummy, dig: LONGREAL;
 
   BEGIN
+	r := arg;
 	DEC(width);
 	IF digits < 0 THEN
 		ecvtflag := TRUE;
@@ -27,9 +36,10 @@ IMPLEMENTATION MODULE RealConversions;
 	END;
 	IF HIGH(str) < ndigits + 3 THEN str[0] := 0C; ok := FALSE; RETURN END;
 	pointpos := 0;
-	sign := r < 0.0;
+	sign := r < 0.0D;
 	IF sign THEN r := -r END;
-	r := FIF(r, 1.0, intpart);
+	r := FIF(r, 1.0D, intpart);
+	fractpart := r;
 	pointpos := 0;
 	ind1 := 0;
 	ok := TRUE;
@@ -37,9 +47,9 @@ IMPLEMENTATION MODULE RealConversions;
 	  Do integer part, which is now in "intpart". "r" now contains the
 	  fraction part.
 	*)
-	IF intpart # 0.0 THEN
+	IF intpart # 0.0D THEN
 		ind2 := 0;
-		WHILE intpart # 0.0 DO
+		WHILE intpart # 0.0D DO
 			IF ind2 > HIGH(str) THEN
 				IF NOT ecvtflag THEN
 					str[0] := 0C;
@@ -51,11 +61,11 @@ IMPLEMENTATION MODULE RealConversions;
 				END;
 				DEC(ind2);
 			END;
-			dummy := FIF(FIF(intpart, 0.1, intpart),10.0, dig);
-			IF (dummy > 0.5) AND (dig < 9.0) THEN
-				dig := dig + 1.0;
+			dummy := FIF(FIF(intpart, 0.1D, intpart),10.0D, dig);
+			IF (dummy > 0.5D) AND (dig < 9.0D) THEN
+				dig := dig + 1.0D;
 			END;
-			str[ind2] := CHR(TRUNC(dig+0.5) + ORD('0'));
+			str[ind2] := CHR(TRUNC(dig+0.5D) + ORD('0'));
 			INC(ind2);
 			INC(pointpos);
 		END;
@@ -70,10 +80,10 @@ IMPLEMENTATION MODULE RealConversions;
 		END;
 	ELSE
 		INC(pointpos);
-		IF r > 0.0 THEN
-			WHILE r < 1.0 DO
+		IF r > 0.0D THEN
+			WHILE r < 1.0D DO
 				fractpart := r;
-				r := r * 10.0;
+				r := r * 10.0D;
 				DEC(pointpos);
 			END;
 		END;
@@ -94,7 +104,7 @@ IMPLEMENTATION MODULE RealConversions;
 		RETURN;
 	END;
 	WHILE ind1 <= ind2 DO
-		fractpart := FIF(fractpart, 10.0, r);
+		fractpart := FIF(fractpart, 10.0D, r);
 		str[ind1] := CHR(TRUNC(r)+ORD('0'));
 		INC(ind1);
 	END;
@@ -191,17 +201,26 @@ IMPLEMENTATION MODULE RealConversions;
 	END;
 	IF (ind1+1) <= HIGH(str) THEN str[ind1+1] := 0C; END;
 
-  END RealToString;
+  END LongRealToString;
 
 	
   PROCEDURE StringToReal(str: ARRAY OF CHAR;
 			 VAR r: REAL; VAR ok: BOOLEAN);
+    VAR x: LONGREAL;
+  BEGIN
+	StringToLongReal(str, x, ok);
+	IF ok THEN
+		r := x;
+	END;
+  END StringToReal;
 
-    CONST	BIG = 1.0E17;
+  PROCEDURE StringToLongReal(str: ARRAY OF CHAR;
+			 VAR r: LONGREAL; VAR ok: BOOLEAN);
+    CONST	BIG = 1.0D17;
     TYPE	SETOFCHAR = SET OF CHAR;
     VAR		pow10 : INTEGER;
 		i : INTEGER;
-		e : REAL;
+		e : LONGREAL;
 		ch : CHAR;
 		signed: BOOLEAN;
 		signedexp: BOOLEAN;
@@ -209,11 +228,11 @@ IMPLEMENTATION MODULE RealConversions;
 
     PROCEDURE dig(ch: CARDINAL);
     BEGIN
-	IF r>BIG THEN INC(pow10) ELSE r:= 10.0*r + FLOAT(ch) END;
+	IF r>BIG THEN INC(pow10) ELSE r:= 10.0D*r + FLOATD(ch) END;
     END dig;
 
   BEGIN
-	r := 0.0;
+	r := 0.0D;
 	pow10 := 0;
 	iB := 0;
 	ok := TRUE;
@@ -276,10 +295,10 @@ IMPLEMENTATION MODULE RealConversions;
 		pow10 := pow10 + i;
 	END;
 	IF pow10 < 0 THEN i := -pow10; ELSE i := pow10; END;
-	e := 1.0;
+	e := 1.0D;
 	DEC(i);
 	WHILE i >= 0 DO
-		e := e * 10.0;
+		e := e * 10.0D;
 		DEC(i)
 	END;
 	IF pow10<0 THEN
@@ -289,6 +308,6 @@ IMPLEMENTATION MODULE RealConversions;
 	END;
 	IF signed THEN r := -r; END;
 	IF (iB <= HIGH(str)) AND (ORD(ch) > ORD(' ')) THEN ok := FALSE; END
-  END StringToReal;
+  END StringToLongReal;
 
 END RealConversions.