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