1 // arith.h Copyright (C) Codemist, 1990-2021
2
3
4 /**************************************************************************
5 * Copyright (C) 2021, Codemist. A C Norman *
6 * *
7 * Redistribution and use in source and binary forms, with or without *
8 * modification, are permitted provided that the following conditions are *
9 * met: *
10 * *
11 * * Redistributions of source code must retain the relevant *
12 * copyright notice, this list of conditions and the following *
13 * disclaimer. *
14 * * Redistributions in binary form must reproduce the above *
15 * copyright notice, this list of conditions and the following *
16 * disclaimer in the documentation and/or other materials provided *
17 * with the distribution. *
18 * *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
20 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT *
21 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS *
22 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE *
23 * COPYRIGHT OWNERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, *
24 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, *
25 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS *
26 * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
27 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR *
28 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF *
29 * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH *
30 * DAMAGE. *
31 *************************************************************************/
32
33 // $Id: arith.h 5609 2021-01-23 22:02:30Z arthurcnorman $
34
35 #ifndef header_arith_h
36 #define header_arith_h 1
37
38 #define TWO_32 4294967296.0 // 2^32
39 #define TWO_31 2147483648.0 // 2^31
40 #define TWO_24 16777216.0 // 2^24
41 #define TWO_22 4194304.0 // 2^22
42 #define TWO_21 2097152.0 // 2^21
43 #define TWO_20 1048576.0 // 2^20
44
45 //
46 // I am going to need to rely on my C compiler turning these strings
47 // into well-rounded versions of the relevant numbers. To increase the
48 // changes of that I will write the constants as
49 // ((A + B)/C)
50 // where A will be an integer and hence representable exactly, B will be
51 // be much smaller and C will be a power of two, so dividing by it should
52 // not introduce any additional error at all.
53 //
54 // #define _pi 3.14159265358979323846
55 // #define _half_pi 1.57079632679489661923
56
57 #define _pi ((12868.0 - 0.036490896206895257)/4096.0)
58 #define _half_pi ((12868.0 - 0.036490896206895257)/8192.0)
59
60 #define boole_clr 0
61 #define boole_and 1
62 #define boole_andc2 2
63 #define boole_1 3
64 #define boole_andc1 4
65 #define boole_2 5
66 #define boole_xor 6
67 #define boole_ior 7
68 #define boole_nor 8
69 #define boole_eqv 9
70 #define boole_c2 10
71 #define boole_orc2 11
72 #define boole_c1 12
73 #define boole_orc1 13
74 #define boole_nand 14
75 #define boole_set 15
76
77 extern unsigned char msd_table[256], lsd_table[256];
78
79 //
80 // Bignums are represented as vectors of digits, where each digit
81 // uses 31 bits, and all but the most significant digit are unsigned
82 // (and thus do not use the 0x80000000L bit). The most significant
83 // digit of a bignum is a signed 2-s complement value in 31 bits that
84 // has been sign extended into the 0x80000000L bit, and thus its top
85 // two bits (in the 32 bit word) will be either '00' or '11'.
86 // NOTE that even on a 64-bit machine I will work with 32-bit values
87 // as digits in bignums.
88 //
89
90 // NOTE please that these are all using 32-bit arthmetic. They are only
91 // intended for use on values that are digits making up parts of bignums.
92 // Note that some C++ compilers try to be VERY clever about assuming that
93 // integer arithmetic will never overflow, so I need to do things in
94 // unsigned mode if the top bit is liable to arise as a carry-bit.
95
96 #define top_bit_set(n) (((int32_t)(n)) < 0)
97 #define top_bit(n) ((int32_t)(((uint32_t)(n)) >> 31))
98 #define set_top_bit(n) ((int32_t)((uint32_t)(n) | (uint32_t)0x80000000U))
99 #define clear_top_bit(n) ((int32_t)((uint32_t)(n) & 0x7fffffff))
100
101 // As with fixnum_of_int I need to take care here that the arithmetic I do
102 // here can not overflow, since if it did the behaviour would be undefined
103 // and a modern optimising compiler would be entitled to do pretty much
104 // whatever it felt like. I hope that good compilers will observe that the
105 // mask operation as written here is not really needed on most machines.
106
107 #define signed_overflow(n) \
108 top_bit_set(static_cast<uint32_t>(n) ^ (static_cast<uint32_t>(n) << 1))
109
110 // ADD32 forces and addition to be done as unsigned arithmetic, and may be
111 // useful when this avoids risk of a carry into the most significant bit
112 // being interpreted as overflow and hence undefined behaviour.
113
114 #define ADD32(a, b) ((uint32_t)(a) + (uint32_t)(b))
115
116
117 // The following all take an intptr_t or int64_t values (as per their name)
118 // and check if their argument would fit in a 29 or 31-bit signed value. I
119 // try fairly hard to avoid overflow and keep the code here so that the
120 // only way in which it goes beyond what the C/C++ standards guarantee is
121 // an assumption that numbers are 2s complement and that casts between
122 // signed and unsigned values leave the bitwise representations unchanged.
123
124 #define signed29_in_64(n) \
125 ((static_cast<int64_t>((static_cast<uint64_t>(n) & 0x1fffffffU) << 35) / (static_cast<int64_t>(1) << 35)) == \
126 static_cast<int64_t>(n))
127
128 #define signed31_in_64(n) \
129 ((static_cast<int64_t>((static_cast<uint64_t>(n) & 0x7fffffffU) << 33) / (static_cast<int64_t>(1) << 33)) == \
130 static_cast<int64_t>(n))
131
132 #define signed31_in_ptr(n) \
133 ((static_cast<intptr_t>((static_cast<uintptr_t>(n)&0x7fffffffU) << (8*sizeof(intptr_t) - 31)) / \
134 (static_cast<intptr_t>(1) << (8*sizeof(intptr_t) - 31))) == static_cast<intptr_t>(n))
135
136 #ifdef HAVE_SOFTFLOAT
137 // I use rounding-direction specifiers from SoftFloat because I can then
138 // pass them down to the float128_t functions when I need to
139 #define FIX_TRUNCATE softfloat_round_minMag
140 #define FIX_ROUND softfloat_round_near_even
141 #define FIX_FLOOR softfloat_round_min
142 #define FIX_CEILING softfloat_round_max
143 #else
144 // If I do not have SoftFloat I will use values that in fact match those I
145 // would be using if I did!
146 #define FIX_TRUNCATE 1
147 #define FIX_ROUND 0
148 #define FIX_FLOOR 2
149 #define FIX_CEILING 3
150 #endif // HAVE_SOFTFLOAT
151
152 extern LispObject lisp_fix(LispObject a, int roundmode);
153 extern LispObject lisp_ifix(LispObject a, LispObject b,
154 int roundmode);
155
156
157 // The following tests for IEEE infinities and NaNs.
158
floating_edge_case(double r)159 inline bool floating_edge_case(double r)
160 { return !std::isfinite(r);
161 }
162
floating_clear_flags()163 inline void floating_clear_flags()
164 {}
165
166 // An alternative possible implementation could go
167 // #include <fenv.h>
168 // {
169 // #pragma STDC FENV_ACCESS ON
170 // <arithmetic>
171 // if (fetestexcept(FE_OVERFLOW | FE_INVALID | FE_DIVBYZERO))
172 // { feclearexcept(FE_ALL_EXCEPT);
173 // ..
174 // }
175 // }
176 // The STDC_FENV_ACCESS pragma came in with C99 and C++11, but at the time
177 // of writing gcc and g++ do not support it on quite enought platforms.
178 // When and if they do I may move to exploit it. Meanwhile I will use the
179 // version that is probably more portable.
180
181 //
182 // Here I do some arithmetic in-line. In the following macros I need to
183 // take care that the names used for local variables do not clash with
184 // those used in the body of the code. Hence the names r64 and c64, which
185 // I must agree not to use elsewhere. Note also the "do {} while (0)" idiom
186 // to avoid nasty problems with C syntax and the need for semicolons.
187 // I should make a transition to use of inline functions!
188
189 #define Dmultiply(hi, lo, a, b, c) \
190 do { uint64_t r64 = static_cast<uint64_t>(a) * \
191 static_cast<uint64_t>(b) + \
192 static_cast<uint32_t>(c); \
193 (lo) = 0x7fffffffu & static_cast<uint32_t>(r64); \
194 (hi) = static_cast<uint32_t>(r64 >> 31); } while (0)
195
196 #define Ddivide(r, q, a, b, c) \
197 do { uint64_t r64 = (static_cast<uint64_t>(a) << 31) | \
198 static_cast<uint64_t>(b); \
199 uint64_t c64 = static_cast<uint64_t>( \
200 static_cast<uint32_t>(c)); \
201 q = static_cast<uint32_t>(r64 / c64); \
202 r = static_cast<uint32_t>(r64 % c64); } while (0)
203
204 #define Ddiv10_9(r, q, a, b) Ddivide(r, q, a, b, 1000000000u)
205
206 #define Ddivideq(q, a, b, c) \
207 do { uint64_t r64 = ((static_cast<uint64_t>(a)) << 31) | static_cast<uint64_t>(b); \
208 uint64_t c64 = static_cast<uint64_t>(static_cast<uint32_t>(c)); \
209 q = static_cast<uint32_t>(r64 / c64); } while (0)
210
211 #define Ddiv10_9q(r, q, a, b) Ddivideq(q, a, b, 1000000000u)
212
213 #define Ddivider(r, a, b, c) \
214 do { uint64_t r64 = ((static_cast<uint64_t>(a)) << 31) | static_cast<uint64_t>(b); \
215 uint64_t c64 = static_cast<uint64_t>(static_cast<uint32_t>(c)); \
216 r = static_cast<uint32_t>(r64 % c64); } while (0)
217
218 #define Ddiv10_9r(r, q, a, b) Ddivider(r, a, b, 1000000000u)
219
220 #define fixnum_minusp(a) (static_cast<intptr_t>(a) < 0)
221
222 #define bignum_minusp(a) \
223 (static_cast<int32_t>(bignum_digits(a)[((bignum_length(a)-CELL)/4)-1])<0)
224
value_of_immediate_float(LispObject a)225 inline double value_of_immediate_float(LispObject a)
226 { Float_union aa;
227 // Worry about strict aliasing here, at least maybe. With GCC I believe I am
228 // safe, but as per the standards I think I am not.
229 if (SIXTY_FOUR_BIT) aa.i = static_cast<int32_t>(static_cast<uint64_t>
230 (a)>>32);
231 else aa.i = static_cast<int32_t>(a - XTAG_SFLOAT);
232 return aa.f;
233 }
234
235 extern LispObject make_boxfloat(double a, int type=TYPE_DOUBLE_FLOAT);
236 #ifdef HAVE_SOFTFLOAT
237 extern LispObject make_boxfloat128(float128_t a);
238 #endif // HAVE_SOFTFLOAT
239
240 // Pack something as a short (28-bit) float
241
pack_short_float(double d)242 inline LispObject pack_short_float(double d)
243 { Float_union aa;
244 aa.f = d;
245 if (trap_floating_overflow &&
246 floating_edge_case(aa.f))
247 { floating_clear_flags();
248 return aerror("exception with short float");
249 }
250 aa.i &= ~0xf;
251 if (SIXTY_FOUR_BIT)
252 return static_cast<LispObject>((static_cast<uint64_t>
253 (aa.i)) << 32) + XTAG_SFLOAT;
254 else return aa.i + XTAG_SFLOAT;
255 }
256
257 // Create a single (32-bit) float. Just like make_single_float but inlined
258 // in the 64-bit case
259
pack_single_float(double d)260 inline LispObject pack_single_float(double d)
261 { if (SIXTY_FOUR_BIT)
262 { Float_union aa;
263 aa.f = d;
264 if (trap_floating_overflow &&
265 floating_edge_case(aa.f))
266 { floating_clear_flags();
267 return aerror("exception with single float");
268 }
269 return static_cast<LispObject>(static_cast<uint64_t>
270 (aa.i) << 32) + XTAG_SFLOAT + XTAG_FLOAT32;
271 }
272 else
273 { LispObject r = get_basic_vector(TAG_BOXFLOAT,
274 TYPE_SINGLE_FLOAT, sizeof(Single_Float));
275 single_float_val(r) = static_cast<float>(d);
276 if (trap_floating_overflow &&
277 floating_edge_case(single_float_val(r)))
278 { floating_clear_flags();
279 return aerror("exception with single float");
280 }
281 return r;
282 }
283 }
284
285 // Pack either a 28 or 32-bit float with type Lisp value "l1" indicating
286 // whether 28 or 32 bits are relevant. On a 32-bit machine like will always
287 // be a 28-bit value. If a second argument l2 is provided the width of the
288 // result will match the wider of the two.
289
290 inline LispObject pack_immediate_float(double d,
291 LispObject l1, LispObject l2=0)
292 { Float_union aa;
293 aa.f = d;
294 if (trap_floating_overflow &&
295 floating_edge_case(aa.f))
296 { floating_clear_flags();
297 if (((l1 | l2) & XTAG_FLOAT32) != 0)
298 return aerror("exception with single float");
299 else return aerror("exception with short float");
300 }
301 if (SIXTY_FOUR_BIT)
302 { if (((l1 | l2) & XTAG_FLOAT32) == 0) aa.i &= ~0xf;
303 return static_cast<LispObject>((static_cast<uint64_t>
304 (aa.i)) << 32) + XTAG_SFLOAT +
305 ((l1 | l2) & XTAG_FLOAT32);
306 }
307 aa.i &= ~0xf;
308 return aa.i + XTAG_SFLOAT;
309 }
310
311 // comparing 64-bit integers against (double precision) is perhaps
312 // unexpectedly delicate. Here is some code to help. You can find two
313 // sources of extra commentary about this. One is by Andrew Koenig in
314 // a Dr Dobbs article in 2013, the other is in (MIT Licensed) Julia and
315 // a discussion at https://github.com/JuliaLang/julia/issues/257.
316
eq_i64d(int64_t a,double b)317 inline bool eq_i64d(int64_t a, double b)
318 {
319 // The integer can always be converted to a double, but of course
320 // sometimes there will be rounding involved. But if the value does not
321 // match the double even after rounding then the two values are certainly
322 // different. Also if the double happens to be a NaN this will lead to
323 // a returned value of false (as required).
324 if (b != static_cast<double>(a)) return false;
325 // Now the two values differ by at most the rounding that happened when
326 // the integer was converted to a double. This ALMOST means that the double
327 // has a value that fits in the range of integers. However if a has a value
328 // just less than 2^63 and b is (double)(2^63) then b can not be cast to
329 // an integer safely. In C++ the consequence of trying to cast a double to
330 // and int where the result would not fit is undefined, and so could
331 // include arbitrary bad behaviour. So I have to filter that case out.
332 if (b == static_cast<double>(static_cast<uint64_t>
333 (1)<<63)) return false;
334 // With the special case out of the way I can afford to case from double to
335 // int64_t. The negative end of the range is safe!
336 return a == static_cast<int64_t>(b);
337 }
338
lessp_i64d(int64_t a,double b)339 inline bool lessp_i64d(int64_t a, double b)
340 {
341 // If the integer is <= 2^53 then converting it to a double does not
342 // introduce any error at all, so I can perform the comparison reliably
343 // on doubles. If d ia a NaN this is still OK.
344 if (a <= (static_cast<int64_t>(1)<<53) &&
345 a >= -(static_cast<int64_t>(1)<<53)) return static_cast<double>
346 (a) < b;
347 // If the float is outside the range of int64_t I can tell how the
348 // comparison must play out. Note that near the value 2^63 the next
349 // double value lower than 2^63 is in integer, as we can not have any
350 // floating point value larger than the largest positive int64_t value
351 // and less then 2^63. I make these tests of the form "if (!xxx)" because
352 // then if b is a NaN the comparison returns false and I end up exiting.
353 if (!(b >= -static_cast<double>(static_cast<uint64_t>
354 (1)<<63))) return false;
355 if (!(b < static_cast<double>(static_cast<uint64_t>
356 (1)<<63))) return true;
357 // Now we know that a is large and b is not huge. I will just discuss the
358 // case of two positive numbers here, but mixed signs and negative values
359 // follow the same.
360 // I am going to convert b to an integer and then compare. Because I have
361 // ensures that b is not too big (including knowing it is not an infinity
362 // or a NaN) I will not get overflow or failure in that conversion. So the
363 // only concern is the effect of rounding in the conversion.
364 // Well if b >= 2^52 it has an exact integer as its value so the conversion
365 // will be exact and the comparison reliable.
366 // if b < 2^52 but a > 2^53 then rounding of b that leaves a fractional part
367 // less than 1 does not matter and again the comparison is reliable.
368 return a < static_cast<int64_t>(b);
369 }
370
lessp_di64(double a,int64_t b)371 inline bool lessp_di64(double a, int64_t b)
372 {
373 // The logic here is much as above - by omitting all the commentary
374 // you can see much more clearly just how long the code is.
375 if (b <= (static_cast<int64_t>(1)<<53) &&
376 b >= -(static_cast<int64_t>(1)<<53)) return a < static_cast<double>
377 (b);
378 if (!(a < static_cast<double>(static_cast<uint64_t>
379 (1)<<63))) return false;
380 if (!(a >= -static_cast<double>(static_cast<uint64_t>
381 (1)<<63))) return true;
382 return static_cast<int64_t>(a) < b;
383 }
384
385 extern LispObject negateb(LispObject);
386 extern LispObject copyb(LispObject);
387 extern LispObject negate(LispObject);
388 extern LispObject plus2(LispObject a, LispObject b);
389 extern LispObject difference2(LispObject a, LispObject b);
390 extern LispObject times2(LispObject a, LispObject b);
391 extern LispObject quot2(LispObject a, LispObject b);
392 extern LispObject CLquot2(LispObject a, LispObject b);
393 extern LispObject quotbn(LispObject a, int32_t n);
394 extern LispObject quotbn1(LispObject a, int32_t n);
395 #define QUOTBB_QUOTIENT_NEEDED 1
396 #define QUOTBB_REMAINDER_NEEDED 2
397 extern LispObject quotbb(LispObject a, LispObject b, int needs);
398 extern LispObject Cremainder(LispObject a, LispObject b);
399 extern LispObject rembi(LispObject a, LispObject b);
400 extern LispObject rembb(LispObject a, LispObject b);
401 extern LispObject shrink_bignum(LispObject a, size_t lena);
402 extern LispObject modulus(LispObject a, LispObject b);
403 extern LispObject rational(LispObject a);
404 extern LispObject rationalize(LispObject a);
405 extern LispObject lcm(LispObject a, LispObject b);
406 extern LispObject lengthen_by_one_bit(LispObject a, int32_t msd);
407 extern bool numeq2(LispObject a, LispObject b);
408 extern bool SL_numeq2(LispObject a, LispObject b);
409 extern bool zerop(LispObject a);
410 extern bool onep(LispObject a);
411 extern bool minusp(LispObject a);
412 extern bool plusp(LispObject a);
413
414 extern LispObject integer_decode_long_float(LispObject a);
415 extern LispObject Linteger_decode_float(LispObject env, LispObject a);
416
417 extern LispObject validate_number(const char *s, LispObject a,
418 LispObject b, LispObject c);
419
420 extern LispObject make_fake_bignum(intptr_t n);
421 extern LispObject make_one_word_bignum(int32_t n);
422 extern LispObject make_two_word_bignum(int32_t a, uint32_t b);
423 extern LispObject make_three_word_bignum(int32_t a, uint32_t b,
424 uint32_t c);
425 extern LispObject make_four_word_bignum(int32_t a, uint32_t b,
426 uint32_t c, uint32_t d);
427 extern LispObject make_five_word_bignum(int32_t a, uint32_t b,
428 uint32_t c, uint32_t d, uint32_t e);
429 extern LispObject make_n_word_bignum(int32_t a2, uint32_t a1,
430 uint32_t a0, size_t n);
431 extern LispObject make_n4_word_bignum(int32_t a3, uint32_t a2,
432 uint32_t a1, uint32_t a0, size_t n);
433 extern LispObject make_n5_word_bignum(int32_t a4, uint32_t a3,
434 uint32_t a2, uint32_t a1,
435 uint32_t a0, size_t n);
436
437 extern LispObject make_power_of_two(size_t n);
438
439 extern LispObject make_lisp_integer32_fn(int32_t n);
make_lisp_integer32(int32_t n)440 inline LispObject make_lisp_integer32(int32_t n)
441 { if (SIXTY_FOUR_BIT ||
442 valid_as_fixnum(n)) return fixnum_of_int(static_cast<intptr_t>(n));
443 else return make_lisp_integer32_fn(n);
444 }
445
446 extern LispObject make_lisp_integer64_fn(int64_t n);
make_lisp_integer64(int64_t n)447 inline LispObject make_lisp_integer64(int64_t n)
448 { if (valid_as_fixnum(n)) return fixnum_of_int(
449 static_cast<intptr_t>(n));
450 else return make_lisp_integer64_fn(n);
451 }
452
453 extern LispObject make_lisp_unsigned64_fn(uint64_t n);
make_lisp_unsigned64(uint64_t n)454 inline LispObject make_lisp_unsigned64(uint64_t n)
455 { if (n < (static_cast<uint64_t>(1))<<(8*sizeof(intptr_t)-5))
456 return fixnum_of_int(static_cast<intptr_t>(n));
457 else return make_lisp_unsigned64_fn(n);
458 }
459
460 extern LispObject make_lisp_integerptr_fn(intptr_t n);
461
462 // There is a little C++ joke here. I had originally used valid_as_fixnum()
463 // to test if the intptr_t value was suitably in range. That function provides
464 // overloads for int32_t and int64_t. And for int128_t if I have a 128-bit
465 // integer type available. It does not provide a direct overload for intptr_t
466 // because that is liable to be one of the previous types and so would clash
467 // and be rejected. Anyway on clang (on a Macintosh) is deemed ambiguous and
468 // the code failed to compile. clang would allow me to add an extra overload
469 // to valid_as_integer so there it is treating intptr_t as a new integer type
470 // not exactly the same as int32_t or int64_t, rather as if it only has
471 // preprocessor levels of understanding so it does not yet know which it is.
472 // All this means that I need to go and read the C++ standard carefully, but
473 // in the meanwhile having conservative code feels safer... so this comment is
474 // the main cost.
475
make_lisp_integerptr(intptr_t n)476 inline LispObject make_lisp_integerptr(intptr_t n)
477 { if (intptr_valid_as_fixnum(n)) return fixnum_of_int(n);
478 else return make_lisp_integerptr_fn(n);
479 }
480
481 extern LispObject make_lisp_unsignedptr_fn(uintptr_t n);
make_lisp_unsignedptr(uintptr_t n)482 inline LispObject make_lisp_unsignedptr(uintptr_t n)
483 { if (n < (static_cast<uintptr_t>(1))<<(8*sizeof(intptr_t)-5))
484 return fixnum_of_int(static_cast<intptr_t>(n));
485 else return make_lisp_unsignedptr_fn(n);
486 }
487
488 extern LispObject make_lisp_integer128_fn(int128_t n);
489
make_lisp_integer128(int128_t n)490 inline LispObject make_lisp_integer128(int128_t n)
491 { if (valid_as_fixnum(n)) return fixnum_of_int(static_cast<int64_t>(n));
492 else return make_lisp_integer128_fn(n);
493 }
494
495 extern LispObject make_lisp_unsigned128_fn(uint128_t n);
496
make_lisp_unsigned128(uint128_t n)497 inline LispObject make_lisp_unsigned128(uint128_t n)
498 { if (uint128_valid_as_fixnum(n))
499 return fixnum_of_int(static_cast<uint64_t>(static_cast<uint64_t>(n)));
500 else return make_lisp_unsigned128_fn(n);
501 }
502
validate_number(LispObject n)503 inline void validate_number(LispObject n)
504 {
505 // Can be used while debugging!
506 }
507
508 extern double float_of_integer(LispObject a);
509 extern LispObject add1(LispObject p);
510 extern LispObject sub1(LispObject p);
511 extern LispObject integerp(LispObject p);
512 extern double float_of_number(LispObject a);
513 #ifdef HAVE_SOFTFLOAT
514 extern float128_t float128_of_number(LispObject a);
515 #endif // HAVE_SOFTFLOAT
516 extern LispObject make_complex(LispObject r, LispObject i);
517 extern LispObject make_ratio(LispObject p, LispObject q);
518 extern LispObject make_approx_ratio(LispObject p, LispObject q,
519 int bits);
520 extern LispObject ash(LispObject a, LispObject b);
521 extern LispObject lognot(LispObject a);
522 extern LispObject logior2(LispObject a, LispObject b);
523 extern LispObject logxor2(LispObject a, LispObject b);
524 extern LispObject logand2(LispObject a, LispObject b);
525 extern LispObject logeqv2(LispObject a, LispObject b);
526 extern LispObject rationalf(double d);
527 #ifdef HAVE_SOFTFLOAT
528 extern LispObject rationalf128(float128_t *d);
529 #endif // HAVE_SOFTFLOAT
530
531 extern int _reduced_exp(double, double *);
532 extern bool lesspbi(LispObject a, LispObject b);
533 extern bool lesspib(LispObject a, LispObject b);
534
535 //
536 // This is going to be a bit of a mess because I will want to use the C
537 // data-type "complex double" when that is available... what happens
538 // here will survive even without that.
539 //
540
541 typedef struct Complex
542 { double real;
543 double imag;
544 } Complex;
545
546 extern Complex Cln(Complex a);
547 extern Complex Ccos(Complex a);
548 extern Complex Cexp(Complex a);
549 extern Complex Cpow(Complex a, Complex b);
550 extern double Cabs(Complex a);
551
552 // For the parallel variant on Karatsuba I need thread support...
553
554 #ifndef HAVE_CILK
555
556 extern std::thread kara_thread[2];
557 #define KARA_0 (1<<0)
558 #define KARA_1 (1<<1)
559 #define KARA_QUIT (1<<2)
560 extern void kara_worker(int id);
561 extern std::mutex kara_mutex;
562 extern std::condition_variable cv_kara_ready,
563 cv_kara_done;
564 extern unsigned int kara_ready;
565 extern int kara_done;
566
567 #endif
568
569 extern size_t kparallel, karatsuba_parallel;
570
571 //
572 // Tests on my Intel i7 4770K system running Windows 7 (64-bit) I find that
573 // on numbers over 10^1000 (or so) the use of a multi-threaded Karatsuba
574 // starts to pay off. By 2000 decimal digits it is quite useful.
575 // The threshold set here of 120 31-bit words represents about the
576 // break-even point.
577 //
578 #ifndef KARATSUBA_PARALLEL_CUTOFF
579 # define KARATSUBA_PARALLEL_CUTOFF 120
580 #endif
581
582 #ifndef KARATSUBA_CUTOFF
583 //
584 // I have conducted some experiments on one machine to find out what the
585 // best cut-off value here is. The exact value chosen is not very
586 // critical, and the fancy techniques do not pay off until numbers get
587 // a lot bigger than this length (which is expressed in units of 31-bit
588 // words, i.e. about 10 decimals). Anyone who wants may recompile with
589 // alternative values to try to get the system fine-tuned for their
590 // own computer - but I do not expect it to be possible to achieve much
591 // by so doing.
592 //
593 #define KARATSUBA_CUTOFF 12
594
595 #endif
596
597
598 #ifdef HAVE_SOFTFLOAT
599 // Now some stuff relating to the float128_t type
600
601 static int f128M_exponent(const float128_t *p);
602 static void f128M_set_exponent(float128_t *p, int n);
603 extern void f128M_ldexp(float128_t *p, int n);
604 extern void f128M_frexp(float128_t *p, float128_t *r, int *x);
605 static bool f128M_infinite(const float128_t *p);
606 static bool f128M_finite(const float128_t *p);
607 static bool f128M_nan(const float128_t *x);
608 static bool f128M_subnorm(const float128_t *x);
609 static bool f128M_negative(const float128_t *x);
610 static void f128M_negate(float128_t *x);
611 extern void f128M_split(
612 const float128_t *x, float128_t *yhi, float128_t *ylo);
613
614 #ifdef LITTLEENDIAN
615 #define HIPART 1
616 #define LOPART 0
617 #else
618 #define HIPART 0
619 #define LOPART 1
620 #endif
621
f128M_zero(const float128_t * p)622 inline bool f128M_zero(const float128_t *p)
623 { return ((p->v[HIPART] & INT64_C(0x7fffffffffffffff)) == 0) &&
624 (p->v[LOPART] == 0);
625 }
626
f128M_infinite(const float128_t * p)627 inline bool f128M_infinite(const float128_t *p)
628 { return (((p->v[HIPART] >> 48) & 0x7fff) == 0x7fff) &&
629 ((p->v[HIPART] & INT64_C(0xffffffffffff)) == 0) &&
630 (p->v[LOPART] == 0);
631 }
632
f128M_finite(const float128_t * p)633 inline bool f128M_finite(const float128_t *p)
634 { return (((p->v[HIPART] >> 48) & 0x7fff) != 0x7fff);
635 }
636
f128M_make_infinite(float128_t * p)637 inline void f128M_make_infinite(float128_t *p)
638 { p->v[HIPART] |= UINT64_C(0x7fff000000000000);
639 p->v[HIPART] &= UINT64_C(0xffff000000000000);
640 p->v[LOPART] = 0;
641 }
642
f128M_make_zero(float128_t * p)643 inline void f128M_make_zero(float128_t *p)
644 { p->v[HIPART] &= UINT64_C(0x8000000000000000);
645 p->v[LOPART] = 0;
646 }
647
648 // Here I do not count 0.0 (or -0.0) as sub-normal.
649
f128M_subnorm(const float128_t * p)650 inline bool f128M_subnorm(const float128_t *p)
651 { return (((p->v[HIPART] >> 48) & 0x7fff) == 0) &&
652 (((p->v[HIPART] & INT64_C(0xffffffffffff)) != 0) ||
653 (p->v[LOPART] != 0));
654 }
655
f128M_nan(const float128_t * p)656 inline bool f128M_nan(const float128_t *p)
657 { return (((p->v[HIPART] >> 48) & 0x7fff) == 0x7fff) &&
658 (((p->v[HIPART] & INT64_C(0xffffffffffff)) != 0) ||
659 (p->v[LOPART] != 0));
660 }
661
f128M_negative(const float128_t * x)662 inline bool f128M_negative(const float128_t *x)
663 { if (f128M_nan(x)) return false;
664 return (static_cast<int64_t>(x->v[HIPART])) < 0;
665 }
666
f128M_exponent(const float128_t * p)667 inline int f128M_exponent(const float128_t *p)
668 { return ((p->v[HIPART] >> 48) & 0x7fff) - 0x3fff;
669 }
670
f128M_set_exponent(float128_t * p,int n)671 inline void f128M_set_exponent(float128_t *p, int n)
672 { p->v[HIPART] = (p->v[HIPART] & INT64_C(0x8000ffffffffffff)) |
673 ((static_cast<uint64_t>(n) + 0x3fff) << 48);
674 }
675
f128M_negate(float128_t * x)676 inline void f128M_negate(float128_t *x)
677 { x->v[HIPART] ^= UINT64_C(0x8000000000000000);
678 }
679
floating_edge_case128(float128_t * r)680 inline bool floating_edge_case128(float128_t *r)
681 { return f128M_infinite(r) || f128M_nan(r);
682 }
683
684 #endif // HAVE_SOFTFLOAT
685
686 extern int double_to_binary(double d, int64_t &m);
687 #ifdef HAVE_SOFTFLOAT
688 extern int float128_to_binary(const float128_t *d,
689 int64_t &mhi, uint64_t &mlo);
690 #endif // HAVE_SOFTFLOAT
691 extern intptr_t double_to_3_digits(double d,
692 int32_t &a2, uint32_t &a1, uint32_t &a0);
693 #ifdef HAVE_SOFTFLOAT
694 extern intptr_t float128_to_5_digits(float128_t *d,
695 int32_t &a4, uint32_t &a3, uint32_t &a2, uint32_t &a1, uint32_t &a0);
696
697 extern float128_t f128_0, // 0.0L0
698 f128_half, // 0.5L0
699 f128_mhalf, // -0.5L0
700 f128_1, // 1.0L0
701 f128_10, // 10.0L0
702 f128_10_17, // 1.0L17
703 f128_10_18, // 1.0L18
704 f128_r10, // 0.1L0
705 f128_N1; // 2.0L0^4096
706
707 typedef struct _float256_t
708 {
709 #ifdef LITTLEENDIAN
710 float128_t lo;
711 float128_t hi;
712 #else
713 float128_t hi;
714 float128_t lo;
715 #endif
716 } float256_t;
717
f128M_to_f256M(const float128_t * a,float256_t * b)718 inline void f128M_to_f256M(const float128_t *a, float256_t *b)
719 { b->hi = *a;
720 b->lo = f128_0;
721 }
722
723 extern void f256M_add(
724 const float256_t *x, const float256_t *y, float256_t *z);
725
726 extern void f256M_mul(
727 const float256_t *x, const float256_t *y, float256_t *z);
728
729 extern void f256M_pow(const float256_t *x, unsigned int y,
730 float256_t *z);
731
732 extern float256_t f256_1, f256_10, f256_r10, f256_5, f256_r5;
733
734 // These print 128-bit floats in the various standard styles, returning the
735 // number of characters used. The "sprint" versions put their result in
736 // a buffer, while "print" goes to stdout.
737
738 extern int f128M_sprint_E(char *s, int width, int precision,
739 float128_t *p);
740 extern int f128M_sprint_F(char *s, int width, int precision,
741 float128_t *p);
742 extern int f128M_sprint_G(char *s, int width, int precision,
743 float128_t *p);
744 extern int f128M_print_E(int width, int precision, float128_t *p);
745 extern int f128M_print_F(int width, int precision, float128_t *p);
746 extern int f128M_print_G(int width, int precision, float128_t *p);
747
748 extern float128_t atof128(const char *s);
749
750 #endif // HAVE_SOFTFLOAT
751
752 // Now let me provide a macro to do full type-dispatch for the implementation
753 // of a generic arithmetic operator. This is a pretty hideously large macro
754 // and it does a full dispatch in every single case, but still writing it
755 // just once is probably a good idea. Well if HAVE_SOFTFLOAT is not defined
756 // I will want the "long float" cases to be supressed!
757
758 // This checks for
759 // i fixnum
760 // b bignum
761 // r ratio
762 // c complex
763 // s short float
764 // f single float
765 // d double float
766 // l long float
767 //
768
769
770 #ifdef HAVE_SOFTFLOAT
771
772 // arith_dispatch_1 switches into one of 8 sub-functions whose names are
773 // got be suffixing the top-level name with one of the above letters.
774 // arith_dispatch_2 ends up with 64 sub-functions to call.
775
776 // First for 1-arg functions
777
778 #define arith_dispatch_1(stgclass, type, name) \
779 stgclass type name(LispObject a1) \
780 { if (is_fixnum(a1)) return name##_i(a1); \
781 switch (a1 & TAG_BITS) \
782 { \
783 case TAG_NUMBERS: \
784 switch (type_of_header(numhdr(a1))) \
785 { \
786 case TYPE_BIGNUM: \
787 return name##_b(a1); \
788 case TYPE_RATNUM: \
789 return name##_r(a1); \
790 case TYPE_COMPLEX_NUM: \
791 return name##_c(a1); \
792 default: \
793 return static_cast<type>(aerror1("bad arg for " #name, a1)); \
794 } \
795 case TAG_BOXFLOAT: \
796 switch (type_of_header(flthdr(a1))) \
797 { \
798 case TYPE_SINGLE_FLOAT: \
799 return name##_f(a1); \
800 case TYPE_DOUBLE_FLOAT: \
801 return name##_d(a1); \
802 case TYPE_LONG_FLOAT: \
803 return name##_l(a1); \
804 default: \
805 return static_cast<type>(aerror1("bad arg for " #name, a1)); \
806 } \
807 default: \
808 return static_cast<type>(aerror1("bad arg for " #name, a1)); \
809 case (XTAG_SFLOAT & TAG_BITS): \
810 if (SIXTY_FOUR_BIT && ((a1 & XTAG_FLOAT32) != 0) \
811 return name##_f(a1); \
812 else return name##_s(a1); \
813 } \
814 }
815
816 // Now a helper macro for 2-arg functions
817
818 #define arith_dispatch_1a(stgclass, type, name, rawname) \
819 stgclass type name(LispObject a1, LispObject a2) \
820 { if (is_fixnum(a2)) return name##_i(a1, a2); \
821 switch (a2 & TAG_BITS) \
822 { \
823 case TAG_NUMBERS: \
824 switch (type_of_header(numhdr(a2))) \
825 { \
826 case TYPE_BIGNUM: \
827 return name##_b(a1, a2); \
828 case TYPE_RATNUM: \
829 return name##_r(a1, a2); \
830 case TYPE_COMPLEX_NUM: \
831 return name##_c(a1, a2); \
832 default: \
833 return static_cast<type>(aerror2("bad arg for " #rawname, a1, a2)); \
834 } \
835 case TAG_BOXFLOAT: \
836 switch (type_of_header(flthdr(a2))) \
837 { \
838 case TYPE_SINGLE_FLOAT: \
839 return name##_f(a1, a2); \
840 case TYPE_DOUBLE_FLOAT: \
841 return name##_d(a1, a2); \
842 case TYPE_LONG_FLOAT: \
843 return name##_l(a1, a2); \
844 default: \
845 return static_cast<type>(aerror2("bad arg for " #rawname, a1, a2)); \
846 } \
847 default: \
848 return static_cast<type>(aerror2("bad arg for " #rawname, a1, a2)); \
849 case (XTAG_SFLOAT & TAG_BITS): \
850 if (SIXTY_FOUR_BIT && ((a2 & XTAG_FLOAT32) != 0)) \
851 return name##_f(a1, a2); \
852 else return name##_s(a1, a2); \
853 } \
854 }
855
856 #define arith_dispatch_2(stgclass, type, name) \
857 arith_dispatch_1a(inline, type, name##_i, name) \
858 \
859 arith_dispatch_1a(inline, type, name##_b, name) \
860 \
861 arith_dispatch_1a(inline, type, name##_r, name) \
862 \
863 arith_dispatch_1a(inline, type, name##_c, name) \
864 \
865 arith_dispatch_1a(inline, type, name##_s, name) \
866 \
867 arith_dispatch_1a(inline, type, name##_f, name) \
868 \
869 arith_dispatch_1a(inline, type, name##_d, name) \
870 \
871 arith_dispatch_1a(inline, type, name##_l, name) \
872 \
873 stgclass type name(LispObject a1, LispObject a2) \
874 { if (is_fixnum(a1)) return name##_i(a1, a2); \
875 switch (a1 & TAG_BITS) \
876 { \
877 case TAG_NUMBERS: \
878 switch (type_of_header(numhdr(a1))) \
879 { \
880 case TYPE_BIGNUM: \
881 return name##_b(a1, a2); \
882 case TYPE_RATNUM: \
883 return name##_r(a1, a2); \
884 case TYPE_COMPLEX_NUM: \
885 return name##_c(a1, a2); \
886 default: \
887 return static_cast<type>(aerror2("bad arg for " #name, a1, a2)); \
888 } \
889 case TAG_BOXFLOAT: \
890 switch (type_of_header(flthdr(a1))) \
891 { \
892 case TYPE_SINGLE_FLOAT: \
893 return name##_f(a1, a2); \
894 case TYPE_DOUBLE_FLOAT: \
895 return name##_d(a1, a2); \
896 case TYPE_LONG_FLOAT: \
897 return name##_l(a1, a2); \
898 default: \
899 return static_cast<type>(aerror2("bad arg for " #name, a1, a2)); \
900 } \
901 default: \
902 return static_cast<type>(aerror2("bad arg for " #name, a1, a2)); \
903 case (XTAG_SFLOAT & TAG_BITS): \
904 if (SIXTY_FOUR_BIT && ((a1 & XTAG_FLOAT32) != 0)) \
905 return name##_f(a1, a2); \
906 else return name##_s(a1, a2); \
907 } \
908 }
909
910 #else // HAVE_SOFTFLOAT
911
912 // arith_dispatch_1 switches into one of 7 sub-functions whose names are
913 // got be suffixing the top-level name with one of the above letters.
914 // arith_dispatch_2 ends up with 49 sub-functions to call.
915
916 // First for 1-arg functions
917
918 #define arith_dispatch_1(stgclass, type, name) \
919 stgclass type name(LispObject a1) \
920 { if (is_fixnum(a1)) return name##_i(a1); \
921 switch (a1 & TAG_BITS) \
922 { \
923 case TAG_NUMBERS: \
924 switch (type_of_header(numhdr(a1))) \
925 { \
926 case TYPE_BIGNUM: \
927 return name##_b(a1); \
928 case TYPE_RATNUM: \
929 return name##_r(a1); \
930 case TYPE_COMPLEX_NUM: \
931 return name##_c(a1); \
932 default: \
933 return static_cast<type>(aerror1("bad arg for " #name, a1)); \
934 } \
935 case TAG_BOXFLOAT: \
936 switch (type_of_header(flthdr(a1))) \
937 { \
938 case TYPE_SINGLE_FLOAT: \
939 return name##_f(a1); \
940 case TYPE_DOUBLE_FLOAT: \
941 return name##_d(a1); \
942 default: \
943 return static_cast<type>(aerror1("bad arg for " #name, a1)); \
944 } \
945 default: \
946 return static_cast<type>(aerror1("bad arg for " #name, a1)); \
947 case (XTAG_SFLOAT & TAG_BITS): \
948 if (SIXTY_FOUR_BIT && ((a1 & XTAG_FLOAT32) != 0)) \
949 return name##_f(a1); \
950 else return name##_s(a1); \
951 } \
952 }
953
954 // Now a helper macro for 2-arg functions
955
956 #define arith_dispatch_1a(stgclass, type, name, rawname) \
957 stgclass type name(LispObject a1, LispObject a2) \
958 { if (is_fixnum(a2)) return name##_i(a1, a2); \
959 switch (a2 & TAG_BITS) \
960 { \
961 case TAG_NUMBERS: \
962 switch (type_of_header(numhdr(a2))) \
963 { \
964 case TYPE_BIGNUM: \
965 return name##_b(a1, a2); \
966 case TYPE_RATNUM: \
967 return name##_r(a1, a2); \
968 case TYPE_COMPLEX_NUM: \
969 return name##_c(a1, a2); \
970 default: \
971 return static_cast<type>(aerror2("bad arg for " #rawname, a1, a2)); \
972 } \
973 case TAG_BOXFLOAT: \
974 switch (type_of_header(flthdr(a2))) \
975 { \
976 case TYPE_SINGLE_FLOAT: \
977 return name##_f(a1, a2); \
978 case TYPE_DOUBLE_FLOAT: \
979 return name##_d(a1, a2); \
980 default: \
981 return static_cast<type>(aerror2("bad arg for " #rawname, a1, a2)); \
982 } \
983 default: \
984 return static_cast<type>(aerror2("bad arg for " #rawname, a1, a2)); \
985 case (XTAG_SFLOAT & TAG_BITS): \
986 if (SIXTY_FOUR_BIT && ((a2 & XTAG_FLOAT32) != 0)) \
987 return name##_f(a1, a2); \
988 else return name##_s(a1, a2); \
989 } \
990 }
991
992 #define arith_dispatch_2(stgclass, type, name) \
993 arith_dispatch_1a(inline, type, name##_i, name) \
994 \
995 arith_dispatch_1a(inline, type, name##_b, name) \
996 \
997 arith_dispatch_1a(inline, type, name##_r, name) \
998 \
999 arith_dispatch_1a(inline, type, name##_c, name) \
1000 \
1001 arith_dispatch_1a(inline, type, name##_s, name) \
1002 \
1003 arith_dispatch_1a(inline, type, name##_f, name) \
1004 \
1005 arith_dispatch_1a(inline, type, name##_d, name) \
1006 \
1007 stgclass type name(LispObject a1, LispObject a2) \
1008 { if (is_fixnum(a1)) return name##_i(a1, a2); \
1009 switch (a1 & TAG_BITS) \
1010 { \
1011 case TAG_NUMBERS: \
1012 switch (type_of_header(numhdr(a1))) \
1013 { \
1014 case TYPE_BIGNUM: \
1015 return name##_b(a1, a2); \
1016 case TYPE_RATNUM: \
1017 return name##_r(a1, a2); \
1018 case TYPE_COMPLEX_NUM: \
1019 return name##_c(a1, a2); \
1020 default: \
1021 return static_cast<type>(aerror2("bad arg for " #name, a1, a2)); \
1022 } \
1023 case TAG_BOXFLOAT: \
1024 switch (type_of_header(flthdr(a1))) \
1025 { \
1026 case TYPE_SINGLE_FLOAT: \
1027 return name##_f(a1, a2); \
1028 case TYPE_DOUBLE_FLOAT: \
1029 return name##_d(a1, a2); \
1030 default: \
1031 return static_cast<type>(aerror2("bad arg for " #name, a1, a2)); \
1032 } \
1033 default: \
1034 return static_cast<type>(aerror2("bad arg for " #name, a1, a2)); \
1035 case (XTAG_SFLOAT & TAG_BITS): \
1036 if (SIXTY_FOUR_BIT && ((a1 & XTAG_FLOAT32) != 0)) \
1037 return name##_f(a1, a2); \
1038 else return name##_s(a1, a2); \
1039 } \
1040 }
1041
1042 #endif // HAVE_SOFTFLOAT
1043
1044 #endif // header_arith_h
1045
1046 // end of arith.h
1047