// CALC.sr

get "ST.DF";
get "CH.DF";
get "CALC.DF";

// Outgoing Procedures

external	[
	RealAdd;
	RealSub;
	RealMult;
	RealDiv;
	RealRadixConvert;
	RealIntExponent;
	Normalize;
	RealFromInt;
	RealFromSt;
	StFromReal;
	WordSizeReal;
	]

// Outgoing Statics

external	[
	cdigitpoint;
	frealmode;
	vintradix;
	];

// Incoming Procedures

external	[
	mult;
	divmod;
	umax;
	umin;
	move;
	stget;
	stsize;
	stput;
	stcopy;
	sbwsize;
	SbSetSize;
	stappend;
	stnum;
	stcompare;
	errhlt;
	]

// Local Statics

static	[
	cdigitpoint = 2;
	frealmode = frealfixed;
	vintradix = 10;
	];

// C H G E T R E A L

let ChGetReal(r, cp) = (cp ls 0 % cp ge r>>REAL.cch ? $0,
	r>>REAL.ch ^ cp);
// end ChGetReal


// C H M A P

and ChMap(ch) = valof
[
if ch ge $: & ch le $? then resultis ch + 7;
if ch ge $A & ch le $F then resultis ch - 7;
resultis ch;
] // end ChMap


// D I G I T A D D

and DigitAdd(a, b, cin, pcout) = valof
[
let sum = a + b + cin - 2 * $0;
if sum ge vintradix+$0 then
	[
	rv pcout = $1;
	resultis sum - vintradix;
	];
rv pcout = $0;
resultis sum;
] // end DigitAdd


// D I G I T M U L T

let Digitmult(a, b, cin, pcout) = valof
[
let prod = mult(a-$0, b-$0) + cin - $0;
rv pcout = divmod(prod, vintradix, lv prod) + $0;
resultis prod + $0;
] // end Digitmult


// M A G N I T U D E A D D

and MagnitudeAdd(rsum, ra, rb) = valof
[
let expa = ra>>REAL.exponent;
let expb = rb>>REAL.exponent;
let dcpa = 1;
let dcpb = 1;
test expa ls expb ifso dcpa = (expb - expa) + 1;
ifnot dcpb = (expa - expb) + 1;
let ccha = ra>>REAL.cch;
let cchb = rb>>REAL.cch;
let cchsum = umin(umax(ccha+dcpa, cchb+dcpb), cchmaxreal);
let carry = $0;
for cp = cchsum-1 to 0 by -1 do
	[
	rsum>>REAL.ch ^ cp = DigitAdd(
		ChGetReal(ra, cp - dcpa), ChGetReal(rb, cp - dcpb),
		carry, lv carry);
	];
rsum>>REAL.exponent = expa + dcpa;
rsum>>REAL.cch = cchsum;
resultis Normalize(rsum);
] // end MagnitudeAdd


// M A G N I T U D E S U B

and MagnitudeSub(rdiff, ra, rb) = valof
[
let expa = ra>>REAL.exponent;
let expb = rb>>REAL.exponent;
let dcpa = 0;
let dcpb = 0;
test expa ls expb ifso dcpa = expb - expa;
ifnot dcpb = expa - expb;
let ccha = ra>>REAL.cch;
let cchb = rb>>REAL.cch;
let cchdiff = umin(umax(ccha+dcpa, cchb+dcpb), cchmaxreal);
let chradix = vintradix + $0 - 1;
let carry = $1;
for cp = cchdiff-1 to 0 by -1 do
	[
	rdiff>>REAL.ch ^ cp = DigitAdd(
		ChGetReal(ra, cp - dcpa),
		chradix + $0 - ChGetReal(rb, cp - dcpb),
		carry, lv carry);
	];
rdiff>>REAL.exponent = expa + dcpa;
rdiff>>REAL.cch = cchdiff;
resultis Normalize(rdiff);
] // end MagnitudeSub


// F M A G N I T U D E C O M P A R E

and FMagnitudeCompare(ra, rb) = valof
[
let expa = ra>>REAL.exponent;
let expb = rb>>REAL.exponent;
if expa gr expb then resultis 1;
if expa ls expb then resultis -1;
resultis stcompare(lv ra>>REAL.asb, lv rb>>REAL.asb);
] // end FMagnitudeCompare


// R E A L A D D

