3f9637545e4a711d5356abc2a80e6e8d424ba963
[mes.git] / lib / mes / div.c
1 /* -*-comment-start: "//";comment-end:""-*-
2  * GNU Mes --- Maxwell Equations of Software
3  * Copyright © 2016,2017,2018,2019 Jan (janneke) Nieuwenhuizen <janneke@gnu.org>
4  * Copyright © 2019 Danny Milosavljevic <dannym@scratchpost.org>
5  *
6  * This file is part of GNU Mes.
7  *
8  * GNU Mes is free software; you can redistribute it and/or modify it
9  * under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 3 of the License, or (at
11  * your option) any later version.
12  *
13  * GNU Mes is distributed in the hope that it will be useful, but
14  * WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with GNU Mes.  If not, see <http://www.gnu.org/licenses/>.
20  */
21
22 #include <mes/lib.h>
23 #include <assert.h>
24 #include <stdint.h>
25 #include <limits.h>
26
27 /*#define LONG_MIN (-(~0UL>>1)-1)*/
28
29 typedef struct
30 {
31   long quot;
32   long rem;
33 } ldiv_t;
34
35 void __mesabi_div0(void)
36 {
37   eputs(" ***MES C LIB*** divide by zero numerator=");
38   eputs("\n");
39   assert(0);
40 }
41
42 /* Compare gcc: __udivdi3 */
43 unsigned long __mesabi_uldiv(unsigned long a, unsigned long b, unsigned long* remainder)
44 {
45   unsigned long tmp;
46   if (!remainder)
47     remainder = &tmp;
48   *remainder = 0;
49   if (b == 1)
50     return a;
51   else if (b == 0)
52     __mesabi_div0();
53   else
54     {
55       unsigned long x;
56       for (x = 0; a >= b; a -= b)
57         ++x;
58       *remainder = a;
59       return x;
60     }
61 }
62
63 /* Note: Rounds towards zero */
64 ldiv_t ldiv(long a, long b)
65 {
66   ldiv_t result;
67   int negate_result = (a < 0) ^ (b < 0);
68   assert(b != LONG_MIN);
69   if (a != LONG_MIN)
70     {
71       if (a < 0)
72         a = -a;
73       if (b < 0)
74         b = -b;
75       result.quot = __mesabi_uldiv(a, b, &result.rem);
76       if (negate_result)
77         result.quot = -result.quot;
78       return result;
79     }
80   else
81     {
82       result.rem = 0;
83       if (b < 0)
84         b = -b;
85       if (b == 1)
86         {
87           result.quot = a;
88           /* Since result.quot is already negative, don't negate it again. */
89           negate_result = !negate_result;
90         }
91       else if (b == 0)
92         __mesabi_div0();
93       else
94         {
95           long x;
96           for (x = 0; a <= -b; a += b)
97             ++x;
98           result.rem = -a;
99           result.quot = x;
100         }
101       if (negate_result)
102         result.quot = -result.quot;
103       return result;
104     }
105 }