1	// SPDX-License-Identifier: GPL-2.0
2	/*---------------------------------------------------------------------------+
3	| fpu_trig.c |
4	| |
5	| Implementation of the FPU "transcendental" functions. |
6	| |
7	| Copyright (C) 1992,1993,1994,1997,1999 |
8	| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
9	| Australia. E-mail billm@melbpc.org.au |
10	| |
11	| |
12	+---------------------------------------------------------------------------*/
13
14	#include "fpu_system.h"
15	#include "exception.h"
16	#include "fpu_emu.h"
17	#include "status_w.h"
18	#include "control_w.h"
19	#include "reg_constant.h"
20
21	static void rem_kernel(unsigned long long st0, unsigned long long *y,
22	unsigned long long st1, unsigned long long q, int n);
23
24	#define BETTER_THAN_486
25
26	#define FCOS 4
27
28	/* Used only by fptan, fsin, fcos, and fsincos. */
29	/* This routine produces very accurate results, similar to
30	using a value of pi with more than 128 bits precision. */
31	/* Limited measurements show no results worse than 64 bit precision
32	except for the results for arguments close to 2^63, where the
33	precision of the result sometimes degrades to about 63.9 bits */
34	static int trig_arg(FPU_REG *st0_ptr, int even)
35	{
36	FPU_REG tmp;
37	u_char tmptag;
38	unsigned long long q;
39	int old_cw = control_word, saved_status = partial_status;
40	int tag, st0_tag = TAG_Valid;
41
42	if (exponent(st0_ptr) >= 63) {
43	partial_status |= SW_C2; /* Reduction incomplete. */
44	return -1;
45	}
46
47	control_word &= ~CW_RC;
48	control_word |= RC_CHOP;
49
50	setpositive(st0_ptr);
51	tag = FPU_u_div(arg1: st0_ptr, arg2: &CONST_PI2, answ: &tmp, PR_64_BITS | RC_CHOP | 0x3f,
52	SIGN_POS);
53
54	FPU_round_to_int(r: &tmp, tag); /* Fortunately, this can't overflow
55	to 2^64 */
56	q = significand(&tmp);
57	if (q) {
58	rem_kernel(significand(st0_ptr),
59	y: &significand(&tmp),
60	significand(&CONST_PI2),
61	q, exponent(st0_ptr) - exponent(&CONST_PI2));
62	setexponent16(&tmp, exponent(&CONST_PI2));
63	st0_tag = FPU_normalize(x: &tmp);
64	FPU_copy_to_reg0(r: &tmp, tag: st0_tag);
65	}
66
67	if ((even && !(q & 1)) || (!even && (q & 1))) {
68	st0_tag =
69	FPU_sub(REV | LOADED | TAG_Valid, rm: (int)&CONST_PI2,
70	FULL_PRECISION);
71
72	#ifdef BETTER_THAN_486
73	/* So far, the results are exact but based upon a 64 bit
74	precision approximation to pi/2. The technique used
75	now is equivalent to using an approximation to pi/2 which
76	is accurate to about 128 bits. */
77	if ((exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64)
78	|| (q > 1)) {
79	/* This code gives the effect of having pi/2 to better than
80	128 bits precision. */
81
82	significand(&tmp) = q + 1;
83	setexponent16(&tmp, 63);
84	FPU_normalize(x: &tmp);
85	tmptag =
86	FPU_u_mul(arg1: &CONST_PI2extra, arg2: &tmp, answ: &tmp,
87	FULL_PRECISION, SIGN_POS,
88	exponent(&CONST_PI2extra) +
89	exponent(&tmp));
90	setsign(&tmp, getsign(&CONST_PI2extra));
91	st0_tag = FPU_add(b: &tmp, tagb: tmptag, destrnr: 0, FULL_PRECISION);
92	if (signnegative(st0_ptr)) {
93	/* CONST_PI2extra is negative, so the result of the addition
94	can be negative. This means that the argument is actually
95	in a different quadrant. The correction is always < pi/2,
96	so it can't overflow into yet another quadrant. */
97	setpositive(st0_ptr);
98	q++;
99	}
100	}
101	#endif /* BETTER_THAN_486 */
102	}
103	#ifdef BETTER_THAN_486
104	else {
105	/* So far, the results are exact but based upon a 64 bit
106	precision approximation to pi/2. The technique used
107	now is equivalent to using an approximation to pi/2 which
108	is accurate to about 128 bits. */
109	if (((q > 0)
110	&& (exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64))
111	|| (q > 1)) {
112	/* This code gives the effect of having p/2 to better than
113	128 bits precision. */
114
115	significand(&tmp) = q;
116	setexponent16(&tmp, 63);
117	FPU_normalize(x: &tmp); /* This must return TAG_Valid */
118	tmptag =
119	FPU_u_mul(arg1: &CONST_PI2extra, arg2: &tmp, answ: &tmp,
120	FULL_PRECISION, SIGN_POS,
121	exponent(&CONST_PI2extra) +
122	exponent(&tmp));
123	setsign(&tmp, getsign(&CONST_PI2extra));
124	st0_tag = FPU_sub(LOADED | (tmptag & 0x0f), rm: (int)&tmp,
125	FULL_PRECISION);
126	if ((exponent(st0_ptr) == exponent(&CONST_PI2)) &&
127	((st0_ptr->sigh > CONST_PI2.sigh)
128	|| ((st0_ptr->sigh == CONST_PI2.sigh)
129	&& (st0_ptr->sigl > CONST_PI2.sigl)))) {
130	/* CONST_PI2extra is negative, so the result of the
131	subtraction can be larger than pi/2. This means
132	that the argument is actually in a different quadrant.
133	The correction is always < pi/2, so it can't overflow
134	into yet another quadrant. */
135	st0_tag =
136	FPU_sub(REV | LOADED | TAG_Valid,
137	rm: (int)&CONST_PI2, FULL_PRECISION);
138	q++;
139	}
140	}
141	}
142	#endif /* BETTER_THAN_486 */
143
144	FPU_settag0(tag: st0_tag);
145	control_word = old_cw;
146	partial_status = saved_status & ~SW_C2; /* Reduction complete. */
147
148	return (q & 3) | even;
149	}
150
151	/* Convert a long to register */
152	static void convert_l2reg(long const *arg, int deststnr)
153	{
154	int tag;
155	long num = *arg;
156	u_char sign;
157	FPU_REG *dest = &st(deststnr);
158
159	if (num == 0) {
160	FPU_copy_to_regi(r: &CONST_Z, TAG_Zero, stnr: deststnr);
161	return;
162	}
163
164	if (num > 0) {
165	sign = SIGN_POS;
166	} else {
167	num = -num;
168	sign = SIGN_NEG;
169	}
170
171	dest->sigh = num;
172	dest->sigl = 0;
173	setexponent16(dest, 31);
174	tag = FPU_normalize(x: dest);
175	FPU_settagi(stnr: deststnr, tag);
176	setsign(dest, sign);
177	return;
178	}
179
180	static void single_arg_error(FPU_REG *st0_ptr, u_char st0_tag)
181	{
182	if (st0_tag == TAG_Empty)
183	FPU_stack_underflow(); /* Puts a QNaN in st(0) */
184	else if (st0_tag == TW_NaN)
185	real_1op_NaN(a: st0_ptr); /* return with a NaN in st(0) */
186	#ifdef PARANOID
187	else
188	EXCEPTION(EX_INTERNAL | 0x0112);
189	#endif /* PARANOID */
190	}
191
192	static void single_arg_2_error(FPU_REG *st0_ptr, u_char st0_tag)
193	{
194	int isNaN;
195
196	switch (st0_tag) {
197	case TW_NaN:
198	isNaN = (exponent(st0_ptr) == EXP_OVER)
199	&& (st0_ptr->sigh & 0x80000000);
200	if (isNaN && !(st0_ptr->sigh & 0x40000000)) { /* Signaling ? */
201	EXCEPTION(EX_Invalid);
202	if (control_word & CW_Invalid) {
203	/* The masked response */
204	/* Convert to a QNaN */
205	st0_ptr->sigh |= 0x40000000;
206	push();
207	FPU_copy_to_reg0(r: st0_ptr, TAG_Special);
208	}
209	} else if (isNaN) {
210	/* A QNaN */
211	push();
212	FPU_copy_to_reg0(r: st0_ptr, TAG_Special);
213	} else {
214	/* pseudoNaN or other unsupported */
215	EXCEPTION(EX_Invalid);
216	if (control_word & CW_Invalid) {
217	/* The masked response */
218	FPU_copy_to_reg0(r: &CONST_QNaN, TAG_Special);
219	push();
220	FPU_copy_to_reg0(r: &CONST_QNaN, TAG_Special);
221	}
222	}
223	break; /* return with a NaN in st(0) */
224	#ifdef PARANOID
225	default:
226	EXCEPTION(EX_INTERNAL | 0x0112);
227	#endif /* PARANOID */
228	}
229	}
230
231	/*---------------------------------------------------------------------------*/
232
233	static void f2xm1(FPU_REG *st0_ptr, u_char tag)
234	{
235	FPU_REG a;
236
237	clear_C1();
238
239	if (tag == TAG_Valid) {
240	/* For an 80486 FPU, the result is undefined if the arg is >= 1.0 */
241	if (exponent(st0_ptr) < 0) {
242	denormal_arg:
243
244	FPU_to_exp16(a: st0_ptr, x: &a);
245
246	/* poly_2xm1(x) requires 0 < st(0) < 1. */
247	poly_2xm1(getsign(st0_ptr), arg: &a, result: st0_ptr);
248	}
249	set_precision_flag_up(); /* 80486 appears to always do this */
250	return;
251	}
252
253	if (tag == TAG_Zero)
254	return;
255
256	if (tag == TAG_Special)
257	tag = FPU_Special(ptr: st0_ptr);
258
259	switch (tag) {
260	case TW_Denormal:
261	if (denormal_operand() < 0)
262	return;
263	goto denormal_arg;
264	case TW_Infinity:
265	if (signnegative(st0_ptr)) {
266	/* -infinity gives -1 (p16-10) */
267	FPU_copy_to_reg0(r: &CONST_1, TAG_Valid);
268	setnegative(st0_ptr);
269	}
270	return;
271	default:
272	single_arg_error(st0_ptr, st0_tag: tag);
273	}
274	}
275
276	static void fptan(FPU_REG *st0_ptr, u_char st0_tag)
277	{
278	FPU_REG *st_new_ptr;
279	int q;
280	u_char arg_sign = getsign(st0_ptr);
281
282	/* Stack underflow has higher priority */
283	if (st0_tag == TAG_Empty) {
284	FPU_stack_underflow(); /* Puts a QNaN in st(0) */
285	if (control_word & CW_Invalid) {
286	st_new_ptr = &st(-1);
287	push();
288	FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
289	}
290	return;
291	}
292
293	if (STACK_OVERFLOW) {
294	FPU_stack_overflow();
295	return;
296	}
297
298	if (st0_tag == TAG_Valid) {
299	if (exponent(st0_ptr) > -40) {
300	if ((q = trig_arg(st0_ptr, even: 0)) == -1) {
301	/* Operand is out of range */
302	return;
303	}
304
305	poly_tan(st0_ptr);
306	setsign(st0_ptr, (q & 1) ^ (arg_sign != 0));
307	set_precision_flag_up(); /* We do not really know if up or down */
308	} else {
309	/* For a small arg, the result == the argument */
310	/* Underflow may happen */
311
312	denormal_arg:
313
314	FPU_to_exp16(a: st0_ptr, x: st0_ptr);
315
316	st0_tag =
317	FPU_round(arg: st0_ptr, extent: 1, dummy: 0, FULL_PRECISION, sign: arg_sign);
318	FPU_settag0(tag: st0_tag);
319	}
320	push();
321	FPU_copy_to_reg0(r: &CONST_1, TAG_Valid);
322	return;
323	}
324
325	if (st0_tag == TAG_Zero) {
326	push();
327	FPU_copy_to_reg0(r: &CONST_1, TAG_Valid);
328	setcc(0);
329	return;
330	}
331
332	if (st0_tag == TAG_Special)
333	st0_tag = FPU_Special(ptr: st0_ptr);
334
335	if (st0_tag == TW_Denormal) {
336	if (denormal_operand() < 0)
337	return;
338
339	goto denormal_arg;
340	}
341
342	if (st0_tag == TW_Infinity) {
343	/* The 80486 treats infinity as an invalid operand */
344	if (arith_invalid(deststnr: 0) >= 0) {
345	st_new_ptr = &st(-1);
346	push();
347	arith_invalid(deststnr: 0);
348	}
349	return;
350	}
351
352	single_arg_2_error(st0_ptr, st0_tag);
353	}
354
355	static void fxtract(FPU_REG *st0_ptr, u_char st0_tag)
356	{
357	FPU_REG *st_new_ptr;
358	u_char sign;
359	register FPU_REG *st1_ptr = st0_ptr; /* anticipate */
360
361	if (STACK_OVERFLOW) {
362	FPU_stack_overflow();
363	return;
364	}
365
366	clear_C1();
367
368	if (st0_tag == TAG_Valid) {
369	long e;
370
371	push();
372	sign = getsign(st1_ptr);
373	reg_copy(x: st1_ptr, y: st_new_ptr);
374	setexponent16(st_new_ptr, exponent(st_new_ptr));
375
376	denormal_arg:
377
378	e = exponent16(st_new_ptr);
379	convert_l2reg(arg: &e, deststnr: 1);
380	setexponentpos(st_new_ptr, 0);
381	setsign(st_new_ptr, sign);
382	FPU_settag0(TAG_Valid); /* Needed if arg was a denormal */
383	return;
384	} else if (st0_tag == TAG_Zero) {
385	sign = getsign(st0_ptr);
386
387	if (FPU_divide_by_zero(deststnr: 0, SIGN_NEG) < 0)
388	return;
389
390	push();
391	FPU_copy_to_reg0(r: &CONST_Z, TAG_Zero);
392	setsign(st_new_ptr, sign);
393	return;
394	}
395
396	if (st0_tag == TAG_Special)
397	st0_tag = FPU_Special(ptr: st0_ptr);
398
399	if (st0_tag == TW_Denormal) {
400	if (denormal_operand() < 0)
401	return;
402
403	push();
404	sign = getsign(st1_ptr);
405	FPU_to_exp16(a: st1_ptr, x: st_new_ptr);
406	goto denormal_arg;
407	} else if (st0_tag == TW_Infinity) {
408	sign = getsign(st0_ptr);
409	setpositive(st0_ptr);
410	push();
411	FPU_copy_to_reg0(r: &CONST_INF, TAG_Special);
412	setsign(st_new_ptr, sign);
413	return;
414	} else if (st0_tag == TW_NaN) {
415	if (real_1op_NaN(a: st0_ptr) < 0)
416	return;
417
418	push();
419	FPU_copy_to_reg0(r: st0_ptr, TAG_Special);
420	return;
421	} else if (st0_tag == TAG_Empty) {
422	/* Is this the correct behaviour? */
423	if (control_word & EX_Invalid) {
424	FPU_stack_underflow();
425	push();
426	FPU_stack_underflow();
427	} else
428	EXCEPTION(EX_StackUnder);
429	}
430	#ifdef PARANOID
431	else
432	EXCEPTION(EX_INTERNAL | 0x119);
433	#endif /* PARANOID */
434	}
435
436	static void fdecstp(void)
437	{
438	clear_C1();
439	top--;
440	}
441
442	static void fincstp(void)
443	{
444	clear_C1();
445	top++;
446	}
447
448	static void fsqrt_(FPU_REG *st0_ptr, u_char st0_tag)
449	{
450	int expon;
451
452	clear_C1();
453
454	if (st0_tag == TAG_Valid) {
455	u_char tag;
456
457	if (signnegative(st0_ptr)) {
458	arith_invalid(deststnr: 0); /* sqrt(negative) is invalid */
459	return;
460	}
461
462	/* make st(0) in [1.0 .. 4.0) */
463	expon = exponent(st0_ptr);
464
465	denormal_arg:
466
467	setexponent16(st0_ptr, (expon & 1));
468
469	/* Do the computation, the sign of the result will be positive. */
470	tag = wm_sqrt(n: st0_ptr, dummy1: 0, dummy2: 0, control_word, SIGN_POS);
471	addexponent(st0_ptr, expon >> 1);
472	FPU_settag0(tag);
473	return;
474	}
475
476	if (st0_tag == TAG_Zero)
477	return;
478
479	if (st0_tag == TAG_Special)
480	st0_tag = FPU_Special(ptr: st0_ptr);
481
482	if (st0_tag == TW_Infinity) {
483	if (signnegative(st0_ptr))
484	arith_invalid(deststnr: 0); /* sqrt(-Infinity) is invalid */
485	return;
486	} else if (st0_tag == TW_Denormal) {
487	if (signnegative(st0_ptr)) {
488	arith_invalid(deststnr: 0); /* sqrt(negative) is invalid */
489	return;
490	}
491
492	if (denormal_operand() < 0)
493	return;
494
495	FPU_to_exp16(a: st0_ptr, x: st0_ptr);
496
497	expon = exponent16(st0_ptr);
498
499	goto denormal_arg;
500	}
501
502	single_arg_error(st0_ptr, st0_tag);
503
504	}
505
506	static void frndint_(FPU_REG *st0_ptr, u_char st0_tag)
507	{
508	int flags, tag;
509
510	if (st0_tag == TAG_Valid) {
511	u_char sign;
512
513	denormal_arg:
514
515	sign = getsign(st0_ptr);
516
517	if (exponent(st0_ptr) > 63)
518	return;
519
520	if (st0_tag == TW_Denormal) {
521	if (denormal_operand() < 0)
522	return;
523	}
524
525	/* Fortunately, this can't overflow to 2^64 */
526	if ((flags = FPU_round_to_int(r: st0_ptr, tag: st0_tag)))
527	set_precision_flag(flags);
528
529	setexponent16(st0_ptr, 63);
530	tag = FPU_normalize(x: st0_ptr);
531	setsign(st0_ptr, sign);
532	FPU_settag0(tag);
533	return;
534	}
535
536	if (st0_tag == TAG_Zero)
537	return;
538
539	if (st0_tag == TAG_Special)
540	st0_tag = FPU_Special(ptr: st0_ptr);
541
542	if (st0_tag == TW_Denormal)
543	goto denormal_arg;
544	else if (st0_tag == TW_Infinity)
545	return;
546	else
547	single_arg_error(st0_ptr, st0_tag);
548	}
549
550	static int f_sin(FPU_REG *st0_ptr, u_char tag)
551	{
552	u_char arg_sign = getsign(st0_ptr);
553
554	if (tag == TAG_Valid) {
555	int q;
556
557	if (exponent(st0_ptr) > -40) {
558	if ((q = trig_arg(st0_ptr, even: 0)) == -1) {
559	/* Operand is out of range */
560	return 1;
561	}
562
563	poly_sine(st0_ptr);
564
565	if (q & 2)
566	changesign(st0_ptr);
567
568	setsign(st0_ptr, getsign(st0_ptr) ^ arg_sign);
569
570	/* We do not really know if up or down */
571	set_precision_flag_up();
572	return 0;
573	} else {
574	/* For a small arg, the result == the argument */
575	set_precision_flag_up(); /* Must be up. */
576	return 0;
577	}
578	}
579
580	if (tag == TAG_Zero) {
581	setcc(0);
582	return 0;
583	}
584
585	if (tag == TAG_Special)
586	tag = FPU_Special(ptr: st0_ptr);
587
588	if (tag == TW_Denormal) {
589	if (denormal_operand() < 0)
590	return 1;
591
592	/* For a small arg, the result == the argument */
593	/* Underflow may happen */
594	FPU_to_exp16(a: st0_ptr, x: st0_ptr);
595
596	tag = FPU_round(arg: st0_ptr, extent: 1, dummy: 0, FULL_PRECISION, sign: arg_sign);
597
598	FPU_settag0(tag);
599
600	return 0;
601	} else if (tag == TW_Infinity) {
602	/* The 80486 treats infinity as an invalid operand */
603	arith_invalid(deststnr: 0);
604	return 1;
605	} else {
606	single_arg_error(st0_ptr, st0_tag: tag);
607	return 1;
608	}
609	}
610
611	static void fsin(FPU_REG *st0_ptr, u_char tag)
612	{
613	f_sin(st0_ptr, tag);
614	}
615
616	static int f_cos(FPU_REG *st0_ptr, u_char tag)
617	{
618	u_char st0_sign;
619
620	st0_sign = getsign(st0_ptr);
621
622	if (tag == TAG_Valid) {
623	int q;
624
625	if (exponent(st0_ptr) > -40) {
626	if ((exponent(st0_ptr) < 0)
627	|| ((exponent(st0_ptr) == 0)
628	&& (significand(st0_ptr) <=
629	0xc90fdaa22168c234LL))) {
630	poly_cos(st0_ptr);
631
632	/* We do not really know if up or down */
633	set_precision_flag_down();
634
635	return 0;
636	} else if ((q = trig_arg(st0_ptr, FCOS)) != -1) {
637	poly_sine(st0_ptr);
638
639	if ((q + 1) & 2)
640	changesign(st0_ptr);
641
642	/* We do not really know if up or down */
643	set_precision_flag_down();
644
645	return 0;
646	} else {
647	/* Operand is out of range */
648	return 1;
649	}
650	} else {
651	denormal_arg:
652
653	setcc(0);
654	FPU_copy_to_reg0(r: &CONST_1, TAG_Valid);
655	#ifdef PECULIAR_486
656	set_precision_flag_down(); /* 80486 appears to do this. */
657	#else
658	set_precision_flag_up(); /* Must be up. */
659	#endif /* PECULIAR_486 */
660	return 0;
661	}
662	} else if (tag == TAG_Zero) {
663	FPU_copy_to_reg0(r: &CONST_1, TAG_Valid);
664	setcc(0);
665	return 0;
666	}
667
668	if (tag == TAG_Special)
669	tag = FPU_Special(ptr: st0_ptr);
670
671	if (tag == TW_Denormal) {
672	if (denormal_operand() < 0)
673	return 1;
674
675	goto denormal_arg;
676	} else if (tag == TW_Infinity) {
677	/* The 80486 treats infinity as an invalid operand */
678	arith_invalid(deststnr: 0);
679	return 1;
680	} else {
681	single_arg_error(st0_ptr, st0_tag: tag); /* requires st0_ptr == &st(0) */
682	return 1;
683	}
684	}
685
686	static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
687	{
688	f_cos(st0_ptr, tag: st0_tag);
689	}
690
691	static void fsincos(FPU_REG *st0_ptr, u_char st0_tag)
692	{
693	FPU_REG *st_new_ptr;
694	FPU_REG arg;
695	u_char tag;
696
697	/* Stack underflow has higher priority */
698	if (st0_tag == TAG_Empty) {
699	FPU_stack_underflow(); /* Puts a QNaN in st(0) */
700	if (control_word & CW_Invalid) {
701	st_new_ptr = &st(-1);
702	push();
703	FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
704	}
705	return;
706	}
707
708	if (STACK_OVERFLOW) {
709	FPU_stack_overflow();
710	return;
711	}
712
713	if (st0_tag == TAG_Special)
714	tag = FPU_Special(ptr: st0_ptr);
715	else
716	tag = st0_tag;
717
718	if (tag == TW_NaN) {
719	single_arg_2_error(st0_ptr, TW_NaN);
720	return;
721	} else if (tag == TW_Infinity) {
722	/* The 80486 treats infinity as an invalid operand */
723	if (arith_invalid(deststnr: 0) >= 0) {
724	/* Masked response */
725	push();
726	arith_invalid(deststnr: 0);
727	}
728	return;
729	}
730
731	reg_copy(x: st0_ptr, y: &arg);
732	if (!f_sin(st0_ptr, tag: st0_tag)) {
733	push();
734	FPU_copy_to_reg0(r: &arg, tag: st0_tag);
735	f_cos(st0_ptr: &st(0), tag: st0_tag);
736	} else {
737	/* An error, so restore st(0) */
738	FPU_copy_to_reg0(r: &arg, tag: st0_tag);
739	}
740	}
741
742	/*---------------------------------------------------------------------------*/
743	/* The following all require two arguments: st(0) and st(1) */
744
745	/* A lean, mean kernel for the fprem instructions. This relies upon
746	the division and rounding to an integer in do_fprem giving an
747	exact result. Because of this, rem_kernel() needs to deal only with
748	the least significant 64 bits, the more significant bits of the
749	result must be zero.
750	*/
751	static void rem_kernel(unsigned long long st0, unsigned long long *y,
752	unsigned long long st1, unsigned long long q, int n)
753	{
754	int dummy;
755	unsigned long long x;
756
757	x = st0 << n;
758
759	/* Do the required multiplication and subtraction in the one operation */
760
761	/* lsw x -= lsw st1 * lsw q */
762	asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
763	(((unsigned *)&x)[0]), "=m"(((unsigned *)&x)[1]),
764	"=a"(dummy)
765	:"2"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[0])
766	:"%dx");
767	/* msw x -= msw st1 * lsw q */
768	asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
769	"=a"(dummy)
770	:"1"(((unsigned *)&st1)[1]), "m"(((unsigned *)&q)[0])
771	:"%dx");
772	/* msw x -= lsw st1 * msw q */
773	asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
774	"=a"(dummy)
775	:"1"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[1])
776	:"%dx");
777
778	*y = x;
779	}
780
781	/* Remainder of st(0) / st(1) */
782	/* This routine produces exact results, i.e. there is never any
783	rounding or truncation, etc of the result. */
784	static void do_fprem(FPU_REG *st0_ptr, u_char st0_tag, int round)
785	{
786	FPU_REG *st1_ptr = &st(1);
787	u_char st1_tag = FPU_gettagi(stnr: 1);
788
789	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
790	FPU_REG tmp, st0, st1;
791	u_char st0_sign, st1_sign;
792	u_char tmptag;
793	int tag;
794	int old_cw;
795	int expdif;
796	long long q;
797	unsigned short saved_status;
798	int cc;
799
800	fprem_valid:
801	/* Convert registers for internal use. */
802	st0_sign = FPU_to_exp16(a: st0_ptr, x: &st0);
803	st1_sign = FPU_to_exp16(a: st1_ptr, x: &st1);
804	expdif = exponent16(&st0) - exponent16(&st1);
805
806	old_cw = control_word;
807	cc = 0;
808
809	/* We want the status following the denorm tests, but don't want
810	the status changed by the arithmetic operations. */
811	saved_status = partial_status;
812	control_word &= ~CW_RC;
813	control_word |= RC_CHOP;
814
815	if (expdif < 64) {
816	/* This should be the most common case */
817
818	if (expdif > -2) {
819	u_char sign = st0_sign ^ st1_sign;
820	tag = FPU_u_div(arg1: &st0, arg2: &st1, answ: &tmp,
821	PR_64_BITS | RC_CHOP | 0x3f,
822	sign);
823	setsign(&tmp, sign);
824
825	if (exponent(&tmp) >= 0) {
826	FPU_round_to_int(r: &tmp, tag); /* Fortunately, this can't
827	overflow to 2^64 */
828	q = significand(&tmp);
829
830	rem_kernel(significand(&st0),
831	y: &significand(&tmp),
832	significand(&st1),
833	q, n: expdif);
834
835	setexponent16(&tmp, exponent16(&st1));
836	} else {
837	reg_copy(x: &st0, y: &tmp);
838	q = 0;
839	}
840
841	if ((round == RC_RND)
842	&& (tmp.sigh & 0xc0000000)) {
843	/* We may need to subtract st(1) once more,
844	to get a result <= 1/2 of st(1). */
845	unsigned long long x;
846	expdif =
847	exponent16(&st1) - exponent16(&tmp);
848	if (expdif <= 1) {
849	if (expdif == 0)
850	x = significand(&st1) -
851	significand(&tmp);
852	else /* expdif is 1 */
853	x = (significand(&st1)
854	<< 1) -
855	significand(&tmp);
856	if ((x < significand(&tmp)) ||
857	/* or equi-distant (from 0 & st(1)) and q is odd */
858	((x == significand(&tmp))
859	&& (q & 1))) {
860	st0_sign = !st0_sign;
861	significand(&tmp) = x;
862	q++;
863	}
864	}
865	}
866
867	if (q & 4)
868	cc |= SW_C0;
869	if (q & 2)
870	cc |= SW_C3;
871	if (q & 1)
872	cc |= SW_C1;
873	} else {
874	control_word = old_cw;
875	setcc(0);
876	return;
877	}
878	} else {
879	/* There is a large exponent difference ( >= 64 ) */
880	/* To make much sense, the code in this section should
881	be done at high precision. */
882	int exp_1, N;
883	u_char sign;
884
885	/* prevent overflow here */
886	/* N is 'a number between 32 and 63' (p26-113) */
887	reg_copy(x: &st0, y: &tmp);
888	tmptag = st0_tag;
889	N = (expdif & 0x0000001f) + 32; /* This choice gives results
890	identical to an AMD 486 */
891	setexponent16(&tmp, N);
892	exp_1 = exponent16(&st1);
893	setexponent16(&st1, 0);
894	expdif -= N;
895
896	sign = getsign(&tmp) ^ st1_sign;
897	tag =
898	FPU_u_div(arg1: &tmp, arg2: &st1, answ: &tmp,
899	PR_64_BITS | RC_CHOP | 0x3f, sign);
900	setsign(&tmp, sign);
901
902	FPU_round_to_int(r: &tmp, tag); /* Fortunately, this can't
903	overflow to 2^64 */
904
905	rem_kernel(significand(&st0),
906	y: &significand(&tmp),
907	significand(&st1),
908	significand(&tmp), exponent(&tmp)
909	);
910	setexponent16(&tmp, exp_1 + expdif);
911
912	/* It is possible for the operation to be complete here.
913	What does the IEEE standard say? The Intel 80486 manual
914	implies that the operation will never be completed at this
915	point, and the behaviour of a real 80486 confirms this.
916	*/
917	if (!(tmp.sigh | tmp.sigl)) {
918	/* The result is zero */
919	control_word = old_cw;
920	partial_status = saved_status;
921	FPU_copy_to_reg0(r: &CONST_Z, TAG_Zero);
922	setsign(&st0, st0_sign);
923	#ifdef PECULIAR_486
924	setcc(SW_C2);
925	#else
926	setcc(0);
927	#endif /* PECULIAR_486 */
928	return;
929	}
930	cc = SW_C2;
931	}
932
933	control_word = old_cw;
934	partial_status = saved_status;
935	tag = FPU_normalize_nuo(x: &tmp);
936	reg_copy(x: &tmp, y: st0_ptr);
937
938	/* The only condition to be looked for is underflow,
939	and it can occur here only if underflow is unmasked. */
940	if ((exponent16(&tmp) <= EXP_UNDER) && (tag != TAG_Zero)
941	&& !(control_word & CW_Underflow)) {
942	setcc(cc);
943	tag = arith_underflow(dest: st0_ptr);
944	setsign(st0_ptr, st0_sign);
945	FPU_settag0(tag);
946	return;
947	} else if ((exponent16(&tmp) > EXP_UNDER) || (tag == TAG_Zero)) {
948	stdexp(st0_ptr);
949	setsign(st0_ptr, st0_sign);
950	} else {
951	tag =
952	FPU_round(arg: st0_ptr, extent: 0, dummy: 0, FULL_PRECISION, sign: st0_sign);
953	}
954	FPU_settag0(tag);
955	setcc(cc);
956
957	return;
958	}
959
960	if (st0_tag == TAG_Special)
961	st0_tag = FPU_Special(ptr: st0_ptr);
962	if (st1_tag == TAG_Special)
963	st1_tag = FPU_Special(ptr: st1_ptr);
964
965	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
966	|| ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
967	|| ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
968	if (denormal_operand() < 0)
969	return;
970	goto fprem_valid;
971	} else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
972	FPU_stack_underflow();
973	return;
974	} else if (st0_tag == TAG_Zero) {
975	if (st1_tag == TAG_Valid) {
976	setcc(0);
977	return;
978	} else if (st1_tag == TW_Denormal) {
979	if (denormal_operand() < 0)
980	return;
981	setcc(0);
982	return;
983	} else if (st1_tag == TAG_Zero) {
984	arith_invalid(deststnr: 0);
985	return;
986	} /* fprem(?,0) always invalid */
987	else if (st1_tag == TW_Infinity) {
988	setcc(0);
989	return;
990	}
991	} else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
992	if (st1_tag == TAG_Zero) {
993	arith_invalid(deststnr: 0); /* fprem(Valid,Zero) is invalid */
994	return;
995	} else if (st1_tag != TW_NaN) {
996	if (((st0_tag == TW_Denormal)
997	|| (st1_tag == TW_Denormal))
998	&& (denormal_operand() < 0))
999	return;
1000
1001	if (st1_tag == TW_Infinity) {
1002	/* fprem(Valid,Infinity) is o.k. */
1003	setcc(0);
1004	return;
1005	}
1006	}
1007	} else if (st0_tag == TW_Infinity) {
1008	if (st1_tag != TW_NaN) {
1009	arith_invalid(deststnr: 0); /* fprem(Infinity,?) is invalid */
1010	return;
1011	}
1012	}
1013
1014	/* One of the registers must contain a NaN if we got here. */
1015
1016	#ifdef PARANOID
1017	if ((st0_tag != TW_NaN) && (st1_tag != TW_NaN))
1018	EXCEPTION(EX_INTERNAL | 0x118);
1019	#endif /* PARANOID */
1020
1021	real_2op_NaN(b: st1_ptr, tagb: st1_tag, deststnr: 0, defaultNaN: st1_ptr);
1022
1023	}
1024
1025	/* ST(1) <- ST(1) * log ST; pop ST */
1026	static void fyl2x(FPU_REG *st0_ptr, u_char st0_tag)
1027	{
1028	FPU_REG *st1_ptr = &st(1), exponent;
1029	u_char st1_tag = FPU_gettagi(stnr: 1);
1030	u_char sign;
1031	int e, tag;
1032
1033	clear_C1();
1034
1035	if ((st0_tag == TAG_Valid) && (st1_tag == TAG_Valid)) {
1036	both_valid:
1037	/* Both regs are Valid or Denormal */
1038	if (signpositive(st0_ptr)) {
1039	if (st0_tag == TW_Denormal)
1040	FPU_to_exp16(a: st0_ptr, x: st0_ptr);
1041	else
1042	/* Convert st(0) for internal use. */
1043	setexponent16(st0_ptr, exponent(st0_ptr));
1044
1045	if ((st0_ptr->sigh == 0x80000000)
1046	&& (st0_ptr->sigl == 0)) {
1047	/* Special case. The result can be precise. */
1048	u_char esign;
1049	e = exponent16(st0_ptr);
1050	if (e >= 0) {
1051	exponent.sigh = e;
1052	esign = SIGN_POS;
1053	} else {
1054	exponent.sigh = -e;
1055	esign = SIGN_NEG;
1056	}
1057	exponent.sigl = 0;
1058	setexponent16(&exponent, 31);
1059	tag = FPU_normalize_nuo(x: &exponent);
1060	stdexp(&exponent);
1061	setsign(&exponent, esign);
1062	tag =
1063	FPU_mul(b: &exponent, tagb: tag, deststnr: 1, FULL_PRECISION);
1064	if (tag >= 0)
1065	FPU_settagi(stnr: 1, tag);
1066	} else {
1067	/* The usual case */
1068	sign = getsign(st1_ptr);
1069	if (st1_tag == TW_Denormal)
1070	FPU_to_exp16(a: st1_ptr, x: st1_ptr);
1071	else
1072	/* Convert st(1) for internal use. */
1073	setexponent16(st1_ptr,
1074	exponent(st1_ptr));
1075	poly_l2(st0_ptr, st1_ptr, st1_sign: sign);
1076	}
1077	} else {
1078	/* negative */
1079	if (arith_invalid(deststnr: 1) < 0)
1080	return;
1081	}
1082
1083	FPU_pop();
1084
1085	return;
1086	}
1087
1088	if (st0_tag == TAG_Special)
1089	st0_tag = FPU_Special(ptr: st0_ptr);
1090	if (st1_tag == TAG_Special)
1091	st1_tag = FPU_Special(ptr: st1_ptr);
1092
1093	if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1094	FPU_stack_underflow_pop(i: 1);
1095	return;
1096	} else if ((st0_tag <= TW_Denormal) && (st1_tag <= TW_Denormal)) {
1097	if (st0_tag == TAG_Zero) {
1098	if (st1_tag == TAG_Zero) {
1099	/* Both args zero is invalid */
1100	if (arith_invalid(deststnr: 1) < 0)
1101	return;
1102	} else {
1103	u_char sign;
1104	sign = getsign(st1_ptr) ^ SIGN_NEG;
1105	if (FPU_divide_by_zero(deststnr: 1, sign) < 0)
1106	return;
1107
1108	setsign(st1_ptr, sign);
1109	}
1110	} else if (st1_tag == TAG_Zero) {
1111	/* st(1) contains zero, st(0) valid <> 0 */
1112	/* Zero is the valid answer */
1113	sign = getsign(st1_ptr);
1114
1115	if (signnegative(st0_ptr)) {
1116	/* log(negative) */
1117	if (arith_invalid(deststnr: 1) < 0)
1118	return;
1119	} else if ((st0_tag == TW_Denormal)
1120	&& (denormal_operand() < 0))
1121	return;
1122	else {
1123	if (exponent(st0_ptr) < 0)
1124	sign ^= SIGN_NEG;
1125
1126	FPU_copy_to_reg1(r: &CONST_Z, TAG_Zero);
1127	setsign(st1_ptr, sign);
1128	}
1129	} else {
1130	/* One or both operands are denormals. */
1131	if (denormal_operand() < 0)
1132	return;
1133	goto both_valid;
1134	}
1135	} else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1136	if (real_2op_NaN(b: st0_ptr, tagb: st0_tag, deststnr: 1, defaultNaN: st0_ptr) < 0)
1137	return;
1138	}
1139	/* One or both arg must be an infinity */
1140	else if (st0_tag == TW_Infinity) {
1141	if ((signnegative(st0_ptr)) || (st1_tag == TAG_Zero)) {
1142	/* log(-infinity) or 0*log(infinity) */
1143	if (arith_invalid(deststnr: 1) < 0)
1144	return;
1145	} else {
1146	u_char sign = getsign(st1_ptr);
1147
1148	if ((st1_tag == TW_Denormal)
1149	&& (denormal_operand() < 0))
1150	return;
1151
1152	FPU_copy_to_reg1(r: &CONST_INF, TAG_Special);
1153	setsign(st1_ptr, sign);
1154	}
1155	}
1156	/* st(1) must be infinity here */
1157	else if (((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal))
1158	&& (signpositive(st0_ptr))) {
1159	if (exponent(st0_ptr) >= 0) {
1160	if ((exponent(st0_ptr) == 0) &&
1161	(st0_ptr->sigh == 0x80000000) &&
1162	(st0_ptr->sigl == 0)) {
1163	/* st(0) holds 1.0 */
1164	/* infinity*log(1) */
1165	if (arith_invalid(deststnr: 1) < 0)
1166	return;
1167	}
1168	/* else st(0) is positive and > 1.0 */
1169	} else {
1170	/* st(0) is positive and < 1.0 */
1171
1172	if ((st0_tag == TW_Denormal)
1173	&& (denormal_operand() < 0))
1174	return;
1175
1176	changesign(st1_ptr);
1177	}
1178	} else {
1179	/* st(0) must be zero or negative */
1180	if (st0_tag == TAG_Zero) {
1181	/* This should be invalid, but a real 80486 is happy with it. */
1182
1183	#ifndef PECULIAR_486
1184	sign = getsign(st1_ptr);
1185	if (FPU_divide_by_zero(1, sign) < 0)
1186	return;
1187	#endif /* PECULIAR_486 */
1188
1189	changesign(st1_ptr);
1190	} else if (arith_invalid(deststnr: 1) < 0) /* log(negative) */
1191	return;
1192	}
1193
1194	FPU_pop();
1195	}
1196
1197	static void fpatan(FPU_REG *st0_ptr, u_char st0_tag)
1198	{
1199	FPU_REG *st1_ptr = &st(1);
1200	u_char st1_tag = FPU_gettagi(stnr: 1);
1201	int tag;
1202
1203	clear_C1();
1204	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1205	valid_atan:
1206
1207	poly_atan(st0_ptr, st0_tag, st1_ptr, st1_tag);
1208
1209	FPU_pop();
1210
1211	return;
1212	}
1213
1214	if (st0_tag == TAG_Special)
1215	st0_tag = FPU_Special(ptr: st0_ptr);
1216	if (st1_tag == TAG_Special)
1217	st1_tag = FPU_Special(ptr: st1_ptr);
1218
1219	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1220	|| ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1221	|| ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1222	if (denormal_operand() < 0)
1223	return;
1224
1225	goto valid_atan;
1226	} else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1227	FPU_stack_underflow_pop(i: 1);
1228	return;
1229	} else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1230	if (real_2op_NaN(b: st0_ptr, tagb: st0_tag, deststnr: 1, defaultNaN: st0_ptr) >= 0)
1231	FPU_pop();
1232	return;
1233	} else if ((st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
1234	u_char sign = getsign(st1_ptr);
1235	if (st0_tag == TW_Infinity) {
1236	if (st1_tag == TW_Infinity) {
1237	if (signpositive(st0_ptr)) {
1238	FPU_copy_to_reg1(r: &CONST_PI4, TAG_Valid);
1239	} else {
1240	setpositive(st1_ptr);
1241	tag =
1242	FPU_u_add(arg1: &CONST_PI4, arg2: &CONST_PI2,
1243	answ: st1_ptr, FULL_PRECISION,
1244	SIGN_POS,
1245	exponent(&CONST_PI4),
1246	exponent(&CONST_PI2));
1247	if (tag >= 0)
1248	FPU_settagi(stnr: 1, tag);
1249	}
1250	} else {
1251	if ((st1_tag == TW_Denormal)
1252	&& (denormal_operand() < 0))
1253	return;
1254
1255	if (signpositive(st0_ptr)) {
1256	FPU_copy_to_reg1(r: &CONST_Z, TAG_Zero);
1257	setsign(st1_ptr, sign); /* An 80486 preserves the sign */
1258	FPU_pop();
1259	return;
1260	} else {
1261	FPU_copy_to_reg1(r: &CONST_PI, TAG_Valid);
1262	}
1263	}
1264	} else {
1265	/* st(1) is infinity, st(0) not infinity */
1266	if ((st0_tag == TW_Denormal)
1267	&& (denormal_operand() < 0))
1268	return;
1269
1270	FPU_copy_to_reg1(r: &CONST_PI2, TAG_Valid);
1271	}
1272	setsign(st1_ptr, sign);
1273	} else if (st1_tag == TAG_Zero) {
1274	/* st(0) must be valid or zero */
1275	u_char sign = getsign(st1_ptr);
1276
1277	if ((st0_tag == TW_Denormal) && (denormal_operand() < 0))
1278	return;
1279
1280	if (signpositive(st0_ptr)) {
1281	/* An 80486 preserves the sign */
1282	FPU_pop();
1283	return;
1284	}
1285
1286	FPU_copy_to_reg1(r: &CONST_PI, TAG_Valid);
1287	setsign(st1_ptr, sign);
1288	} else if (st0_tag == TAG_Zero) {
1289	/* st(1) must be TAG_Valid here */
1290	u_char sign = getsign(st1_ptr);
1291
1292	if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1293	return;
1294
1295	FPU_copy_to_reg1(r: &CONST_PI2, TAG_Valid);
1296	setsign(st1_ptr, sign);
1297	}
1298	#ifdef PARANOID
1299	else
1300	EXCEPTION(EX_INTERNAL | 0x125);
1301	#endif /* PARANOID */
1302
1303	FPU_pop();
1304	set_precision_flag_up(); /* We do not really know if up or down */
1305	}
1306
1307	static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1308	{
1309	do_fprem(st0_ptr, st0_tag, RC_CHOP);
1310	}
1311
1312	static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1313	{
1314	do_fprem(st0_ptr, st0_tag, RC_RND);
1315	}
1316
1317	static void fyl2xp1(FPU_REG *st0_ptr, u_char st0_tag)
1318	{
1319	u_char sign, sign1;
1320	FPU_REG *st1_ptr = &st(1), a, b;
1321	u_char st1_tag = FPU_gettagi(stnr: 1);
1322
1323	clear_C1();
1324	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1325	valid_yl2xp1:
1326
1327	sign = getsign(st0_ptr);
1328	sign1 = getsign(st1_ptr);
1329
1330	FPU_to_exp16(a: st0_ptr, x: &a);
1331	FPU_to_exp16(a: st1_ptr, x: &b);
1332
1333	if (poly_l2p1(s0: sign, s1: sign1, r0: &a, r1: &b, d: st1_ptr))
1334	return;
1335
1336	FPU_pop();
1337	return;
1338	}
1339
1340	if (st0_tag == TAG_Special)
1341	st0_tag = FPU_Special(ptr: st0_ptr);
1342	if (st1_tag == TAG_Special)
1343	st1_tag = FPU_Special(ptr: st1_ptr);
1344
1345	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1346	|| ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1347	|| ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1348	if (denormal_operand() < 0)
1349	return;
1350
1351	goto valid_yl2xp1;
1352	} else if ((st0_tag == TAG_Empty) | (st1_tag == TAG_Empty)) {
1353	FPU_stack_underflow_pop(i: 1);
1354	return;
1355	} else if (st0_tag == TAG_Zero) {
1356	switch (st1_tag) {
1357	case TW_Denormal:
1358	if (denormal_operand() < 0)
1359	return;
1360	fallthrough;
1361	case TAG_Zero:
1362	case TAG_Valid:
1363	setsign(st0_ptr, getsign(st0_ptr) ^ getsign(st1_ptr));
1364	FPU_copy_to_reg1(r: st0_ptr, tag: st0_tag);
1365	break;
1366
1367	case TW_Infinity:
1368	/* Infinity*log(1) */
1369	if (arith_invalid(deststnr: 1) < 0)
1370	return;
1371	break;
1372
1373	case TW_NaN:
1374	if (real_2op_NaN(b: st0_ptr, tagb: st0_tag, deststnr: 1, defaultNaN: st0_ptr) < 0)
1375	return;
1376	break;
1377
1378	default:
1379	#ifdef PARANOID
1380	EXCEPTION(EX_INTERNAL | 0x116);
1381	return;
1382	#endif /* PARANOID */
1383	break;
1384	}
1385	} else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1386	switch (st1_tag) {
1387	case TAG_Zero:
1388	if (signnegative(st0_ptr)) {
1389	if (exponent(st0_ptr) >= 0) {
1390	/* st(0) holds <= -1.0 */
1391	#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1392	changesign(st1_ptr);
1393	#else
1394	if (arith_invalid(1) < 0)
1395	return;
1396	#endif /* PECULIAR_486 */
1397	} else if ((st0_tag == TW_Denormal)
1398	&& (denormal_operand() < 0))
1399	return;
1400	else
1401	changesign(st1_ptr);
1402	} else if ((st0_tag == TW_Denormal)
1403	&& (denormal_operand() < 0))
1404	return;
1405	break;
1406
1407	case TW_Infinity:
1408	if (signnegative(st0_ptr)) {
1409	if ((exponent(st0_ptr) >= 0) &&
1410	!((st0_ptr->sigh == 0x80000000) &&
1411	(st0_ptr->sigl == 0))) {
1412	/* st(0) holds < -1.0 */
1413	#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1414	changesign(st1_ptr);
1415	#else
1416	if (arith_invalid(1) < 0)
1417	return;
1418	#endif /* PECULIAR_486 */
1419	} else if ((st0_tag == TW_Denormal)
1420	&& (denormal_operand() < 0))
1421	return;
1422	else
1423	changesign(st1_ptr);
1424	} else if ((st0_tag == TW_Denormal)
1425	&& (denormal_operand() < 0))
1426	return;
1427	break;
1428
1429	case TW_NaN:
1430	if (real_2op_NaN(b: st0_ptr, tagb: st0_tag, deststnr: 1, defaultNaN: st0_ptr) < 0)
1431	return;
1432	}
1433
1434	} else if (st0_tag == TW_NaN) {
1435	if (real_2op_NaN(b: st0_ptr, tagb: st0_tag, deststnr: 1, defaultNaN: st0_ptr) < 0)
1436	return;
1437	} else if (st0_tag == TW_Infinity) {
1438	if (st1_tag == TW_NaN) {
1439	if (real_2op_NaN(b: st0_ptr, tagb: st0_tag, deststnr: 1, defaultNaN: st0_ptr) < 0)
1440	return;
1441	} else if (signnegative(st0_ptr)) {
1442	#ifndef PECULIAR_486
1443	/* This should have higher priority than denormals, but... */
1444	if (arith_invalid(1) < 0) /* log(-infinity) */
1445	return;
1446	#endif /* PECULIAR_486 */
1447	if ((st1_tag == TW_Denormal)
1448	&& (denormal_operand() < 0))
1449	return;
1450	#ifdef PECULIAR_486
1451	/* Denormal operands actually get higher priority */
1452	if (arith_invalid(deststnr: 1) < 0) /* log(-infinity) */
1453	return;
1454	#endif /* PECULIAR_486 */
1455	} else if (st1_tag == TAG_Zero) {
1456	/* log(infinity) */
1457	if (arith_invalid(deststnr: 1) < 0)
1458	return;
1459	}
1460
1461	/* st(1) must be valid here. */
1462
1463	else if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1464	return;
1465
1466	/* The Manual says that log(Infinity) is invalid, but a real
1467	80486 sensibly says that it is o.k. */
1468	else {
1469	u_char sign = getsign(st1_ptr);
1470	FPU_copy_to_reg1(r: &CONST_INF, TAG_Special);
1471	setsign(st1_ptr, sign);
1472	}
1473	}
1474	#ifdef PARANOID
1475	else {
1476	EXCEPTION(EX_INTERNAL | 0x117);
1477	return;
1478	}
1479	#endif /* PARANOID */
1480
1481	FPU_pop();
1482	return;
1483
1484	}
1485
1486	static void fscale(FPU_REG *st0_ptr, u_char st0_tag)
1487	{
1488	FPU_REG *st1_ptr = &st(1);
1489	u_char st1_tag = FPU_gettagi(stnr: 1);
1490	int old_cw = control_word;
1491	u_char sign = getsign(st0_ptr);
1492
1493	clear_C1();
1494	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1495	long scale;
1496	FPU_REG tmp;
1497
1498	/* Convert register for internal use. */
1499	setexponent16(st0_ptr, exponent(st0_ptr));
1500
1501	valid_scale:
1502
1503	if (exponent(st1_ptr) > 30) {
1504	/* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1505
1506	if (signpositive(st1_ptr)) {
1507	EXCEPTION(EX_Overflow);
1508	FPU_copy_to_reg0(r: &CONST_INF, TAG_Special);
1509	} else {
1510	EXCEPTION(EX_Underflow);
1511	FPU_copy_to_reg0(r: &CONST_Z, TAG_Zero);
1512	}
1513	setsign(st0_ptr, sign);
1514	return;
1515	}
1516
1517	control_word &= ~CW_RC;
1518	control_word |= RC_CHOP;
1519	reg_copy(x: st1_ptr, y: &tmp);
1520	FPU_round_to_int(r: &tmp, tag: st1_tag); /* This can never overflow here */
1521	control_word = old_cw;
1522	scale = signnegative(st1_ptr) ? -tmp.sigl : tmp.sigl;
1523	scale += exponent16(st0_ptr);
1524
1525	setexponent16(st0_ptr, scale);
1526
1527	/* Use FPU_round() to properly detect under/overflow etc */
1528	FPU_round(arg: st0_ptr, extent: 0, dummy: 0, control_word, sign);
1529
1530	return;
1531	}
1532
1533	if (st0_tag == TAG_Special)
1534	st0_tag = FPU_Special(ptr: st0_ptr);
1535	if (st1_tag == TAG_Special)
1536	st1_tag = FPU_Special(ptr: st1_ptr);
1537
1538	if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1539	switch (st1_tag) {
1540	case TAG_Valid:
1541	/* st(0) must be a denormal */
1542	if ((st0_tag == TW_Denormal)
1543	&& (denormal_operand() < 0))
1544	return;
1545
1546	FPU_to_exp16(a: st0_ptr, x: st0_ptr); /* Will not be left on stack */
1547	goto valid_scale;
1548
1549	case TAG_Zero:
1550	if (st0_tag == TW_Denormal)
1551	denormal_operand();
1552	return;
1553
1554	case TW_Denormal:
1555	denormal_operand();
1556	return;
1557
1558	case TW_Infinity:
1559	if ((st0_tag == TW_Denormal)
1560	&& (denormal_operand() < 0))
1561	return;
1562
1563	if (signpositive(st1_ptr))
1564	FPU_copy_to_reg0(r: &CONST_INF, TAG_Special);
1565	else
1566	FPU_copy_to_reg0(r: &CONST_Z, TAG_Zero);
1567	setsign(st0_ptr, sign);
1568	return;
1569
1570	case TW_NaN:
1571	real_2op_NaN(b: st1_ptr, tagb: st1_tag, deststnr: 0, defaultNaN: st0_ptr);
1572	return;
1573	}
1574	} else if (st0_tag == TAG_Zero) {
1575	switch (st1_tag) {
1576	case TAG_Valid:
1577	case TAG_Zero:
1578	return;
1579
1580	case TW_Denormal:
1581	denormal_operand();
1582	return;
1583
1584	case TW_Infinity:
1585	if (signpositive(st1_ptr))
1586	arith_invalid(deststnr: 0); /* Zero scaled by +Infinity */
1587	return;
1588
1589	case TW_NaN:
1590	real_2op_NaN(b: st1_ptr, tagb: st1_tag, deststnr: 0, defaultNaN: st0_ptr);
1591	return;
1592	}
1593	} else if (st0_tag == TW_Infinity) {
1594	switch (st1_tag) {
1595	case TAG_Valid:
1596	case TAG_Zero:
1597	return;
1598
1599	case TW_Denormal:
1600	denormal_operand();
1601	return;
1602
1603	case TW_Infinity:
1604	if (signnegative(st1_ptr))
1605	arith_invalid(deststnr: 0); /* Infinity scaled by -Infinity */
1606	return;
1607
1608	case TW_NaN:
1609	real_2op_NaN(b: st1_ptr, tagb: st1_tag, deststnr: 0, defaultNaN: st0_ptr);
1610	return;
1611	}
1612	} else if (st0_tag == TW_NaN) {
1613	if (st1_tag != TAG_Empty) {
1614	real_2op_NaN(b: st1_ptr, tagb: st1_tag, deststnr: 0, defaultNaN: st0_ptr);
1615	return;
1616	}
1617	}
1618	#ifdef PARANOID
1619	if (!((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty))) {
1620	EXCEPTION(EX_INTERNAL | 0x115);
1621	return;
1622	}
1623	#endif
1624
1625	/* At least one of st(0), st(1) must be empty */
1626	FPU_stack_underflow();
1627
1628	}
1629
1630	/*---------------------------------------------------------------------------*/
1631
1632	static FUNC_ST0 const trig_table_a[] = {
1633	f2xm1, fyl2x, fptan, fpatan,
1634	fxtract, fprem1, (FUNC_ST0) fdecstp, (FUNC_ST0) fincstp
1635	};
1636
1637	void FPU_triga(void)
1638	{
1639	(trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1640	}
1641
1642	static FUNC_ST0 const trig_table_b[] = {
1643	fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, fsin, fcos
1644	};
1645
1646	void FPU_trigb(void)
1647	{
1648	(trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1649	}
1650