aboutsummaryrefslogtreecommitdiffstats
path: root/arch/mips/math-emu
diff options
context:
space:
mode:
Diffstat (limited to 'arch/mips/math-emu')
-rw-r--r--arch/mips/math-emu/Makefile11
-rw-r--r--arch/mips/math-emu/cp1emu.c1322
-rw-r--r--arch/mips/math-emu/dp_add.c183
-rw-r--r--arch/mips/math-emu/dp_cmp.c67
-rw-r--r--arch/mips/math-emu/dp_div.c157
-rw-r--r--arch/mips/math-emu/dp_fint.c80
-rw-r--r--arch/mips/math-emu/dp_flong.c78
-rw-r--r--arch/mips/math-emu/dp_frexp.c53
-rw-r--r--arch/mips/math-emu/dp_fsp.c74
-rw-r--r--arch/mips/math-emu/dp_logb.c54
-rw-r--r--arch/mips/math-emu/dp_modf.c80
-rw-r--r--arch/mips/math-emu/dp_mul.c177
-rw-r--r--arch/mips/math-emu/dp_scalb.c58
-rw-r--r--arch/mips/math-emu/dp_simple.c84
-rw-r--r--arch/mips/math-emu/dp_sqrt.c165
-rw-r--r--arch/mips/math-emu/dp_sub.c191
-rw-r--r--arch/mips/math-emu/dp_tint.c124
-rw-r--r--arch/mips/math-emu/dp_tlong.c127
-rw-r--r--arch/mips/math-emu/dsemul.c172
-rw-r--r--arch/mips/math-emu/dsemul.h23
-rw-r--r--arch/mips/math-emu/ieee754.c138
-rw-r--r--arch/mips/math-emu/ieee754.h489
-rw-r--r--arch/mips/math-emu/ieee754d.c138
-rw-r--r--arch/mips/math-emu/ieee754dp.c243
-rw-r--r--arch/mips/math-emu/ieee754dp.h83
-rw-r--r--arch/mips/math-emu/ieee754int.h165
-rw-r--r--arch/mips/math-emu/ieee754m.c56
-rw-r--r--arch/mips/math-emu/ieee754sp.c243
-rw-r--r--arch/mips/math-emu/ieee754sp.h89
-rw-r--r--arch/mips/math-emu/ieee754xcpt.c49
-rw-r--r--arch/mips/math-emu/kernel_linkage.c125
-rw-r--r--arch/mips/math-emu/sp_add.c177
-rw-r--r--arch/mips/math-emu/sp_cmp.c67
-rw-r--r--arch/mips/math-emu/sp_div.c157
-rw-r--r--arch/mips/math-emu/sp_fdp.c77
-rw-r--r--arch/mips/math-emu/sp_fint.c80
-rw-r--r--arch/mips/math-emu/sp_flong.c79
-rw-r--r--arch/mips/math-emu/sp_frexp.c53
-rw-r--r--arch/mips/math-emu/sp_logb.c54
-rw-r--r--arch/mips/math-emu/sp_modf.c80
-rw-r--r--arch/mips/math-emu/sp_mul.c171
-rw-r--r--arch/mips/math-emu/sp_scalb.c58
-rw-r--r--arch/mips/math-emu/sp_simple.c84
-rw-r--r--arch/mips/math-emu/sp_sqrt.c117
-rw-r--r--arch/mips/math-emu/sp_sub.c184
-rw-r--r--arch/mips/math-emu/sp_tint.c128
-rw-r--r--arch/mips/math-emu/sp_tlong.c123
47 files changed, 6787 insertions, 0 deletions
diff --git a/arch/mips/math-emu/Makefile b/arch/mips/math-emu/Makefile
new file mode 100644
index 000000000000..121a848a3594
--- /dev/null
+++ b/arch/mips/math-emu/Makefile
@@ -0,0 +1,11 @@
1#
2# Makefile for the Linux/MIPS kernel FPU emulation.
3#
4
5obj-y := cp1emu.o ieee754m.o ieee754d.o ieee754dp.o ieee754sp.o ieee754.o \
6 ieee754xcpt.o dp_frexp.o dp_modf.o dp_div.o dp_mul.o dp_sub.o \
7 dp_add.o dp_fsp.o dp_cmp.o dp_logb.o dp_scalb.o dp_simple.o \
8 dp_tint.o dp_fint.o dp_tlong.o dp_flong.o sp_frexp.o sp_modf.o \
9 sp_div.o sp_mul.o sp_sub.o sp_add.o sp_fdp.o sp_cmp.o sp_logb.o \
10 sp_scalb.o sp_simple.o sp_tint.o sp_fint.o sp_tlong.o sp_flong.o \
11 dp_sqrt.o sp_sqrt.o kernel_linkage.o dsemul.o
diff --git a/arch/mips/math-emu/cp1emu.c b/arch/mips/math-emu/cp1emu.c
new file mode 100644
index 000000000000..20a552be02ee
--- /dev/null
+++ b/arch/mips/math-emu/cp1emu.c
@@ -0,0 +1,1322 @@
1/*
2 * cp1emu.c: a MIPS coprocessor 1 (fpu) instruction emulator
3 *
4 * MIPS floating point support
5 * Copyright (C) 1994-2000 Algorithmics Ltd.
6 * http://www.algor.co.uk
7 *
8 * Kevin D. Kissell, kevink@mips.com and Carsten Langgaard, carstenl@mips.com
9 * Copyright (C) 2000 MIPS Technologies, Inc.
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * A complete emulator for MIPS coprocessor 1 instructions. This is
25 * required for #float(switch) or #float(trap), where it catches all
26 * COP1 instructions via the "CoProcessor Unusable" exception.
27 *
28 * More surprisingly it is also required for #float(ieee), to help out
29 * the hardware fpu at the boundaries of the IEEE-754 representation
30 * (denormalised values, infinities, underflow, etc). It is made
31 * quite nasty because emulation of some non-COP1 instructions is
32 * required, e.g. in branch delay slots.
33 *
34 * Note if you know that you won't have an fpu, then you'll get much
35 * better performance by compiling with -msoft-float!
36 */
37#include <linux/sched.h>
38
39#include <asm/inst.h>
40#include <asm/bootinfo.h>
41#include <asm/cpu.h>
42#include <asm/cpu-features.h>
43#include <asm/processor.h>
44#include <asm/ptrace.h>
45#include <asm/signal.h>
46#include <asm/mipsregs.h>
47#include <asm/fpu_emulator.h>
48#include <asm/uaccess.h>
49#include <asm/branch.h>
50
51#include "ieee754.h"
52#include "dsemul.h"
53
54/* Strap kernel emulator for full MIPS IV emulation */
55
56#ifdef __mips
57#undef __mips
58#endif
59#define __mips 4
60
61/* Function which emulates a floating point instruction. */
62
63static int fpu_emu(struct pt_regs *, struct mips_fpu_soft_struct *,
64 mips_instruction);
65
66#if __mips >= 4 && __mips != 32
67static int fpux_emu(struct pt_regs *,
68 struct mips_fpu_soft_struct *, mips_instruction);
69#endif
70
71/* Further private data for which no space exists in mips_fpu_soft_struct */
72
73struct mips_fpu_emulator_private fpuemuprivate;
74
75/* Control registers */
76
77#define FPCREG_RID 0 /* $0 = revision id */
78#define FPCREG_CSR 31 /* $31 = csr */
79
80/* Convert Mips rounding mode (0..3) to IEEE library modes. */
81static const unsigned char ieee_rm[4] = {
82 IEEE754_RN, IEEE754_RZ, IEEE754_RU, IEEE754_RD
83};
84
85#if __mips >= 4
86/* convert condition code register number to csr bit */
87static const unsigned int fpucondbit[8] = {
88 FPU_CSR_COND0,
89 FPU_CSR_COND1,
90 FPU_CSR_COND2,
91 FPU_CSR_COND3,
92 FPU_CSR_COND4,
93 FPU_CSR_COND5,
94 FPU_CSR_COND6,
95 FPU_CSR_COND7
96};
97#endif
98
99
100/*
101 * Redundant with logic already in kernel/branch.c,
102 * embedded in compute_return_epc. At some point,
103 * a single subroutine should be used across both
104 * modules.
105 */
106static int isBranchInstr(mips_instruction * i)
107{
108 switch (MIPSInst_OPCODE(*i)) {
109 case spec_op:
110 switch (MIPSInst_FUNC(*i)) {
111 case jalr_op:
112 case jr_op:
113 return 1;
114 }
115 break;
116
117 case bcond_op:
118 switch (MIPSInst_RT(*i)) {
119 case bltz_op:
120 case bgez_op:
121 case bltzl_op:
122 case bgezl_op:
123 case bltzal_op:
124 case bgezal_op:
125 case bltzall_op:
126 case bgezall_op:
127 return 1;
128 }
129 break;
130
131 case j_op:
132 case jal_op:
133 case jalx_op:
134 case beq_op:
135 case bne_op:
136 case blez_op:
137 case bgtz_op:
138 case beql_op:
139 case bnel_op:
140 case blezl_op:
141 case bgtzl_op:
142 return 1;
143
144 case cop0_op:
145 case cop1_op:
146 case cop2_op:
147 case cop1x_op:
148 if (MIPSInst_RS(*i) == bc_op)
149 return 1;
150 break;
151 }
152
153 return 0;
154}
155
156/*
157 * In the Linux kernel, we support selection of FPR format on the
158 * basis of the Status.FR bit. This does imply that, if a full 32
159 * FPRs are desired, there needs to be a flip-flop that can be written
160 * to one at that bit position. In any case, O32 MIPS ABI uses
161 * only the even FPRs (Status.FR = 0).
162 */
163
164#define CP0_STATUS_FR_SUPPORT
165
166#ifdef CP0_STATUS_FR_SUPPORT
167#define FR_BIT ST0_FR
168#else
169#define FR_BIT 0
170#endif
171
172#define SIFROMREG(si,x) ((si) = \
173 (xcp->cp0_status & FR_BIT) || !(x & 1) ? \
174 (int)ctx->fpr[x] : \
175 (int)(ctx->fpr[x & ~1] >> 32 ))
176#define SITOREG(si,x) (ctx->fpr[x & ~((xcp->cp0_status & FR_BIT) == 0)] = \
177 (xcp->cp0_status & FR_BIT) || !(x & 1) ? \
178 ctx->fpr[x & ~1] >> 32 << 32 | (u32)(si) : \
179 ctx->fpr[x & ~1] << 32 >> 32 | (u64)(si) << 32)
180
181#define DIFROMREG(di,x) ((di) = \
182 ctx->fpr[x & ~((xcp->cp0_status & FR_BIT) == 0)])
183#define DITOREG(di,x) (ctx->fpr[x & ~((xcp->cp0_status & FR_BIT) == 0)] \
184 = (di))
185
186#define SPFROMREG(sp,x) SIFROMREG((sp).bits,x)
187#define SPTOREG(sp,x) SITOREG((sp).bits,x)
188#define DPFROMREG(dp,x) DIFROMREG((dp).bits,x)
189#define DPTOREG(dp,x) DITOREG((dp).bits,x)
190
191/*
192 * Emulate the single floating point instruction pointed at by EPC.
193 * Two instructions if the instruction is in a branch delay slot.
194 */
195
196static int cop1Emulate(struct pt_regs *xcp, struct mips_fpu_soft_struct *ctx)
197{
198 mips_instruction ir;
199 vaddr_t emulpc, contpc;
200 unsigned int cond;
201
202 if (get_user(ir, (mips_instruction *) xcp->cp0_epc)) {
203 fpuemuprivate.stats.errors++;
204 return SIGBUS;
205 }
206
207 /* XXX NEC Vr54xx bug workaround */
208 if ((xcp->cp0_cause & CAUSEF_BD) && !isBranchInstr(&ir))
209 xcp->cp0_cause &= ~CAUSEF_BD;
210
211 if (xcp->cp0_cause & CAUSEF_BD) {
212 /*
213 * The instruction to be emulated is in a branch delay slot
214 * which means that we have to emulate the branch instruction
215 * BEFORE we do the cop1 instruction.
216 *
217 * This branch could be a COP1 branch, but in that case we
218 * would have had a trap for that instruction, and would not
219 * come through this route.
220 *
221 * Linux MIPS branch emulator operates on context, updating the
222 * cp0_epc.
223 */
224 emulpc = REG_TO_VA(xcp->cp0_epc + 4); /* Snapshot emulation target */
225
226 if (__compute_return_epc(xcp)) {
227#ifdef CP1DBG
228 printk("failed to emulate branch at %p\n",
229 REG_TO_VA(xcp->cp0_epc));
230#endif
231 return SIGILL;
232 }
233 if (get_user(ir, (mips_instruction *) emulpc)) {
234 fpuemuprivate.stats.errors++;
235 return SIGBUS;
236 }
237 /* __compute_return_epc() will have updated cp0_epc */
238 contpc = REG_TO_VA xcp->cp0_epc;
239 /* In order not to confuse ptrace() et al, tweak context */
240 xcp->cp0_epc = VA_TO_REG emulpc - 4;
241 }
242 else {
243 emulpc = REG_TO_VA xcp->cp0_epc;
244 contpc = REG_TO_VA(xcp->cp0_epc + 4);
245 }
246
247 emul:
248 fpuemuprivate.stats.emulated++;
249 switch (MIPSInst_OPCODE(ir)) {
250#ifndef SINGLE_ONLY_FPU
251 case ldc1_op:{
252 u64 *va = REG_TO_VA(xcp->regs[MIPSInst_RS(ir)] +
253 MIPSInst_SIMM(ir));
254 u64 val;
255
256 fpuemuprivate.stats.loads++;
257 if (get_user(val, va)) {
258 fpuemuprivate.stats.errors++;
259 return SIGBUS;
260 }
261 DITOREG(val, MIPSInst_RT(ir));
262 break;
263 }
264
265 case sdc1_op:{
266 u64 *va = REG_TO_VA(xcp->regs[MIPSInst_RS(ir)] +
267 MIPSInst_SIMM(ir));
268 u64 val;
269
270 fpuemuprivate.stats.stores++;
271 DIFROMREG(val, MIPSInst_RT(ir));
272 if (put_user(val, va)) {
273 fpuemuprivate.stats.errors++;
274 return SIGBUS;
275 }
276 break;
277 }
278#endif
279
280 case lwc1_op:{
281 u32 *va = REG_TO_VA(xcp->regs[MIPSInst_RS(ir)] +
282 MIPSInst_SIMM(ir));
283 u32 val;
284
285 fpuemuprivate.stats.loads++;
286 if (get_user(val, va)) {
287 fpuemuprivate.stats.errors++;
288 return SIGBUS;
289 }
290#ifdef SINGLE_ONLY_FPU
291 if (MIPSInst_RT(ir) & 1) {
292 /* illegal register in single-float mode */
293 return SIGILL;
294 }
295#endif
296 SITOREG(val, MIPSInst_RT(ir));
297 break;
298 }
299
300 case swc1_op:{
301 u32 *va = REG_TO_VA(xcp->regs[MIPSInst_RS(ir)] +
302 MIPSInst_SIMM(ir));
303 u32 val;
304
305 fpuemuprivate.stats.stores++;
306#ifdef SINGLE_ONLY_FPU
307 if (MIPSInst_RT(ir) & 1) {
308 /* illegal register in single-float mode */
309 return SIGILL;
310 }
311#endif
312 SIFROMREG(val, MIPSInst_RT(ir));
313 if (put_user(val, va)) {
314 fpuemuprivate.stats.errors++;
315 return SIGBUS;
316 }
317 break;
318 }
319
320 case cop1_op:
321 switch (MIPSInst_RS(ir)) {
322
323#if __mips64 && !defined(SINGLE_ONLY_FPU)
324 case dmfc_op:
325 /* copregister fs -> gpr[rt] */
326 if (MIPSInst_RT(ir) != 0) {
327 DIFROMREG(xcp->regs[MIPSInst_RT(ir)],
328 MIPSInst_RD(ir));
329 }
330 break;
331
332 case dmtc_op:
333 /* copregister fs <- rt */
334 DITOREG(xcp->regs[MIPSInst_RT(ir)], MIPSInst_RD(ir));
335 break;
336#endif
337
338 case mfc_op:
339 /* copregister rd -> gpr[rt] */
340#ifdef SINGLE_ONLY_FPU
341 if (MIPSInst_RD(ir) & 1) {
342 /* illegal register in single-float mode */
343 return SIGILL;
344 }
345#endif
346 if (MIPSInst_RT(ir) != 0) {
347 SIFROMREG(xcp->regs[MIPSInst_RT(ir)],
348 MIPSInst_RD(ir));
349 }
350 break;
351
352 case mtc_op:
353 /* copregister rd <- rt */
354#ifdef SINGLE_ONLY_FPU
355 if (MIPSInst_RD(ir) & 1) {
356 /* illegal register in single-float mode */
357 return SIGILL;
358 }
359#endif
360 SITOREG(xcp->regs[MIPSInst_RT(ir)], MIPSInst_RD(ir));
361 break;
362
363 case cfc_op:{
364 /* cop control register rd -> gpr[rt] */
365 u32 value;
366
367 if (ir == CP1UNDEF) {
368 return do_dsemulret(xcp);
369 }
370 if (MIPSInst_RD(ir) == FPCREG_CSR) {
371 value = ctx->fcr31;
372#ifdef CSRTRACE
373 printk("%p gpr[%d]<-csr=%08x\n",
374 REG_TO_VA(xcp->cp0_epc),
375 MIPSInst_RT(ir), value);
376#endif
377 }
378 else if (MIPSInst_RD(ir) == FPCREG_RID)
379 value = 0;
380 else
381 value = 0;
382 if (MIPSInst_RT(ir))
383 xcp->regs[MIPSInst_RT(ir)] = value;
384 break;
385 }
386
387 case ctc_op:{
388 /* copregister rd <- rt */
389 u32 value;
390
391 if (MIPSInst_RT(ir) == 0)
392 value = 0;
393 else
394 value = xcp->regs[MIPSInst_RT(ir)];
395
396 /* we only have one writable control reg
397 */
398 if (MIPSInst_RD(ir) == FPCREG_CSR) {
399#ifdef CSRTRACE
400 printk("%p gpr[%d]->csr=%08x\n",
401 REG_TO_VA(xcp->cp0_epc),
402 MIPSInst_RT(ir), value);
403#endif
404 ctx->fcr31 = value;
405 /* copy new rounding mode and
406 flush bit to ieee library state! */
407 ieee754_csr.nod = (ctx->fcr31 & 0x1000000) != 0;
408 ieee754_csr.rm = ieee_rm[value & 0x3];
409 }
410 if ((ctx->fcr31 >> 5) & ctx->fcr31 & FPU_CSR_ALL_E) {
411 return SIGFPE;
412 }
413 break;
414 }
415
416 case bc_op:{
417 int likely = 0;
418
419 if (xcp->cp0_cause & CAUSEF_BD)
420 return SIGILL;
421
422#if __mips >= 4
423 cond = ctx->fcr31 & fpucondbit[MIPSInst_RT(ir) >> 2];
424#else
425 cond = ctx->fcr31 & FPU_CSR_COND;
426#endif
427 switch (MIPSInst_RT(ir) & 3) {
428 case bcfl_op:
429 likely = 1;
430 case bcf_op:
431 cond = !cond;
432 break;
433 case bctl_op:
434 likely = 1;
435 case bct_op:
436 break;
437 default:
438 /* thats an illegal instruction */
439 return SIGILL;
440 }
441
442 xcp->cp0_cause |= CAUSEF_BD;
443 if (cond) {
444 /* branch taken: emulate dslot
445 * instruction
446 */
447 xcp->cp0_epc += 4;
448 contpc = REG_TO_VA
449 (xcp->cp0_epc +
450 (MIPSInst_SIMM(ir) << 2));
451
452 if (get_user(ir, (mips_instruction *)
453 REG_TO_VA xcp->cp0_epc)) {
454 fpuemuprivate.stats.errors++;
455 return SIGBUS;
456 }
457
458 switch (MIPSInst_OPCODE(ir)) {
459 case lwc1_op:
460 case swc1_op:
461#if (__mips >= 2 || __mips64) && !defined(SINGLE_ONLY_FPU)
462 case ldc1_op:
463 case sdc1_op:
464#endif
465 case cop1_op:
466#if __mips >= 4 && __mips != 32
467 case cop1x_op:
468#endif
469 /* its one of ours */
470 goto emul;
471#if __mips >= 4
472 case spec_op:
473 if (MIPSInst_FUNC(ir) == movc_op)
474 goto emul;
475 break;
476#endif
477 }
478
479 /*
480 * Single step the non-cp1
481 * instruction in the dslot
482 */
483 return mips_dsemul(xcp, ir, VA_TO_REG contpc);
484 }
485 else {
486 /* branch not taken */
487 if (likely) {
488 /*
489 * branch likely nullifies
490 * dslot if not taken
491 */
492 xcp->cp0_epc += 4;
493 contpc += 4;
494 /*
495 * else continue & execute
496 * dslot as normal insn
497 */
498 }
499 }
500 break;
501 }
502
503 default:
504 if (!(MIPSInst_RS(ir) & 0x10))
505 return SIGILL;
506 {
507 int sig;
508
509 /* a real fpu computation instruction */
510 if ((sig = fpu_emu(xcp, ctx, ir)))
511 return sig;
512 }
513 }
514 break;
515
516#if __mips >= 4 && __mips != 32
517 case cop1x_op:{
518 int sig;
519
520 if ((sig = fpux_emu(xcp, ctx, ir)))
521 return sig;
522 break;
523 }
524#endif
525
526#if __mips >= 4
527 case spec_op:
528 if (MIPSInst_FUNC(ir) != movc_op)
529 return SIGILL;
530 cond = fpucondbit[MIPSInst_RT(ir) >> 2];
531 if (((ctx->fcr31 & cond) != 0) == ((MIPSInst_RT(ir) & 1) != 0))
532 xcp->regs[MIPSInst_RD(ir)] =
533 xcp->regs[MIPSInst_RS(ir)];
534 break;
535#endif
536
537 default:
538 return SIGILL;
539 }
540
541 /* we did it !! */
542 xcp->cp0_epc = VA_TO_REG(contpc);
543 xcp->cp0_cause &= ~CAUSEF_BD;
544 return 0;
545}
546
547/*
548 * Conversion table from MIPS compare ops 48-63
549 * cond = ieee754dp_cmp(x,y,IEEE754_UN,sig);
550 */
551static const unsigned char cmptab[8] = {
552 0, /* cmp_0 (sig) cmp_sf */
553 IEEE754_CUN, /* cmp_un (sig) cmp_ngle */
554 IEEE754_CEQ, /* cmp_eq (sig) cmp_seq */
555 IEEE754_CEQ | IEEE754_CUN, /* cmp_ueq (sig) cmp_ngl */
556 IEEE754_CLT, /* cmp_olt (sig) cmp_lt */
557 IEEE754_CLT | IEEE754_CUN, /* cmp_ult (sig) cmp_nge */
558 IEEE754_CLT | IEEE754_CEQ, /* cmp_ole (sig) cmp_le */
559 IEEE754_CLT | IEEE754_CEQ | IEEE754_CUN, /* cmp_ule (sig) cmp_ngt */
560};
561
562
563#if __mips >= 4 && __mips != 32
564
565/*
566 * Additional MIPS4 instructions
567 */
568
569#define DEF3OP(name, p, f1, f2, f3) \
570static ieee754##p fpemu_##p##_##name (ieee754##p r, ieee754##p s, \
571 ieee754##p t) \
572{ \
573 struct ieee754_csr ieee754_csr_save; \
574 s = f1 (s, t); \
575 ieee754_csr_save = ieee754_csr; \
576 s = f2 (s, r); \
577 ieee754_csr_save.cx |= ieee754_csr.cx; \
578 ieee754_csr_save.sx |= ieee754_csr.sx; \
579 s = f3 (s); \
580 ieee754_csr.cx |= ieee754_csr_save.cx; \
581 ieee754_csr.sx |= ieee754_csr_save.sx; \
582 return s; \
583}
584
585static ieee754dp fpemu_dp_recip(ieee754dp d)
586{
587 return ieee754dp_div(ieee754dp_one(0), d);
588}
589
590static ieee754dp fpemu_dp_rsqrt(ieee754dp d)
591{
592 return ieee754dp_div(ieee754dp_one(0), ieee754dp_sqrt(d));
593}
594
595static ieee754sp fpemu_sp_recip(ieee754sp s)
596{
597 return ieee754sp_div(ieee754sp_one(0), s);
598}
599
600static ieee754sp fpemu_sp_rsqrt(ieee754sp s)
601{
602 return ieee754sp_div(ieee754sp_one(0), ieee754sp_sqrt(s));
603}
604
605DEF3OP(madd, sp, ieee754sp_mul, ieee754sp_add,);
606DEF3OP(msub, sp, ieee754sp_mul, ieee754sp_sub,);
607DEF3OP(nmadd, sp, ieee754sp_mul, ieee754sp_add, ieee754sp_neg);
608DEF3OP(nmsub, sp, ieee754sp_mul, ieee754sp_sub, ieee754sp_neg);
609DEF3OP(madd, dp, ieee754dp_mul, ieee754dp_add,);
610DEF3OP(msub, dp, ieee754dp_mul, ieee754dp_sub,);
611DEF3OP(nmadd, dp, ieee754dp_mul, ieee754dp_add, ieee754dp_neg);
612DEF3OP(nmsub, dp, ieee754dp_mul, ieee754dp_sub, ieee754dp_neg);
613
614static int fpux_emu(struct pt_regs *xcp, struct mips_fpu_soft_struct *ctx,
615 mips_instruction ir)
616{
617 unsigned rcsr = 0; /* resulting csr */
618
619 fpuemuprivate.stats.cp1xops++;
620
621 switch (MIPSInst_FMA_FFMT(ir)) {
622 case s_fmt:{ /* 0 */
623
624 ieee754sp(*handler) (ieee754sp, ieee754sp, ieee754sp);
625 ieee754sp fd, fr, fs, ft;
626 u32 *va;
627 u32 val;
628
629 switch (MIPSInst_FUNC(ir)) {
630 case lwxc1_op:
631 va = REG_TO_VA(xcp->regs[MIPSInst_FR(ir)] +
632 xcp->regs[MIPSInst_FT(ir)]);
633
634 fpuemuprivate.stats.loads++;
635 if (get_user(val, va)) {
636 fpuemuprivate.stats.errors++;
637 return SIGBUS;
638 }
639#ifdef SINGLE_ONLY_FPU
640 if (MIPSInst_FD(ir) & 1) {
641 /* illegal register in single-float
642 * mode
643 */
644 return SIGILL;
645 }
646#endif
647 SITOREG(val, MIPSInst_FD(ir));
648 break;
649
650 case swxc1_op:
651 va = REG_TO_VA(xcp->regs[MIPSInst_FR(ir)] +
652 xcp->regs[MIPSInst_FT(ir)]);
653
654 fpuemuprivate.stats.stores++;
655#ifdef SINGLE_ONLY_FPU
656 if (MIPSInst_FS(ir) & 1) {
657 /* illegal register in single-float
658 * mode
659 */
660 return SIGILL;
661 }
662#endif
663
664 SIFROMREG(val, MIPSInst_FS(ir));
665 if (put_user(val, va)) {
666 fpuemuprivate.stats.errors++;
667 return SIGBUS;
668 }
669 break;
670
671 case madd_s_op:
672 handler = fpemu_sp_madd;
673 goto scoptop;
674 case msub_s_op:
675 handler = fpemu_sp_msub;
676 goto scoptop;
677 case nmadd_s_op:
678 handler = fpemu_sp_nmadd;
679 goto scoptop;
680 case nmsub_s_op:
681 handler = fpemu_sp_nmsub;
682 goto scoptop;
683
684 scoptop:
685 SPFROMREG(fr, MIPSInst_FR(ir));
686 SPFROMREG(fs, MIPSInst_FS(ir));
687 SPFROMREG(ft, MIPSInst_FT(ir));
688 fd = (*handler) (fr, fs, ft);
689 SPTOREG(fd, MIPSInst_FD(ir));
690
691 copcsr:
692 if (ieee754_cxtest(IEEE754_INEXACT))
693 rcsr |= FPU_CSR_INE_X | FPU_CSR_INE_S;
694 if (ieee754_cxtest(IEEE754_UNDERFLOW))
695 rcsr |= FPU_CSR_UDF_X | FPU_CSR_UDF_S;
696 if (ieee754_cxtest(IEEE754_OVERFLOW))
697 rcsr |= FPU_CSR_OVF_X | FPU_CSR_OVF_S;
698 if (ieee754_cxtest(IEEE754_INVALID_OPERATION))
699 rcsr |= FPU_CSR_INV_X | FPU_CSR_INV_S;
700
701 ctx->fcr31 = (ctx->fcr31 & ~FPU_CSR_ALL_X) | rcsr;
702 if (ieee754_csr.nod)
703 ctx->fcr31 |= 0x1000000;
704 if ((ctx->fcr31 >> 5) & ctx->fcr31 & FPU_CSR_ALL_E) {
705 /*printk ("SIGFPE: fpu csr = %08x\n",
706 ctx->fcr31); */
707 return SIGFPE;
708 }
709
710 break;
711
712 default:
713 return SIGILL;
714 }
715 break;
716 }
717
718#ifndef SINGLE_ONLY_FPU
719 case d_fmt:{ /* 1 */
720 ieee754dp(*handler) (ieee754dp, ieee754dp, ieee754dp);
721 ieee754dp fd, fr, fs, ft;
722 u64 *va;
723 u64 val;
724
725 switch (MIPSInst_FUNC(ir)) {
726 case ldxc1_op:
727 va = REG_TO_VA(xcp->regs[MIPSInst_FR(ir)] +
728 xcp->regs[MIPSInst_FT(ir)]);
729
730 fpuemuprivate.stats.loads++;
731 if (get_user(val, va)) {
732 fpuemuprivate.stats.errors++;
733 return SIGBUS;
734 }
735 DITOREG(val, MIPSInst_FD(ir));
736 break;
737
738 case sdxc1_op:
739 va = REG_TO_VA(xcp->regs[MIPSInst_FR(ir)] +
740 xcp->regs[MIPSInst_FT(ir)]);
741
742 fpuemuprivate.stats.stores++;
743 DIFROMREG(val, MIPSInst_FS(ir));
744 if (put_user(val, va)) {
745 fpuemuprivate.stats.errors++;
746 return SIGBUS;
747 }
748 break;
749
750 case madd_d_op:
751 handler = fpemu_dp_madd;
752 goto dcoptop;
753 case msub_d_op:
754 handler = fpemu_dp_msub;
755 goto dcoptop;
756 case nmadd_d_op:
757 handler = fpemu_dp_nmadd;
758 goto dcoptop;
759 case nmsub_d_op:
760 handler = fpemu_dp_nmsub;
761 goto dcoptop;
762
763 dcoptop:
764 DPFROMREG(fr, MIPSInst_FR(ir));
765 DPFROMREG(fs, MIPSInst_FS(ir));
766 DPFROMREG(ft, MIPSInst_FT(ir));
767 fd = (*handler) (fr, fs, ft);
768 DPTOREG(fd, MIPSInst_FD(ir));
769 goto copcsr;
770
771 default:
772 return SIGILL;
773 }
774 break;
775 }
776#endif
777
778 case 0x7: /* 7 */
779 if (MIPSInst_FUNC(ir) != pfetch_op) {
780 return SIGILL;
781 }
782 /* ignore prefx operation */
783 break;
784
785 default:
786 return SIGILL;
787 }
788
789 return 0;
790}
791#endif
792
793
794
795/*
796 * Emulate a single COP1 arithmetic instruction.
797 */
798static int fpu_emu(struct pt_regs *xcp, struct mips_fpu_soft_struct *ctx,
799 mips_instruction ir)
800{
801 int rfmt; /* resulting format */
802 unsigned rcsr = 0; /* resulting csr */
803 unsigned cond;
804 union {
805 ieee754dp d;
806 ieee754sp s;
807 int w;
808#if __mips64
809 s64 l;
810#endif
811 } rv; /* resulting value */
812
813 fpuemuprivate.stats.cp1ops++;
814 switch (rfmt = (MIPSInst_FFMT(ir) & 0xf)) {
815 case s_fmt:{ /* 0 */
816 union {
817 ieee754sp(*b) (ieee754sp, ieee754sp);
818 ieee754sp(*u) (ieee754sp);
819 } handler;
820
821 switch (MIPSInst_FUNC(ir)) {
822 /* binary ops */
823 case fadd_op:
824 handler.b = ieee754sp_add;
825 goto scopbop;
826 case fsub_op:
827 handler.b = ieee754sp_sub;
828 goto scopbop;
829 case fmul_op:
830 handler.b = ieee754sp_mul;
831 goto scopbop;
832 case fdiv_op:
833 handler.b = ieee754sp_div;
834 goto scopbop;
835
836 /* unary ops */
837#if __mips >= 2 || __mips64
838 case fsqrt_op:
839 handler.u = ieee754sp_sqrt;
840 goto scopuop;
841#endif
842#if __mips >= 4 && __mips != 32
843 case frsqrt_op:
844 handler.u = fpemu_sp_rsqrt;
845 goto scopuop;
846 case frecip_op:
847 handler.u = fpemu_sp_recip;
848 goto scopuop;
849#endif
850#if __mips >= 4
851 case fmovc_op:
852 cond = fpucondbit[MIPSInst_FT(ir) >> 2];
853 if (((ctx->fcr31 & cond) != 0) !=
854 ((MIPSInst_FT(ir) & 1) != 0))
855 return 0;
856 SPFROMREG(rv.s, MIPSInst_FS(ir));
857 break;
858 case fmovz_op:
859 if (xcp->regs[MIPSInst_FT(ir)] != 0)
860 return 0;
861 SPFROMREG(rv.s, MIPSInst_FS(ir));
862 break;
863 case fmovn_op:
864 if (xcp->regs[MIPSInst_FT(ir)] == 0)
865 return 0;
866 SPFROMREG(rv.s, MIPSInst_FS(ir));
867 break;
868#endif
869 case fabs_op:
870 handler.u = ieee754sp_abs;
871 goto scopuop;
872 case fneg_op:
873 handler.u = ieee754sp_neg;
874 goto scopuop;
875 case fmov_op:
876 /* an easy one */
877 SPFROMREG(rv.s, MIPSInst_FS(ir));
878 goto copcsr;
879
880 /* binary op on handler */
881 scopbop:
882 {
883 ieee754sp fs, ft;
884
885 SPFROMREG(fs, MIPSInst_FS(ir));
886 SPFROMREG(ft, MIPSInst_FT(ir));
887
888 rv.s = (*handler.b) (fs, ft);
889 goto copcsr;
890 }
891 scopuop:
892 {
893 ieee754sp fs;
894
895 SPFROMREG(fs, MIPSInst_FS(ir));
896 rv.s = (*handler.u) (fs);
897 goto copcsr;
898 }
899 copcsr:
900 if (ieee754_cxtest(IEEE754_INEXACT))
901 rcsr |= FPU_CSR_INE_X | FPU_CSR_INE_S;
902 if (ieee754_cxtest(IEEE754_UNDERFLOW))
903 rcsr |= FPU_CSR_UDF_X | FPU_CSR_UDF_S;
904 if (ieee754_cxtest(IEEE754_OVERFLOW))
905 rcsr |= FPU_CSR_OVF_X | FPU_CSR_OVF_S;
906 if (ieee754_cxtest(IEEE754_ZERO_DIVIDE))
907 rcsr |= FPU_CSR_DIV_X | FPU_CSR_DIV_S;
908 if (ieee754_cxtest(IEEE754_INVALID_OPERATION))
909 rcsr |= FPU_CSR_INV_X | FPU_CSR_INV_S;
910 break;
911
912 /* unary conv ops */
913 case fcvts_op:
914 return SIGILL; /* not defined */
915 case fcvtd_op:{
916#ifdef SINGLE_ONLY_FPU
917 return SIGILL; /* not defined */
918#else
919 ieee754sp fs;
920
921 SPFROMREG(fs, MIPSInst_FS(ir));
922 rv.d = ieee754dp_fsp(fs);
923 rfmt = d_fmt;
924 goto copcsr;
925 }
926#endif
927 case fcvtw_op:{
928 ieee754sp fs;
929
930 SPFROMREG(fs, MIPSInst_FS(ir));
931 rv.w = ieee754sp_tint(fs);
932 rfmt = w_fmt;
933 goto copcsr;
934 }
935
936#if __mips >= 2 || __mips64
937 case fround_op:
938 case ftrunc_op:
939 case fceil_op:
940 case ffloor_op:{
941 unsigned int oldrm = ieee754_csr.rm;
942 ieee754sp fs;
943
944 SPFROMREG(fs, MIPSInst_FS(ir));
945 ieee754_csr.rm = ieee_rm[MIPSInst_FUNC(ir) & 0x3];
946 rv.w = ieee754sp_tint(fs);
947 ieee754_csr.rm = oldrm;
948 rfmt = w_fmt;
949 goto copcsr;
950 }
951#endif /* __mips >= 2 */
952
953#if __mips64 && !defined(SINGLE_ONLY_FPU)
954 case fcvtl_op:{
955 ieee754sp fs;
956
957 SPFROMREG(fs, MIPSInst_FS(ir));
958 rv.l = ieee754sp_tlong(fs);
959 rfmt = l_fmt;
960 goto copcsr;
961 }
962
963 case froundl_op:
964 case ftruncl_op:
965 case fceill_op:
966 case ffloorl_op:{
967 unsigned int oldrm = ieee754_csr.rm;
968 ieee754sp fs;
969
970 SPFROMREG(fs, MIPSInst_FS(ir));
971 ieee754_csr.rm = ieee_rm[MIPSInst_FUNC(ir) & 0x3];
972 rv.l = ieee754sp_tlong(fs);
973 ieee754_csr.rm = oldrm;
974 rfmt = l_fmt;
975 goto copcsr;
976 }
977#endif /* __mips64 && !fpu(single) */
978
979 default:
980 if (MIPSInst_FUNC(ir) >= fcmp_op) {
981 unsigned cmpop = MIPSInst_FUNC(ir) - fcmp_op;
982 ieee754sp fs, ft;
983
984 SPFROMREG(fs, MIPSInst_FS(ir));
985 SPFROMREG(ft, MIPSInst_FT(ir));
986 rv.w = ieee754sp_cmp(fs, ft,
987 cmptab[cmpop & 0x7], cmpop & 0x8);
988 rfmt = -1;
989 if ((cmpop & 0x8) && ieee754_cxtest
990 (IEEE754_INVALID_OPERATION))
991 rcsr = FPU_CSR_INV_X | FPU_CSR_INV_S;
992 else
993 goto copcsr;
994
995 }
996 else {
997 return SIGILL;
998 }
999 break;
1000 }
1001 break;
1002 }
1003
1004#ifndef SINGLE_ONLY_FPU
1005 case d_fmt:{
1006 union {
1007 ieee754dp(*b) (ieee754dp, ieee754dp);
1008 ieee754dp(*u) (ieee754dp);
1009 } handler;
1010
1011 switch (MIPSInst_FUNC(ir)) {
1012 /* binary ops */
1013 case fadd_op:
1014 handler.b = ieee754dp_add;
1015 goto dcopbop;
1016 case fsub_op:
1017 handler.b = ieee754dp_sub;
1018 goto dcopbop;
1019 case fmul_op:
1020 handler.b = ieee754dp_mul;
1021 goto dcopbop;
1022 case fdiv_op:
1023 handler.b = ieee754dp_div;
1024 goto dcopbop;
1025
1026 /* unary ops */
1027#if __mips >= 2 || __mips64
1028 case fsqrt_op:
1029 handler.u = ieee754dp_sqrt;
1030 goto dcopuop;
1031#endif
1032#if __mips >= 4 && __mips != 32
1033 case frsqrt_op:
1034 handler.u = fpemu_dp_rsqrt;
1035 goto dcopuop;
1036 case frecip_op:
1037 handler.u = fpemu_dp_recip;
1038 goto dcopuop;
1039#endif
1040#if __mips >= 4
1041 case fmovc_op:
1042 cond = fpucondbit[MIPSInst_FT(ir) >> 2];
1043 if (((ctx->fcr31 & cond) != 0) !=
1044 ((MIPSInst_FT(ir) & 1) != 0))
1045 return 0;
1046 DPFROMREG(rv.d, MIPSInst_FS(ir));
1047 break;
1048 case fmovz_op:
1049 if (xcp->regs[MIPSInst_FT(ir)] != 0)
1050 return 0;
1051 DPFROMREG(rv.d, MIPSInst_FS(ir));
1052 break;
1053 case fmovn_op:
1054 if (xcp->regs[MIPSInst_FT(ir)] == 0)
1055 return 0;
1056 DPFROMREG(rv.d, MIPSInst_FS(ir));
1057 break;
1058#endif
1059 case fabs_op:
1060 handler.u = ieee754dp_abs;
1061 goto dcopuop;
1062
1063 case fneg_op:
1064 handler.u = ieee754dp_neg;
1065 goto dcopuop;
1066
1067 case fmov_op:
1068 /* an easy one */
1069 DPFROMREG(rv.d, MIPSInst_FS(ir));
1070 goto copcsr;
1071
1072 /* binary op on handler */
1073 dcopbop:{
1074 ieee754dp fs, ft;
1075
1076 DPFROMREG(fs, MIPSInst_FS(ir));
1077 DPFROMREG(ft, MIPSInst_FT(ir));
1078
1079 rv.d = (*handler.b) (fs, ft);
1080 goto copcsr;
1081 }
1082 dcopuop:{
1083 ieee754dp fs;
1084
1085 DPFROMREG(fs, MIPSInst_FS(ir));
1086 rv.d = (*handler.u) (fs);
1087 goto copcsr;
1088 }
1089
1090 /* unary conv ops */
1091 case fcvts_op:{
1092 ieee754dp fs;
1093
1094 DPFROMREG(fs, MIPSInst_FS(ir));
1095 rv.s = ieee754sp_fdp(fs);
1096 rfmt = s_fmt;
1097 goto copcsr;
1098 }
1099 case fcvtd_op:
1100 return SIGILL; /* not defined */
1101
1102 case fcvtw_op:{
1103 ieee754dp fs;
1104
1105 DPFROMREG(fs, MIPSInst_FS(ir));
1106 rv.w = ieee754dp_tint(fs); /* wrong */
1107 rfmt = w_fmt;
1108 goto copcsr;
1109 }
1110
1111#if __mips >= 2 || __mips64
1112 case fround_op:
1113 case ftrunc_op:
1114 case fceil_op:
1115 case ffloor_op:{
1116 unsigned int oldrm = ieee754_csr.rm;
1117 ieee754dp fs;
1118
1119 DPFROMREG(fs, MIPSInst_FS(ir));
1120 ieee754_csr.rm = ieee_rm[MIPSInst_FUNC(ir) & 0x3];
1121 rv.w = ieee754dp_tint(fs);
1122 ieee754_csr.rm = oldrm;
1123 rfmt = w_fmt;
1124 goto copcsr;
1125 }
1126#endif
1127
1128#if __mips64 && !defined(SINGLE_ONLY_FPU)
1129 case fcvtl_op:{
1130 ieee754dp fs;
1131
1132 DPFROMREG(fs, MIPSInst_FS(ir));
1133 rv.l = ieee754dp_tlong(fs);
1134 rfmt = l_fmt;
1135 goto copcsr;
1136 }
1137
1138 case froundl_op:
1139 case ftruncl_op:
1140 case fceill_op:
1141 case ffloorl_op:{
1142 unsigned int oldrm = ieee754_csr.rm;
1143 ieee754dp fs;
1144
1145 DPFROMREG(fs, MIPSInst_FS(ir));
1146 ieee754_csr.rm = ieee_rm[MIPSInst_FUNC(ir) & 0x3];
1147 rv.l = ieee754dp_tlong(fs);
1148 ieee754_csr.rm = oldrm;
1149 rfmt = l_fmt;
1150 goto copcsr;
1151 }
1152#endif /* __mips >= 3 && !fpu(single) */
1153
1154 default:
1155 if (MIPSInst_FUNC(ir) >= fcmp_op) {
1156 unsigned cmpop = MIPSInst_FUNC(ir) - fcmp_op;
1157 ieee754dp fs, ft;
1158
1159 DPFROMREG(fs, MIPSInst_FS(ir));
1160 DPFROMREG(ft, MIPSInst_FT(ir));
1161 rv.w = ieee754dp_cmp(fs, ft,
1162 cmptab[cmpop & 0x7], cmpop & 0x8);
1163 rfmt = -1;
1164 if ((cmpop & 0x8)
1165 &&
1166 ieee754_cxtest
1167 (IEEE754_INVALID_OPERATION))
1168 rcsr = FPU_CSR_INV_X | FPU_CSR_INV_S;
1169 else
1170 goto copcsr;
1171
1172 }
1173 else {
1174 return SIGILL;
1175 }
1176 break;
1177 }
1178 break;
1179 }
1180#endif /* ifndef SINGLE_ONLY_FPU */
1181
1182 case w_fmt:{
1183 ieee754sp fs;
1184
1185 switch (MIPSInst_FUNC(ir)) {
1186 case fcvts_op:
1187 /* convert word to single precision real */
1188 SPFROMREG(fs, MIPSInst_FS(ir));
1189 rv.s = ieee754sp_fint(fs.bits);
1190 rfmt = s_fmt;
1191 goto copcsr;
1192#ifndef SINGLE_ONLY_FPU
1193 case fcvtd_op:
1194 /* convert word to double precision real */
1195 SPFROMREG(fs, MIPSInst_FS(ir));
1196 rv.d = ieee754dp_fint(fs.bits);
1197 rfmt = d_fmt;
1198 goto copcsr;
1199#endif
1200 default:
1201 return SIGILL;
1202 }
1203 break;
1204 }
1205
1206#if __mips64 && !defined(SINGLE_ONLY_FPU)
1207 case l_fmt:{
1208 switch (MIPSInst_FUNC(ir)) {
1209 case fcvts_op:
1210 /* convert long to single precision real */
1211 rv.s = ieee754sp_flong(ctx->fpr[MIPSInst_FS(ir)]);
1212 rfmt = s_fmt;
1213 goto copcsr;
1214 case fcvtd_op:
1215 /* convert long to double precision real */
1216 rv.d = ieee754dp_flong(ctx->fpr[MIPSInst_FS(ir)]);
1217 rfmt = d_fmt;
1218 goto copcsr;
1219 default:
1220 return SIGILL;
1221 }
1222 break;
1223 }
1224#endif
1225
1226 default:
1227 return SIGILL;
1228 }
1229
1230 /*
1231 * Update the fpu CSR register for this operation.
1232 * If an exception is required, generate a tidy SIGFPE exception,
1233 * without updating the result register.
1234 * Note: cause exception bits do not accumulate, they are rewritten
1235 * for each op; only the flag/sticky bits accumulate.
1236 */
1237 ctx->fcr31 = (ctx->fcr31 & ~FPU_CSR_ALL_X) | rcsr;
1238 if ((ctx->fcr31 >> 5) & ctx->fcr31 & FPU_CSR_ALL_E) {
1239 /*printk ("SIGFPE: fpu csr = %08x\n",ctx->fcr31); */
1240 return SIGFPE;
1241 }
1242
1243 /*
1244 * Now we can safely write the result back to the register file.
1245 */
1246 switch (rfmt) {
1247 case -1:{
1248#if __mips >= 4
1249 cond = fpucondbit[MIPSInst_FD(ir) >> 2];
1250#else
1251 cond = FPU_CSR_COND;
1252#endif
1253 if (rv.w)
1254 ctx->fcr31 |= cond;
1255 else
1256 ctx->fcr31 &= ~cond;
1257 break;
1258 }
1259#ifndef SINGLE_ONLY_FPU
1260 case d_fmt:
1261 DPTOREG(rv.d, MIPSInst_FD(ir));
1262 break;
1263#endif
1264 case s_fmt:
1265 SPTOREG(rv.s, MIPSInst_FD(ir));
1266 break;
1267 case w_fmt:
1268 SITOREG(rv.w, MIPSInst_FD(ir));
1269 break;
1270#if __mips64 && !defined(SINGLE_ONLY_FPU)
1271 case l_fmt:
1272 DITOREG(rv.l, MIPSInst_FD(ir));
1273 break;
1274#endif
1275 default:
1276 return SIGILL;
1277 }
1278
1279 return 0;
1280}
1281
1282int fpu_emulator_cop1Handler(int xcptno, struct pt_regs *xcp,
1283 struct mips_fpu_soft_struct *ctx)
1284{
1285 gpreg_t oldepc, prevepc;
1286 mips_instruction insn;
1287 int sig = 0;
1288
1289 oldepc = xcp->cp0_epc;
1290 do {
1291 prevepc = xcp->cp0_epc;
1292
1293 if (get_user(insn, (mips_instruction *) xcp->cp0_epc)) {
1294 fpuemuprivate.stats.errors++;
1295 return SIGBUS;
1296 }
1297 if (insn == 0)
1298 xcp->cp0_epc += 4; /* skip nops */
1299 else {
1300 /* Update ieee754_csr. Only relevant if we have a
1301 h/w FPU */
1302 ieee754_csr.nod = (ctx->fcr31 & 0x1000000) != 0;
1303 ieee754_csr.rm = ieee_rm[ctx->fcr31 & 0x3];
1304 ieee754_csr.cx = (ctx->fcr31 >> 12) & 0x1f;
1305 sig = cop1Emulate(xcp, ctx);
1306 }
1307
1308 if (cpu_has_fpu)
1309 break;
1310 if (sig)
1311 break;
1312
1313 cond_resched();
1314 } while (xcp->cp0_epc > prevepc);
1315
1316 /* SIGILL indicates a non-fpu instruction */
1317 if (sig == SIGILL && xcp->cp0_epc != oldepc)
1318 /* but if epc has advanced, then ignore it */
1319 sig = 0;
1320
1321 return sig;
1322}
diff --git a/arch/mips/math-emu/dp_add.c b/arch/mips/math-emu/dp_add.c
new file mode 100644
index 000000000000..bcf73bb5c33a
--- /dev/null
+++ b/arch/mips/math-emu/dp_add.c
@@ -0,0 +1,183 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 *
26 */
27
28
29#include "ieee754dp.h"
30
31ieee754dp ieee754dp_add(ieee754dp x, ieee754dp y)
32{
33 COMPXDP;
34 COMPYDP;
35
36 EXPLODEXDP;
37 EXPLODEYDP;
38
39 CLEARCX;
40
41 FLUSHXDP;
42 FLUSHYDP;
43
44 switch (CLPAIR(xc, yc)) {
45 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
46 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
47 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
48 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
49 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
50 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
51 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
52 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
53 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
54 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
55 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
56 SETCX(IEEE754_INVALID_OPERATION);
57 return ieee754dp_nanxcpt(ieee754dp_indef(), "add", x, y);
58
59 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
60 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
61 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
62 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
63 return y;
64
65 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
66 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
67 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
68 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
69 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
70 return x;
71
72
73 /* Infinity handling
74 */
75
76 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
77 if (xs == ys)
78 return x;
79 SETCX(IEEE754_INVALID_OPERATION);
80 return ieee754dp_xcpt(ieee754dp_indef(), "add", x, y);
81
82 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
83 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
84 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
85 return y;
86
87 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
88 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
89 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
90 return x;
91
92 /* Zero handling
93 */
94
95 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
96 if (xs == ys)
97 return x;
98 else
99 return ieee754dp_zero(ieee754_csr.rm ==
100 IEEE754_RD);
101
102 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
103 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
104 return x;
105
106 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
107 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
108 return y;
109
110 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
111 DPDNORMX;
112
113 /* FALL THROUGH */
114
115 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
116 DPDNORMY;
117 break;
118
119 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
120 DPDNORMX;
121 break;
122
123 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
124 break;
125 }
126 assert(xm & DP_HIDDEN_BIT);
127 assert(ym & DP_HIDDEN_BIT);
128
129 /* provide guard,round and stick bit space */
130 xm <<= 3;
131 ym <<= 3;
132
133 if (xe > ye) {
134 /* have to shift y fraction right to align
135 */
136 int s = xe - ye;
137 ym = XDPSRS(ym, s);
138 ye += s;
139 } else if (ye > xe) {
140 /* have to shift x fraction right to align
141 */
142 int s = ye - xe;
143 xm = XDPSRS(xm, s);
144 xe += s;
145 }
146 assert(xe == ye);
147 assert(xe <= DP_EMAX);
148
149 if (xs == ys) {
150 /* generate 28 bit result of adding two 27 bit numbers
151 * leaving result in xm,xs,xe
152 */
153 xm = xm + ym;
154 xe = xe;
155 xs = xs;
156
157 if (xm >> (DP_MBITS + 1 + 3)) { /* carry out */
158 xm = XDPSRS1(xm);
159 xe++;
160 }
161 } else {
162 if (xm >= ym) {
163 xm = xm - ym;
164 xe = xe;
165 xs = xs;
166 } else {
167 xm = ym - xm;
168 xe = xe;
169 xs = ys;
170 }
171 if (xm == 0)
172 return ieee754dp_zero(ieee754_csr.rm ==
173 IEEE754_RD);
174
175 /* normalize to rounding precision */
176 while ((xm >> (DP_MBITS + 3)) == 0) {
177 xm <<= 1;
178 xe--;
179 }
180
181 }
182 DPNORMRET2(xs, xe, xm, "add", x, y);
183}
diff --git a/arch/mips/math-emu/dp_cmp.c b/arch/mips/math-emu/dp_cmp.c
new file mode 100644
index 000000000000..8ab4f320a478
--- /dev/null
+++ b/arch/mips/math-emu/dp_cmp.c
@@ -0,0 +1,67 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30int ieee754dp_cmp(ieee754dp x, ieee754dp y, int cmp, int sig)
31{
32 COMPXDP;
33 COMPYDP;
34
35 EXPLODEXDP;
36 EXPLODEYDP;
37 FLUSHXDP;
38 FLUSHYDP;
39 CLEARCX; /* Even clear inexact flag here */
40
41 if (ieee754dp_isnan(x) || ieee754dp_isnan(y)) {
42 if (sig || xc == IEEE754_CLASS_SNAN || yc == IEEE754_CLASS_SNAN)
43 SETCX(IEEE754_INVALID_OPERATION);
44 if (cmp & IEEE754_CUN)
45 return 1;
46 if (cmp & (IEEE754_CLT | IEEE754_CGT)) {
47 if (sig && SETANDTESTCX(IEEE754_INVALID_OPERATION))
48 return ieee754si_xcpt(0, "fcmpf", x);
49 }
50 return 0;
51 } else {
52 s64 vx = x.bits;
53 s64 vy = y.bits;
54
55 if (vx < 0)
56 vx = -vx ^ DP_SIGN_BIT;
57 if (vy < 0)
58 vy = -vy ^ DP_SIGN_BIT;
59
60 if (vx < vy)
61 return (cmp & IEEE754_CLT) != 0;
62 else if (vx == vy)
63 return (cmp & IEEE754_CEQ) != 0;
64 else
65 return (cmp & IEEE754_CGT) != 0;
66 }
67}
diff --git a/arch/mips/math-emu/dp_div.c b/arch/mips/math-emu/dp_div.c
new file mode 100644
index 000000000000..6acedce3b32d
--- /dev/null
+++ b/arch/mips/math-emu/dp_div.c
@@ -0,0 +1,157 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30ieee754dp ieee754dp_div(ieee754dp x, ieee754dp y)
31{
32 COMPXDP;
33 COMPYDP;
34
35 EXPLODEXDP;
36 EXPLODEYDP;
37
38 CLEARCX;
39
40 FLUSHXDP;
41 FLUSHYDP;
42
43 switch (CLPAIR(xc, yc)) {
44 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
45 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
46 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
47 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
48 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
49 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
50 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
51 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
52 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
53 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
54 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754dp_nanxcpt(ieee754dp_indef(), "div", x, y);
57
58 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
59 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
60 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
61 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
62 return y;
63
64 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
65 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
66 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
67 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
68 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
69 return x;
70
71
72 /* Infinity handling
73 */
74
75 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
76 SETCX(IEEE754_INVALID_OPERATION);
77 return ieee754dp_xcpt(ieee754dp_indef(), "div", x, y);
78
79 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
80 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
81 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
82 return ieee754dp_zero(xs ^ ys);
83
84 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
85 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
86 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
87 return ieee754dp_inf(xs ^ ys);
88
89 /* Zero handling
90 */
91
92 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
93 SETCX(IEEE754_INVALID_OPERATION);
94 return ieee754dp_xcpt(ieee754dp_indef(), "div", x, y);
95
96 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
97 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
98 SETCX(IEEE754_ZERO_DIVIDE);
99 return ieee754dp_xcpt(ieee754dp_inf(xs ^ ys), "div", x, y);
100
101 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
102 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
103 return ieee754dp_zero(xs == ys ? 0 : 1);
104
105 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
106 DPDNORMX;
107
108 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
109 DPDNORMY;
110 break;
111
112 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
113 DPDNORMX;
114 break;
115
116 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
117 break;
118 }
119 assert(xm & DP_HIDDEN_BIT);
120 assert(ym & DP_HIDDEN_BIT);
121
122 /* provide rounding space */
123 xm <<= 3;
124 ym <<= 3;
125
126 {
127 /* now the dirty work */
128
129 u64 rm = 0;
130 int re = xe - ye;
131 u64 bm;
132
133 for (bm = DP_MBIT(DP_MBITS + 2); bm; bm >>= 1) {
134 if (xm >= ym) {
135 xm -= ym;
136 rm |= bm;
137 if (xm == 0)
138 break;
139 }
140 xm <<= 1;
141 }
142 rm <<= 1;
143 if (xm)
144 rm |= 1; /* have remainder, set sticky */
145
146 assert(rm);
147
148 /* normalise rm to rounding precision ?
149 */
150 while ((rm >> (DP_MBITS + 3)) == 0) {
151 rm <<= 1;
152 re--;
153 }
154
155 DPNORMRET2(xs == ys ? 0 : 1, re, rm, "div", x, y);
156 }
157}
diff --git a/arch/mips/math-emu/dp_fint.c b/arch/mips/math-emu/dp_fint.c
new file mode 100644
index 000000000000..0065deaee24b
--- /dev/null
+++ b/arch/mips/math-emu/dp_fint.c
@@ -0,0 +1,80 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30ieee754dp ieee754dp_fint(int x)
31{
32 COMPXDP;
33
34 CLEARCX;
35
36 xc = ( 0 ? xc : xc );
37
38 if (x == 0)
39 return ieee754dp_zero(0);
40 if (x == 1 || x == -1)
41 return ieee754dp_one(x < 0);
42 if (x == 10 || x == -10)
43 return ieee754dp_ten(x < 0);
44
45 xs = (x < 0);
46 if (xs) {
47 if (x == (1 << 31))
48 xm = ((unsigned) 1 << 31); /* max neg can't be safely negated */
49 else
50 xm = -x;
51 } else {
52 xm = x;
53 }
54
55#if 1
56 /* normalize - result can never be inexact or overflow */
57 xe = DP_MBITS;
58 while ((xm >> DP_MBITS) == 0) {
59 xm <<= 1;
60 xe--;
61 }
62 return builddp(xs, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
63#else
64 /* normalize */
65 xe = DP_MBITS + 3;
66 while ((xm >> (DP_MBITS + 3)) == 0) {
67 xm <<= 1;
68 xe--;
69 }
70 DPNORMRET1(xs, xe, xm, "fint", x);
71#endif
72}
73
74ieee754dp ieee754dp_funs(unsigned int u)
75{
76 if ((int) u < 0)
77 return ieee754dp_add(ieee754dp_1e31(),
78 ieee754dp_fint(u & ~(1 << 31)));
79 return ieee754dp_fint(u);
80}
diff --git a/arch/mips/math-emu/dp_flong.c b/arch/mips/math-emu/dp_flong.c
new file mode 100644
index 000000000000..cb105b1dd860
--- /dev/null
+++ b/arch/mips/math-emu/dp_flong.c
@@ -0,0 +1,78 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30ieee754dp ieee754dp_flong(s64 x)
31{
32 COMPXDP;
33
34 CLEARCX;
35
36 xc = ( 0 ? xc : xc );
37
38 if (x == 0)
39 return ieee754dp_zero(0);
40 if (x == 1 || x == -1)
41 return ieee754dp_one(x < 0);
42 if (x == 10 || x == -10)
43 return ieee754dp_ten(x < 0);
44
45 xs = (x < 0);
46 if (xs) {
47 if (x == (1ULL << 63))
48 xm = (1ULL << 63); /* max neg can't be safely negated */
49 else
50 xm = -x;
51 } else {
52 xm = x;
53 }
54
55 /* normalize */
56 xe = DP_MBITS + 3;
57 if (xm >> (DP_MBITS + 1 + 3)) {
58 /* shunt out overflow bits */
59 while (xm >> (DP_MBITS + 1 + 3)) {
60 XDPSRSX1();
61 }
62 } else {
63 /* normalize in grs extended double precision */
64 while ((xm >> (DP_MBITS + 3)) == 0) {
65 xm <<= 1;
66 xe--;
67 }
68 }
69 DPNORMRET1(xs, xe, xm, "dp_flong", x);
70}
71
72ieee754dp ieee754dp_fulong(u64 u)
73{
74 if ((s64) u < 0)
75 return ieee754dp_add(ieee754dp_1e63(),
76 ieee754dp_flong(u & ~(1ULL << 63)));
77 return ieee754dp_flong(u);
78}
diff --git a/arch/mips/math-emu/dp_frexp.c b/arch/mips/math-emu/dp_frexp.c
new file mode 100644
index 000000000000..e650cb10c947
--- /dev/null
+++ b/arch/mips/math-emu/dp_frexp.c
@@ -0,0 +1,53 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30/* close to ieeep754dp_logb
31*/
32ieee754dp ieee754dp_frexp(ieee754dp x, int *eptr)
33{
34 COMPXDP;
35 CLEARCX;
36 EXPLODEXDP;
37
38 switch (xc) {
39 case IEEE754_CLASS_SNAN:
40 case IEEE754_CLASS_QNAN:
41 case IEEE754_CLASS_INF:
42 case IEEE754_CLASS_ZERO:
43 *eptr = 0;
44 return x;
45 case IEEE754_CLASS_DNORM:
46 DPDNORMX;
47 break;
48 case IEEE754_CLASS_NORM:
49 break;
50 }
51 *eptr = xe + 1;
52 return builddp(xs, -1 + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
53}
diff --git a/arch/mips/math-emu/dp_fsp.c b/arch/mips/math-emu/dp_fsp.c
new file mode 100644
index 000000000000..494d19ac7049
--- /dev/null
+++ b/arch/mips/math-emu/dp_fsp.c
@@ -0,0 +1,74 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30ieee754dp ieee754dp_fsp(ieee754sp x)
31{
32 COMPXSP;
33
34 EXPLODEXSP;
35
36 CLEARCX;
37
38 FLUSHXSP;
39
40 switch (xc) {
41 case IEEE754_CLASS_SNAN:
42 SETCX(IEEE754_INVALID_OPERATION);
43 return ieee754dp_nanxcpt(ieee754dp_indef(), "fsp");
44 case IEEE754_CLASS_QNAN:
45 return ieee754dp_nanxcpt(builddp(xs,
46 DP_EMAX + 1 + DP_EBIAS,
47 ((u64) xm
48 << (DP_MBITS -
49 SP_MBITS))), "fsp",
50 x);
51 case IEEE754_CLASS_INF:
52 return ieee754dp_inf(xs);
53 case IEEE754_CLASS_ZERO:
54 return ieee754dp_zero(xs);
55 case IEEE754_CLASS_DNORM:
56 /* normalize */
57 while ((xm >> SP_MBITS) == 0) {
58 xm <<= 1;
59 xe--;
60 }
61 break;
62 case IEEE754_CLASS_NORM:
63 break;
64 }
65
66 /* CANT possibly overflow,underflow, or need rounding
67 */
68
69 /* drop the hidden bit */
70 xm &= ~SP_HIDDEN_BIT;
71
72 return builddp(xs, xe + DP_EBIAS,
73 (u64) xm << (DP_MBITS - SP_MBITS));
74}
diff --git a/arch/mips/math-emu/dp_logb.c b/arch/mips/math-emu/dp_logb.c
new file mode 100644
index 000000000000..603388621ca5
--- /dev/null
+++ b/arch/mips/math-emu/dp_logb.c
@@ -0,0 +1,54 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30ieee754dp ieee754dp_logb(ieee754dp x)
31{
32 COMPXDP;
33
34 CLEARCX;
35
36 EXPLODEXDP;
37
38 switch (xc) {
39 case IEEE754_CLASS_SNAN:
40 return ieee754dp_nanxcpt(x, "logb", x);
41 case IEEE754_CLASS_QNAN:
42 return x;
43 case IEEE754_CLASS_INF:
44 return ieee754dp_inf(0);
45 case IEEE754_CLASS_ZERO:
46 return ieee754dp_inf(1);
47 case IEEE754_CLASS_DNORM:
48 DPDNORMX;
49 break;
50 case IEEE754_CLASS_NORM:
51 break;
52 }
53 return ieee754dp_fint(xe);
54}
diff --git a/arch/mips/math-emu/dp_modf.c b/arch/mips/math-emu/dp_modf.c
new file mode 100644
index 000000000000..25861a42c36f
--- /dev/null
+++ b/arch/mips/math-emu/dp_modf.c
@@ -0,0 +1,80 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30/* modf function is always exact for a finite number
31*/
32ieee754dp ieee754dp_modf(ieee754dp x, ieee754dp * ip)
33{
34 COMPXDP;
35
36 CLEARCX;
37
38 EXPLODEXDP;
39
40 switch (xc) {
41 case IEEE754_CLASS_SNAN:
42 case IEEE754_CLASS_QNAN:
43 case IEEE754_CLASS_INF:
44 case IEEE754_CLASS_ZERO:
45 *ip = x;
46 return x;
47 case IEEE754_CLASS_DNORM:
48 /* far to small */
49 *ip = ieee754dp_zero(xs);
50 return x;
51 case IEEE754_CLASS_NORM:
52 break;
53 }
54 if (xe < 0) {
55 *ip = ieee754dp_zero(xs);
56 return x;
57 }
58 if (xe >= DP_MBITS) {
59 *ip = x;
60 return ieee754dp_zero(xs);
61 }
62 /* generate ipart mantissa by clearing bottom bits
63 */
64 *ip = builddp(xs, xe + DP_EBIAS,
65 ((xm >> (DP_MBITS - xe)) << (DP_MBITS - xe)) &
66 ~DP_HIDDEN_BIT);
67
68 /* generate fpart mantissa by clearing top bits
69 * and normalizing (must be able to normalize)
70 */
71 xm = (xm << (64 - (DP_MBITS - xe))) >> (64 - (DP_MBITS - xe));
72 if (xm == 0)
73 return ieee754dp_zero(xs);
74
75 while ((xm >> DP_MBITS) == 0) {
76 xm <<= 1;
77 xe--;
78 }
79 return builddp(xs, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
80}
diff --git a/arch/mips/math-emu/dp_mul.c b/arch/mips/math-emu/dp_mul.c
new file mode 100644
index 000000000000..f2373902f524
--- /dev/null
+++ b/arch/mips/math-emu/dp_mul.c
@@ -0,0 +1,177 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y)
31{
32 COMPXDP;
33 COMPYDP;
34
35 EXPLODEXDP;
36 EXPLODEYDP;
37
38 CLEARCX;
39
40 FLUSHXDP;
41 FLUSHYDP;
42
43 switch (CLPAIR(xc, yc)) {
44 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
45 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
46 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
47 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
48 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
49 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
50 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
51 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
52 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
53 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
54 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754dp_nanxcpt(ieee754dp_indef(), "mul", x, y);
57
58 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
59 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
60 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
61 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
62 return y;
63
64 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
65 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
66 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
67 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
68 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
69 return x;
70
71
72 /* Infinity handling */
73
74 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
75 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
76 SETCX(IEEE754_INVALID_OPERATION);
77 return ieee754dp_xcpt(ieee754dp_indef(), "mul", x, y);
78
79 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
80 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
81 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
82 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
83 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
84 return ieee754dp_inf(xs ^ ys);
85
86 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
87 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
88 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
89 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
90 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
91 return ieee754dp_zero(xs ^ ys);
92
93
94 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
95 DPDNORMX;
96
97 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
98 DPDNORMY;
99 break;
100
101 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
102 DPDNORMX;
103 break;
104
105 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
106 break;
107 }
108 /* rm = xm * ym, re = xe+ye basicly */
109 assert(xm & DP_HIDDEN_BIT);
110 assert(ym & DP_HIDDEN_BIT);
111 {
112 int re = xe + ye;
113 int rs = xs ^ ys;
114 u64 rm;
115
116 /* shunt to top of word */
117 xm <<= 64 - (DP_MBITS + 1);
118 ym <<= 64 - (DP_MBITS + 1);
119
120 /* multiply 32bits xm,ym to give high 32bits rm with stickness
121 */
122
123 /* 32 * 32 => 64 */
124#define DPXMULT(x,y) ((u64)(x) * (u64)y)
125
126 {
127 unsigned lxm = xm;
128 unsigned hxm = xm >> 32;
129 unsigned lym = ym;
130 unsigned hym = ym >> 32;
131 u64 lrm;
132 u64 hrm;
133
134 lrm = DPXMULT(lxm, lym);
135 hrm = DPXMULT(hxm, hym);
136
137 {
138 u64 t = DPXMULT(lxm, hym);
139 {
140 u64 at =
141 lrm + (t << 32);
142 hrm += at < lrm;
143 lrm = at;
144 }
145 hrm = hrm + (t >> 32);
146 }
147
148 {
149 u64 t = DPXMULT(hxm, lym);
150 {
151 u64 at =
152 lrm + (t << 32);
153 hrm += at < lrm;
154 lrm = at;
155 }
156 hrm = hrm + (t >> 32);
157 }
158 rm = hrm | (lrm != 0);
159 }
160
161 /*
162 * sticky shift down to normal rounding precision
163 */
164 if ((s64) rm < 0) {
165 rm =
166 (rm >> (64 - (DP_MBITS + 1 + 3))) |
167 ((rm << (DP_MBITS + 1 + 3)) != 0);
168 re++;
169 } else {
170 rm =
171 (rm >> (64 - (DP_MBITS + 1 + 3 + 1))) |
172 ((rm << (DP_MBITS + 1 + 3 + 1)) != 0);
173 }
174 assert(rm & (DP_HIDDEN_BIT << 3));
175 DPNORMRET2(rs, re, rm, "mul", x, y);
176 }
177}
diff --git a/arch/mips/math-emu/dp_scalb.c b/arch/mips/math-emu/dp_scalb.c
new file mode 100644
index 000000000000..b84e6338330e
--- /dev/null
+++ b/arch/mips/math-emu/dp_scalb.c
@@ -0,0 +1,58 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30ieee754dp ieee754dp_scalb(ieee754dp x, int n)
31{
32 COMPXDP;
33
34 CLEARCX;
35
36 EXPLODEXDP;
37
38 switch (xc) {
39 case IEEE754_CLASS_SNAN:
40 return ieee754dp_nanxcpt(x, "scalb", x, n);
41 case IEEE754_CLASS_QNAN:
42 case IEEE754_CLASS_INF:
43 case IEEE754_CLASS_ZERO:
44 return x;
45 case IEEE754_CLASS_DNORM:
46 DPDNORMX;
47 break;
48 case IEEE754_CLASS_NORM:
49 break;
50 }
51 DPNORMRET2(xs, xe + n, xm << 3, "scalb", x, n);
52}
53
54
55ieee754dp ieee754dp_ldexp(ieee754dp x, int n)
56{
57 return ieee754dp_scalb(x, n);
58}
diff --git a/arch/mips/math-emu/dp_simple.c b/arch/mips/math-emu/dp_simple.c
new file mode 100644
index 000000000000..495c1ac94298
--- /dev/null
+++ b/arch/mips/math-emu/dp_simple.c
@@ -0,0 +1,84 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30int ieee754dp_finite(ieee754dp x)
31{
32 return DPBEXP(x) != DP_EMAX + 1 + DP_EBIAS;
33}
34
35ieee754dp ieee754dp_copysign(ieee754dp x, ieee754dp y)
36{
37 CLEARCX;
38 DPSIGN(x) = DPSIGN(y);
39 return x;
40}
41
42
43ieee754dp ieee754dp_neg(ieee754dp x)
44{
45 COMPXDP;
46
47 EXPLODEXDP;
48 CLEARCX;
49 FLUSHXDP;
50
51 if (xc == IEEE754_CLASS_SNAN) {
52 SETCX(IEEE754_INVALID_OPERATION);
53 return ieee754dp_nanxcpt(ieee754dp_indef(), "neg");
54 }
55
56 if (ieee754dp_isnan(x)) /* but not infinity */
57 return ieee754dp_nanxcpt(x, "neg", x);
58
59 /* quick fix up */
60 DPSIGN(x) ^= 1;
61 return x;
62}
63
64
65ieee754dp ieee754dp_abs(ieee754dp x)
66{
67 COMPXDP;
68
69 EXPLODEXDP;
70 CLEARCX;
71 FLUSHXDP;
72
73 if (xc == IEEE754_CLASS_SNAN) {
74 SETCX(IEEE754_INVALID_OPERATION);
75 return ieee754dp_nanxcpt(ieee754dp_indef(), "neg");
76 }
77
78 if (ieee754dp_isnan(x)) /* but not infinity */
79 return ieee754dp_nanxcpt(x, "abs", x);
80
81 /* quick fix up */
82 DPSIGN(x) = 0;
83 return x;
84}
diff --git a/arch/mips/math-emu/dp_sqrt.c b/arch/mips/math-emu/dp_sqrt.c
new file mode 100644
index 000000000000..c35e871ae975
--- /dev/null
+++ b/arch/mips/math-emu/dp_sqrt.c
@@ -0,0 +1,165 @@
1/* IEEE754 floating point arithmetic
2 * double precision square root
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30static const unsigned table[] = {
31 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592,
32 29598, 36145, 43202, 50740, 58733, 67158, 75992,
33 85215, 83599, 71378, 60428, 50647, 41945, 34246,
34 27478, 21581, 16499, 12183, 8588, 5674, 3403,
35 1742, 661, 130
36};
37
38ieee754dp ieee754dp_sqrt(ieee754dp x)
39{
40 struct ieee754_csr oldcsr;
41 ieee754dp y, z, t;
42 unsigned scalx, yh;
43 COMPXDP;
44
45 EXPLODEXDP;
46 CLEARCX;
47 FLUSHXDP;
48
49 /* x == INF or NAN? */
50 switch (xc) {
51 case IEEE754_CLASS_QNAN:
52 /* sqrt(Nan) = Nan */
53 return ieee754dp_nanxcpt(x, "sqrt");
54 case IEEE754_CLASS_SNAN:
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt");
57 case IEEE754_CLASS_ZERO:
58 /* sqrt(0) = 0 */
59 return x;
60 case IEEE754_CLASS_INF:
61 if (xs) {
62 /* sqrt(-Inf) = Nan */
63 SETCX(IEEE754_INVALID_OPERATION);
64 return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt");
65 }
66 /* sqrt(+Inf) = Inf */
67 return x;
68 case IEEE754_CLASS_DNORM:
69 DPDNORMX;
70 /* fall through */
71 case IEEE754_CLASS_NORM:
72 if (xs) {
73 /* sqrt(-x) = Nan */
74 SETCX(IEEE754_INVALID_OPERATION);
75 return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt");
76 }
77 break;
78 }
79
80 /* save old csr; switch off INX enable & flag; set RN rounding */
81 oldcsr = ieee754_csr;
82 ieee754_csr.mx &= ~IEEE754_INEXACT;
83 ieee754_csr.sx &= ~IEEE754_INEXACT;
84 ieee754_csr.rm = IEEE754_RN;
85
86 /* adjust exponent to prevent overflow */
87 scalx = 0;
88 if (xe > 512) { /* x > 2**-512? */
89 xe -= 512; /* x = x / 2**512 */
90 scalx += 256;
91 } else if (xe < -512) { /* x < 2**-512? */
92 xe += 512; /* x = x * 2**512 */
93 scalx -= 256;
94 }
95
96 y = x = builddp(0, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
97
98 /* magic initial approximation to almost 8 sig. bits */
99 yh = y.bits >> 32;
100 yh = (yh >> 1) + 0x1ff80000;
101 yh = yh - table[(yh >> 15) & 31];
102 y.bits = ((u64) yh << 32) | (y.bits & 0xffffffff);
103
104 /* Heron's rule once with correction to improve to ~18 sig. bits */
105 /* t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; */
106 t = ieee754dp_div(x, y);
107 y = ieee754dp_add(y, t);
108 y.bits -= 0x0010000600000000LL;
109 y.bits &= 0xffffffff00000000LL;
110
111 /* triple to almost 56 sig. bits: y ~= sqrt(x) to within 1 ulp */
112 /* t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; */
113 z = t = ieee754dp_mul(y, y);
114 t.parts.bexp += 0x001;
115 t = ieee754dp_add(t, z);
116 z = ieee754dp_mul(ieee754dp_sub(x, z), y);
117
118 /* t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; */
119 t = ieee754dp_div(z, ieee754dp_add(t, x));
120 t.parts.bexp += 0x001;
121 y = ieee754dp_add(y, t);
122
123 /* twiddle last bit to force y correctly rounded */
124
125 /* set RZ, clear INEX flag */
126 ieee754_csr.rm = IEEE754_RZ;
127 ieee754_csr.sx &= ~IEEE754_INEXACT;
128
129 /* t=x/y; ...chopped quotient, possibly inexact */
130 t = ieee754dp_div(x, y);
131
132 if (ieee754_csr.sx & IEEE754_INEXACT || t.bits != y.bits) {
133
134 if (!(ieee754_csr.sx & IEEE754_INEXACT))
135 /* t = t-ulp */
136 t.bits -= 1;
137
138 /* add inexact to result status */
139 oldcsr.cx |= IEEE754_INEXACT;
140 oldcsr.sx |= IEEE754_INEXACT;
141
142 switch (oldcsr.rm) {
143 case IEEE754_RP:
144 y.bits += 1;
145 /* drop through */
146 case IEEE754_RN:
147 t.bits += 1;
148 break;
149 }
150
151 /* y=y+t; ...chopped sum */
152 y = ieee754dp_add(y, t);
153
154 /* adjust scalx for correctly rounded sqrt(x) */
155 scalx -= 1;
156 }
157
158 /* py[n0]=py[n0]+scalx; ...scale back y */
159 y.parts.bexp += scalx;
160
161 /* restore rounding mode, possibly set inexact */
162 ieee754_csr = oldcsr;
163
164 return y;
165}
diff --git a/arch/mips/math-emu/dp_sub.c b/arch/mips/math-emu/dp_sub.c
new file mode 100644
index 000000000000..b30c5b1f1a2c
--- /dev/null
+++ b/arch/mips/math-emu/dp_sub.c
@@ -0,0 +1,191 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30ieee754dp ieee754dp_sub(ieee754dp x, ieee754dp y)
31{
32 COMPXDP;
33 COMPYDP;
34
35 EXPLODEXDP;
36 EXPLODEYDP;
37
38 CLEARCX;
39
40 FLUSHXDP;
41 FLUSHYDP;
42
43 switch (CLPAIR(xc, yc)) {
44 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
45 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
46 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
47 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
48 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
49 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
50 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
51 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
52 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
53 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
54 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754dp_nanxcpt(ieee754dp_indef(), "sub", x, y);
57
58 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
59 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
60 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
61 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
62 return y;
63
64 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
65 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
66 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
67 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
68 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
69 return x;
70
71
72 /* Infinity handling
73 */
74
75 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
76 if (xs != ys)
77 return x;
78 SETCX(IEEE754_INVALID_OPERATION);
79 return ieee754dp_xcpt(ieee754dp_indef(), "sub", x, y);
80
81 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
82 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
83 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
84 return ieee754dp_inf(ys ^ 1);
85
86 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
87 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
88 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
89 return x;
90
91 /* Zero handling
92 */
93
94 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
95 if (xs != ys)
96 return x;
97 else
98 return ieee754dp_zero(ieee754_csr.rm ==
99 IEEE754_RD);
100
101 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
102 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
103 return x;
104
105 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
106 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
107 /* quick fix up */
108 DPSIGN(y) ^= 1;
109 return y;
110
111 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
112 DPDNORMX;
113 /* FAAL THOROUGH */
114
115 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
116 /* normalize ym,ye */
117 DPDNORMY;
118 break;
119
120 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
121 /* normalize xm,xe */
122 DPDNORMX;
123 break;
124
125 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
126 break;
127 }
128 /* flip sign of y and handle as add */
129 ys ^= 1;
130
131 assert(xm & DP_HIDDEN_BIT);
132 assert(ym & DP_HIDDEN_BIT);
133
134
135 /* provide guard,round and stick bit dpace */
136 xm <<= 3;
137 ym <<= 3;
138
139 if (xe > ye) {
140 /* have to shift y fraction right to align
141 */
142 int s = xe - ye;
143 ym = XDPSRS(ym, s);
144 ye += s;
145 } else if (ye > xe) {
146 /* have to shift x fraction right to align
147 */
148 int s = ye - xe;
149 xm = XDPSRS(xm, s);
150 xe += s;
151 }
152 assert(xe == ye);
153 assert(xe <= DP_EMAX);
154
155 if (xs == ys) {
156 /* generate 28 bit result of adding two 27 bit numbers
157 */
158 xm = xm + ym;
159 xe = xe;
160 xs = xs;
161
162 if (xm >> (DP_MBITS + 1 + 3)) { /* carry out */
163 xm = XDPSRS1(xm); /* shift preserving sticky */
164 xe++;
165 }
166 } else {
167 if (xm >= ym) {
168 xm = xm - ym;
169 xe = xe;
170 xs = xs;
171 } else {
172 xm = ym - xm;
173 xe = xe;
174 xs = ys;
175 }
176 if (xm == 0) {
177 if (ieee754_csr.rm == IEEE754_RD)
178 return ieee754dp_zero(1); /* round negative inf. => sign = -1 */
179 else
180 return ieee754dp_zero(0); /* other round modes => sign = 1 */
181 }
182
183 /* normalize to rounding precision
184 */
185 while ((xm >> (DP_MBITS + 3)) == 0) {
186 xm <<= 1;
187 xe--;
188 }
189 }
190 DPNORMRET2(xs, xe, xm, "sub", x, y);
191}
diff --git a/arch/mips/math-emu/dp_tint.c b/arch/mips/math-emu/dp_tint.c
new file mode 100644
index 000000000000..77b2b7ccf28a
--- /dev/null
+++ b/arch/mips/math-emu/dp_tint.c
@@ -0,0 +1,124 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include <linux/kernel.h>
29#include "ieee754dp.h"
30
31int ieee754dp_tint(ieee754dp x)
32{
33 COMPXDP;
34
35 CLEARCX;
36
37 EXPLODEXDP;
38 FLUSHXDP;
39
40 switch (xc) {
41 case IEEE754_CLASS_SNAN:
42 case IEEE754_CLASS_QNAN:
43 case IEEE754_CLASS_INF:
44 SETCX(IEEE754_INVALID_OPERATION);
45 return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x);
46 case IEEE754_CLASS_ZERO:
47 return 0;
48 case IEEE754_CLASS_DNORM:
49 case IEEE754_CLASS_NORM:
50 break;
51 }
52 if (xe > 31) {
53 /* Set invalid. We will only use overflow for floating
54 point overflow */
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x);
57 }
58 /* oh gawd */
59 if (xe > DP_MBITS) {
60 xm <<= xe - DP_MBITS;
61 } else if (xe < DP_MBITS) {
62 u64 residue;
63 int round;
64 int sticky;
65 int odd;
66
67 if (xe < -1) {
68 residue = xm;
69 round = 0;
70 sticky = residue != 0;
71 xm = 0;
72 }
73 else {
74 residue = xm << (64 - DP_MBITS + xe);
75 round = (residue >> 63) != 0;
76 sticky = (residue << 1) != 0;
77 xm >>= DP_MBITS - xe;
78 }
79 /* Note: At this point upper 32 bits of xm are guaranteed
80 to be zero */
81 odd = (xm & 0x1) != 0x0;
82 switch (ieee754_csr.rm) {
83 case IEEE754_RN:
84 if (round && (sticky || odd))
85 xm++;
86 break;
87 case IEEE754_RZ:
88 break;
89 case IEEE754_RU: /* toward +Infinity */
90 if ((round || sticky) && !xs)
91 xm++;
92 break;
93 case IEEE754_RD: /* toward -Infinity */
94 if ((round || sticky) && xs)
95 xm++;
96 break;
97 }
98 /* look for valid corner case 0x80000000 */
99 if ((xm >> 31) != 0 && (xs == 0 || xm != 0x80000000)) {
100 /* This can happen after rounding */
101 SETCX(IEEE754_INVALID_OPERATION);
102 return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x);
103 }
104 if (round || sticky)
105 SETCX(IEEE754_INEXACT);
106 }
107 if (xs)
108 return -xm;
109 else
110 return xm;
111}
112
113
114unsigned int ieee754dp_tuns(ieee754dp x)
115{
116 ieee754dp hb = ieee754dp_1e31();
117
118 /* what if x < 0 ?? */
119 if (ieee754dp_lt(x, hb))
120 return (unsigned) ieee754dp_tint(x);
121
122 return (unsigned) ieee754dp_tint(ieee754dp_sub(x, hb)) |
123 ((unsigned) 1 << 31);
124}
diff --git a/arch/mips/math-emu/dp_tlong.c b/arch/mips/math-emu/dp_tlong.c
new file mode 100644
index 000000000000..d71113e07164
--- /dev/null
+++ b/arch/mips/math-emu/dp_tlong.c
@@ -0,0 +1,127 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30s64 ieee754dp_tlong(ieee754dp x)
31{
32 COMPXDP;
33
34 CLEARCX;
35
36 EXPLODEXDP;
37 FLUSHXDP;
38
39 switch (xc) {
40 case IEEE754_CLASS_SNAN:
41 case IEEE754_CLASS_QNAN:
42 case IEEE754_CLASS_INF:
43 SETCX(IEEE754_INVALID_OPERATION);
44 return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x);
45 case IEEE754_CLASS_ZERO:
46 return 0;
47 case IEEE754_CLASS_DNORM:
48 case IEEE754_CLASS_NORM:
49 break;
50 }
51 if (xe >= 63) {
52 /* look for valid corner case */
53 if (xe == 63 && xs && xm == DP_HIDDEN_BIT)
54 return -0x8000000000000000LL;
55 /* Set invalid. We will only use overflow for floating
56 point overflow */
57 SETCX(IEEE754_INVALID_OPERATION);
58 return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x);
59 }
60 /* oh gawd */
61 if (xe > DP_MBITS) {
62 xm <<= xe - DP_MBITS;
63 } else if (xe < DP_MBITS) {
64 u64 residue;
65 int round;
66 int sticky;
67 int odd;
68
69 if (xe < -1) {
70 residue = xm;
71 round = 0;
72 sticky = residue != 0;
73 xm = 0;
74 }
75 else {
76 /* Shifting a u64 64 times does not work,
77 * so we do it in two steps. Be aware that xe
78 * may be -1 */
79 residue = xm << (xe + 1);
80 residue <<= 63 - DP_MBITS;
81 round = (residue >> 63) != 0;
82 sticky = (residue << 1) != 0;
83 xm >>= DP_MBITS - xe;
84 }
85 odd = (xm & 0x1) != 0x0;
86 switch (ieee754_csr.rm) {
87 case IEEE754_RN:
88 if (round && (sticky || odd))
89 xm++;
90 break;
91 case IEEE754_RZ:
92 break;
93 case IEEE754_RU: /* toward +Infinity */
94 if ((round || sticky) && !xs)
95 xm++;
96 break;
97 case IEEE754_RD: /* toward -Infinity */
98 if ((round || sticky) && xs)
99 xm++;
100 break;
101 }
102 if ((xm >> 63) != 0) {
103 /* This can happen after rounding */
104 SETCX(IEEE754_INVALID_OPERATION);
105 return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x);
106 }
107 if (round || sticky)
108 SETCX(IEEE754_INEXACT);
109 }
110 if (xs)
111 return -xm;
112 else
113 return xm;
114}
115
116
117u64 ieee754dp_tulong(ieee754dp x)
118{
119 ieee754dp hb = ieee754dp_1e63();
120
121 /* what if x < 0 ?? */
122 if (ieee754dp_lt(x, hb))
123 return (u64) ieee754dp_tlong(x);
124
125 return (u64) ieee754dp_tlong(ieee754dp_sub(x, hb)) |
126 (1ULL << 63);
127}
diff --git a/arch/mips/math-emu/dsemul.c b/arch/mips/math-emu/dsemul.c
new file mode 100644
index 000000000000..aa989c2246da
--- /dev/null
+++ b/arch/mips/math-emu/dsemul.c
@@ -0,0 +1,172 @@
1#include <linux/compiler.h>
2#include <linux/mm.h>
3#include <linux/signal.h>
4#include <linux/smp.h>
5#include <linux/smp_lock.h>
6
7#include <asm/asm.h>
8#include <asm/bootinfo.h>
9#include <asm/byteorder.h>
10#include <asm/cpu.h>
11#include <asm/inst.h>
12#include <asm/processor.h>
13#include <asm/uaccess.h>
14#include <asm/branch.h>
15#include <asm/mipsregs.h>
16#include <asm/system.h>
17#include <asm/cacheflush.h>
18
19#include <asm/fpu_emulator.h>
20
21#include "ieee754.h"
22#include "dsemul.h"
23
24/* Strap kernel emulator for full MIPS IV emulation */
25
26#ifdef __mips
27#undef __mips
28#endif
29#define __mips 4
30
31extern struct mips_fpu_emulator_private fpuemuprivate;
32
33
34/*
35 * Emulate the arbritrary instruction ir at xcp->cp0_epc. Required when
36 * we have to emulate the instruction in a COP1 branch delay slot. Do
37 * not change cp0_epc due to the instruction
38 *
39 * According to the spec:
40 * 1) it shouldnt be a branch :-)
41 * 2) it can be a COP instruction :-(
42 * 3) if we are tring to run a protected memory space we must take
43 * special care on memory access instructions :-(
44 */
45
46/*
47 * "Trampoline" return routine to catch exception following
48 * execution of delay-slot instruction execution.
49 */
50
51struct emuframe {
52 mips_instruction emul;
53 mips_instruction badinst;
54 mips_instruction cookie;
55 gpreg_t epc;
56};
57
58int mips_dsemul(struct pt_regs *regs, mips_instruction ir, gpreg_t cpc)
59{
60 extern asmlinkage void handle_dsemulret(void);
61 mips_instruction *dsemul_insns;
62 struct emuframe *fr;
63 int err;
64
65 if (ir == 0) { /* a nop is easy */
66 regs->cp0_epc = cpc;
67 regs->cp0_cause &= ~CAUSEF_BD;
68 return 0;
69 }
70#ifdef DSEMUL_TRACE
71 printk("dsemul %lx %lx\n", regs->cp0_epc, cpc);
72
73#endif
74
75 /*
76 * The strategy is to push the instruction onto the user stack
77 * and put a trap after it which we can catch and jump to
78 * the required address any alternative apart from full
79 * instruction emulation!!.
80 *
81 * Algorithmics used a system call instruction, and
82 * borrowed that vector. MIPS/Linux version is a bit
83 * more heavyweight in the interests of portability and
84 * multiprocessor support. For Linux we generate a
85 * an unaligned access and force an address error exception.
86 *
87 * For embedded systems (stand-alone) we prefer to use a
88 * non-existing CP1 instruction. This prevents us from emulating
89 * branches, but gives us a cleaner interface to the exception
90 * handler (single entry point).
91 */
92
93 /* Ensure that the two instructions are in the same cache line */
94 dsemul_insns = (mips_instruction *) REG_TO_VA ((regs->regs[29] - sizeof(struct emuframe)) & ~0x7);
95 fr = (struct emuframe *) dsemul_insns;
96
97 /* Verify that the stack pointer is not competely insane */
98 if (unlikely(!access_ok(VERIFY_WRITE, fr, sizeof(struct emuframe))))
99 return SIGBUS;
100
101 err = __put_user(ir, &fr->emul);
102 err |= __put_user((mips_instruction)BADINST, &fr->badinst);
103 err |= __put_user((mips_instruction)BD_COOKIE, &fr->cookie);
104 err |= __put_user(cpc, &fr->epc);
105
106 if (unlikely(err)) {
107 fpuemuprivate.stats.errors++;
108 return SIGBUS;
109 }
110
111 regs->cp0_epc = VA_TO_REG & fr->emul;
112
113 flush_cache_sigtramp((unsigned long)&fr->badinst);
114
115 return SIGILL; /* force out of emulation loop */
116}
117
118int do_dsemulret(struct pt_regs *xcp)
119{
120 struct emuframe *fr;
121 gpreg_t epc;
122 u32 insn, cookie;
123 int err = 0;
124
125 fr = (struct emuframe *) (xcp->cp0_epc - sizeof(mips_instruction));
126
127 /*
128 * If we can't even access the area, something is very wrong, but we'll
129 * leave that to the default handling
130 */
131 if (!access_ok(VERIFY_READ, fr, sizeof(struct emuframe)))
132 return 0;
133
134 /*
135 * Do some sanity checking on the stackframe:
136 *
137 * - Is the instruction pointed to by the EPC an BADINST?
138 * - Is the following memory word the BD_COOKIE?
139 */
140 err = __get_user(insn, &fr->badinst);
141 err |= __get_user(cookie, &fr->cookie);
142
143 if (unlikely(err || (insn != BADINST) || (cookie != BD_COOKIE))) {
144 fpuemuprivate.stats.errors++;
145 return 0;
146 }
147
148 /*
149 * At this point, we are satisfied that it's a BD emulation trap. Yes,
150 * a user might have deliberately put two malformed and useless
151 * instructions in a row in his program, in which case he's in for a
152 * nasty surprise - the next instruction will be treated as a
153 * continuation address! Alas, this seems to be the only way that we
154 * can handle signals, recursion, and longjmps() in the context of
155 * emulating the branch delay instruction.
156 */
157
158#ifdef DSEMUL_TRACE
159 printk("dsemulret\n");
160#endif
161 if (__get_user(epc, &fr->epc)) { /* Saved EPC */
162 /* This is not a good situation to be in */
163 force_sig(SIGBUS, current);
164
165 return 0;
166 }
167
168 /* Set EPC to return to post-branch instruction */
169 xcp->cp0_epc = epc;
170
171 return 1;
172}
diff --git a/arch/mips/math-emu/dsemul.h b/arch/mips/math-emu/dsemul.h
new file mode 100644
index 000000000000..dbd85f95268d
--- /dev/null
+++ b/arch/mips/math-emu/dsemul.h
@@ -0,0 +1,23 @@
1typedef long gpreg_t;
2typedef void *vaddr_t;
3
4#define REG_TO_VA (vaddr_t)
5#define VA_TO_REG (gpreg_t)
6
7int mips_dsemul(struct pt_regs *regs, mips_instruction ir, gpreg_t cpc);
8int do_dsemulret(struct pt_regs *xcp);
9
10/* Instruction which will always cause an address error */
11#define AdELOAD 0x8c000001 /* lw $0,1($0) */
12/* Instruction which will plainly cause a CP1 exception when FPU is disabled */
13#define CP1UNDEF 0x44400001 /* cfc1 $0,$0 undef */
14
15/* Instruction inserted following the badinst to further tag the sequence */
16#define BD_COOKIE 0x0000bd36 /* tne $0,$0 with baggage */
17
18/* Setup which instruction to use for trampoline */
19#ifdef STANDALONE_EMULATOR
20#define BADINST CP1UNDEF
21#else
22#define BADINST AdELOAD
23#endif
diff --git a/arch/mips/math-emu/ieee754.c b/arch/mips/math-emu/ieee754.c
new file mode 100644
index 000000000000..f0a364adbf34
--- /dev/null
+++ b/arch/mips/math-emu/ieee754.c
@@ -0,0 +1,138 @@
1/* ieee754 floating point arithmetic
2 * single and double precision
3 *
4 * BUGS
5 * not much dp done
6 * doesn't generate IEEE754_INEXACT
7 *
8 */
9/*
10 * MIPS floating point support
11 * Copyright (C) 1994-2000 Algorithmics Ltd.
12 * http://www.algor.co.uk
13 *
14 * ########################################################################
15 *
16 * This program is free software; you can distribute it and/or modify it
17 * under the terms of the GNU General Public License (Version 2) as
18 * published by the Free Software Foundation.
19 *
20 * This program is distributed in the hope it will be useful, but WITHOUT
21 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23 * for more details.
24 *
25 * You should have received a copy of the GNU General Public License along
26 * with this program; if not, write to the Free Software Foundation, Inc.,
27 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
28 *
29 * ########################################################################
30 */
31
32
33#include "ieee754int.h"
34
35#define DP_EBIAS 1023
36#define DP_EMIN (-1022)
37#define DP_EMAX 1023
38
39#define SP_EBIAS 127
40#define SP_EMIN (-126)
41#define SP_EMAX 127
42
43/* indexed by class */
44const char *const ieee754_cname[] = {
45 "Normal",
46 "Zero",
47 "Denormal",
48 "Infinity",
49 "QNaN",
50 "SNaN",
51};
52
53/* the control status register
54*/
55struct ieee754_csr ieee754_csr;
56
57/* special constants
58*/
59
60
61#if (defined(BYTE_ORDER) && BYTE_ORDER == LITTLE_ENDIAN) || defined(__MIPSEL__)
62#define SPSTR(s,b,m) {m,b,s}
63#define DPSTR(s,b,mh,ml) {ml,mh,b,s}
64#endif
65
66#ifdef __MIPSEB__
67#define SPSTR(s,b,m) {s,b,m}
68#define DPSTR(s,b,mh,ml) {s,b,mh,ml}
69#endif
70
71const struct ieee754dp_konst __ieee754dp_spcvals[] = {
72 DPSTR(0, DP_EMIN - 1 + DP_EBIAS, 0, 0), /* + zero */
73 DPSTR(1, DP_EMIN - 1 + DP_EBIAS, 0, 0), /* - zero */
74 DPSTR(0, DP_EBIAS, 0, 0), /* + 1.0 */
75 DPSTR(1, DP_EBIAS, 0, 0), /* - 1.0 */
76 DPSTR(0, 3 + DP_EBIAS, 0x40000, 0), /* + 10.0 */
77 DPSTR(1, 3 + DP_EBIAS, 0x40000, 0), /* - 10.0 */
78 DPSTR(0, DP_EMAX + 1 + DP_EBIAS, 0, 0), /* + infinity */
79 DPSTR(1, DP_EMAX + 1 + DP_EBIAS, 0, 0), /* - infinity */
80 DPSTR(0,DP_EMAX+1+DP_EBIAS,0x7FFFF,0xFFFFFFFF), /* + indef quiet Nan */
81 DPSTR(0, DP_EMAX + DP_EBIAS, 0xFFFFF, 0xFFFFFFFF), /* + max */
82 DPSTR(1, DP_EMAX + DP_EBIAS, 0xFFFFF, 0xFFFFFFFF), /* - max */
83 DPSTR(0, DP_EMIN + DP_EBIAS, 0, 0), /* + min normal */
84 DPSTR(1, DP_EMIN + DP_EBIAS, 0, 0), /* - min normal */
85 DPSTR(0, DP_EMIN - 1 + DP_EBIAS, 0, 1), /* + min denormal */
86 DPSTR(1, DP_EMIN - 1 + DP_EBIAS, 0, 1), /* - min denormal */
87 DPSTR(0, 31 + DP_EBIAS, 0, 0), /* + 1.0e31 */
88 DPSTR(0, 63 + DP_EBIAS, 0, 0), /* + 1.0e63 */
89};
90
91const struct ieee754sp_konst __ieee754sp_spcvals[] = {
92 SPSTR(0, SP_EMIN - 1 + SP_EBIAS, 0), /* + zero */
93 SPSTR(1, SP_EMIN - 1 + SP_EBIAS, 0), /* - zero */
94 SPSTR(0, SP_EBIAS, 0), /* + 1.0 */
95 SPSTR(1, SP_EBIAS, 0), /* - 1.0 */
96 SPSTR(0, 3 + SP_EBIAS, 0x200000), /* + 10.0 */
97 SPSTR(1, 3 + SP_EBIAS, 0x200000), /* - 10.0 */
98 SPSTR(0, SP_EMAX + 1 + SP_EBIAS, 0), /* + infinity */
99 SPSTR(1, SP_EMAX + 1 + SP_EBIAS, 0), /* - infinity */
100 SPSTR(0,SP_EMAX+1+SP_EBIAS,0x3FFFFF), /* + indef quiet Nan */
101 SPSTR(0, SP_EMAX + SP_EBIAS, 0x7FFFFF), /* + max normal */
102 SPSTR(1, SP_EMAX + SP_EBIAS, 0x7FFFFF), /* - max normal */
103 SPSTR(0, SP_EMIN + SP_EBIAS, 0), /* + min normal */
104 SPSTR(1, SP_EMIN + SP_EBIAS, 0), /* - min normal */
105 SPSTR(0, SP_EMIN - 1 + SP_EBIAS, 1), /* + min denormal */
106 SPSTR(1, SP_EMIN - 1 + SP_EBIAS, 1), /* - min denormal */
107 SPSTR(0, 31 + SP_EBIAS, 0), /* + 1.0e31 */
108 SPSTR(0, 63 + SP_EBIAS, 0), /* + 1.0e63 */
109};
110
111
112int ieee754si_xcpt(int r, const char *op, ...)
113{
114 struct ieee754xctx ax;
115
116 if (!TSTX())
117 return r;
118 ax.op = op;
119 ax.rt = IEEE754_RT_SI;
120 ax.rv.si = r;
121 va_start(ax.ap, op);
122 ieee754_xcpt(&ax);
123 return ax.rv.si;
124}
125
126s64 ieee754di_xcpt(s64 r, const char *op, ...)
127{
128 struct ieee754xctx ax;
129
130 if (!TSTX())
131 return r;
132 ax.op = op;
133 ax.rt = IEEE754_RT_DI;
134 ax.rv.di = r;
135 va_start(ax.ap, op);
136 ieee754_xcpt(&ax);
137 return ax.rv.di;
138}
diff --git a/arch/mips/math-emu/ieee754.h b/arch/mips/math-emu/ieee754.h
new file mode 100644
index 000000000000..b8772f46972d
--- /dev/null
+++ b/arch/mips/math-emu/ieee754.h
@@ -0,0 +1,489 @@
1/* single and double precision fp ops
2 * missing extended precision.
3*/
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27/**************************************************************************
28 * Nov 7, 2000
29 * Modification to allow integration with Linux kernel
30 *
31 * Kevin D. Kissell, kevink@mips.com and Carsten Langgard, carstenl@mips.com
32 * Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
33 *************************************************************************/
34
35#ifdef __KERNEL__
36/* Going from Algorithmics to Linux native environment, add this */
37#include <linux/types.h>
38
39/*
40 * Not very pretty, but the Linux kernel's normal va_list definition
41 * does not allow it to be used as a structure element, as it is here.
42 */
43#ifndef _STDARG_H
44#include <stdarg.h>
45#endif
46
47#else
48
49/* Note that __KERNEL__ is taken to mean Linux kernel */
50
51#if #system(OpenBSD)
52#include <machine/types.h>
53#endif
54#include <machine/endian.h>
55
56#endif /* __KERNEL__ */
57
58#if (defined(BYTE_ORDER) && BYTE_ORDER == LITTLE_ENDIAN) || defined(__MIPSEL__)
59struct ieee754dp_konst {
60 unsigned mantlo:32;
61 unsigned manthi:20;
62 unsigned bexp:11;
63 unsigned sign:1;
64};
65struct ieee754sp_konst {
66 unsigned mant:23;
67 unsigned bexp:8;
68 unsigned sign:1;
69};
70
71typedef union _ieee754dp {
72 struct ieee754dp_konst oparts;
73 struct {
74 u64 mant:52;
75 unsigned int bexp:11;
76 unsigned int sign:1;
77 } parts;
78 u64 bits;
79 double d;
80} ieee754dp;
81
82typedef union _ieee754sp {
83 struct ieee754sp_konst parts;
84 float f;
85 u32 bits;
86} ieee754sp;
87#endif
88
89#if (defined(BYTE_ORDER) && BYTE_ORDER == BIG_ENDIAN) || defined(__MIPSEB__)
90struct ieee754dp_konst {
91 unsigned sign:1;
92 unsigned bexp:11;
93 unsigned manthi:20;
94 unsigned mantlo:32;
95};
96typedef union _ieee754dp {
97 struct ieee754dp_konst oparts;
98 struct {
99 unsigned int sign:1;
100 unsigned int bexp:11;
101 u64 mant:52;
102 } parts;
103 double d;
104 u64 bits;
105} ieee754dp;
106
107struct ieee754sp_konst {
108 unsigned sign:1;
109 unsigned bexp:8;
110 unsigned mant:23;
111};
112
113typedef union _ieee754sp {
114 struct ieee754sp_konst parts;
115 float f;
116 u32 bits;
117} ieee754sp;
118#endif
119
120/*
121 * single precision (often aka float)
122*/
123int ieee754sp_finite(ieee754sp x);
124int ieee754sp_class(ieee754sp x);
125
126ieee754sp ieee754sp_abs(ieee754sp x);
127ieee754sp ieee754sp_neg(ieee754sp x);
128ieee754sp ieee754sp_scalb(ieee754sp x, int);
129ieee754sp ieee754sp_logb(ieee754sp x);
130
131/* x with sign of y */
132ieee754sp ieee754sp_copysign(ieee754sp x, ieee754sp y);
133
134ieee754sp ieee754sp_add(ieee754sp x, ieee754sp y);
135ieee754sp ieee754sp_sub(ieee754sp x, ieee754sp y);
136ieee754sp ieee754sp_mul(ieee754sp x, ieee754sp y);
137ieee754sp ieee754sp_div(ieee754sp x, ieee754sp y);
138
139ieee754sp ieee754sp_fint(int x);
140ieee754sp ieee754sp_funs(unsigned x);
141ieee754sp ieee754sp_flong(s64 x);
142ieee754sp ieee754sp_fulong(u64 x);
143ieee754sp ieee754sp_fdp(ieee754dp x);
144
145int ieee754sp_tint(ieee754sp x);
146unsigned int ieee754sp_tuns(ieee754sp x);
147s64 ieee754sp_tlong(ieee754sp x);
148u64 ieee754sp_tulong(ieee754sp x);
149
150int ieee754sp_cmp(ieee754sp x, ieee754sp y, int cop, int sig);
151/*
152 * basic sp math
153 */
154ieee754sp ieee754sp_modf(ieee754sp x, ieee754sp * ip);
155ieee754sp ieee754sp_frexp(ieee754sp x, int *exp);
156ieee754sp ieee754sp_ldexp(ieee754sp x, int exp);
157
158ieee754sp ieee754sp_ceil(ieee754sp x);
159ieee754sp ieee754sp_floor(ieee754sp x);
160ieee754sp ieee754sp_trunc(ieee754sp x);
161
162ieee754sp ieee754sp_sqrt(ieee754sp x);
163
164/*
165 * double precision (often aka double)
166*/
167int ieee754dp_finite(ieee754dp x);
168int ieee754dp_class(ieee754dp x);
169
170/* x with sign of y */
171ieee754dp ieee754dp_copysign(ieee754dp x, ieee754dp y);
172
173ieee754dp ieee754dp_add(ieee754dp x, ieee754dp y);
174ieee754dp ieee754dp_sub(ieee754dp x, ieee754dp y);
175ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y);
176ieee754dp ieee754dp_div(ieee754dp x, ieee754dp y);
177
178ieee754dp ieee754dp_abs(ieee754dp x);
179ieee754dp ieee754dp_neg(ieee754dp x);
180ieee754dp ieee754dp_scalb(ieee754dp x, int);
181
182/* return exponent as integer in floating point format
183 */
184ieee754dp ieee754dp_logb(ieee754dp x);
185
186ieee754dp ieee754dp_fint(int x);
187ieee754dp ieee754dp_funs(unsigned x);
188ieee754dp ieee754dp_flong(s64 x);
189ieee754dp ieee754dp_fulong(u64 x);
190ieee754dp ieee754dp_fsp(ieee754sp x);
191
192ieee754dp ieee754dp_ceil(ieee754dp x);
193ieee754dp ieee754dp_floor(ieee754dp x);
194ieee754dp ieee754dp_trunc(ieee754dp x);
195
196int ieee754dp_tint(ieee754dp x);
197unsigned int ieee754dp_tuns(ieee754dp x);
198s64 ieee754dp_tlong(ieee754dp x);
199u64 ieee754dp_tulong(ieee754dp x);
200
201int ieee754dp_cmp(ieee754dp x, ieee754dp y, int cop, int sig);
202/*
203 * basic sp math
204 */
205ieee754dp ieee754dp_modf(ieee754dp x, ieee754dp * ip);
206ieee754dp ieee754dp_frexp(ieee754dp x, int *exp);
207ieee754dp ieee754dp_ldexp(ieee754dp x, int exp);
208
209ieee754dp ieee754dp_ceil(ieee754dp x);
210ieee754dp ieee754dp_floor(ieee754dp x);
211ieee754dp ieee754dp_trunc(ieee754dp x);
212
213ieee754dp ieee754dp_sqrt(ieee754dp x);
214
215
216
217/* 5 types of floating point number
218*/
219#define IEEE754_CLASS_NORM 0x00
220#define IEEE754_CLASS_ZERO 0x01
221#define IEEE754_CLASS_DNORM 0x02
222#define IEEE754_CLASS_INF 0x03
223#define IEEE754_CLASS_SNAN 0x04
224#define IEEE754_CLASS_QNAN 0x05
225extern const char *const ieee754_cname[];
226
227/* exception numbers */
228#define IEEE754_INEXACT 0x01
229#define IEEE754_UNDERFLOW 0x02
230#define IEEE754_OVERFLOW 0x04
231#define IEEE754_ZERO_DIVIDE 0x08
232#define IEEE754_INVALID_OPERATION 0x10
233
234/* cmp operators
235*/
236#define IEEE754_CLT 0x01
237#define IEEE754_CEQ 0x02
238#define IEEE754_CGT 0x04
239#define IEEE754_CUN 0x08
240
241/* rounding mode
242*/
243#define IEEE754_RN 0 /* round to nearest */
244#define IEEE754_RZ 1 /* round toward zero */
245#define IEEE754_RD 2 /* round toward -Infinity */
246#define IEEE754_RU 3 /* round toward +Infinity */
247
248/* other naming */
249#define IEEE754_RM IEEE754_RD
250#define IEEE754_RP IEEE754_RU
251
252/* "normal" comparisons
253*/
254static __inline int ieee754sp_eq(ieee754sp x, ieee754sp y)
255{
256 return ieee754sp_cmp(x, y, IEEE754_CEQ, 0);
257}
258
259static __inline int ieee754sp_ne(ieee754sp x, ieee754sp y)
260{
261 return ieee754sp_cmp(x, y,
262 IEEE754_CLT | IEEE754_CGT | IEEE754_CUN, 0);
263}
264
265static __inline int ieee754sp_lt(ieee754sp x, ieee754sp y)
266{
267 return ieee754sp_cmp(x, y, IEEE754_CLT, 0);
268}
269
270static __inline int ieee754sp_le(ieee754sp x, ieee754sp y)
271{
272 return ieee754sp_cmp(x, y, IEEE754_CLT | IEEE754_CEQ, 0);
273}
274
275static __inline int ieee754sp_gt(ieee754sp x, ieee754sp y)
276{
277 return ieee754sp_cmp(x, y, IEEE754_CGT, 0);
278}
279
280
281static __inline int ieee754sp_ge(ieee754sp x, ieee754sp y)
282{
283 return ieee754sp_cmp(x, y, IEEE754_CGT | IEEE754_CEQ, 0);
284}
285
286static __inline int ieee754dp_eq(ieee754dp x, ieee754dp y)
287{
288 return ieee754dp_cmp(x, y, IEEE754_CEQ, 0);
289}
290
291static __inline int ieee754dp_ne(ieee754dp x, ieee754dp y)
292{
293 return ieee754dp_cmp(x, y,
294 IEEE754_CLT | IEEE754_CGT | IEEE754_CUN, 0);
295}
296
297static __inline int ieee754dp_lt(ieee754dp x, ieee754dp y)
298{
299 return ieee754dp_cmp(x, y, IEEE754_CLT, 0);
300}
301
302static __inline int ieee754dp_le(ieee754dp x, ieee754dp y)
303{
304 return ieee754dp_cmp(x, y, IEEE754_CLT | IEEE754_CEQ, 0);
305}
306
307static __inline int ieee754dp_gt(ieee754dp x, ieee754dp y)
308{
309 return ieee754dp_cmp(x, y, IEEE754_CGT, 0);
310}
311
312static __inline int ieee754dp_ge(ieee754dp x, ieee754dp y)
313{
314 return ieee754dp_cmp(x, y, IEEE754_CGT | IEEE754_CEQ, 0);
315}
316
317
318/* like strtod
319*/
320ieee754dp ieee754dp_fstr(const char *s, char **endp);
321char *ieee754dp_tstr(ieee754dp x, int prec, int fmt, int af);
322
323
324/* the control status register
325*/
326struct ieee754_csr {
327 unsigned pad:13;
328 unsigned nod:1; /* set 1 for no denormalised numbers */
329 unsigned cx:5; /* exceptions this operation */
330 unsigned mx:5; /* exception enable mask */
331 unsigned sx:5; /* exceptions total */
332 unsigned rm:2; /* current rounding mode */
333};
334extern struct ieee754_csr ieee754_csr;
335
336static __inline unsigned ieee754_getrm(void)
337{
338 return (ieee754_csr.rm);
339}
340static __inline unsigned ieee754_setrm(unsigned rm)
341{
342 return (ieee754_csr.rm = rm);
343}
344
345/*
346 * get current exceptions
347 */
348static __inline unsigned ieee754_getcx(void)
349{
350 return (ieee754_csr.cx);
351}
352
353/* test for current exception condition
354 */
355static __inline int ieee754_cxtest(unsigned n)
356{
357 return (ieee754_csr.cx & n);
358}
359
360/*
361 * get sticky exceptions
362 */
363static __inline unsigned ieee754_getsx(void)
364{
365 return (ieee754_csr.sx);
366}
367
368/* clear sticky conditions
369*/
370static __inline unsigned ieee754_clrsx(void)
371{
372 return (ieee754_csr.sx = 0);
373}
374
375/* test for sticky exception condition
376 */
377static __inline int ieee754_sxtest(unsigned n)
378{
379 return (ieee754_csr.sx & n);
380}
381
382/* debugging */
383ieee754sp ieee754sp_dump(char *s, ieee754sp x);
384ieee754dp ieee754dp_dump(char *s, ieee754dp x);
385
386#define IEEE754_SPCVAL_PZERO 0
387#define IEEE754_SPCVAL_NZERO 1
388#define IEEE754_SPCVAL_PONE 2
389#define IEEE754_SPCVAL_NONE 3
390#define IEEE754_SPCVAL_PTEN 4
391#define IEEE754_SPCVAL_NTEN 5
392#define IEEE754_SPCVAL_PINFINITY 6
393#define IEEE754_SPCVAL_NINFINITY 7
394#define IEEE754_SPCVAL_INDEF 8
395#define IEEE754_SPCVAL_PMAX 9 /* +max norm */
396#define IEEE754_SPCVAL_NMAX 10 /* -max norm */
397#define IEEE754_SPCVAL_PMIN 11 /* +min norm */
398#define IEEE754_SPCVAL_NMIN 12 /* +min norm */
399#define IEEE754_SPCVAL_PMIND 13 /* +min denorm */
400#define IEEE754_SPCVAL_NMIND 14 /* +min denorm */
401#define IEEE754_SPCVAL_P1E31 15 /* + 1.0e31 */
402#define IEEE754_SPCVAL_P1E63 16 /* + 1.0e63 */
403
404extern const struct ieee754dp_konst __ieee754dp_spcvals[];
405extern const struct ieee754sp_konst __ieee754sp_spcvals[];
406#define ieee754dp_spcvals ((const ieee754dp *)__ieee754dp_spcvals)
407#define ieee754sp_spcvals ((const ieee754sp *)__ieee754sp_spcvals)
408
409/* return infinity with given sign
410*/
411#define ieee754dp_inf(sn) \
412 (ieee754dp_spcvals[IEEE754_SPCVAL_PINFINITY+(sn)])
413#define ieee754dp_zero(sn) \
414 (ieee754dp_spcvals[IEEE754_SPCVAL_PZERO+(sn)])
415#define ieee754dp_one(sn) \
416 (ieee754dp_spcvals[IEEE754_SPCVAL_PONE+(sn)])
417#define ieee754dp_ten(sn) \
418 (ieee754dp_spcvals[IEEE754_SPCVAL_PTEN+(sn)])
419#define ieee754dp_indef() \
420 (ieee754dp_spcvals[IEEE754_SPCVAL_INDEF])
421#define ieee754dp_max(sn) \
422 (ieee754dp_spcvals[IEEE754_SPCVAL_PMAX+(sn)])
423#define ieee754dp_min(sn) \
424 (ieee754dp_spcvals[IEEE754_SPCVAL_PMIN+(sn)])
425#define ieee754dp_mind(sn) \
426 (ieee754dp_spcvals[IEEE754_SPCVAL_PMIND+(sn)])
427#define ieee754dp_1e31() \
428 (ieee754dp_spcvals[IEEE754_SPCVAL_P1E31])
429#define ieee754dp_1e63() \
430 (ieee754dp_spcvals[IEEE754_SPCVAL_P1E63])
431
432#define ieee754sp_inf(sn) \
433 (ieee754sp_spcvals[IEEE754_SPCVAL_PINFINITY+(sn)])
434#define ieee754sp_zero(sn) \
435 (ieee754sp_spcvals[IEEE754_SPCVAL_PZERO+(sn)])
436#define ieee754sp_one(sn) \
437 (ieee754sp_spcvals[IEEE754_SPCVAL_PONE+(sn)])
438#define ieee754sp_ten(sn) \
439 (ieee754sp_spcvals[IEEE754_SPCVAL_PTEN+(sn)])
440#define ieee754sp_indef() \
441 (ieee754sp_spcvals[IEEE754_SPCVAL_INDEF])
442#define ieee754sp_max(sn) \
443 (ieee754sp_spcvals[IEEE754_SPCVAL_PMAX+(sn)])
444#define ieee754sp_min(sn) \
445 (ieee754sp_spcvals[IEEE754_SPCVAL_PMIN+(sn)])
446#define ieee754sp_mind(sn) \
447 (ieee754sp_spcvals[IEEE754_SPCVAL_PMIND+(sn)])
448#define ieee754sp_1e31() \
449 (ieee754sp_spcvals[IEEE754_SPCVAL_P1E31])
450#define ieee754sp_1e63() \
451 (ieee754sp_spcvals[IEEE754_SPCVAL_P1E63])
452
453/* indefinite integer value
454*/
455#define ieee754si_indef() INT_MAX
456#ifdef LONG_LONG_MAX
457#define ieee754di_indef() LONG_LONG_MAX
458#else
459#define ieee754di_indef() ((s64)(~0ULL>>1))
460#endif
461
462/* IEEE exception context, passed to handler */
463struct ieee754xctx {
464 const char *op; /* operation name */
465 int rt; /* result type */
466 union {
467 ieee754sp sp; /* single precision */
468 ieee754dp dp; /* double precision */
469#ifdef IEEE854_XP
470 ieee754xp xp; /* extended precision */
471#endif
472 int si; /* standard signed integer (32bits) */
473 s64 di; /* extended signed integer (64bits) */
474 } rv; /* default result format implied by op */
475 va_list ap;
476};
477
478/* result types for xctx.rt */
479#define IEEE754_RT_SP 0
480#define IEEE754_RT_DP 1
481#define IEEE754_RT_XP 2
482#define IEEE754_RT_SI 3
483#define IEEE754_RT_DI 4
484
485extern void ieee754_xcpt(struct ieee754xctx *xcp);
486
487/* compat */
488#define ieee754dp_fix(x) ieee754dp_tint(x)
489#define ieee754sp_fix(x) ieee754sp_tint(x)
diff --git a/arch/mips/math-emu/ieee754d.c b/arch/mips/math-emu/ieee754d.c
new file mode 100644
index 000000000000..7e900f30987e
--- /dev/null
+++ b/arch/mips/math-emu/ieee754d.c
@@ -0,0 +1,138 @@
1/*
2 * Some debug functions
3 *
4 * MIPS floating point support
5 *
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * This program is free software; you can distribute it and/or modify it
10 * under the terms of the GNU General Public License (Version 2) as
11 * published by the Free Software Foundation.
12 *
13 * This program is distributed in the hope it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 * for more details.
17 *
18 * You should have received a copy of the GNU General Public License along
19 * with this program; if not, write to the Free Software Foundation, Inc.,
20 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
21 *
22 * Nov 7, 2000
23 * Modified to build and operate in Linux kernel environment.
24 *
25 * Kevin D. Kissell, kevink@mips.com and Carsten Langgaard, carstenl@mips.com
26 * Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
27 */
28
29#include <linux/kernel.h>
30#include "ieee754.h"
31
32#define DP_EBIAS 1023
33#define DP_EMIN (-1022)
34#define DP_EMAX 1023
35#define DP_FBITS 52
36
37#define SP_EBIAS 127
38#define SP_EMIN (-126)
39#define SP_EMAX 127
40#define SP_FBITS 23
41
42#define DP_MBIT(x) ((u64)1 << (x))
43#define DP_HIDDEN_BIT DP_MBIT(DP_FBITS)
44#define DP_SIGN_BIT DP_MBIT(63)
45
46
47#define SP_MBIT(x) ((u32)1 << (x))
48#define SP_HIDDEN_BIT SP_MBIT(SP_FBITS)
49#define SP_SIGN_BIT SP_MBIT(31)
50
51
52#define SPSIGN(sp) (sp.parts.sign)
53#define SPBEXP(sp) (sp.parts.bexp)
54#define SPMANT(sp) (sp.parts.mant)
55
56#define DPSIGN(dp) (dp.parts.sign)
57#define DPBEXP(dp) (dp.parts.bexp)
58#define DPMANT(dp) (dp.parts.mant)
59
60ieee754dp ieee754dp_dump(char *m, ieee754dp x)
61{
62 int i;
63
64 printk("%s", m);
65 printk("<%08x,%08x>\n", (unsigned) (x.bits >> 32),
66 (unsigned) x.bits);
67 printk("\t=");
68 switch (ieee754dp_class(x)) {
69 case IEEE754_CLASS_QNAN:
70 case IEEE754_CLASS_SNAN:
71 printk("Nan %c", DPSIGN(x) ? '-' : '+');
72 for (i = DP_FBITS - 1; i >= 0; i--)
73 printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
74 break;
75 case IEEE754_CLASS_INF:
76 printk("%cInfinity", DPSIGN(x) ? '-' : '+');
77 break;
78 case IEEE754_CLASS_ZERO:
79 printk("%cZero", DPSIGN(x) ? '-' : '+');
80 break;
81 case IEEE754_CLASS_DNORM:
82 printk("%c0.", DPSIGN(x) ? '-' : '+');
83 for (i = DP_FBITS - 1; i >= 0; i--)
84 printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
85 printk("e%d", DPBEXP(x) - DP_EBIAS);
86 break;
87 case IEEE754_CLASS_NORM:
88 printk("%c1.", DPSIGN(x) ? '-' : '+');
89 for (i = DP_FBITS - 1; i >= 0; i--)
90 printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
91 printk("e%d", DPBEXP(x) - DP_EBIAS);
92 break;
93 default:
94 printk("Illegal/Unknown IEEE754 value class");
95 }
96 printk("\n");
97 return x;
98}
99
100ieee754sp ieee754sp_dump(char *m, ieee754sp x)
101{
102 int i;
103
104 printk("%s=", m);
105 printk("<%08x>\n", (unsigned) x.bits);
106 printk("\t=");
107 switch (ieee754sp_class(x)) {
108 case IEEE754_CLASS_QNAN:
109 case IEEE754_CLASS_SNAN:
110 printk("Nan %c", SPSIGN(x) ? '-' : '+');
111 for (i = SP_FBITS - 1; i >= 0; i--)
112 printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
113 break;
114 case IEEE754_CLASS_INF:
115 printk("%cInfinity", SPSIGN(x) ? '-' : '+');
116 break;
117 case IEEE754_CLASS_ZERO:
118 printk("%cZero", SPSIGN(x) ? '-' : '+');
119 break;
120 case IEEE754_CLASS_DNORM:
121 printk("%c0.", SPSIGN(x) ? '-' : '+');
122 for (i = SP_FBITS - 1; i >= 0; i--)
123 printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
124 printk("e%d", SPBEXP(x) - SP_EBIAS);
125 break;
126 case IEEE754_CLASS_NORM:
127 printk("%c1.", SPSIGN(x) ? '-' : '+');
128 for (i = SP_FBITS - 1; i >= 0; i--)
129 printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
130 printk("e%d", SPBEXP(x) - SP_EBIAS);
131 break;
132 default:
133 printk("Illegal/Unknown IEEE754 value class");
134 }
135 printk("\n");
136 return x;
137}
138
diff --git a/arch/mips/math-emu/ieee754dp.c b/arch/mips/math-emu/ieee754dp.c
new file mode 100644
index 000000000000..3e214aac4b12
--- /dev/null
+++ b/arch/mips/math-emu/ieee754dp.c
@@ -0,0 +1,243 @@
1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30int ieee754dp_class(ieee754dp x)
31{
32 COMPXDP;
33 EXPLODEXDP;
34 return xc;
35}
36
37int ieee754dp_isnan(ieee754dp x)
38{
39 return ieee754dp_class(x) >= IEEE754_CLASS_SNAN;
40}
41
42int ieee754dp_issnan(ieee754dp x)
43{
44 assert(ieee754dp_isnan(x));
45 return ((DPMANT(x) & DP_MBIT(DP_MBITS-1)) == DP_MBIT(DP_MBITS-1));
46}
47
48
49ieee754dp ieee754dp_xcpt(ieee754dp r, const char *op, ...)
50{
51 struct ieee754xctx ax;
52 if (!TSTX())
53 return r;
54
55 ax.op = op;
56 ax.rt = IEEE754_RT_DP;
57 ax.rv.dp = r;
58 va_start(ax.ap, op);
59 ieee754_xcpt(&ax);
60 return ax.rv.dp;
61}
62
63ieee754dp ieee754dp_nanxcpt(ieee754dp r, const char *op, ...)
64{
65 struct ieee754xctx ax;
66
67 assert(ieee754dp_isnan(r));
68
69 if (!ieee754dp_issnan(r)) /* QNAN does not cause invalid op !! */
70 return r;
71
72 if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) {
73 /* not enabled convert to a quiet NaN */
74 DPMANT(r) &= (~DP_MBIT(DP_MBITS-1));
75 if (ieee754dp_isnan(r))
76 return r;
77 else
78 return ieee754dp_indef();
79 }
80
81 ax.op = op;
82 ax.rt = 0;
83 ax.rv.dp = r;
84 va_start(ax.ap, op);
85 ieee754_xcpt(&ax);
86 return ax.rv.dp;
87}
88
89ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y)
90{
91 assert(ieee754dp_isnan(x));
92 assert(ieee754dp_isnan(y));
93
94 if (DPMANT(x) > DPMANT(y))
95 return x;
96 else
97 return y;
98}
99
100
101static u64 get_rounding(int sn, u64 xm)
102{
103 /* inexact must round of 3 bits
104 */
105 if (xm & (DP_MBIT(3) - 1)) {
106 switch (ieee754_csr.rm) {
107 case IEEE754_RZ:
108 break;
109 case IEEE754_RN:
110 xm += 0x3 + ((xm >> 3) & 1);
111 /* xm += (xm&0x8)?0x4:0x3 */
112 break;
113 case IEEE754_RU: /* toward +Infinity */
114 if (!sn) /* ?? */
115 xm += 0x8;
116 break;
117 case IEEE754_RD: /* toward -Infinity */
118 if (sn) /* ?? */
119 xm += 0x8;
120 break;
121 }
122 }
123 return xm;
124}
125
126
127/* generate a normal/denormal number with over,under handling
128 * sn is sign
129 * xe is an unbiased exponent
130 * xm is 3bit extended precision value.
131 */
132ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
133{
134 assert(xm); /* we don't gen exact zeros (probably should) */
135
136 assert((xm >> (DP_MBITS + 1 + 3)) == 0); /* no execess */
137 assert(xm & (DP_HIDDEN_BIT << 3));
138
139 if (xe < DP_EMIN) {
140 /* strip lower bits */
141 int es = DP_EMIN - xe;
142
143 if (ieee754_csr.nod) {
144 SETCX(IEEE754_UNDERFLOW);
145 SETCX(IEEE754_INEXACT);
146
147 switch(ieee754_csr.rm) {
148 case IEEE754_RN:
149 return ieee754dp_zero(sn);
150 case IEEE754_RZ:
151 return ieee754dp_zero(sn);
152 case IEEE754_RU: /* toward +Infinity */
153 if(sn == 0)
154 return ieee754dp_min(0);
155 else
156 return ieee754dp_zero(1);
157 case IEEE754_RD: /* toward -Infinity */
158 if(sn == 0)
159 return ieee754dp_zero(0);
160 else
161 return ieee754dp_min(1);
162 }
163 }
164
165 if (xe == DP_EMIN - 1
166 && get_rounding(sn, xm) >> (DP_MBITS + 1 + 3))
167 {
168 /* Not tiny after rounding */
169 SETCX(IEEE754_INEXACT);
170 xm = get_rounding(sn, xm);
171 xm >>= 1;
172 /* Clear grs bits */
173 xm &= ~(DP_MBIT(3) - 1);
174 xe++;
175 }
176 else {
177 /* sticky right shift es bits
178 */
179 xm = XDPSRS(xm, es);
180 xe += es;
181 assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
182 assert(xe == DP_EMIN);
183 }
184 }
185 if (xm & (DP_MBIT(3) - 1)) {
186 SETCX(IEEE754_INEXACT);
187 if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
188 SETCX(IEEE754_UNDERFLOW);
189 }
190
191 /* inexact must round of 3 bits
192 */
193 xm = get_rounding(sn, xm);
194 /* adjust exponent for rounding add overflowing
195 */
196 if (xm >> (DP_MBITS + 3 + 1)) {
197 /* add causes mantissa overflow */
198 xm >>= 1;
199 xe++;
200 }
201 }
202 /* strip grs bits */
203 xm >>= 3;
204
205 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */
206 assert(xe >= DP_EMIN);
207
208 if (xe > DP_EMAX) {
209 SETCX(IEEE754_OVERFLOW);
210 SETCX(IEEE754_INEXACT);
211 /* -O can be table indexed by (rm,sn) */
212 switch (ieee754_csr.rm) {
213 case IEEE754_RN:
214 return ieee754dp_inf(sn);
215 case IEEE754_RZ:
216 return ieee754dp_max(sn);
217 case IEEE754_RU: /* toward +Infinity */
218 if (sn == 0)
219 return ieee754dp_inf(0);
220 else
221 return ieee754dp_max(1);
222 case IEEE754_RD: /* toward -Infinity */
223 if (sn == 0)
224 return ieee754dp_max(0);
225 else
226 return ieee754dp_inf(1);
227 }
228 }
229 /* gen norm/denorm/zero */
230
231 if ((xm & DP_HIDDEN_BIT) == 0) {
232 /* we underflow (tiny/zero) */
233 assert(xe == DP_EMIN);
234 if (ieee754_csr.mx & IEEE754_UNDERFLOW)
235 SETCX(IEEE754_UNDERFLOW);
236 return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
237 } else {
238 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */
239 assert(xm & DP_HIDDEN_BIT);
240
241 return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
242 }
243}
diff --git a/arch/mips/math-emu/ieee754dp.h b/arch/mips/math-emu/ieee754dp.h
new file mode 100644
index 000000000000..a37370dae232
--- /dev/null
+++ b/arch/mips/math-emu/ieee754dp.h
@@ -0,0 +1,83 @@
1/*
2 * IEEE754 floating point
3 * double precision internal header file
4 */
5/*
6 * MIPS floating point support
7 * Copyright (C) 1994-2000 Algorithmics Ltd.
8 * http://www.algor.co.uk
9 *
10 * ########################################################################
11 *
12 * This program is free software; you can distribute it and/or modify it
13 * under the terms of the GNU General Public License (Version 2) as
14 * published by the Free Software Foundation.
15 *
16 * This program is distributed in the hope it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19 * for more details.
20 *
21 * You should have received a copy of the GNU General Public License along
22 * with this program; if not, write to the Free Software Foundation, Inc.,
23 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
24 *
25 * ########################################################################
26 */
27
28
29#include "ieee754int.h"
30
31#define assert(expr) ((void)0)
32
33/* 3bit extended double precision sticky right shift */
34#define XDPSRS(v,rs) \
35 ((rs > (DP_MBITS+3))?1:((v) >> (rs)) | ((v) << (64-(rs)) != 0))
36
37#define XDPSRSX1() \
38 (xe++, (xm = (xm >> 1) | (xm & 1)))
39
40#define XDPSRS1(v) \
41 (((v) >> 1) | ((v) & 1))
42
43/* convert denormal to normalized with extended exponent */
44#define DPDNORMx(m,e) \
45 while( (m >> DP_MBITS) == 0) { m <<= 1; e--; }
46#define DPDNORMX DPDNORMx(xm,xe)
47#define DPDNORMY DPDNORMx(ym,ye)
48
49static __inline ieee754dp builddp(int s, int bx, u64 m)
50{
51 ieee754dp r;
52
53 assert((s) == 0 || (s) == 1);
54 assert((bx) >= DP_EMIN - 1 + DP_EBIAS
55 && (bx) <= DP_EMAX + 1 + DP_EBIAS);
56 assert(((m) >> DP_MBITS) == 0);
57
58 r.parts.sign = s;
59 r.parts.bexp = bx;
60 r.parts.mant = m;
61 return r;
62}
63
64extern int ieee754dp_isnan(ieee754dp);
65extern int ieee754dp_issnan(ieee754dp);
66extern int ieee754si_xcpt(int, const char *, ...);
67extern s64 ieee754di_xcpt(s64, const char *, ...);
68extern ieee754dp ieee754dp_xcpt(ieee754dp, const char *, ...);
69extern ieee754dp ieee754dp_nanxcpt(ieee754dp, const char *, ...);
70extern ieee754dp ieee754dp_bestnan(ieee754dp, ieee754dp);
71extern ieee754dp ieee754dp_format(int, int, u64);
72
73
74#define DPNORMRET2(s,e,m,name,a0,a1) \
75{ \
76 ieee754dp V = ieee754dp_format(s,e,m); \
77 if(TSTX()) \
78 return ieee754dp_xcpt(V,name,a0,a1); \
79 else \
80 return V; \
81}
82
83#define DPNORMRET1(s,e,m,name,a0) DPNORMRET2(s,e,m,name,a0,a0)
diff --git a/arch/mips/math-emu/ieee754int.h b/arch/mips/math-emu/ieee754int.h
new file mode 100644
index 000000000000..4a5a81d6b893
--- /dev/null
+++ b/arch/mips/math-emu/ieee754int.h
@@ -0,0 +1,165 @@
1/*
2 * IEEE754 floating point
3 * common internal header file
4 */
5/*
6 * MIPS floating point support
7 * Copyright (C) 1994-2000 Algorithmics Ltd.
8 * http://www.algor.co.uk
9 *
10 * ########################################################################
11 *
12 * This program is free software; you can distribute it and/or modify it
13 * under the terms of the GNU General Public License (Version 2) as
14 * published by the Free Software Foundation.
15 *
16 * This program is distributed in the hope it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19 * for more details.
20 *
21 * You should have received a copy of the GNU General Public License along
22 * with this program; if not, write to the Free Software Foundation, Inc.,
23 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
24 *
25 * ########################################################################
26 */
27
28
29#include "ieee754.h"
30
31#define DP_EBIAS 1023
32#define DP_EMIN (-1022)
33#define DP_EMAX 1023
34#define DP_MBITS 52
35
36#define SP_EBIAS 127
37#define SP_EMIN (-126)
38#define SP_EMAX 127
39#define SP_MBITS 23
40
41#define DP_MBIT(x) ((u64)1 << (x))
42#define DP_HIDDEN_BIT DP_MBIT(DP_MBITS)
43#define DP_SIGN_BIT DP_MBIT(63)
44
45#define SP_MBIT(x) ((u32)1 << (x))
46#define SP_HIDDEN_BIT SP_MBIT(SP_MBITS)
47#define SP_SIGN_BIT SP_MBIT(31)
48
49
50#define SPSIGN(sp) (sp.parts.sign)
51#define SPBEXP(sp) (sp.parts.bexp)
52#define SPMANT(sp) (sp.parts.mant)
53
54#define DPSIGN(dp) (dp.parts.sign)
55#define DPBEXP(dp) (dp.parts.bexp)
56#define DPMANT(dp) (dp.parts.mant)
57
58#define CLPAIR(x,y) ((x)*6+(y))
59
60#define CLEARCX \
61 (ieee754_csr.cx = 0)
62
63#define SETCX(x) \
64 (ieee754_csr.cx |= (x),ieee754_csr.sx |= (x))
65
66#define SETANDTESTCX(x) \
67 (SETCX(x),ieee754_csr.mx & (x))
68
69#define TSTX() \
70 (ieee754_csr.cx & ieee754_csr.mx)
71
72
73#define COMPXSP \
74 unsigned xm; int xe; int xs; int xc
75
76#define COMPYSP \
77 unsigned ym; int ye; int ys; int yc
78
79#define EXPLODESP(v,vc,vs,ve,vm) \
80{\
81 vs = SPSIGN(v);\
82 ve = SPBEXP(v);\
83 vm = SPMANT(v);\
84 if(ve == SP_EMAX+1+SP_EBIAS){\
85 if(vm == 0)\
86 vc = IEEE754_CLASS_INF;\
87 else if(vm & SP_MBIT(SP_MBITS-1)) \
88 vc = IEEE754_CLASS_SNAN;\
89 else \
90 vc = IEEE754_CLASS_QNAN;\
91 } else if(ve == SP_EMIN-1+SP_EBIAS) {\
92 if(vm) {\
93 ve = SP_EMIN;\
94 vc = IEEE754_CLASS_DNORM;\
95 } else\
96 vc = IEEE754_CLASS_ZERO;\
97 } else {\
98 ve -= SP_EBIAS;\
99 vm |= SP_HIDDEN_BIT;\
100 vc = IEEE754_CLASS_NORM;\
101 }\
102}
103#define EXPLODEXSP EXPLODESP(x,xc,xs,xe,xm)
104#define EXPLODEYSP EXPLODESP(y,yc,ys,ye,ym)
105
106
107#define COMPXDP \
108u64 xm; int xe; int xs; int xc
109
110#define COMPYDP \
111u64 ym; int ye; int ys; int yc
112
113#define EXPLODEDP(v,vc,vs,ve,vm) \
114{\
115 vm = DPMANT(v);\
116 vs = DPSIGN(v);\
117 ve = DPBEXP(v);\
118 if(ve == DP_EMAX+1+DP_EBIAS){\
119 if(vm == 0)\
120 vc = IEEE754_CLASS_INF;\
121 else if(vm & DP_MBIT(DP_MBITS-1)) \
122 vc = IEEE754_CLASS_SNAN;\
123 else \
124 vc = IEEE754_CLASS_QNAN;\
125 } else if(ve == DP_EMIN-1+DP_EBIAS) {\
126 if(vm) {\
127 ve = DP_EMIN;\
128 vc = IEEE754_CLASS_DNORM;\
129 } else\
130 vc = IEEE754_CLASS_ZERO;\
131 } else {\
132 ve -= DP_EBIAS;\
133 vm |= DP_HIDDEN_BIT;\
134 vc = IEEE754_CLASS_NORM;\
135 }\
136}
137#define EXPLODEXDP EXPLODEDP(x,xc,xs,xe,xm)
138#define EXPLODEYDP EXPLODEDP(y,yc,ys,ye,ym)
139
140#define FLUSHDP(v,vc,vs,ve,vm) \
141 if(vc==IEEE754_CLASS_DNORM) {\
142 if(ieee754_csr.nod) {\
143 SETCX(IEEE754_INEXACT);\
144 vc = IEEE754_CLASS_ZERO;\
145 ve = DP_EMIN-1+DP_EBIAS;\
146 vm = 0;\
147 v = ieee754dp_zero(vs);\
148 }\
149 }
150
151#define FLUSHSP(v,vc,vs,ve,vm) \
152 if(vc==IEEE754_CLASS_DNORM) {\
153 if(ieee754_csr.nod) {\
154 SETCX(IEEE754_INEXACT);\
155 vc = IEEE754_CLASS_ZERO;\
156 ve = SP_EMIN-1+SP_EBIAS;\
157 vm = 0;\
158 v = ieee754sp_zero(vs);\
159 }\
160 }
161
162#define FLUSHXDP FLUSHDP(x,xc,xs,xe,xm)
163#define FLUSHYDP FLUSHDP(y,yc,ys,ye,ym)
164#define FLUSHXSP FLUSHSP(x,xc,xs,xe,xm)
165#define FLUSHYSP FLUSHSP(y,yc,ys,ye,ym)
diff --git a/arch/mips/math-emu/ieee754m.c b/arch/mips/math-emu/ieee754m.c
new file mode 100644
index 000000000000..d66896cd8f21
--- /dev/null
+++ b/arch/mips/math-emu/ieee754m.c
@@ -0,0 +1,56 @@
1/*
2 * floor, trunc, ceil
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754.h"
29
30ieee754dp ieee754dp_floor(ieee754dp x)
31{
32 ieee754dp i;
33
34 if (ieee754dp_lt(ieee754dp_modf(x, &i), ieee754dp_zero(0)))
35 return ieee754dp_sub(i, ieee754dp_one(0));
36 else
37 return i;
38}
39
40ieee754dp ieee754dp_ceil(ieee754dp x)
41{
42 ieee754dp i;
43
44 if (ieee754dp_gt(ieee754dp_modf(x, &i), ieee754dp_zero(0)))
45 return ieee754dp_add(i, ieee754dp_one(0));
46 else
47 return i;
48}
49
50ieee754dp ieee754dp_trunc(ieee754dp x)
51{
52 ieee754dp i;
53
54 (void) ieee754dp_modf(x, &i);
55 return i;
56}
diff --git a/arch/mips/math-emu/ieee754sp.c b/arch/mips/math-emu/ieee754sp.c
new file mode 100644
index 000000000000..adda851cd04f
--- /dev/null
+++ b/arch/mips/math-emu/ieee754sp.c
@@ -0,0 +1,243 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30int ieee754sp_class(ieee754sp x)
31{
32 COMPXSP;
33 EXPLODEXSP;
34 return xc;
35}
36
37int ieee754sp_isnan(ieee754sp x)
38{
39 return ieee754sp_class(x) >= IEEE754_CLASS_SNAN;
40}
41
42int ieee754sp_issnan(ieee754sp x)
43{
44 assert(ieee754sp_isnan(x));
45 return (SPMANT(x) & SP_MBIT(SP_MBITS-1));
46}
47
48
49ieee754sp ieee754sp_xcpt(ieee754sp r, const char *op, ...)
50{
51 struct ieee754xctx ax;
52
53 if (!TSTX())
54 return r;
55
56 ax.op = op;
57 ax.rt = IEEE754_RT_SP;
58 ax.rv.sp = r;
59 va_start(ax.ap, op);
60 ieee754_xcpt(&ax);
61 return ax.rv.sp;
62}
63
64ieee754sp ieee754sp_nanxcpt(ieee754sp r, const char *op, ...)
65{
66 struct ieee754xctx ax;
67
68 assert(ieee754sp_isnan(r));
69
70 if (!ieee754sp_issnan(r)) /* QNAN does not cause invalid op !! */
71 return r;
72
73 if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) {
74 /* not enabled convert to a quiet NaN */
75 SPMANT(r) &= (~SP_MBIT(SP_MBITS-1));
76 if (ieee754sp_isnan(r))
77 return r;
78 else
79 return ieee754sp_indef();
80 }
81
82 ax.op = op;
83 ax.rt = 0;
84 ax.rv.sp = r;
85 va_start(ax.ap, op);
86 ieee754_xcpt(&ax);
87 return ax.rv.sp;
88}
89
90ieee754sp ieee754sp_bestnan(ieee754sp x, ieee754sp y)
91{
92 assert(ieee754sp_isnan(x));
93 assert(ieee754sp_isnan(y));
94
95 if (SPMANT(x) > SPMANT(y))
96 return x;
97 else
98 return y;
99}
100
101
102static unsigned get_rounding(int sn, unsigned xm)
103{
104 /* inexact must round of 3 bits
105 */
106 if (xm & (SP_MBIT(3) - 1)) {
107 switch (ieee754_csr.rm) {
108 case IEEE754_RZ:
109 break;
110 case IEEE754_RN:
111 xm += 0x3 + ((xm >> 3) & 1);
112 /* xm += (xm&0x8)?0x4:0x3 */
113 break;
114 case IEEE754_RU: /* toward +Infinity */
115 if (!sn) /* ?? */
116 xm += 0x8;
117 break;
118 case IEEE754_RD: /* toward -Infinity */
119 if (sn) /* ?? */
120 xm += 0x8;
121 break;
122 }
123 }
124 return xm;
125}
126
127
128/* generate a normal/denormal number with over,under handling
129 * sn is sign
130 * xe is an unbiased exponent
131 * xm is 3bit extended precision value.
132 */
133ieee754sp ieee754sp_format(int sn, int xe, unsigned xm)
134{
135 assert(xm); /* we don't gen exact zeros (probably should) */
136
137 assert((xm >> (SP_MBITS + 1 + 3)) == 0); /* no execess */
138 assert(xm & (SP_HIDDEN_BIT << 3));
139
140 if (xe < SP_EMIN) {
141 /* strip lower bits */
142 int es = SP_EMIN - xe;
143
144 if (ieee754_csr.nod) {
145 SETCX(IEEE754_UNDERFLOW);
146 SETCX(IEEE754_INEXACT);
147
148 switch(ieee754_csr.rm) {
149 case IEEE754_RN:
150 return ieee754sp_zero(sn);
151 case IEEE754_RZ:
152 return ieee754sp_zero(sn);
153 case IEEE754_RU: /* toward +Infinity */
154 if(sn == 0)
155 return ieee754sp_min(0);
156 else
157 return ieee754sp_zero(1);
158 case IEEE754_RD: /* toward -Infinity */
159 if(sn == 0)
160 return ieee754sp_zero(0);
161 else
162 return ieee754sp_min(1);
163 }
164 }
165
166 if (xe == SP_EMIN - 1
167 && get_rounding(sn, xm) >> (SP_MBITS + 1 + 3))
168 {
169 /* Not tiny after rounding */
170 SETCX(IEEE754_INEXACT);
171 xm = get_rounding(sn, xm);
172 xm >>= 1;
173 /* Clear grs bits */
174 xm &= ~(SP_MBIT(3) - 1);
175 xe++;
176 }
177 else {
178 /* sticky right shift es bits
179 */
180 SPXSRSXn(es);
181 assert((xm & (SP_HIDDEN_BIT << 3)) == 0);
182 assert(xe == SP_EMIN);
183 }
184 }
185 if (xm & (SP_MBIT(3) - 1)) {
186 SETCX(IEEE754_INEXACT);
187 if ((xm & (SP_HIDDEN_BIT << 3)) == 0) {
188 SETCX(IEEE754_UNDERFLOW);
189 }
190
191 /* inexact must round of 3 bits
192 */
193 xm = get_rounding(sn, xm);
194 /* adjust exponent for rounding add overflowing
195 */
196 if (xm >> (SP_MBITS + 1 + 3)) {
197 /* add causes mantissa overflow */
198 xm >>= 1;
199 xe++;
200 }
201 }
202 /* strip grs bits */
203 xm >>= 3;
204
205 assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */
206 assert(xe >= SP_EMIN);
207
208 if (xe > SP_EMAX) {
209 SETCX(IEEE754_OVERFLOW);
210 SETCX(IEEE754_INEXACT);
211 /* -O can be table indexed by (rm,sn) */
212 switch (ieee754_csr.rm) {
213 case IEEE754_RN:
214 return ieee754sp_inf(sn);
215 case IEEE754_RZ:
216 return ieee754sp_max(sn);
217 case IEEE754_RU: /* toward +Infinity */
218 if (sn == 0)
219 return ieee754sp_inf(0);
220 else
221 return ieee754sp_max(1);
222 case IEEE754_RD: /* toward -Infinity */
223 if (sn == 0)
224 return ieee754sp_max(0);
225 else
226 return ieee754sp_inf(1);
227 }
228 }
229 /* gen norm/denorm/zero */
230
231 if ((xm & SP_HIDDEN_BIT) == 0) {
232 /* we underflow (tiny/zero) */
233 assert(xe == SP_EMIN);
234 if (ieee754_csr.mx & IEEE754_UNDERFLOW)
235 SETCX(IEEE754_UNDERFLOW);
236 return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm);
237 } else {
238 assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */
239 assert(xm & SP_HIDDEN_BIT);
240
241 return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
242 }
243}
diff --git a/arch/mips/math-emu/ieee754sp.h b/arch/mips/math-emu/ieee754sp.h
new file mode 100644
index 000000000000..ae82f51297e5
--- /dev/null
+++ b/arch/mips/math-emu/ieee754sp.h
@@ -0,0 +1,89 @@
1/*
2 * IEEE754 floating point
3 * double precision internal header file
4 */
5/*
6 * MIPS floating point support
7 * Copyright (C) 1994-2000 Algorithmics Ltd.
8 * http://www.algor.co.uk
9 *
10 * ########################################################################
11 *
12 * This program is free software; you can distribute it and/or modify it
13 * under the terms of the GNU General Public License (Version 2) as
14 * published by the Free Software Foundation.
15 *
16 * This program is distributed in the hope it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19 * for more details.
20 *
21 * You should have received a copy of the GNU General Public License along
22 * with this program; if not, write to the Free Software Foundation, Inc.,
23 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
24 *
25 * ########################################################################
26 */
27
28
29#include "ieee754int.h"
30
31#define assert(expr) ((void)0)
32
33/* 3bit extended single precision sticky right shift */
34#define SPXSRSXn(rs) \
35 (xe += rs, \
36 xm = (rs > (SP_MBITS+3))?1:((xm) >> (rs)) | ((xm) << (32-(rs)) != 0))
37
38#define SPXSRSX1() \
39 (xe++, (xm = (xm >> 1) | (xm & 1)))
40
41#define SPXSRSYn(rs) \
42 (ye+=rs, \
43 ym = (rs > (SP_MBITS+3))?1:((ym) >> (rs)) | ((ym) << (32-(rs)) != 0))
44
45#define SPXSRSY1() \
46 (ye++, (ym = (ym >> 1) | (ym & 1)))
47
48/* convert denormal to normalized with extended exponent */
49#define SPDNORMx(m,e) \
50 while( (m >> SP_MBITS) == 0) { m <<= 1; e--; }
51#define SPDNORMX SPDNORMx(xm,xe)
52#define SPDNORMY SPDNORMx(ym,ye)
53
54static __inline ieee754sp buildsp(int s, int bx, unsigned m)
55{
56 ieee754sp r;
57
58 assert((s) == 0 || (s) == 1);
59 assert((bx) >= SP_EMIN - 1 + SP_EBIAS
60 && (bx) <= SP_EMAX + 1 + SP_EBIAS);
61 assert(((m) >> SP_MBITS) == 0);
62
63 r.parts.sign = s;
64 r.parts.bexp = bx;
65 r.parts.mant = m;
66
67 return r;
68}
69
70extern int ieee754sp_isnan(ieee754sp);
71extern int ieee754sp_issnan(ieee754sp);
72extern int ieee754si_xcpt(int, const char *, ...);
73extern s64 ieee754di_xcpt(s64, const char *, ...);
74extern ieee754sp ieee754sp_xcpt(ieee754sp, const char *, ...);
75extern ieee754sp ieee754sp_nanxcpt(ieee754sp, const char *, ...);
76extern ieee754sp ieee754sp_bestnan(ieee754sp, ieee754sp);
77extern ieee754sp ieee754sp_format(int, int, unsigned);
78
79
80#define SPNORMRET2(s,e,m,name,a0,a1) \
81{ \
82 ieee754sp V = ieee754sp_format(s,e,m); \
83 if(TSTX()) \
84 return ieee754sp_xcpt(V,name,a0,a1); \
85 else \
86 return V; \
87}
88
89#define SPNORMRET1(s,e,m,name,a0) SPNORMRET2(s,e,m,name,a0,a0)
diff --git a/arch/mips/math-emu/ieee754xcpt.c b/arch/mips/math-emu/ieee754xcpt.c
new file mode 100644
index 000000000000..7d8ef8965067
--- /dev/null
+++ b/arch/mips/math-emu/ieee754xcpt.c
@@ -0,0 +1,49 @@
1/*
2 * MIPS floating point support
3 * Copyright (C) 1994-2000 Algorithmics Ltd.
4 * http://www.algor.co.uk
5 *
6 * ########################################################################
7 *
8 * This program is free software; you can distribute it and/or modify it
9 * under the terms of the GNU General Public License (Version 2) as
10 * published by the Free Software Foundation.
11 *
12 * This program is distributed in the hope it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 * for more details.
16 *
17 * You should have received a copy of the GNU General Public License along
18 * with this program; if not, write to the Free Software Foundation, Inc.,
19 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
20 *
21 * ########################################################################
22 */
23
24/**************************************************************************
25 * Nov 7, 2000
26 * Added preprocessor hacks to map to Linux kernel diagnostics.
27 *
28 * Kevin D. Kissell, kevink@mips.com and Carsten Langgaard, carstenl@mips.com
29 * Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
30 *************************************************************************/
31
32#include <linux/kernel.h>
33#include "ieee754.h"
34
35/*
36 * Very naff exception handler (you can plug in your own and
37 * override this).
38 */
39
40static const char *const rtnames[] = {
41 "sp", "dp", "xp", "si", "di"
42};
43
44void ieee754_xcpt(struct ieee754xctx *xcp)
45{
46 printk(KERN_DEBUG "floating point exception in \"%s\", type=%s\n",
47 xcp->op, rtnames[xcp->rt]);
48}
49
diff --git a/arch/mips/math-emu/kernel_linkage.c b/arch/mips/math-emu/kernel_linkage.c
new file mode 100644
index 000000000000..04397fec30fc
--- /dev/null
+++ b/arch/mips/math-emu/kernel_linkage.c
@@ -0,0 +1,125 @@
1/*
2 * Kevin D. Kissell, kevink@mips and Carsten Langgaard, carstenl@mips.com
3 * Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
4 *
5 * This program is free software; you can distribute it and/or modify it
6 * under the terms of the GNU General Public License (Version 2) as
7 * published by the Free Software Foundation.
8 *
9 * This program is distributed in the hope it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
12 * for more details.
13 *
14 * You should have received a copy of the GNU General Public License along
15 * with this program; if not, write to the Free Software Foundation, Inc.,
16 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
17 *
18 * Routines corresponding to Linux kernel FP context
19 * manipulation primitives for the Algorithmics MIPS
20 * FPU Emulator
21 */
22#include <linux/config.h>
23#include <linux/sched.h>
24#include <asm/processor.h>
25#include <asm/signal.h>
26#include <asm/uaccess.h>
27
28#include <asm/fpu_emulator.h>
29
30extern struct mips_fpu_emulator_private fpuemuprivate;
31
32#define SIGNALLING_NAN 0x7ff800007ff80000LL
33
34void fpu_emulator_init_fpu(void)
35{
36 static int first = 1;
37 int i;
38
39 if (first) {
40 first = 0;
41 printk("Algorithmics/MIPS FPU Emulator v1.5\n");
42 }
43
44 current->thread.fpu.soft.fcr31 = 0;
45 for (i = 0; i < 32; i++) {
46 current->thread.fpu.soft.fpr[i] = SIGNALLING_NAN;
47 }
48}
49
50
51/*
52 * Emulator context save/restore to/from a signal context
53 * presumed to be on the user stack, and therefore accessed
54 * with appropriate macros from uaccess.h
55 */
56
57int fpu_emulator_save_context(struct sigcontext *sc)
58{
59 int i;
60 int err = 0;
61
62 for (i = 0; i < 32; i++) {
63 err |=
64 __put_user(current->thread.fpu.soft.fpr[i],
65 &sc->sc_fpregs[i]);
66 }
67 err |= __put_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr);
68 err |= __put_user(fpuemuprivate.eir, &sc->sc_fpc_eir);
69
70 return err;
71}
72
73int fpu_emulator_restore_context(struct sigcontext *sc)
74{
75 int i;
76 int err = 0;
77
78 for (i = 0; i < 32; i++) {
79 err |=
80 __get_user(current->thread.fpu.soft.fpr[i],
81 &sc->sc_fpregs[i]);
82 }
83 err |= __get_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr);
84 err |= __get_user(fpuemuprivate.eir, &sc->sc_fpc_eir);
85
86 return err;
87}
88
89#ifdef CONFIG_MIPS64
90/*
91 * This is the o32 version
92 */
93
94int fpu_emulator_save_context32(struct sigcontext32 *sc)
95{
96 int i;
97 int err = 0;
98
99 for (i = 0; i < 32; i+=2) {
100 err |=
101 __put_user(current->thread.fpu.soft.fpr[i],
102 &sc->sc_fpregs[i]);
103 }
104 err |= __put_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr);
105 err |= __put_user(fpuemuprivate.eir, &sc->sc_fpc_eir);
106
107 return err;
108}
109
110int fpu_emulator_restore_context32(struct sigcontext32 *sc)
111{
112 int i;
113 int err = 0;
114
115 for (i = 0; i < 32; i+=2) {
116 err |=
117 __get_user(current->thread.fpu.soft.fpr[i],
118 &sc->sc_fpregs[i]);
119 }
120 err |= __get_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr);
121 err |= __get_user(fpuemuprivate.eir, &sc->sc_fpc_eir);
122
123 return err;
124}
125#endif
diff --git a/arch/mips/math-emu/sp_add.c b/arch/mips/math-emu/sp_add.c
new file mode 100644
index 000000000000..d8c4211bcfbe
--- /dev/null
+++ b/arch/mips/math-emu/sp_add.c
@@ -0,0 +1,177 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_add(ieee754sp x, ieee754sp y)
31{
32 COMPXSP;
33 COMPYSP;
34
35 EXPLODEXSP;
36 EXPLODEYSP;
37
38 CLEARCX;
39
40 FLUSHXSP;
41 FLUSHYSP;
42
43 switch (CLPAIR(xc, yc)) {
44 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
45 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
46 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
47 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
48 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
49 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
50 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
51 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
52 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
53 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
54 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754sp_nanxcpt(ieee754sp_indef(), "add", x, y);
57
58 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
59 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
60 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
61 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
62 return y;
63
64 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
65 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
66 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
67 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
68 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
69 return x;
70
71
72 /* Infinity handling
73 */
74
75 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
76 if (xs == ys)
77 return x;
78 SETCX(IEEE754_INVALID_OPERATION);
79 return ieee754sp_xcpt(ieee754sp_indef(), "add", x, y);
80
81 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
82 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
83 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
84 return y;
85
86 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
87 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
88 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
89 return x;
90
91 /* Zero handling
92 */
93
94 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
95 if (xs == ys)
96 return x;
97 else
98 return ieee754sp_zero(ieee754_csr.rm ==
99 IEEE754_RD);
100
101 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
102 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
103 return x;
104
105 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
106 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
107 return y;
108
109 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
110 SPDNORMX;
111
112 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
113 SPDNORMY;
114 break;
115
116 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
117 SPDNORMX;
118 break;
119
120 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
121 break;
122 }
123 assert(xm & SP_HIDDEN_BIT);
124 assert(ym & SP_HIDDEN_BIT);
125
126 /* provide guard,round and stick bit space */
127 xm <<= 3;
128 ym <<= 3;
129
130 if (xe > ye) {
131 /* have to shift y fraction right to align
132 */
133 int s = xe - ye;
134 SPXSRSYn(s);
135 } else if (ye > xe) {
136 /* have to shift x fraction right to align
137 */
138 int s = ye - xe;
139 SPXSRSXn(s);
140 }
141 assert(xe == ye);
142 assert(xe <= SP_EMAX);
143
144 if (xs == ys) {
145 /* generate 28 bit result of adding two 27 bit numbers
146 * leaving result in xm,xs,xe
147 */
148 xm = xm + ym;
149 xe = xe;
150 xs = xs;
151
152 if (xm >> (SP_MBITS + 1 + 3)) { /* carry out */
153 SPXSRSX1();
154 }
155 } else {
156 if (xm >= ym) {
157 xm = xm - ym;
158 xe = xe;
159 xs = xs;
160 } else {
161 xm = ym - xm;
162 xe = xe;
163 xs = ys;
164 }
165 if (xm == 0)
166 return ieee754sp_zero(ieee754_csr.rm ==
167 IEEE754_RD);
168
169 /* normalize in extended single precision */
170 while ((xm >> (SP_MBITS + 3)) == 0) {
171 xm <<= 1;
172 xe--;
173 }
174
175 }
176 SPNORMRET2(xs, xe, xm, "add", x, y);
177}
diff --git a/arch/mips/math-emu/sp_cmp.c b/arch/mips/math-emu/sp_cmp.c
new file mode 100644
index 000000000000..d3eff6b04b5a
--- /dev/null
+++ b/arch/mips/math-emu/sp_cmp.c
@@ -0,0 +1,67 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30int ieee754sp_cmp(ieee754sp x, ieee754sp y, int cmp, int sig)
31{
32 COMPXSP;
33 COMPYSP;
34
35 EXPLODEXSP;
36 EXPLODEYSP;
37 FLUSHXSP;
38 FLUSHYSP;
39 CLEARCX; /* Even clear inexact flag here */
40
41 if (ieee754sp_isnan(x) || ieee754sp_isnan(y)) {
42 if (sig || xc == IEEE754_CLASS_SNAN || yc == IEEE754_CLASS_SNAN)
43 SETCX(IEEE754_INVALID_OPERATION);
44 if (cmp & IEEE754_CUN)
45 return 1;
46 if (cmp & (IEEE754_CLT | IEEE754_CGT)) {
47 if (sig && SETANDTESTCX(IEEE754_INVALID_OPERATION))
48 return ieee754si_xcpt(0, "fcmpf", x);
49 }
50 return 0;
51 } else {
52 int vx = x.bits;
53 int vy = y.bits;
54
55 if (vx < 0)
56 vx = -vx ^ SP_SIGN_BIT;
57 if (vy < 0)
58 vy = -vy ^ SP_SIGN_BIT;
59
60 if (vx < vy)
61 return (cmp & IEEE754_CLT) != 0;
62 else if (vx == vy)
63 return (cmp & IEEE754_CEQ) != 0;
64 else
65 return (cmp & IEEE754_CGT) != 0;
66 }
67}
diff --git a/arch/mips/math-emu/sp_div.c b/arch/mips/math-emu/sp_div.c
new file mode 100644
index 000000000000..2b437fcfdad9
--- /dev/null
+++ b/arch/mips/math-emu/sp_div.c
@@ -0,0 +1,157 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_div(ieee754sp x, ieee754sp y)
31{
32 COMPXSP;
33 COMPYSP;
34
35 EXPLODEXSP;
36 EXPLODEYSP;
37
38 CLEARCX;
39
40 FLUSHXSP;
41 FLUSHYSP;
42
43 switch (CLPAIR(xc, yc)) {
44 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
45 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
46 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
47 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
48 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
49 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
50 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
51 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
52 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
53 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
54 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754sp_nanxcpt(ieee754sp_indef(), "div", x, y);
57
58 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
59 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
60 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
61 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
62 return y;
63
64 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
65 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
66 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
67 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
68 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
69 return x;
70
71
72 /* Infinity handling
73 */
74
75 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
76 SETCX(IEEE754_INVALID_OPERATION);
77 return ieee754sp_xcpt(ieee754sp_indef(), "div", x, y);
78
79 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
80 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
81 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
82 return ieee754sp_zero(xs ^ ys);
83
84 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
85 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
86 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
87 return ieee754sp_inf(xs ^ ys);
88
89 /* Zero handling
90 */
91
92 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
93 SETCX(IEEE754_INVALID_OPERATION);
94 return ieee754sp_xcpt(ieee754sp_indef(), "div", x, y);
95
96 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
97 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
98 SETCX(IEEE754_ZERO_DIVIDE);
99 return ieee754sp_xcpt(ieee754sp_inf(xs ^ ys), "div", x, y);
100
101 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
102 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
103 return ieee754sp_zero(xs == ys ? 0 : 1);
104
105 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
106 SPDNORMX;
107
108 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
109 SPDNORMY;
110 break;
111
112 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
113 SPDNORMX;
114 break;
115
116 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
117 break;
118 }
119 assert(xm & SP_HIDDEN_BIT);
120 assert(ym & SP_HIDDEN_BIT);
121
122 /* provide rounding space */
123 xm <<= 3;
124 ym <<= 3;
125
126 {
127 /* now the dirty work */
128
129 unsigned rm = 0;
130 int re = xe - ye;
131 unsigned bm;
132
133 for (bm = SP_MBIT(SP_MBITS + 2); bm; bm >>= 1) {
134 if (xm >= ym) {
135 xm -= ym;
136 rm |= bm;
137 if (xm == 0)
138 break;
139 }
140 xm <<= 1;
141 }
142 rm <<= 1;
143 if (xm)
144 rm |= 1; /* have remainder, set sticky */
145
146 assert(rm);
147
148 /* normalise rm to rounding precision ?
149 */
150 while ((rm >> (SP_MBITS + 3)) == 0) {
151 rm <<= 1;
152 re--;
153 }
154
155 SPNORMRET2(xs == ys ? 0 : 1, re, rm, "div", x, y);
156 }
157}
diff --git a/arch/mips/math-emu/sp_fdp.c b/arch/mips/math-emu/sp_fdp.c
new file mode 100644
index 000000000000..4093723d1aa5
--- /dev/null
+++ b/arch/mips/math-emu/sp_fdp.c
@@ -0,0 +1,77 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_fdp(ieee754dp x)
31{
32 COMPXDP;
33 ieee754sp nan;
34
35 EXPLODEXDP;
36
37 CLEARCX;
38
39 FLUSHXDP;
40
41 switch (xc) {
42 case IEEE754_CLASS_SNAN:
43 SETCX(IEEE754_INVALID_OPERATION);
44 return ieee754sp_nanxcpt(ieee754sp_indef(), "fdp");
45 case IEEE754_CLASS_QNAN:
46 nan = buildsp(xs, SP_EMAX + 1 + SP_EBIAS, (u32)
47 (xm >> (DP_MBITS - SP_MBITS)));
48 if (!ieee754sp_isnan(nan))
49 nan = ieee754sp_indef();
50 return ieee754sp_nanxcpt(nan, "fdp", x);
51 case IEEE754_CLASS_INF:
52 return ieee754sp_inf(xs);
53 case IEEE754_CLASS_ZERO:
54 return ieee754sp_zero(xs);
55 case IEEE754_CLASS_DNORM:
56 /* can't possibly be sp representable */
57 SETCX(IEEE754_UNDERFLOW);
58 SETCX(IEEE754_INEXACT);
59 if ((ieee754_csr.rm == IEEE754_RU && !xs) ||
60 (ieee754_csr.rm == IEEE754_RD && xs))
61 return ieee754sp_xcpt(ieee754sp_mind(xs), "fdp", x);
62 return ieee754sp_xcpt(ieee754sp_zero(xs), "fdp", x);
63 case IEEE754_CLASS_NORM:
64 break;
65 }
66
67 {
68 u32 rm;
69
70 /* convert from DP_MBITS to SP_MBITS+3 with sticky right shift
71 */
72 rm = (xm >> (DP_MBITS - (SP_MBITS + 3))) |
73 ((xm << (64 - (DP_MBITS - (SP_MBITS + 3)))) != 0);
74
75 SPNORMRET1(xs, xe, rm, "fdp", x);
76 }
77}
diff --git a/arch/mips/math-emu/sp_fint.c b/arch/mips/math-emu/sp_fint.c
new file mode 100644
index 000000000000..42d9ed4b9a94
--- /dev/null
+++ b/arch/mips/math-emu/sp_fint.c
@@ -0,0 +1,80 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_fint(int x)
31{
32 COMPXSP;
33
34 CLEARCX;
35
36 xc = ( 0 ? xc : xc );
37
38 if (x == 0)
39 return ieee754sp_zero(0);
40 if (x == 1 || x == -1)
41 return ieee754sp_one(x < 0);
42 if (x == 10 || x == -10)
43 return ieee754sp_ten(x < 0);
44
45 xs = (x < 0);
46 if (xs) {
47 if (x == (1 << 31))
48 xm = ((unsigned) 1 << 31); /* max neg can't be safely negated */
49 else
50 xm = -x;
51 } else {
52 xm = x;
53 }
54 xe = SP_MBITS + 3;
55
56 if (xm >> (SP_MBITS + 1 + 3)) {
57 /* shunt out overflow bits
58 */
59 while (xm >> (SP_MBITS + 1 + 3)) {
60 SPXSRSX1();
61 }
62 } else {
63 /* normalize in grs extended single precision
64 */
65 while ((xm >> (SP_MBITS + 3)) == 0) {
66 xm <<= 1;
67 xe--;
68 }
69 }
70 SPNORMRET1(xs, xe, xm, "fint", x);
71}
72
73
74ieee754sp ieee754sp_funs(unsigned int u)
75{
76 if ((int) u < 0)
77 return ieee754sp_add(ieee754sp_1e31(),
78 ieee754sp_fint(u & ~(1 << 31)));
79 return ieee754sp_fint(u);
80}
diff --git a/arch/mips/math-emu/sp_flong.c b/arch/mips/math-emu/sp_flong.c
new file mode 100644
index 000000000000..1e26795ccecb
--- /dev/null
+++ b/arch/mips/math-emu/sp_flong.c
@@ -0,0 +1,79 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_flong(s64 x)
31{
32 COMPXDP; /* <--- need 64-bit mantissa temp */
33
34 CLEARCX;
35
36 xc = ( 0 ? xc : xc );
37
38 if (x == 0)
39 return ieee754sp_zero(0);
40 if (x == 1 || x == -1)
41 return ieee754sp_one(x < 0);
42 if (x == 10 || x == -10)
43 return ieee754sp_ten(x < 0);
44
45 xs = (x < 0);
46 if (xs) {
47 if (x == (1ULL << 63))
48 xm = (1ULL << 63); /* max neg can't be safely negated */
49 else
50 xm = -x;
51 } else {
52 xm = x;
53 }
54 xe = SP_MBITS + 3;
55
56 if (xm >> (SP_MBITS + 1 + 3)) {
57 /* shunt out overflow bits
58 */
59 while (xm >> (SP_MBITS + 1 + 3)) {
60 SPXSRSX1();
61 }
62 } else {
63 /* normalize in grs extended single precision */
64 while ((xm >> (SP_MBITS + 3)) == 0) {
65 xm <<= 1;
66 xe--;
67 }
68 }
69 SPNORMRET1(xs, xe, xm, "sp_flong", x);
70}
71
72
73ieee754sp ieee754sp_fulong(u64 u)
74{
75 if ((s64) u < 0)
76 return ieee754sp_add(ieee754sp_1e63(),
77 ieee754sp_flong(u & ~(1ULL << 63)));
78 return ieee754sp_flong(u);
79}
diff --git a/arch/mips/math-emu/sp_frexp.c b/arch/mips/math-emu/sp_frexp.c
new file mode 100644
index 000000000000..359c6483dbfa
--- /dev/null
+++ b/arch/mips/math-emu/sp_frexp.c
@@ -0,0 +1,53 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30/* close to ieeep754sp_logb
31*/
32ieee754sp ieee754sp_frexp(ieee754sp x, int *eptr)
33{
34 COMPXSP;
35 CLEARCX;
36 EXPLODEXSP;
37
38 switch (xc) {
39 case IEEE754_CLASS_SNAN:
40 case IEEE754_CLASS_QNAN:
41 case IEEE754_CLASS_INF:
42 case IEEE754_CLASS_ZERO:
43 *eptr = 0;
44 return x;
45 case IEEE754_CLASS_DNORM:
46 SPDNORMX;
47 break;
48 case IEEE754_CLASS_NORM:
49 break;
50 }
51 *eptr = xe + 1;
52 return buildsp(xs, -1 + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
53}
diff --git a/arch/mips/math-emu/sp_logb.c b/arch/mips/math-emu/sp_logb.c
new file mode 100644
index 000000000000..3c337219ca32
--- /dev/null
+++ b/arch/mips/math-emu/sp_logb.c
@@ -0,0 +1,54 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_logb(ieee754sp x)
31{
32 COMPXSP;
33
34 CLEARCX;
35
36 EXPLODEXSP;
37
38 switch (xc) {
39 case IEEE754_CLASS_SNAN:
40 return ieee754sp_nanxcpt(x, "logb", x);
41 case IEEE754_CLASS_QNAN:
42 return x;
43 case IEEE754_CLASS_INF:
44 return ieee754sp_inf(0);
45 case IEEE754_CLASS_ZERO:
46 return ieee754sp_inf(1);
47 case IEEE754_CLASS_DNORM:
48 SPDNORMX;
49 break;
50 case IEEE754_CLASS_NORM:
51 break;
52 }
53 return ieee754sp_fint(xe);
54}
diff --git a/arch/mips/math-emu/sp_modf.c b/arch/mips/math-emu/sp_modf.c
new file mode 100644
index 000000000000..4b1dbac796f8
--- /dev/null
+++ b/arch/mips/math-emu/sp_modf.c
@@ -0,0 +1,80 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30/* modf function is always exact for a finite number
31*/
32ieee754sp ieee754sp_modf(ieee754sp x, ieee754sp * ip)
33{
34 COMPXSP;
35
36 CLEARCX;
37
38 EXPLODEXSP;
39
40 switch (xc) {
41 case IEEE754_CLASS_SNAN:
42 case IEEE754_CLASS_QNAN:
43 case IEEE754_CLASS_INF:
44 case IEEE754_CLASS_ZERO:
45 *ip = x;
46 return x;
47 case IEEE754_CLASS_DNORM:
48 /* far to small */
49 *ip = ieee754sp_zero(xs);
50 return x;
51 case IEEE754_CLASS_NORM:
52 break;
53 }
54 if (xe < 0) {
55 *ip = ieee754sp_zero(xs);
56 return x;
57 }
58 if (xe >= SP_MBITS) {
59 *ip = x;
60 return ieee754sp_zero(xs);
61 }
62 /* generate ipart mantissa by clearing bottom bits
63 */
64 *ip = buildsp(xs, xe + SP_EBIAS,
65 ((xm >> (SP_MBITS - xe)) << (SP_MBITS - xe)) &
66 ~SP_HIDDEN_BIT);
67
68 /* generate fpart mantissa by clearing top bits
69 * and normalizing (must be able to normalize)
70 */
71 xm = (xm << (32 - (SP_MBITS - xe))) >> (32 - (SP_MBITS - xe));
72 if (xm == 0)
73 return ieee754sp_zero(xs);
74
75 while ((xm >> SP_MBITS) == 0) {
76 xm <<= 1;
77 xe--;
78 }
79 return buildsp(xs, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
80}
diff --git a/arch/mips/math-emu/sp_mul.c b/arch/mips/math-emu/sp_mul.c
new file mode 100644
index 000000000000..3f070f82212f
--- /dev/null
+++ b/arch/mips/math-emu/sp_mul.c
@@ -0,0 +1,171 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_mul(ieee754sp x, ieee754sp y)
31{
32 COMPXSP;
33 COMPYSP;
34
35 EXPLODEXSP;
36 EXPLODEYSP;
37
38 CLEARCX;
39
40 FLUSHXSP;
41 FLUSHYSP;
42
43 switch (CLPAIR(xc, yc)) {
44 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
45 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
46 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
47 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
48 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
49 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
50 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
51 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
52 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
53 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
54 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754sp_nanxcpt(ieee754sp_indef(), "mul", x, y);
57
58 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
59 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
60 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
61 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
62 return y;
63
64 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
65 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
66 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
67 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
68 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
69 return x;
70
71
72 /* Infinity handling */
73
74 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
75 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
76 SETCX(IEEE754_INVALID_OPERATION);
77 return ieee754sp_xcpt(ieee754sp_indef(), "mul", x, y);
78
79 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
80 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
81 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
82 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
83 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
84 return ieee754sp_inf(xs ^ ys);
85
86 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
87 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
88 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
89 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
90 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
91 return ieee754sp_zero(xs ^ ys);
92
93
94 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
95 SPDNORMX;
96
97 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
98 SPDNORMY;
99 break;
100
101 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
102 SPDNORMX;
103 break;
104
105 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
106 break;
107 }
108 /* rm = xm * ym, re = xe+ye basicly */
109 assert(xm & SP_HIDDEN_BIT);
110 assert(ym & SP_HIDDEN_BIT);
111
112 {
113 int re = xe + ye;
114 int rs = xs ^ ys;
115 unsigned rm;
116
117 /* shunt to top of word */
118 xm <<= 32 - (SP_MBITS + 1);
119 ym <<= 32 - (SP_MBITS + 1);
120
121 /* multiply 32bits xm,ym to give high 32bits rm with stickness
122 */
123 {
124 unsigned short lxm = xm & 0xffff;
125 unsigned short hxm = xm >> 16;
126 unsigned short lym = ym & 0xffff;
127 unsigned short hym = ym >> 16;
128 unsigned lrm;
129 unsigned hrm;
130
131 lrm = lxm * lym; /* 16 * 16 => 32 */
132 hrm = hxm * hym; /* 16 * 16 => 32 */
133
134 {
135 unsigned t = lxm * hym; /* 16 * 16 => 32 */
136 {
137 unsigned at = lrm + (t << 16);
138 hrm += at < lrm;
139 lrm = at;
140 }
141 hrm = hrm + (t >> 16);
142 }
143
144 {
145 unsigned t = hxm * lym; /* 16 * 16 => 32 */
146 {
147 unsigned at = lrm + (t << 16);
148 hrm += at < lrm;
149 lrm = at;
150 }
151 hrm = hrm + (t >> 16);
152 }
153 rm = hrm | (lrm != 0);
154 }
155
156 /*
157 * sticky shift down to normal rounding precision
158 */
159 if ((int) rm < 0) {
160 rm = (rm >> (32 - (SP_MBITS + 1 + 3))) |
161 ((rm << (SP_MBITS + 1 + 3)) != 0);
162 re++;
163 } else {
164 rm = (rm >> (32 - (SP_MBITS + 1 + 3 + 1))) |
165 ((rm << (SP_MBITS + 1 + 3 + 1)) != 0);
166 }
167 assert(rm & (SP_HIDDEN_BIT << 3));
168
169 SPNORMRET2(rs, re, rm, "mul", x, y);
170 }
171}
diff --git a/arch/mips/math-emu/sp_scalb.c b/arch/mips/math-emu/sp_scalb.c
new file mode 100644
index 000000000000..44ceb87ea944
--- /dev/null
+++ b/arch/mips/math-emu/sp_scalb.c
@@ -0,0 +1,58 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_scalb(ieee754sp x, int n)
31{
32 COMPXSP;
33
34 CLEARCX;
35
36 EXPLODEXSP;
37
38 switch (xc) {
39 case IEEE754_CLASS_SNAN:
40 return ieee754sp_nanxcpt(x, "scalb", x, n);
41 case IEEE754_CLASS_QNAN:
42 case IEEE754_CLASS_INF:
43 case IEEE754_CLASS_ZERO:
44 return x;
45 case IEEE754_CLASS_DNORM:
46 SPDNORMX;
47 break;
48 case IEEE754_CLASS_NORM:
49 break;
50 }
51 SPNORMRET2(xs, xe + n, xm << 3, "scalb", x, n);
52}
53
54
55ieee754sp ieee754sp_ldexp(ieee754sp x, int n)
56{
57 return ieee754sp_scalb(x, n);
58}
diff --git a/arch/mips/math-emu/sp_simple.c b/arch/mips/math-emu/sp_simple.c
new file mode 100644
index 000000000000..c809830dffb4
--- /dev/null
+++ b/arch/mips/math-emu/sp_simple.c
@@ -0,0 +1,84 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30int ieee754sp_finite(ieee754sp x)
31{
32 return SPBEXP(x) != SP_EMAX + 1 + SP_EBIAS;
33}
34
35ieee754sp ieee754sp_copysign(ieee754sp x, ieee754sp y)
36{
37 CLEARCX;
38 SPSIGN(x) = SPSIGN(y);
39 return x;
40}
41
42
43ieee754sp ieee754sp_neg(ieee754sp x)
44{
45 COMPXSP;
46
47 EXPLODEXSP;
48 CLEARCX;
49 FLUSHXSP;
50
51 if (xc == IEEE754_CLASS_SNAN) {
52 SETCX(IEEE754_INVALID_OPERATION);
53 return ieee754sp_nanxcpt(ieee754sp_indef(), "neg");
54 }
55
56 if (ieee754sp_isnan(x)) /* but not infinity */
57 return ieee754sp_nanxcpt(x, "neg", x);
58
59 /* quick fix up */
60 SPSIGN(x) ^= 1;
61 return x;
62}
63
64
65ieee754sp ieee754sp_abs(ieee754sp x)
66{
67 COMPXSP;
68
69 EXPLODEXSP;
70 CLEARCX;
71 FLUSHXSP;
72
73 if (xc == IEEE754_CLASS_SNAN) {
74 SETCX(IEEE754_INVALID_OPERATION);
75 return ieee754sp_nanxcpt(ieee754sp_indef(), "abs");
76 }
77
78 if (ieee754sp_isnan(x)) /* but not infinity */
79 return ieee754sp_nanxcpt(x, "abs", x);
80
81 /* quick fix up */
82 SPSIGN(x) = 0;
83 return x;
84}
diff --git a/arch/mips/math-emu/sp_sqrt.c b/arch/mips/math-emu/sp_sqrt.c
new file mode 100644
index 000000000000..8a934b9f7eb8
--- /dev/null
+++ b/arch/mips/math-emu/sp_sqrt.c
@@ -0,0 +1,117 @@
1/* IEEE754 floating point arithmetic
2 * single precision square root
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_sqrt(ieee754sp x)
31{
32 int ix, s, q, m, t, i;
33 unsigned int r;
34 COMPXSP;
35
36 /* take care of Inf and NaN */
37
38 EXPLODEXSP;
39 CLEARCX;
40 FLUSHXSP;
41
42 /* x == INF or NAN? */
43 switch (xc) {
44 case IEEE754_CLASS_QNAN:
45 /* sqrt(Nan) = Nan */
46 return ieee754sp_nanxcpt(x, "sqrt");
47 case IEEE754_CLASS_SNAN:
48 SETCX(IEEE754_INVALID_OPERATION);
49 return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
50 case IEEE754_CLASS_ZERO:
51 /* sqrt(0) = 0 */
52 return x;
53 case IEEE754_CLASS_INF:
54 if (xs) {
55 /* sqrt(-Inf) = Nan */
56 SETCX(IEEE754_INVALID_OPERATION);
57 return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
58 }
59 /* sqrt(+Inf) = Inf */
60 return x;
61 case IEEE754_CLASS_DNORM:
62 case IEEE754_CLASS_NORM:
63 if (xs) {
64 /* sqrt(-x) = Nan */
65 SETCX(IEEE754_INVALID_OPERATION);
66 return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
67 }
68 break;
69 }
70
71 ix = x.bits;
72
73 /* normalize x */
74 m = (ix >> 23);
75 if (m == 0) { /* subnormal x */
76 for (i = 0; (ix & 0x00800000) == 0; i++)
77 ix <<= 1;
78 m -= i - 1;
79 }
80 m -= 127; /* unbias exponent */
81 ix = (ix & 0x007fffff) | 0x00800000;
82 if (m & 1) /* odd m, double x to make it even */
83 ix += ix;
84 m >>= 1; /* m = [m/2] */
85
86 /* generate sqrt(x) bit by bit */
87 ix += ix;
88 q = s = 0; /* q = sqrt(x) */
89 r = 0x01000000; /* r = moving bit from right to left */
90
91 while (r != 0) {
92 t = s + r;
93 if (t <= ix) {
94 s = t + r;
95 ix -= t;
96 q += r;
97 }
98 ix += ix;
99 r >>= 1;
100 }
101
102 if (ix != 0) {
103 SETCX(IEEE754_INEXACT);
104 switch (ieee754_csr.rm) {
105 case IEEE754_RP:
106 q += 2;
107 break;
108 case IEEE754_RN:
109 q += (q & 1);
110 break;
111 }
112 }
113 ix = (q >> 1) + 0x3f000000;
114 ix += (m << 23);
115 x.bits = ix;
116 return x;
117}
diff --git a/arch/mips/math-emu/sp_sub.c b/arch/mips/math-emu/sp_sub.c
new file mode 100644
index 000000000000..dbb802c1a086
--- /dev/null
+++ b/arch/mips/math-emu/sp_sub.c
@@ -0,0 +1,184 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30ieee754sp ieee754sp_sub(ieee754sp x, ieee754sp y)
31{
32 COMPXSP;
33 COMPYSP;
34
35 EXPLODEXSP;
36 EXPLODEYSP;
37
38 CLEARCX;
39
40 FLUSHXSP;
41 FLUSHYSP;
42
43 switch (CLPAIR(xc, yc)) {
44 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
45 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
46 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
47 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
48 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
49 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
50 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
51 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
52 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
53 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
54 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
55 SETCX(IEEE754_INVALID_OPERATION);
56 return ieee754sp_nanxcpt(ieee754sp_indef(), "sub", x, y);
57
58 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
59 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
60 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
61 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
62 return y;
63
64 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
65 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
66 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
67 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
68 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
69 return x;
70
71
72 /* Infinity handling
73 */
74
75 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
76 if (xs != ys)
77 return x;
78 SETCX(IEEE754_INVALID_OPERATION);
79 return ieee754sp_xcpt(ieee754sp_indef(), "sub", x, y);
80
81 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
82 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
83 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
84 return ieee754sp_inf(ys ^ 1);
85
86 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
87 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
88 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
89 return x;
90
91 /* Zero handling
92 */
93
94 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
95 if (xs != ys)
96 return x;
97 else
98 return ieee754sp_zero(ieee754_csr.rm ==
99 IEEE754_RD);
100
101 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
102 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
103 return x;
104
105 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
106 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
107 /* quick fix up */
108 DPSIGN(y) ^= 1;
109 return y;
110
111 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
112 SPDNORMX;
113
114 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
115 SPDNORMY;
116 break;
117
118 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
119 SPDNORMX;
120 break;
121
122 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
123 break;
124 }
125 /* flip sign of y and handle as add */
126 ys ^= 1;
127
128 assert(xm & SP_HIDDEN_BIT);
129 assert(ym & SP_HIDDEN_BIT);
130
131
132 /* provide guard,round and stick bit space */
133 xm <<= 3;
134 ym <<= 3;
135
136 if (xe > ye) {
137 /* have to shift y fraction right to align
138 */
139 int s = xe - ye;
140 SPXSRSYn(s);
141 } else if (ye > xe) {
142 /* have to shift x fraction right to align
143 */
144 int s = ye - xe;
145 SPXSRSXn(s);
146 }
147 assert(xe == ye);
148 assert(xe <= SP_EMAX);
149
150 if (xs == ys) {
151 /* generate 28 bit result of adding two 27 bit numbers
152 */
153 xm = xm + ym;
154 xe = xe;
155 xs = xs;
156
157 if (xm >> (SP_MBITS + 1 + 3)) { /* carry out */
158 SPXSRSX1(); /* shift preserving sticky */
159 }
160 } else {
161 if (xm >= ym) {
162 xm = xm - ym;
163 xe = xe;
164 xs = xs;
165 } else {
166 xm = ym - xm;
167 xe = xe;
168 xs = ys;
169 }
170 if (xm == 0) {
171 if (ieee754_csr.rm == IEEE754_RD)
172 return ieee754sp_zero(1); /* round negative inf. => sign = -1 */
173 else
174 return ieee754sp_zero(0); /* other round modes => sign = 1 */
175 }
176 /* normalize to rounding precision
177 */
178 while ((xm >> (SP_MBITS + 3)) == 0) {
179 xm <<= 1;
180 xe--;
181 }
182 }
183 SPNORMRET2(xs, xe, xm, "sub", x, y);
184}
diff --git a/arch/mips/math-emu/sp_tint.c b/arch/mips/math-emu/sp_tint.c
new file mode 100644
index 000000000000..1d73d2abe0b5
--- /dev/null
+++ b/arch/mips/math-emu/sp_tint.c
@@ -0,0 +1,128 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include <linux/kernel.h>
29#include "ieee754sp.h"
30
31int ieee754sp_tint(ieee754sp x)
32{
33 COMPXSP;
34
35 CLEARCX;
36
37 EXPLODEXSP;
38 FLUSHXSP;
39
40 switch (xc) {
41 case IEEE754_CLASS_SNAN:
42 case IEEE754_CLASS_QNAN:
43 case IEEE754_CLASS_INF:
44 SETCX(IEEE754_INVALID_OPERATION);
45 return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x);
46 case IEEE754_CLASS_ZERO:
47 return 0;
48 case IEEE754_CLASS_DNORM:
49 case IEEE754_CLASS_NORM:
50 break;
51 }
52 if (xe >= 31) {
53 /* look for valid corner case */
54 if (xe == 31 && xs && xm == SP_HIDDEN_BIT)
55 return -0x80000000;
56 /* Set invalid. We will only use overflow for floating
57 point overflow */
58 SETCX(IEEE754_INVALID_OPERATION);
59 return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x);
60 }
61 /* oh gawd */
62 if (xe > SP_MBITS) {
63 xm <<= xe - SP_MBITS;
64 } else {
65 u32 residue;
66 int round;
67 int sticky;
68 int odd;
69
70 if (xe < -1) {
71 residue = xm;
72 round = 0;
73 sticky = residue != 0;
74 xm = 0;
75 }
76 else {
77 /* Shifting a u32 32 times does not work,
78 * so we do it in two steps. Be aware that xe
79 * may be -1 */
80 residue = xm << (xe + 1);
81 residue <<= 31 - SP_MBITS;
82 round = (residue >> 31) != 0;
83 sticky = (residue << 1) != 0;
84 xm >>= SP_MBITS - xe;
85 }
86 odd = (xm & 0x1) != 0x0;
87 switch (ieee754_csr.rm) {
88 case IEEE754_RN:
89 if (round && (sticky || odd))
90 xm++;
91 break;
92 case IEEE754_RZ:
93 break;
94 case IEEE754_RU: /* toward +Infinity */
95 if ((round || sticky) && !xs)
96 xm++;
97 break;
98 case IEEE754_RD: /* toward -Infinity */
99 if ((round || sticky) && xs)
100 xm++;
101 break;
102 }
103 if ((xm >> 31) != 0) {
104 /* This can happen after rounding */
105 SETCX(IEEE754_INVALID_OPERATION);
106 return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x);
107 }
108 if (round || sticky)
109 SETCX(IEEE754_INEXACT);
110 }
111 if (xs)
112 return -xm;
113 else
114 return xm;
115}
116
117
118unsigned int ieee754sp_tuns(ieee754sp x)
119{
120 ieee754sp hb = ieee754sp_1e31();
121
122 /* what if x < 0 ?? */
123 if (ieee754sp_lt(x, hb))
124 return (unsigned) ieee754sp_tint(x);
125
126 return (unsigned) ieee754sp_tint(ieee754sp_sub(x, hb)) |
127 ((unsigned) 1 << 31);
128}
diff --git a/arch/mips/math-emu/sp_tlong.c b/arch/mips/math-emu/sp_tlong.c
new file mode 100644
index 000000000000..4be21aa81fbf
--- /dev/null
+++ b/arch/mips/math-emu/sp_tlong.c
@@ -0,0 +1,123 @@
1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 * This program is free software; you can distribute it and/or modify it
12 * under the terms of the GNU General Public License (Version 2) as
13 * published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 * for more details.
19 *
20 * You should have received a copy of the GNU General Public License along
21 * with this program; if not, write to the Free Software Foundation, Inc.,
22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754sp.h"
29
30s64 ieee754sp_tlong(ieee754sp x)
31{
32 COMPXDP; /* <-- need 64-bit mantissa tmp */
33
34 CLEARCX;
35
36 EXPLODEXSP;
37 FLUSHXSP;
38
39 switch (xc) {
40 case IEEE754_CLASS_SNAN:
41 case IEEE754_CLASS_QNAN:
42 case IEEE754_CLASS_INF:
43 SETCX(IEEE754_INVALID_OPERATION);
44 return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x);
45 case IEEE754_CLASS_ZERO:
46 return 0;
47 case IEEE754_CLASS_DNORM:
48 case IEEE754_CLASS_NORM:
49 break;
50 }
51 if (xe >= 63) {
52 /* look for valid corner case */
53 if (xe == 63 && xs && xm == SP_HIDDEN_BIT)
54 return -0x8000000000000000LL;
55 /* Set invalid. We will only use overflow for floating
56 point overflow */
57 SETCX(IEEE754_INVALID_OPERATION);
58 return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x);
59 }
60 /* oh gawd */
61 if (xe > SP_MBITS) {
62 xm <<= xe - SP_MBITS;
63 } else if (xe < SP_MBITS) {
64 u32 residue;
65 int round;
66 int sticky;
67 int odd;
68
69 if (xe < -1) {
70 residue = xm;
71 round = 0;
72 sticky = residue != 0;
73 xm = 0;
74 }
75 else {
76 residue = xm << (32 - SP_MBITS + xe);
77 round = (residue >> 31) != 0;
78 sticky = (residue << 1) != 0;
79 xm >>= SP_MBITS - xe;
80 }
81 odd = (xm & 0x1) != 0x0;
82 switch (ieee754_csr.rm) {
83 case IEEE754_RN:
84 if (round && (sticky || odd))
85 xm++;
86 break;
87 case IEEE754_RZ:
88 break;
89 case IEEE754_RU: /* toward +Infinity */
90 if ((round || sticky) && !xs)
91 xm++;
92 break;
93 case IEEE754_RD: /* toward -Infinity */
94 if ((round || sticky) && xs)
95 xm++;
96 break;
97 }
98 if ((xm >> 63) != 0) {
99 /* This can happen after rounding */
100 SETCX(IEEE754_INVALID_OPERATION);
101 return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x);
102 }
103 if (round || sticky)
104 SETCX(IEEE754_INEXACT);
105 }
106 if (xs)
107 return -xm;
108 else
109 return xm;
110}
111
112
113u64 ieee754sp_tulong(ieee754sp x)
114{
115 ieee754sp hb = ieee754sp_1e63();
116
117 /* what if x < 0 ?? */
118 if (ieee754sp_lt(x, hb))
119 return (u64) ieee754sp_tlong(x);
120
121 return (u64) ieee754sp_tlong(ieee754sp_sub(x, hb)) |
122 (1ULL << 63);
123}
generated by cgit v1.2.2 (git 2.25.0) at 2024-11-30 16:21:23 -0500