and RealAdd(rsum, ra, rb) = valof
[
let fsigna = ra>>REAL.fsign;
let fsignb = rb>>REAL.fsign;
if not fsigna & not fsignb then
	[
	rsum>>REAL.fsign = false;
	resultis MagnitudeAdd(rsum, ra, rb);
	];
if fsigna & fsignb then
	[
	rsum>>REAL.fsign = true;
	resultis MagnitudeAdd(rsum, ra, rb);
	];
let fcompare = FMagnitudeCompare(ra, rb) ge 0;
rsum>>REAL.fsign = fsigna eqv fcompare;
if fcompare then resultis MagnitudeSub(rsum, ra, rb);
resultis MagnitudeSub(rsum, rb, ra);
] // end RealAdd


// R E A L S U B

and RealSub(rdiff, ra, rb) = valof
[
let fsigna = ra>>REAL.fsign;
let fsignb = rb>>REAL.fsign;
if not fsigna & fsignb then
	[
	rdiff>>REAL.fsign = false;
	resultis MagnitudeAdd(rdiff, ra, rb);
	];
if fsigna & not fsignb then
	[
	rdiff>>REAL.fsign = true;
	resultis MagnitudeAdd(rdiff, ra, rb);
	];
let fcompare = FMagnitudeCompare(ra, rb) ge 0;
rdiff>>REAL.fsign = fsigna eqv fcompare;
if fcompare then resultis MagnitudeSub(rdiff, ra, rb);
resultis MagnitudeSub(rdiff, rb, ra);
] // end RealSub


// R E A L M U L T

and RealMult(rprod, ra, rb) = valof
[
let trprod = vec cwmaxreal;
let ccha = ra>>REAL.cch;
let cchb = rb>>REAL.cch;
let cchprod = umin(ccha + cchb, cchmaxreal);
for cpprod = 0 to cchprod-1 do trprod>>REAL.ch ^ cpprod = $0;
for cpb = cchb-1 to 0 by -1 do
	[
	let addcarry = $0;
	let multcarry = $0;
	for cpa = ccha-1 to 0 by -1 do
		[
		let cpprod = cpa + cpb + 1;
		if cpprod ge cchprod then loop;
		trprod>>REAL.ch ^ cpprod = DigitAdd(
			trprod>>REAL.ch ^ cpprod,
			Digitmult(ra>>REAL.ch ^ cpa,
				rb>>REAL.ch ^ cpb,
				multcarry, lv multcarry),
			addcarry, lv addcarry);
		];
	trprod>>REAL.ch ^ cpb = DigitAdd(
		trprod>>REAL.ch ^ cpb,
		multcarry,
		addcarry, lv addcarry);
	];
trprod>>REAL.cch = cchprod;
trprod>>REAL.exponent = ra>>REAL.exponent + rb>>REAL.exponent + 1;
trprod>>REAL.fsign = ra>>REAL.fsign xor rb>>REAL.fsign;
move(Normalize(trprod), rprod, WordSizeReal(trprod));
resultis rprod;
] // end RealMult


// R E A L D I V

and RealDiv(rquo, ra, rb, cchquo; numargs carg) = valof
[
let trquo = vec cwmaxreal;
let ccha = ra>>REAL.cch;
let cchb = rb>>REAL.cch;
if cchb eq 0 then resultis realnil;
if carg ls 4 % cchquo eq 0 then cchquo = cchmaxreal;
cchquo = umin(cchquo, cchmaxreal);
let chradix = vintradix + $0 - 1;
SbSetSize(lv trquo>>REAL.asb, 0);
stcopy(lv trquo>>REAL.asb, lv ra>>REAL.asb);
for cpquo = ccha to cchquo-1 do trquo>>REAL.ch ^ cpquo = $0;
let chfirst = $0;
for cpquo = 0 to cchquo-1 do
	[
	let chquo = $0;
	let carry = $1;
	while carry eq $1 do
		[
		for cpb = cchb-1 to 0 by -1 do
			[
			let tcp = cpquo + cpb;
			if tcp ge cchquo then loop;
			trquo>>REAL.ch ^ tcp = DigitAdd(
				trquo>>REAL.ch ^ tcp,
				chradix + $0 - rb>>REAL.ch ^ cpb,
				carry, lv carry);
			];
		chfirst = DigitAdd(chfirst, chradix, carry, lv carry);
		chquo = chquo + 1;
		];
	for cpb = cchb-1 to 0 by -1 do
		[
		let tcp = cpquo + cpb;
		if tcp ge cchquo then loop;
		trquo>>REAL.ch ^ tcp = DigitAdd(
			trquo>>REAL.ch ^ tcp,
			rb>>REAL.ch ^ cpb,
			carry, lv carry);
		];
	chfirst = trquo>>REAL.ch ^ cpquo;
	trquo>>REAL.ch ^ cpquo = chquo - 1;
	];
trquo>>REAL.cch = cchquo;
trquo>>REAL.exponent = ra>>REAL.exponent - rb>>REAL.exponent;
trquo>>REAL.fsign = ra>>REAL.fsign xor rb>>REAL.fsign;
move(Normalize(trquo), rquo, WordSizeReal(trquo));
resultis rquo;
] // end RealDiv


