1/* Copyright (c) 2002  Michael Stumpf  <mistumpf@de.pepperl-fuchs.com>
2   Copyright (c) 2006  Dmitry Xmelkov
3   Copyright (c) 2008  Ruud v Gessel
4   All rights reserved.
5
6   Redistribution and use in source and binary forms, with or without
7   modification, are permitted provided that the following conditions are met:
8
9   * Redistributions of source code must retain the above copyright
10     notice, this list of conditions and the following disclaimer.
11   * Redistributions in binary form must reproduce the above copyright
12     notice, this list of conditions and the following disclaimer in
13     the documentation and/or other materials provided with the
14     distribution.
15   * Neither the name of the copyright holders nor the names of
16     contributors may be used to endorse or promote products derived
17     from this software without specific prior written permission.
18
19   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
23   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29   POSSIBILITY OF SUCH DAMAGE. */
30
31/* $Id: fp_rempio2.S 2473 2015-04-09 08:10:22Z pitchumani $ */
32
33#include "fp32def.h"
34#include "asmdef.h"
35
36/* <non_standart> __fp_rempio2 (float x);
37     The __fp_rempio2() function computes the remainder of dividing
38     absolute value of x by Pi/2. The return value is x - n*Pi/2, where
39     n is the quotient of abs(x)/(Pi/2), rounded towards zero to an integer.
40   Output:
41	rA3.rA2.rA1.rA0.rAE	- flt40_t remainder
42	ZL			- low byte of n
43 */
44
45#define	HI40_PIO2	0x3FC90FDA	/* (flt40_t) Pi/2	*/
46#define	LO40_PIO2	0xA2
47
48FUNCTION __fp_rempio2
49
500:	XJMP	_U(__fp_nan)
51
52ENTRY __fp_rempio2
53  ; split and check finite
54	XCALL	_U(__fp_splitA)
55	brcs	0b		; only finite numbers are valid
56	clt			; ignore a sign
57  ; init division result
58	ldi	ZL, 0
59  ; extend A
60	clr	rAE
61  ; check (and modify) exponent
62	subi	rA3, hhi8(HI40_PIO2 << 1)
63	brlo	5f		; fabs(A) < 1.0 radian, C is set
64  ; prepare loop
65	ldi	rB0,  lo8(HI40_PIO2)
66	ldi	rB1,  hi8(HI40_PIO2)
67	ldi	rB2, hlo8(HI40_PIO2 | 0x800000)		; + hidden bit
68	rjmp	1f
69
70.Loop:	lsl	ZL
71	lsl	rAE
72	rol	rA0
73	rol	rA1
74	rol	rA2
75	brcs	2f
761:	cpi	rAE, lo8(LO40_PIO2)
77	cpc	rA0, rB0
78	cpc	rA1, rB1
79	cpc	rA2, rB2
80	brlo	3f
812:	subi	rAE, lo8(LO40_PIO2)
82	sbc	rA0, rB0
83	sbc	rA1, rB1
84	sbc	rA2, rB2
85	inc	ZL
863:	dec	rA3
87	brpl	.Loop
88
89  ; Normalize, we know that rA2.1.0.E >= 0x0E. You can check this with
90  ; a test program below.
91	cpi	rA2,0x80
92	brcc	5f
934:	dec	rA3
94	lsl	rAE
95	rol	rA0
96	rol	rA1
97	rol	rA2			; C := 0
98	brpl	4b
99
1005:	sbci	rA3, hhi8((HI40_PIO2<<1) + 0x01000000)	; undo the subi 0x7f
101	XJMP	_U(__fp_mpack_finite)
102ENDFUNC
103
104#if 0
105/* This is a test program to find the smallest value of rA2.1.0.E after
106   division. The nonzero value gives a garanty that normalization loop
107   is finite.	*/
108
109#include <stdio.h>
110#define	MNT32_PIO2	0xC90FDAA2
111
112int main ()
113{
114    unsigned long rA210;
115    unsigned long rA210E;
116    int rA3;
117    unsigned long c;
118    unsigned long amin = 0xffffffff;
119
120    for (rA210 = 0x800000; rA210 <= 0xffffff; rA210 += 1) {
121	rA210E = rA210 << 8;
122	c = 0;
123	rA3 = 127;	/* this is max for finite number	*/
124	goto m;
125	do {
126	    c = rA210E & 0x80000000;
127	    rA210E <<= 1;
128	  m:
129	    if (c || (rA210E >= MNT32_PIO2))
130		rA210E -= MNT32_PIO2;
131	    if (rA210E < amin) {
132		amin = rA210E;
133	        printf ("min of rA210E: 0x%08lx\r", amin);
134		fflush (stdout);
135	    }
136	} while (--rA3 >= 0);
137    }
138    putchar ('\n');
139    return 0;
140}
141#endif
142