// R E A L R A D I X C O N V E R T

and RealRadixConvert(rresult, r, rradix) = valof
[
let tr = vec cwmaxreal;
rresult>>REAL.fsign = false;
rresult>>REAL.exponent = 0;
SbSetSize(lv rresult>>REAL.asb, 0);
for cp = 0 to r>>REAL.cch - 1 do
	RealAdd(rresult, RealMult(rresult, rresult, rradix),
		RealFromInt(tr, r>>REAL.ch ^ cp - $0));
resultis RealMult(rresult, rresult,
	RealIntExponent(tr, rradix, r>>REAL.exponent + 1 - r>>REAL.cch));
] // end RealRadixConvert


// R E A L I N T E X P O N E N T

and RealIntExponent(rresult, r, int) = valof
[
let fsign = false;
if int ls 0 then
	[
	fsign = true;
	int = -int;
	];
let rone = table
	[
	0;
	0;
	1 lshift 8 + $1;
	];
move(rone, rresult, 3);
let tr = vec cwmaxreal;
move(r, tr, WordSizeReal(r));
	[
	if (int & 1) ne 0 then RealMult(rresult, rresult, tr);
	int = int rshift 1;
	if int eq 0 then break;
	RealMult(tr, tr, tr);
	] repeat;
resultis (not fsign ? rresult,
	RealDiv(rresult, rone, rresult));
] // end RealIntExponent


// N O R M A L I Z E

and Normalize(r) = valof
[
let cch = r>>REAL.cch;
let sb = lv r>>REAL.asb;
let cpfirst = 0;
while cpfirst ls cch do
	[
	if sb>>SB.ch ^ cpfirst ne $0 then break;
	cpfirst = cpfirst + 1;
	];
let cplast = cch-1;
while cplast ge cpfirst do
	[
	if sb>>SB.ch ^ cplast ne $0 then break;
	cplast = cplast - 1;
	];
if cpfirst gr cplast then
	[
	r>>REAL.fsign = false;
	r>>REAL.exponent = 0;
	r>>REAL.asb = 0;
	resultis r;
	];
let cchnew = cplast+1 - cpfirst
if cchnew eq cch then resultis r;
for cp = cpfirst to cplast do
	sb>>SB.ch ^ (cp-cpfirst) = sb>>SB.ch ^ cp
sb>>SB.cch = cchnew
r>>REAL.exponent = r>>REAL.exponent - cpfirst;
resultis r;
] // end Normalize


// R E A L F R O M I N T

and RealFromInt(r, int) = valof
[
r>>REAL.fsign = false;
if int ls 0 then
	[
	r>>REAL.fsign = true;
	int = -int;
	];
SbSetSize(lv r>>REAL.asb, 0);
stnum(lv r>>REAL.asb, int, vintradix);
r>>REAL.exponent = r>>REAL.cch - 1;
resultis Normalize(r);
] // end RealFromInt


// R E A L F R O M S T

and RealFromSt(r, st) = valof
[
let cp = 0;
let state = 0;
let cpreal = 0;
let cppoint = cpnil;
let fexpsign = false;
let exponent = 0;
let chmax = vintradix + $0;
let ch = nil;
while valof [ ch = stget(st, cp); cp = cp+1; resultis ch ] ne chnil do
teststate:
	switchon state into
		[
	case 0:	state = 1;
		r>>REAL.fsign = ch eq $-;
		if ch eq $+  % ch eq $- then loop;
		if ChMap(ch) ge $0 & ChMap(ch) ls chmax
			% ch eq $. then goto teststate;
		resultis realnil;

	case 1:	if ch eq $. then
			[
			cppoint = cpreal;
			state = 2;
			loop;
			];
		if ch eq $^ then
			[
			cppoint = cpreal;
			state = 3;
			loop;
			];
		ch = ChMap(ch);
		if ch ls $0 % ch ge chmax then resultis realnil;
		if cpreal ls cchmaxreal then
			[
			r>>REAL.ch ^ cpreal = ch;
			cpreal = cpreal+1;
			];
		loop;

	case 2:	if ch eq $^ then
			[
			state = 3;
			loop;
			];
		ch = ChMap(ch);
		if ch ls $0 % ch ge chmax then resultis realnil;
		if cpreal ls cchmaxreal then
			[
			r>>REAL.ch ^ cpreal = ch;
			cpreal = cpreal+1;
			];
		loop;

	case 3:	state = 4;
		fexpsign = ch eq $-;
		if ch eq $+ % ch eq $- then loop;
		if ChMap(ch) ge $0 & ChMap(ch) ls chmax then
			goto teststate;
		resultis realnil;

	case 4:	ch = ChMap(ch);
		if ch ls $0 % ch ge chmax then resultis realnil;
		exponent = mult(exponent, 10) + ch-$0;
		if exponent ls 0 then resultis realnil;
		loop;
		];
if state eq 0 % state eq 3 then resultis realnil;
SbSetSize(lv r>>REAL.asb, cpreal);
exponent = (fexpsign ? -exponent, exponent);
r>>REAL.exponent = (cppoint eq cpnil ? cpreal-1, exponent + cppoint-1);
resultis Normalize(r);
] // end RealFromSt


// S T F R O M R E A L

and StFromReal(st, r, cchst; numargs carg) = valof
[
if carg ls 3 then cchst = -1
let tr = vec cwmaxreal;
let tfrealmode = frealmode;
let tcdigitpoint = cdigitpoint;
let cppoint = nil;
let cpmax = nil;
let cpreal = nil;
round:
let rround = vec 3;
rround>>REAL.exponent = selecton tfrealmode into
	[
case frealfixed:	0
case frealsci:	r>>REAL.exponent
case frealeng:	mult(divmod(r>>REAL.exponent + #100001, 3, lv cppoint), 3) - #100001
	] - (tcdigitpoint + 1);
rround>>REAL.asb = (r>>REAL.cch eq 0 ? 0, 1 lshift 8 + vintradix rshift 1 + $0);
MagnitudeAdd(tr, r, rround);
let exp = tr>>REAL.exponent;
if tfrealmode eq frealeng then
	exp = mult(divmod(exp + #100001, 3, lv cppoint), 3) - #100001
let cchreal = tr>>REAL.cch;
let sbexp = vec 4;
sbexp>>SB.ch ^ 0 = $^;
SbSetSize(sbexp, 1);
let tsb = vec 4;
SbSetSize(tsb, 0);
switchon tfrealmode into
	[
case frealfixed:
	sbexp = "";
	cppoint = (exp ls 0 ? 1, exp+1);
	cpmax = cppoint + (tcdigitpoint le 0 ? 0, tcdigitpoint + 1);
	cpreal = (exp ge 0 ? 0, exp);
	endcase;

case frealsci:
	stappend(sbexp, stnum(tsb, exp));
	cppoint = 1;
	cpmax = tcdigitpoint + 2;
	cpreal = 0;
	endcase;

case frealeng:
	stappend(sbexp, stnum(tsb, exp));
	cppoint = cppoint + 1;
	cpmax = cppoint + tcdigitpoint + 1;
	cpreal = 0;
	endcase;
	];
let cchtotal = cpmax + stsize(sbexp) + (r>>REAL.fsign ? 1, 0);
cchst = umin(cchst, cchtotal)
// if rv st eq 1 then cchst = st>>ST.cch;
let dcch = cchst - cchtotal;
if dcch ls 0 then
	[
	if tfrealmode ne frealsci then
		[
		tfrealmode = frealsci;
		goto round;
		];
	tcdigitpoint = tcdigitpoint + dcch;
	if tcdigitpoint ge 0 then goto round;
	for tcp = 0 to cchst-1 do
		stput(st, tcp, $**);
// 	if rv st ne 1 then
		SbSetSize(st, cchst);
	resultis st;
 	];
let cp = 0;
while cp ls dcch do
	[
	stput(st, cp, chblank);
	cp = cp + 1;
	];
if r>>REAL.fsign then
	[
	stput(st, cp, $-);
	cp = cp + 1;
	];
let dcp = 0;
while dcp ls cpmax do
	[
	if dcp eq cppoint then
		[
		stput(st, cp+dcp, $.);
		dcp = dcp + 1;
		loop;
		];
	stput(st, cp+dcp, ChMap(ChGetReal(tr, cpreal)));
	dcp = dcp + 1;
	cpreal = cpreal + 1;
	];
let dcp = cp+cpmax
for tcp = 0 to stsize(sbexp)-1 do
	stput(st, tcp+dcp, stget(sbexp, tcp))
// if rv st ne 1 then
	SbSetSize(st, cchst);
resultis st;
] // end StFromReal


// W O R D S I Z E R E A L

and WordSizeReal(r) = (r>>REAL.cch rshift 1)+offasbreal+1;
// end WordSizeReal