mirror of
https://github.com/lcn2/calc.git
synced 2025-08-16 01:03:29 +03:00
db77e29a23
Converted all ASCII tabs to ASCII spaces using a 8 character tab stop, for all files, except for all Makefiles (plus rpm.mk). The `git diff -w` reports no changes.
2314 lines
75 KiB
C
2314 lines
75 KiB
C
/*
|
|
* zfunc - extended precision integral arithmetic non-primitive routines
|
|
*
|
|
* Copyright (C) 1999-2007,2021-2023 David I. Bell, Landon Curt Noll and Ernest Bowen
|
|
*
|
|
* Primary author: David I. Bell
|
|
*
|
|
* Calc is open software; you can redistribute it and/or modify it under
|
|
* the terms of the version 2.1 of the GNU Lesser General Public License
|
|
* as published by the Free Software Foundation.
|
|
*
|
|
* Calc is distributed in the hope that it will be useful, but WITHOUT
|
|
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
|
|
* Public License for more details.
|
|
*
|
|
* A copy of version 2.1 of the GNU Lesser General Public License is
|
|
* distributed with calc under the filename COPYING-LGPL. You should have
|
|
* received a copy with calc; if not, write to Free Software Foundation, Inc.
|
|
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
*
|
|
* Under source code control: 1990/02/15 01:48:27
|
|
* File existed as early as: before 1990
|
|
*
|
|
* Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
|
|
*/
|
|
|
|
|
|
#include "zmath.h"
|
|
#include "alloc.h"
|
|
|
|
|
|
#include "errtbl.h"
|
|
#include "banned.h" /* include after system header <> includes */
|
|
|
|
|
|
ZVALUE _tenpowers_[TEN_MAX+1]; /* table of 10^2^n */
|
|
|
|
STATIC long *power10 = NULL;
|
|
STATIC int max_power10_exp = 0;
|
|
|
|
/*
|
|
* given:
|
|
*
|
|
* unsigned long x
|
|
* or: unsigned long long x
|
|
* or: long x and x >= 0
|
|
* or: long long x and x >= 0
|
|
*
|
|
* If issq_mod4k[x & 0xfff] == 0, then x cannot be a perfect square
|
|
* else x might be a perfect square.
|
|
*/
|
|
STATIC USB8 issq_mod4k[1<<12] = {
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
|
|
};
|
|
|
|
|
|
/*
|
|
* Compute the factorial of a number.
|
|
*/
|
|
void
|
|
zfact(ZVALUE z, ZVALUE *dest)
|
|
{
|
|
long ptwo; /* count of powers of two */
|
|
long n; /* current multiplication value */
|
|
long m; /* reduced multiplication value */
|
|
long mul; /* collected value to multiply by */
|
|
ZVALUE res, temp;
|
|
|
|
/* firewall */
|
|
if (dest == NULL) {
|
|
math_error("%s: dest NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (zisneg(z)) {
|
|
math_error("Negative argument for factorial");
|
|
not_reached();
|
|
}
|
|
if (zge31b(z)) {
|
|
math_error("Very large factorial");
|
|
not_reached();
|
|
}
|
|
n = ztolong(z);
|
|
ptwo = 0;
|
|
mul = 1;
|
|
res = _one_;
|
|
/*
|
|
* Multiply numbers together, but squeeze out all powers of two.
|
|
* We will put them back in at the end. Also collect multiple
|
|
* numbers together until there is a risk of overflow.
|
|
*/
|
|
for (; n > 1; n--) {
|
|
for (m = n; ((m & 0x1) == 0); m >>= 1)
|
|
ptwo++;
|
|
if (mul <= MAXLONG/m) {
|
|
mul *= m;
|
|
continue;
|
|
}
|
|
zmuli(res, mul, &temp);
|
|
zfree(res);
|
|
res = temp;
|
|
mul = m;
|
|
}
|
|
/*
|
|
* Multiply by the remaining value, then scale result by
|
|
* the proper power of two.
|
|
*/
|
|
if (mul > 1) {
|
|
zmuli(res, mul, &temp);
|
|
zfree(res);
|
|
res = temp;
|
|
}
|
|
zshift(res, ptwo, &temp);
|
|
zfree(res);
|
|
*dest = temp;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the permutation function M! / (M - N)!.
|
|
*/
|
|
void
|
|
zperm(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
SFULL count;
|
|
ZVALUE cur, tmp, ans;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (zisneg(z1) || zisneg(z2)) {
|
|
math_error("Negative argument for permutation");
|
|
not_reached();
|
|
}
|
|
if (zrel(z1, z2) < 0) {
|
|
math_error("Second arg larger than first in permutation");
|
|
not_reached();
|
|
}
|
|
if (zge31b(z2)) {
|
|
math_error("Very large permutation");
|
|
not_reached();
|
|
}
|
|
count = ztolong(z2);
|
|
zcopy(z1, &ans);
|
|
zsub(z1, _one_, &cur);
|
|
while (--count > 0) {
|
|
zmul(ans, cur, &tmp);
|
|
zfree(ans);
|
|
ans = tmp;
|
|
zsub(cur, _one_, &tmp);
|
|
zfree(cur);
|
|
cur = tmp;
|
|
}
|
|
zfree(cur);
|
|
*res = ans;
|
|
}
|
|
|
|
/*
|
|
* docomb evaluates binomial coefficient when z1 >= 0, z2 >= 0
|
|
*/
|
|
S_FUNC int
|
|
docomb(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
ZVALUE ans;
|
|
ZVALUE mul, div, temp;
|
|
FULL count, i;
|
|
#if BASEB == 16
|
|
HALF dh[2];
|
|
#else
|
|
HALF dh[1];
|
|
#endif
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (zrel(z2, z1) > 0)
|
|
return 0;
|
|
zsub(z1, z2, &temp);
|
|
|
|
if (zge31b(z2) && zge31b(temp)) {
|
|
zfree(temp);
|
|
return -2;
|
|
}
|
|
if (zrel(temp, z2) < 0)
|
|
count = ztofull(temp);
|
|
else
|
|
count = ztofull(z2);
|
|
zfree(temp);
|
|
if (count == 0)
|
|
return 1;
|
|
if (count == 1)
|
|
return 2;
|
|
div.sign = 0;
|
|
div.v = dh;
|
|
div.len = 1;
|
|
zcopy(z1, &mul);
|
|
zcopy(z1, &ans);
|
|
for (i = 2; i <= count; i++) {
|
|
#if BASEB == 16
|
|
dh[0] = (HALF)(i & BASE1);
|
|
dh[1] = (HALF)(i >> BASEB);
|
|
div.len = 1 + (dh[1] != 0);
|
|
#else
|
|
dh[0] = (HALF) i;
|
|
#endif
|
|
zsub(mul, _one_, &temp);
|
|
zfree(mul);
|
|
mul = temp;
|
|
zmul(ans, mul, &temp);
|
|
zfree(ans);
|
|
zquo(temp, div, &ans, 0);
|
|
zfree(temp);
|
|
}
|
|
zfree(mul);
|
|
*res = ans;
|
|
return 3;
|
|
}
|
|
|
|
/*
|
|
* Compute the combinatorial function M! / ( N! * (M - N)! ).
|
|
* Returns 0 if result is 0
|
|
* 1 1
|
|
* 2 z1
|
|
* -1 -1
|
|
* -2 if too complicated
|
|
* 3 result stored at res
|
|
*/
|
|
int
|
|
zcomb(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
ZVALUE z3, z4;
|
|
int r;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (z2.sign || (!z1.sign && zrel(z2, z1) > 0))
|
|
return 0;
|
|
if (zisone(z2))
|
|
return 2;
|
|
if (z1.sign) {
|
|
z1.sign = 0;
|
|
zsub(z1, _one_, &z3);
|
|
zadd(z3, z2, &z4);
|
|
zfree(z3);
|
|
r = docomb(z4, z2, res);
|
|
if (r == 2) {
|
|
*res = z4;
|
|
r = 3;
|
|
}
|
|
else
|
|
zfree(z4);
|
|
if (z2.v[0] & 1) {
|
|
if (r == 1)
|
|
r = -1;
|
|
if (r == 3)
|
|
res->sign = 1;
|
|
}
|
|
return r;
|
|
}
|
|
return docomb(z1, z2, res);
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the Jacobi function (m / n) for odd n.
|
|
*
|
|
* The property of the Jacobi function is: If n>2 is prime then
|
|
*
|
|
* the result is 1 if m == x^2 (mod n) for some x.
|
|
* otherwise the result is -1.
|
|
*
|
|
* If n is not prime, then the result does not prove that n is not prime
|
|
* when the value of the Jacobi is 1.
|
|
*
|
|
* Jacobi evaluation of (m / n), where n > 0 is odd AND m > 0 is odd:
|
|
*
|
|
* rule 0: (0 / n) == 0
|
|
* rule 1: (1 / n) == 1
|
|
* rule 2: (m / n) == (a / n) if m == a % n
|
|
* rule 3: (m / n) == (2*m / n) if n == 1 % 8 OR n == 7 % 8
|
|
* rule 4: (m / n) == -(2*m / n) if n != 1 & 8 AND n != 7 % 8
|
|
* rule 5: (m / n) == (n / m) if m == 3 % 4 AND n == 3 % 4
|
|
* rule 6: (m / n) == -(n / m) if m != 3 % 4 OR n != 3 % 4
|
|
*
|
|
* NOTE: This function returns 0 in invalid Jacobi parameters:
|
|
* m < 0 OR n is even OR n < 1.
|
|
*/
|
|
FLAG
|
|
zjacobi(ZVALUE z1, ZVALUE z2)
|
|
{
|
|
ZVALUE p, q, tmp;
|
|
long lowbit;
|
|
int val;
|
|
|
|
/* firewall */
|
|
if (ziszero(z1) || zisneg(z1))
|
|
return 0;
|
|
if (ziseven(z2) || zisneg(z2))
|
|
return 0;
|
|
|
|
/* assume a value of 1 unless we find otherwise */
|
|
if (zisone(z1))
|
|
return 1;
|
|
val = 1;
|
|
zcopy(z1, &p);
|
|
zcopy(z2, &q);
|
|
for (;;) {
|
|
zmod(p, q, &tmp, 0);
|
|
zfree(p);
|
|
p = tmp;
|
|
if (ziszero(p)) {
|
|
zfree(p);
|
|
zfree(q);
|
|
return 0;
|
|
}
|
|
if (ziseven(p)) {
|
|
lowbit = zlowbit(p);
|
|
zshift(p, -lowbit, &tmp);
|
|
zfree(p);
|
|
p = tmp;
|
|
if ((lowbit & 1) && (((*q.v & 0x7) == 3) ||
|
|
((*q.v & 0x7) == 5)))
|
|
val = -val;
|
|
}
|
|
if (zisunit(p)) {
|
|
zfree(p);
|
|
zfree(q);
|
|
return val;
|
|
}
|
|
if ((*p.v & *q.v & 0x3) == 3)
|
|
val = -val;
|
|
tmp = q;
|
|
q = p;
|
|
p = tmp;
|
|
}
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the Fibonacci number F(n).
|
|
* This is evaluated by recursively using the formulas:
|
|
* F(2N+1) = F(N+1)^2 + F(N)^2
|
|
* and
|
|
* F(2N) = F(N+1)^2 - F(N-1)^2
|
|
*/
|
|
void
|
|
zfib(ZVALUE z, ZVALUE *res)
|
|
{
|
|
long n;
|
|
int sign;
|
|
ZVALUE fnm1, fn, fnp1; /* consecutive Fibonacci values */
|
|
ZVALUE t1, t2, t3;
|
|
FULL i;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (zge31b(z)) {
|
|
math_error("Very large Fibonacci number");
|
|
not_reached();
|
|
}
|
|
n = ztolong(z);
|
|
if (n == 0) {
|
|
*res = _zero_;
|
|
return;
|
|
}
|
|
sign = z.sign && ((n & 0x1) == 0);
|
|
if (n <= 2) {
|
|
*res = _one_;
|
|
res->sign = (bool)sign;
|
|
return;
|
|
}
|
|
i = TOPFULL;
|
|
while ((i & n) == 0)
|
|
i >>= (FULL)1;
|
|
i >>= (FULL)1;
|
|
fnm1 = _zero_;
|
|
fn = _one_;
|
|
fnp1 = _one_;
|
|
while (i) {
|
|
zsquare(fnm1, &t1);
|
|
zsquare(fn, &t2);
|
|
zsquare(fnp1, &t3);
|
|
zfree(fnm1);
|
|
zfree(fn);
|
|
zfree(fnp1);
|
|
zadd(t2, t3, &fnp1);
|
|
zsub(t3, t1, &fn);
|
|
zfree(t1);
|
|
zfree(t2);
|
|
zfree(t3);
|
|
if (i & n) {
|
|
fnm1 = fn;
|
|
fn = fnp1;
|
|
zadd(fnm1, fn, &fnp1);
|
|
} else {
|
|
zsub(fnp1, fn, &fnm1);
|
|
}
|
|
i >>= (FULL)1;
|
|
}
|
|
zfree(fnm1);
|
|
zfree(fnp1);
|
|
*res = fn;
|
|
res->sign = (bool)sign;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the result of raising one number to the power of another
|
|
* The second number is assumed to be non-negative.
|
|
* It cannot be too large except for trivial cases.
|
|
*/
|
|
void
|
|
zpowi(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
int sign; /* final sign of number */
|
|
unsigned long power; /* power to raise to */
|
|
FULL bit; /* current bit value */
|
|
long twos; /* count of times 2 is in result */
|
|
ZVALUE ans, temp;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
sign = (z1.sign && zisodd(z2));
|
|
z1.sign = 0;
|
|
z2.sign = 0;
|
|
if (ziszero(z2) && !ziszero(z1)) { /* number raised to power 0 */
|
|
*res = _one_;
|
|
return;
|
|
}
|
|
if (zisabsleone(z1)) { /* 0, 1, or -1 raised to a power */
|
|
ans = _one_;
|
|
ans.sign = (bool)sign;
|
|
if (*z1.v == 0)
|
|
ans = _zero_;
|
|
*res = ans;
|
|
return;
|
|
}
|
|
if (zge31b(z2)) {
|
|
math_error("Raising to very large power");
|
|
not_reached();
|
|
}
|
|
power = ztoulong(z2);
|
|
if (zistwo(z1)) { /* two raised to a power */
|
|
zbitvalue((long) power, res);
|
|
return;
|
|
}
|
|
/*
|
|
* See if this is a power of ten
|
|
*/
|
|
if (zistiny(z1) && (*z1.v == 10)) {
|
|
ztenpow((long) power, res);
|
|
res->sign = (bool)sign;
|
|
return;
|
|
}
|
|
/*
|
|
* Handle low powers specially
|
|
*/
|
|
if (power <= 4) {
|
|
switch ((int) power) {
|
|
case 1:
|
|
ans.len = z1.len;
|
|
ans.v = alloc(ans.len);
|
|
zcopyval(z1, ans);
|
|
ans.sign = (bool)sign;
|
|
*res = ans;
|
|
return;
|
|
case 2:
|
|
zsquare(z1, res);
|
|
return;
|
|
case 3:
|
|
zsquare(z1, &temp);
|
|
zmul(z1, temp, res);
|
|
zfree(temp);
|
|
res->sign = (bool)sign;
|
|
return;
|
|
case 4:
|
|
zsquare(z1, &temp);
|
|
zsquare(temp, res);
|
|
zfree(temp);
|
|
return;
|
|
}
|
|
}
|
|
/*
|
|
* Shift out all powers of twos so the multiplies are smaller.
|
|
* We will shift back the right amount when done.
|
|
*/
|
|
twos = 0;
|
|
if (ziseven(z1)) {
|
|
twos = zlowbit(z1);
|
|
ans.v = alloc(z1.len);
|
|
ans.len = z1.len;
|
|
ans.sign = z1.sign;
|
|
zcopyval(z1, ans);
|
|
zshiftr(ans, twos);
|
|
ztrim(&ans);
|
|
z1 = ans;
|
|
twos *= power;
|
|
}
|
|
/*
|
|
* Compute the power by squaring and multiplying.
|
|
* This uses the left to right method of power raising.
|
|
*/
|
|
bit = TOPFULL;
|
|
while ((bit & power) == 0)
|
|
bit >>= 1;
|
|
bit >>= 1;
|
|
zsquare(z1, &ans);
|
|
if (bit & power) {
|
|
zmul(ans, z1, &temp);
|
|
zfree(ans);
|
|
ans = temp;
|
|
}
|
|
bit >>= 1;
|
|
while (bit) {
|
|
zsquare(ans, &temp);
|
|
zfree(ans);
|
|
ans = temp;
|
|
if (bit & power) {
|
|
zmul(ans, z1, &temp);
|
|
zfree(ans);
|
|
ans = temp;
|
|
}
|
|
bit >>= 1;
|
|
}
|
|
/*
|
|
* Scale back up by proper power of two
|
|
*/
|
|
if (twos) {
|
|
zshift(ans, twos, &temp);
|
|
zfree(ans);
|
|
ans = temp;
|
|
zfree(z1);
|
|
}
|
|
ans.sign = (bool)sign;
|
|
*res = ans;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute ten to the specified power
|
|
* This saves some work since the squares of ten are saved.
|
|
*/
|
|
void
|
|
ztenpow(long power, ZVALUE *res)
|
|
{
|
|
long i;
|
|
ZVALUE ans;
|
|
ZVALUE temp;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (power <= 0) {
|
|
*res = _one_;
|
|
return;
|
|
}
|
|
ans = _one_;
|
|
_tenpowers_[0] = _ten_;
|
|
for (i = 0; power; i++) {
|
|
if (_tenpowers_[i].len == 0) {
|
|
if (i <= TEN_MAX) {
|
|
zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
|
|
} else {
|
|
math_error("cannot compute 10^2^(TEN_MAX+1)");
|
|
not_reached();
|
|
}
|
|
}
|
|
if (power & 0x1) {
|
|
zmul(ans, _tenpowers_[i], &temp);
|
|
zfree(ans);
|
|
ans = temp;
|
|
}
|
|
power /= 2;
|
|
}
|
|
*res = ans;
|
|
}
|
|
|
|
|
|
/*
|
|
* Calculate modular inverse suppressing unnecessary divisions.
|
|
* This is based on the Euclidean algorithm for large numbers.
|
|
* (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
|
|
* Returns true if there is no solution because the numbers
|
|
* are not relatively prime.
|
|
*/
|
|
bool
|
|
zmodinv(ZVALUE u, ZVALUE v, ZVALUE *res)
|
|
{
|
|
FULL q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
|
|
ZVALUE u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
v.sign = 0;
|
|
if (zisneg(u) || (zrel(u, v) >= 0))
|
|
zmod(u, v, &v3, 0);
|
|
else
|
|
zcopy(u, &v3);
|
|
zcopy(v, &u3);
|
|
u2 = _zero_;
|
|
v2 = _one_;
|
|
|
|
/*
|
|
* Loop here while the size of the numbers remain above
|
|
* the size of a HALF. Throughout this loop u3 >= v3.
|
|
*/
|
|
while ((u3.len > 1) && !ziszero(v3)) {
|
|
vh = 0;
|
|
#if LONG_BITS == BASEB
|
|
uh = u3.v[u3.len - 1];
|
|
if (v3.len == u3.len)
|
|
vh = v3.v[v3.len - 1];
|
|
#else
|
|
uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
|
|
if ((v3.len + 1) >= u3.len)
|
|
vh = v3.v[v3.len - 1];
|
|
if (v3.len == u3.len)
|
|
vh = (vh << BASEB) + v3.v[v3.len - 2];
|
|
#endif
|
|
A = 1;
|
|
B = 0;
|
|
C = 0;
|
|
D = 1;
|
|
|
|
/*
|
|
* Calculate successive quotients of the continued fraction
|
|
* expansion using only single precision arithmetic until
|
|
* greater precision is required.
|
|
*/
|
|
while ((vh + C) && (vh + D)) {
|
|
q1 = (uh + A) / (vh + C);
|
|
q2 = (uh + B) / (vh + D);
|
|
if (q1 != q2)
|
|
break;
|
|
T = A - q1 * C;
|
|
A = C;
|
|
C = T;
|
|
T = B - q1 * D;
|
|
B = D;
|
|
D = T;
|
|
T = uh - q1 * vh;
|
|
uh = vh;
|
|
vh = T;
|
|
}
|
|
|
|
/*
|
|
* If B is zero, then we made no progress because
|
|
* the calculation requires a very large quotient.
|
|
* So we must do this step of the calculation in
|
|
* full precision
|
|
*/
|
|
if (B == 0) {
|
|
zquo(u3, v3, &qz, 0);
|
|
zmul(qz, v2, &tmp1);
|
|
zsub(u2, tmp1, &tmp2);
|
|
zfree(tmp1);
|
|
zfree(u2);
|
|
u2 = v2;
|
|
v2 = tmp2;
|
|
zmul(qz, v3, &tmp1);
|
|
zsub(u3, tmp1, &tmp2);
|
|
zfree(tmp1);
|
|
zfree(u3);
|
|
u3 = v3;
|
|
v3 = tmp2;
|
|
zfree(qz);
|
|
continue;
|
|
}
|
|
/*
|
|
* Apply the calculated A,B,C,D numbers to the current
|
|
* values to update them as if the full precision
|
|
* calculations had been carried out.
|
|
*/
|
|
zmuli(u2, (long) A, &tmp1);
|
|
zmuli(v2, (long) B, &tmp2);
|
|
zadd(tmp1, tmp2, &tmp3);
|
|
zfree(tmp1);
|
|
zfree(tmp2);
|
|
zmuli(u2, (long) C, &tmp1);
|
|
zmuli(v2, (long) D, &tmp2);
|
|
zfree(u2);
|
|
zfree(v2);
|
|
u2 = tmp3;
|
|
zadd(tmp1, tmp2, &v2);
|
|
zfree(tmp1);
|
|
zfree(tmp2);
|
|
zmuli(u3, (long) A, &tmp1);
|
|
zmuli(v3, (long) B, &tmp2);
|
|
zadd(tmp1, tmp2, &tmp3);
|
|
zfree(tmp1);
|
|
zfree(tmp2);
|
|
zmuli(u3, (long) C, &tmp1);
|
|
zmuli(v3, (long) D, &tmp2);
|
|
zfree(u3);
|
|
zfree(v3);
|
|
u3 = tmp3;
|
|
zadd(tmp1, tmp2, &v3);
|
|
zfree(tmp1);
|
|
zfree(tmp2);
|
|
}
|
|
|
|
/*
|
|
* Here when the remaining numbers become single precision in size.
|
|
* Finish the procedure using single precision calculations.
|
|
*/
|
|
if (ziszero(v3) && !zisone(u3)) {
|
|
zfree(u3);
|
|
zfree(v3);
|
|
zfree(u2);
|
|
zfree(v2);
|
|
return true;
|
|
}
|
|
ui3 = ztofull(u3);
|
|
vi3 = ztofull(v3);
|
|
zfree(u3);
|
|
zfree(v3);
|
|
while (vi3) {
|
|
q1 = ui3 / vi3;
|
|
zmuli(v2, (long) q1, &tmp1);
|
|
zsub(u2, tmp1, &tmp2);
|
|
zfree(tmp1);
|
|
zfree(u2);
|
|
u2 = v2;
|
|
v2 = tmp2;
|
|
q2 = ui3 - q1 * vi3;
|
|
ui3 = vi3;
|
|
vi3 = q2;
|
|
}
|
|
zfree(v2);
|
|
if (ui3 != 1) {
|
|
zfree(u2);
|
|
return true;
|
|
}
|
|
if (zisneg(u2)) {
|
|
zadd(v, u2, res);
|
|
zfree(u2);
|
|
return false;
|
|
}
|
|
*res = u2;
|
|
return false;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the greatest common divisor of a pair of integers.
|
|
*/
|
|
void
|
|
zgcd(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
int h, i, j, k;
|
|
LEN len, l, m, n, o, p, q;
|
|
HALF u, v, w, x;
|
|
HALF *a, *a0, *A, *b, *b0, *B, *c, *d;
|
|
FULL f, g;
|
|
ZVALUE gcd;
|
|
bool needw;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (zisunit(z1) || zisunit(z2)) {
|
|
*res = _one_;
|
|
return;
|
|
}
|
|
z1.sign = 0;
|
|
z2.sign = 0;
|
|
if (ziszero(z1) || !zcmp(z1, z2)) {
|
|
zcopy(z2, res);
|
|
return;
|
|
}
|
|
if (ziszero(z2)) {
|
|
zcopy(z1, res);
|
|
return;
|
|
}
|
|
|
|
o = 0;
|
|
while (!(z1.v[o] | z2.v[o])) o++; /* Count common zero digits */
|
|
|
|
c = z1.v + o;
|
|
d = z2.v + o;
|
|
|
|
m = z1.len - o;
|
|
n = z2.len - o;
|
|
u = *c | *d; /* Count common zero bits */
|
|
v = 1;
|
|
p = 0;
|
|
while (!(u & v)) {
|
|
v <<= 1;
|
|
p++;
|
|
}
|
|
|
|
while (!*c) { /* Removing zero digits */
|
|
c++;
|
|
m--;
|
|
}
|
|
|
|
while (!*d) {
|
|
d++;
|
|
n--;
|
|
}
|
|
|
|
|
|
u = *d; /* Count zero bits for *d */
|
|
v = 1;
|
|
q = 0;
|
|
while (!(u & v)) {
|
|
v <<= 1;
|
|
q++;
|
|
}
|
|
|
|
a0 = A = alloc(m);
|
|
b0 = B = alloc(n);
|
|
|
|
memcpy(A, c, m * sizeof(HALF)); /* Copy c[] to A[] */
|
|
|
|
/* Copy d[] to B[], shifting if necessary */
|
|
if (q) {
|
|
i = n;
|
|
b = B + n;
|
|
d += n;
|
|
f = 0;
|
|
while (i--) {
|
|
f = f << BASEB | *--d;
|
|
*--b = (HALF) (f >> q);
|
|
}
|
|
if (B[n-1] == 0) n--;
|
|
}
|
|
else memcpy(B, d, n * sizeof(HALF));
|
|
|
|
if (n == 1) { /* One digit case; use Euclid's algorithm */
|
|
n = m;
|
|
b0 = A;
|
|
m = 1;
|
|
a0 = B;
|
|
if (m == 1) { /* a has one digit */
|
|
v = *a0;
|
|
if (v > 1) { /* Euclid's algorithm */
|
|
b = b0 + n;
|
|
i = n;
|
|
u = 0;
|
|
while (i--) {
|
|
f = (FULL) u << BASEB | *--b;
|
|
u = (HALF) (f % v);
|
|
}
|
|
while (u) { w = v % u; v = u; u = w; }
|
|
}
|
|
*b0 = v;
|
|
n = 1;
|
|
}
|
|
len = n + o;
|
|
gcd.v = alloc(len + 1);
|
|
/* Common zero digits */
|
|
if (o) memset(gcd.v, 0, o * sizeof(HALF));
|
|
/* Left shift for common zero bits */
|
|
if (p) {
|
|
i = n;
|
|
f = 0;
|
|
b = b0;
|
|
a = gcd.v + o;
|
|
while (i--) {
|
|
f = f >> BASEB | (FULL) *b++ << p;
|
|
*a++ = (HALF) f;
|
|
}
|
|
if (f >>= BASEB) {len++; *a = (HALF) f;}
|
|
} else {
|
|
memcpy(gcd.v + o, b0, n * sizeof(HALF));
|
|
}
|
|
gcd.len = len;
|
|
gcd.sign = 0;
|
|
freeh(A);
|
|
freeh(B);
|
|
*res = gcd;
|
|
return;
|
|
}
|
|
|
|
u = B[n-1]; /* Bit count for b */
|
|
k = (n - 1) * BASEB;
|
|
while (u >>= 1) k++;
|
|
|
|
needw = true;
|
|
|
|
w = 0;
|
|
j = 0;
|
|
while (m) { /* START OF MAIN LOOP */
|
|
if (m - n < 2 || needw) {
|
|
q = 0;
|
|
u = *a0;
|
|
v = 1;
|
|
while (!(u & v)) { /* count zero bits for *a0 */
|
|
q++;
|
|
v <<= 1;
|
|
}
|
|
|
|
if (q) { /* right-justify a */
|
|
a = a0 + m;
|
|
i = m;
|
|
f = 0;
|
|
while (i--) {
|
|
f = f << BASEB | *--a;
|
|
*a = (HALF) (f >> q);
|
|
}
|
|
if (!a0[m-1]) m--; /* top digit vanishes */
|
|
}
|
|
|
|
if (m == 1) break;
|
|
|
|
u = a0[m-1];
|
|
j = (m - 1) * BASEB;
|
|
while (u >>= 1) j++; /* counting bits for a */
|
|
h = j - k;
|
|
if (h < 0) { /* swapping to get h > 0 */
|
|
l = m;
|
|
m = n;
|
|
n = l;
|
|
a = a0;
|
|
a0 = b0;
|
|
b0 = a;
|
|
k = j;
|
|
h = -h;
|
|
needw = true;
|
|
}
|
|
if (h > 1) {
|
|
if (needw) { /* find w = minv(*b0, h0) */
|
|
u = 1;
|
|
v = *b0;
|
|
w = 0;
|
|
x = 1;
|
|
i = h;
|
|
while (i-- && x) {
|
|
if (u & x) { u -= v * x; w |= x;}
|
|
x <<= 1;
|
|
}
|
|
needw = false;
|
|
}
|
|
g = (FULL) (*a0 * w);
|
|
if (h < BASEB) {
|
|
g &= (FULL)lowhalf[h];
|
|
} else {
|
|
g &= BASE1;
|
|
}
|
|
} else {
|
|
g = 1;
|
|
}
|
|
} else {
|
|
g = (FULL) (*a0 * w);
|
|
}
|
|
a = a0;
|
|
b = b0;
|
|
i = n;
|
|
if (g > 1) { /* a - g * b case */
|
|
f = 0;
|
|
while (i--) {
|
|
f = (FULL) *a - g * *b++ - f;
|
|
*a++ = (HALF) f;
|
|
f >>= BASEB;
|
|
f = -f & BASE1;
|
|
}
|
|
if (f) {
|
|
i = m - n;
|
|
while (i-- && f) {
|
|
f = *a - f;
|
|
*a++ = (HALF) f;
|
|
f >>= BASEB;
|
|
f = -f & BASE1;
|
|
}
|
|
}
|
|
while (m && !*a0) { /* Removing trailing zeros */
|
|
m--;
|
|
a0++;
|
|
}
|
|
if (f) { /* a - g * b < 0 */
|
|
while (m > 1 && a0[m-1] == BASE1) m--;
|
|
*a0 = - *a0;
|
|
a = a0;
|
|
i = m;
|
|
while (--i) {
|
|
a++;
|
|
*a = ~*a;
|
|
}
|
|
}
|
|
} else { /* abs(a - b) case */
|
|
while (i && *a++ == *b++) i--;
|
|
q = n - i;
|
|
if (m == n) { /* a and b same length */
|
|
if (i) { /* a not equal to b */
|
|
while (m && a0[m-1] == b0[m-1]) m--;
|
|
if (a0[m-1] < b0[m-1]) {
|
|
/* Swapping since a < b */
|
|
a = a0;
|
|
a0 = b0;
|
|
b0 = a;
|
|
k = j;
|
|
}
|
|
a = a0 + q;
|
|
b = b0 + q;
|
|
i = m - q;
|
|
f = 0;
|
|
while (i--) {
|
|
f = (FULL) *a - *b++ - f;
|
|
*a++ = (HALF) f;
|
|
f >>= BASEB;
|
|
f = -f & BASE1;
|
|
}
|
|
}
|
|
} else { /* a has more digits than b */
|
|
a = a0 + q;
|
|
b = b0 + q;
|
|
i = n - q;
|
|
f = 0;
|
|
while (i--) {
|
|
f = (FULL) *a - *b++ - f;
|
|
*a++ = (HALF) f;
|
|
f >>= BASEB;
|
|
f = -f & BASE1;
|
|
}
|
|
if (f) { while (!*a) *a++ = BASE1;
|
|
(*a)--;
|
|
}
|
|
}
|
|
a0 += q;
|
|
m -= q;
|
|
while (m && !*a0) { /* Removing trailing zeros */
|
|
m--;
|
|
a0++;
|
|
}
|
|
}
|
|
while (m && !a0[m-1]) m--; /* Removing leading zeros */
|
|
}
|
|
if (m == 1) { /* a has one digit */
|
|
v = *a0;
|
|
if (v > 1) { /* Euclid's algorithm */
|
|
b = b0 + n;
|
|
i = n;
|
|
u = 0;
|
|
while (i--) {
|
|
f = (FULL) u << BASEB | *--b;
|
|
u = (HALF) (f % v);
|
|
}
|
|
while (u) { w = v % u; v = u; u = w; }
|
|
}
|
|
*b0 = v;
|
|
n = 1;
|
|
}
|
|
len = n + o;
|
|
gcd.v = alloc(len + 1);
|
|
if (o) memset(gcd.v, 0, o * sizeof(HALF)); /* Common zero digits */
|
|
if (p) { /* Left shift for common zero bits */
|
|
i = n;
|
|
f = 0;
|
|
b = b0;
|
|
a = gcd.v + o;
|
|
while (i--) {
|
|
f = (FULL) *b++ << p | f;
|
|
*a++ = (HALF) f;
|
|
f >>= BASEB;
|
|
}
|
|
if (f) {
|
|
len++; *a = (HALF) f;
|
|
}
|
|
} else {
|
|
memcpy(gcd.v + o, b0, n * sizeof(HALF));
|
|
}
|
|
gcd.len = len;
|
|
gcd.sign = 0;
|
|
freeh(A);
|
|
freeh(B);
|
|
*res = gcd;
|
|
return;
|
|
}
|
|
|
|
/*
|
|
* Compute the lcm of two integers (least common multiple).
|
|
* This is done using the formula: gcd(a,b) * lcm(a,b) = a * b.
|
|
*/
|
|
void
|
|
zlcm(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
ZVALUE temp1, temp2;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
zgcd(z1, z2, &temp1);
|
|
zequo(z1, temp1, &temp2);
|
|
zfree(temp1);
|
|
zmul(temp2, z2, res);
|
|
zfree(temp2);
|
|
}
|
|
|
|
|
|
/*
|
|
* Return whether or not two numbers are relatively prime to each other.
|
|
*/
|
|
bool
|
|
zrelprime(ZVALUE z1, ZVALUE z2)
|
|
{
|
|
FULL rem1, rem2; /* remainders */
|
|
ZVALUE rem;
|
|
bool result;
|
|
|
|
z1.sign = 0;
|
|
z2.sign = 0;
|
|
if (ziseven(z1) && ziseven(z2)) /* false if both even */
|
|
return false;
|
|
if (zisunit(z1) || zisunit(z2)) /* true if either is a unit */
|
|
return true;
|
|
if (ziszero(z1) || ziszero(z2)) /* false if either is zero */
|
|
return false;
|
|
if (zistwo(z1) || zistwo(z2)) /* true if either is two */
|
|
return true;
|
|
/*
|
|
* Try reducing each number by the product of the first few odd primes
|
|
* to see if any of them are a common factor.
|
|
*/
|
|
rem1 = zmodi(z1, (FULL)3 * 5 * 7 * 11 * 13);
|
|
rem2 = zmodi(z2, (FULL)3 * 5 * 7 * 11 * 13);
|
|
if (((rem1 % 3) == 0) && ((rem2 % 3) == 0))
|
|
return false;
|
|
if (((rem1 % 5) == 0) && ((rem2 % 5) == 0))
|
|
return false;
|
|
if (((rem1 % 7) == 0) && ((rem2 % 7) == 0))
|
|
return false;
|
|
if (((rem1 % 11) == 0) && ((rem2 % 11) == 0))
|
|
return false;
|
|
if (((rem1 % 13) == 0) && ((rem2 % 13) == 0))
|
|
return false;
|
|
/*
|
|
* Try a new batch of primes now
|
|
*/
|
|
rem1 = zmodi(z1, (FULL)17 * 19 * 23);
|
|
rem2 = zmodi(z2, (FULL)17 * 19 * 23);
|
|
if (((rem1 % 17) == 0) && ((rem2 % 17) == 0))
|
|
return false;
|
|
if (((rem1 % 19) == 0) && ((rem2 % 19) == 0))
|
|
return false;
|
|
if (((rem1 % 23) == 0) && ((rem2 % 23) == 0))
|
|
return false;
|
|
/*
|
|
* Yuk, we must actually compute the gcd to know the answer
|
|
*/
|
|
zgcd(z1, z2, &rem);
|
|
result = zisunit(rem);
|
|
zfree(rem);
|
|
return result;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the integer floor of the log of an integer to a specified base.
|
|
* The signs of the integers and base are ignored.
|
|
* Example: zlog(123456, 10) = 5.
|
|
*/
|
|
long
|
|
zlog(ZVALUE z, ZVALUE base)
|
|
{
|
|
ZVALUE *zp; /* current square */
|
|
long power; /* current power */
|
|
ZVALUE temp; /* temporary */
|
|
ZVALUE squares[32]; /* table of squares of base */
|
|
|
|
/* ignore signs */
|
|
|
|
z.sign = 0;
|
|
base.sign = 0;
|
|
|
|
/*
|
|
* Make sure that the numbers are nonzero and the base is > 1
|
|
*/
|
|
if (ziszero(z) || ziszero(base) || zisone(base)) {
|
|
math_error("Zero or too small argument argument for zlog!!!");
|
|
not_reached();
|
|
}
|
|
|
|
/*
|
|
* Some trivial cases.
|
|
*/
|
|
power = zrel(z, base);
|
|
if (power <= 0)
|
|
return (power + 1);
|
|
|
|
/* base - power of two */
|
|
if (zisonebit(base))
|
|
return (zhighbit(z) / zlowbit(base));
|
|
|
|
/* base = 10 */
|
|
if (base.len == 1 && base.v[0] == 10)
|
|
return zlog10(z, NULL);
|
|
/*
|
|
* Now loop by squaring the base each time, and see whether or
|
|
* not each successive square is still smaller than the number.
|
|
*/
|
|
zp = &squares[0];
|
|
*zp = base;
|
|
while (zp->len * 2 - 1 <= z.len && zrel(z, *zp) > 0) {
|
|
/* while square not too large */
|
|
zsquare(*zp, zp + 1);
|
|
zp++;
|
|
}
|
|
|
|
/*
|
|
* Now back down the squares,
|
|
*/
|
|
power = 0;
|
|
for (; zp > squares; zp--) {
|
|
if (zrel(z, *zp) >= 0) {
|
|
zquo(z, *zp, &temp, 0);
|
|
if (power)
|
|
zfree(z);
|
|
z = temp;
|
|
power++;
|
|
}
|
|
zfree(*zp);
|
|
power <<= 1;
|
|
}
|
|
if (zrel(z, *zp) >= 0)
|
|
power++;
|
|
if (power > 1)
|
|
zfree(z);
|
|
return power;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the integral log base 10 of a number.
|
|
*
|
|
* If was_10_power != NULL, then this flag is set to true if the
|
|
* value was a power of 10, false otherwise.
|
|
*/
|
|
long
|
|
zlog10(ZVALUE z, bool *was_10_power)
|
|
{
|
|
ZVALUE *zp; /* current square */
|
|
long power; /* current power */
|
|
ZVALUE temp; /* temporary */
|
|
ZVALUE pow10; /* power of 10 */
|
|
FLAG rel; /* relationship */
|
|
int i;
|
|
|
|
/* NOTE: It is OK if was_10_power == NULL */
|
|
|
|
if (ziszero(z)) {
|
|
math_error("Zero argument argument for zlog10");
|
|
not_reached();
|
|
}
|
|
|
|
/* Ignore sign of z */
|
|
z.sign = 0;
|
|
|
|
/* preload power10 table if missing */
|
|
if (power10 == NULL) {
|
|
long v;
|
|
|
|
/* determine power10 table size */
|
|
for (v=1, max_power10_exp=0;
|
|
v <= (long)(MAXLONG/10L);
|
|
v *= 10L, ++max_power10_exp) {
|
|
}
|
|
|
|
/* create power10 table */
|
|
power10 = calloc(max_power10_exp+1, sizeof(long));
|
|
if (power10 == NULL) {
|
|
math_error("cannot malloc power10 table");
|
|
not_reached();
|
|
}
|
|
|
|
/* load power10 table */
|
|
for (i=0, v = 1L; i < max_power10_exp; ++i, v *= 10L) {
|
|
power10[i] = v;
|
|
}
|
|
}
|
|
|
|
/* assume not a power of ten unless we find out otherwise */
|
|
if (was_10_power != NULL) {
|
|
*was_10_power = false;
|
|
}
|
|
|
|
/* quick exit for small values */
|
|
if (! zgtmaxlong(z)) {
|
|
long value = ztolong(z);
|
|
|
|
for (i=0; i <= max_power10_exp; ++i) {
|
|
if (value == power10[i]) {
|
|
if (was_10_power != NULL) {
|
|
*was_10_power = true;
|
|
}
|
|
return i;
|
|
} else if (value < power10[i]) {
|
|
return i-1;
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Loop by squaring the base each time, and see whether or
|
|
* not each successive square is still smaller than the number.
|
|
*/
|
|
zp = &_tenpowers_[0];
|
|
*zp = _ten_;
|
|
while (((zp->len * 2) - 1) <= z.len) { /* while square not too large */
|
|
if (zp >= &_tenpowers_[TEN_MAX]) {
|
|
math_error("Maximum storable power of 10 reached!");
|
|
not_reached();
|
|
}
|
|
if (zp[1].len == 0)
|
|
zsquare(*zp, zp + 1);
|
|
zp++;
|
|
}
|
|
|
|
/*
|
|
* Now back down the squares, and multiply them together to see
|
|
* exactly how many times the base can be raised by.
|
|
*/
|
|
/* find the tenpower table entry < z */
|
|
do {
|
|
rel = zrel(*zp, z);
|
|
if (rel == 0) {
|
|
/* quick return - we match a tenpower entry */
|
|
if (was_10_power != NULL) {
|
|
*was_10_power = true;
|
|
}
|
|
return (1L << (zp - _tenpowers_));
|
|
}
|
|
} while (rel > 0 && --zp >= _tenpowers_);
|
|
if (zp < _tenpowers_) {
|
|
math_error("fell off bottom of tenpower table!");
|
|
not_reached();
|
|
}
|
|
|
|
/* the tenpower value is now our starting comparison value */
|
|
zcopy(*zp, &pow10);
|
|
power = (1L << (zp - _tenpowers_));
|
|
|
|
/* try to build up a power of 10 from tenpower table entries */
|
|
while (--zp >= _tenpowers_) {
|
|
|
|
/* try the next lower tenpower value */
|
|
zmul(pow10, *zp, &temp);
|
|
rel = zrel(temp, z);
|
|
if (rel == 0) {
|
|
/* exact power of 10 match */
|
|
power += (1L << (zp - _tenpowers_));
|
|
if (was_10_power != NULL) {
|
|
*was_10_power = true;
|
|
}
|
|
zfree(pow10);
|
|
zfree(temp);
|
|
return power;
|
|
|
|
/* ignore this entry if we went too large */
|
|
} else if (rel > 0) {
|
|
zfree(temp);
|
|
|
|
/* otherwise increase power and keep going */
|
|
} else {
|
|
power += (1L << (zp - _tenpowers_));
|
|
zfree(pow10);
|
|
pow10 = temp;
|
|
}
|
|
}
|
|
zfree(pow10);
|
|
return power;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the number of times that one number will divide another.
|
|
* This works similarly to zlog, except that divisions must be exact.
|
|
* For example, zdivcount(540, 3) = 3, since 3^3 divides 540 but 3^4 won't.
|
|
*/
|
|
long
|
|
zdivcount(ZVALUE z1, ZVALUE z2)
|
|
{
|
|
long count; /* number of factors removed */
|
|
ZVALUE tmp; /* ignored return value */
|
|
|
|
if (ziszero(z1) || ziszero(z2) || zisunit(z2))
|
|
return 0;
|
|
count = zfacrem(z1, z2, &tmp);
|
|
zfree(tmp);
|
|
return count;
|
|
}
|
|
|
|
|
|
/*
|
|
* Remove all occurrences of the specified factor from a number.
|
|
* Also returns the number of factors removed as a function return value.
|
|
* Example: zfacrem(540, 3, &x) returns 3 and sets x to 20.
|
|
*/
|
|
long
|
|
zfacrem(ZVALUE z1, ZVALUE z2, ZVALUE *rem)
|
|
{
|
|
register ZVALUE *zp; /* current square */
|
|
long count; /* total count of divisions */
|
|
long worth; /* worth of current square */
|
|
long lowbit; /* for zlowbit(z2) */
|
|
ZVALUE temp1, temp2, temp3; /* temporaries */
|
|
ZVALUE squares[32]; /* table of squares of factor */
|
|
|
|
/* firewall */
|
|
if (rem == NULL) {
|
|
math_error("%s: rem NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
z1.sign = 0;
|
|
z2.sign = 0;
|
|
/*
|
|
* Reject trivial cases.
|
|
*/
|
|
if ((z1.len < z2.len) || (zisodd(z1) && ziseven(z2)) ||
|
|
ziszero(z2) || zisone(z2) ||
|
|
((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))) {
|
|
rem->v = alloc(z1.len);
|
|
rem->len = z1.len;
|
|
rem->sign = 0;
|
|
zcopyval(z1, *rem);
|
|
return 0;
|
|
}
|
|
/*
|
|
* Handle any power of two special.
|
|
*/
|
|
if (zisonebit(z2)) {
|
|
lowbit = zlowbit(z2);
|
|
count = zlowbit(z1) / lowbit;
|
|
rem->v = alloc(z1.len);
|
|
rem->len = z1.len;
|
|
rem->sign = 0;
|
|
zcopyval(z1, *rem);
|
|
zshiftr(*rem, count * lowbit);
|
|
ztrim(rem);
|
|
return count;
|
|
}
|
|
/*
|
|
* See if the factor goes in even once.
|
|
*/
|
|
zdiv(z1, z2, &temp1, &temp2, 0);
|
|
if (!ziszero(temp2)) {
|
|
zfree(temp1);
|
|
zfree(temp2);
|
|
rem->v = alloc(z1.len);
|
|
rem->len = z1.len;
|
|
rem->sign = 0;
|
|
zcopyval(z1, *rem);
|
|
return 0;
|
|
}
|
|
zfree(temp2);
|
|
z1 = temp1;
|
|
/*
|
|
* Now loop by squaring the factor each time, and see whether
|
|
* or not each successive square will still divide the number.
|
|
*/
|
|
count = 1;
|
|
worth = 1;
|
|
zp = &squares[0];
|
|
*zp = z2;
|
|
while (((zp->len * 2) - 1) <= z1.len) { /* while square not too large */
|
|
zsquare(*zp, &temp1);
|
|
zdiv(z1, temp1, &temp2, &temp3, 0);
|
|
if (!ziszero(temp3)) {
|
|
zfree(temp1);
|
|
zfree(temp2);
|
|
zfree(temp3);
|
|
break;
|
|
}
|
|
zfree(temp3);
|
|
zfree(z1);
|
|
z1 = temp2;
|
|
*++zp = temp1;
|
|
worth *= 2;
|
|
count += worth;
|
|
}
|
|
|
|
/*
|
|
* Now back down the list of squares, and see if the lower powers
|
|
* will divide any more times.
|
|
*/
|
|
/*
|
|
* We prevent the zp pointer from walking behind squares
|
|
* by stopping one short of the end and running the loop one
|
|
* more time.
|
|
*
|
|
* We could stop the loop with just zp >= squares, but stopping
|
|
* short and running the loop one last time manually helps make
|
|
* code checkers such as insure happy.
|
|
*/
|
|
for (; zp > squares; zp--, worth /= 2) {
|
|
if (zp->len <= z1.len) {
|
|
zdiv(z1, *zp, &temp1, &temp2, 0);
|
|
if (ziszero(temp2)) {
|
|
temp3 = z1;
|
|
z1 = temp1;
|
|
temp1 = temp3;
|
|
count += worth;
|
|
}
|
|
zfree(temp1);
|
|
zfree(temp2);
|
|
}
|
|
if (zp != squares)
|
|
zfree(*zp);
|
|
}
|
|
/* run the loop manually one last time */
|
|
if (zp == squares) {
|
|
if (zp->len <= z1.len) {
|
|
zdiv(z1, *zp, &temp1, &temp2, 0);
|
|
if (ziszero(temp2)) {
|
|
temp3 = z1;
|
|
z1 = temp1;
|
|
temp1 = temp3;
|
|
count += worth;
|
|
}
|
|
zfree(temp1);
|
|
zfree(temp2);
|
|
}
|
|
if (zp != squares)
|
|
zfree(*zp);
|
|
}
|
|
|
|
*rem = z1;
|
|
return count;
|
|
}
|
|
|
|
|
|
/*
|
|
* Keep dividing a number by the gcd of it with another number until the
|
|
* result is relatively prime to the second number. Returns the number
|
|
* of divisions made, and if this is positive, stores result at res.
|
|
*/
|
|
long
|
|
zgcdrem(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
ZVALUE tmp1, tmp2, tmp3, tmp4;
|
|
long count, onecount;
|
|
long sh;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (ziszero(z1) || ziszero(z2)) {
|
|
math_error("Zero argument in call to zgcdrem!!!");
|
|
not_reached();
|
|
}
|
|
/*
|
|
* Begin by taking the gcd for the first time.
|
|
* If the number is already relatively prime, then we are done.
|
|
*/
|
|
z1.sign = 0;
|
|
z2.sign = 0;
|
|
if (zisone(z2))
|
|
return 0;
|
|
if (zisonebit(z2)) {
|
|
sh = zlowbit(z1);
|
|
if (sh == 0)
|
|
return 0;
|
|
zshift(z1, -sh, res);
|
|
return 1 + (sh - 1)/zlowbit(z2);
|
|
}
|
|
if (zisonebit(z1)) {
|
|
if (zisodd(z2))
|
|
return 0;
|
|
*res = _one_;
|
|
return zlowbit(z1);
|
|
}
|
|
|
|
zgcd(z1, z2, &tmp1);
|
|
if (zisunit(tmp1) || ziszero(tmp1)) {
|
|
zfree(tmp1);
|
|
return 0;
|
|
}
|
|
zequo(z1, tmp1, &tmp2);
|
|
count = 1;
|
|
z1 = tmp2;
|
|
z2 = tmp1;
|
|
/*
|
|
* Now keep alternately taking the gcd and removing factors until
|
|
* the gcd becomes one.
|
|
*/
|
|
while (!zisunit(z2)) {
|
|
onecount = zfacrem(z1, z2, &tmp3);
|
|
if (onecount) {
|
|
count += onecount;
|
|
zfree(z1);
|
|
z1 = tmp3;
|
|
} else {
|
|
zfree(tmp3);
|
|
}
|
|
zgcd(z1, z2, &tmp4);
|
|
zfree(z2);
|
|
z2 = tmp4;
|
|
}
|
|
zfree(z2);
|
|
*res = z1;
|
|
return count;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the number of digits (base 10) in a number, ignoring the sign.
|
|
*/
|
|
long
|
|
zdigits(ZVALUE z1)
|
|
{
|
|
long count, val;
|
|
|
|
z1.sign = 0;
|
|
if (!zge16b(z1)) { /* do small numbers ourself */
|
|
count = 1;
|
|
val = 10;
|
|
while (*z1.v >= (HALF)val) {
|
|
count++;
|
|
val *= 10;
|
|
}
|
|
return count;
|
|
}
|
|
return (zlog10(z1, NULL) + 1);
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the single digit at the specified decimal place of a number,
|
|
* where 0 means the rightmost digit. Example: zdigit(1234, 1) = 3.
|
|
*/
|
|
long
|
|
zdigit(ZVALUE z1, long n)
|
|
{
|
|
ZVALUE tmp1, tmp2;
|
|
long res;
|
|
|
|
z1.sign = 0;
|
|
if (ziszero(z1) || (n < 0) || (n / BASEDIG >= z1.len))
|
|
return 0;
|
|
if (n == 0)
|
|
return zmodi(z1, 10L);
|
|
if (n == 1)
|
|
return zmodi(z1, 100L) / 10;
|
|
if (n == 2)
|
|
return zmodi(z1, 1000L) / 100;
|
|
if (n == 3)
|
|
return zmodi(z1, 10000L) / 1000;
|
|
ztenpow(n, &tmp1);
|
|
zquo(z1, tmp1, &tmp2, 0);
|
|
res = zmodi(tmp2, 10L);
|
|
zfree(tmp1);
|
|
zfree(tmp2);
|
|
return res;
|
|
}
|
|
|
|
|
|
/*
|
|
* z is to be a nonnegative integer
|
|
* If z is the square of a integer stores at dest the square root of z;
|
|
* otherwise stores at z an integer differing from the square root
|
|
* by less than 1. Returns the sign of the true square root minus
|
|
* the calculated integer. Type of rounding is determined by
|
|
* rnd as follows: rnd = 0 gives round down, rnd = 1
|
|
* rounds up, rnd = 8 rounds to even integer, rnd = 9 rounds to odd
|
|
* integer, rnd = 16 rounds to nearest integer.
|
|
*/
|
|
FLAG
|
|
zsqrt(ZVALUE z, ZVALUE *dest, long rnd)
|
|
{
|
|
HALF *a, *A, *b, *a0, u;
|
|
int i, j, j1, j2, k, k1, m, m0, m1, n, n0, o;
|
|
FULL d, e, f, g, h, s, t, x, topbit;
|
|
int remsign;
|
|
bool up, onebit;
|
|
ZVALUE sqrt;
|
|
|
|
/* firewall */
|
|
if (dest == NULL) {
|
|
math_error("%s: dest NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (z.sign) {
|
|
math_error("Square root of negative number");
|
|
not_reached();
|
|
}
|
|
if (ziszero(z)) {
|
|
*dest = _zero_;
|
|
return 0;
|
|
}
|
|
m0 = z.len;
|
|
o = m0 & 1;
|
|
m = m0 + o; /* m is smallest even number >= z.len */
|
|
n0 = n = m / 2;
|
|
f = z.v[z.len - 1];
|
|
k = 1;
|
|
while (f >>= 2)
|
|
k++;
|
|
if (!o)
|
|
k += BASEB/2;
|
|
j = BASEB - k;
|
|
m1 = m;
|
|
if (k == BASEB) {
|
|
m1 += 2;
|
|
n0++;
|
|
}
|
|
A = alloc(m1);
|
|
A[m1] = 0;
|
|
a0 = A + n0;
|
|
memcpy(A, z.v, m0 * sizeof(HALF));
|
|
if (o)
|
|
A[m - 1] = 0;
|
|
if (n == 1) {
|
|
if (j)
|
|
f = (FULL) A[1] << j | A[0] >> k;
|
|
else
|
|
f = A[1];
|
|
g = (FULL) A[0] << (j + BASEB);
|
|
d = e = topbit = (FULL)1 << (k - 1);
|
|
} else {
|
|
if (j)
|
|
f = (FULL) A[m-1] << (j + BASEB) | (FULL) A[m-2] << j |
|
|
A[m-3] >> k;
|
|
else
|
|
f = (FULL) A[m-1] << BASEB | A[m-2];
|
|
g = (FULL) A[m-3] << (j + BASEB) | (FULL) A[m-4] << j;
|
|
d = e = topbit = (FULL)1 << (BASEB + k - 1);
|
|
}
|
|
|
|
s = (f & topbit);
|
|
f <<= 1;
|
|
if (g & TOPFULL)
|
|
f++;
|
|
g <<= 1;
|
|
if (s) {
|
|
f -= 4 * d;
|
|
e = 2 * d - 1;
|
|
}
|
|
else
|
|
f -= d;
|
|
while (d >>= 1) {
|
|
if (!(s | f | g))
|
|
break;
|
|
while (d && (f & topbit) == s) {
|
|
d >>= 1;
|
|
f <<= 1;
|
|
if (g & TOPFULL)
|
|
f++;
|
|
g <<= 1;
|
|
}
|
|
if (d == 0)
|
|
break;
|
|
if (s)
|
|
f += e + 1;
|
|
else
|
|
f -= e;
|
|
t = f & topbit;
|
|
f <<= 1;
|
|
if (g & TOPFULL)
|
|
f++;
|
|
g <<= 1;
|
|
if (t == 0 && f < d)
|
|
t = topbit;
|
|
f -= d;
|
|
if (s)
|
|
e -= d - !t;
|
|
else
|
|
e += d - (t > 0);
|
|
s = t;
|
|
}
|
|
if (n0 == 1) {
|
|
A[1] = (HALF)e;
|
|
A[0] = (HALF)f;
|
|
m = 1;
|
|
goto done;
|
|
}
|
|
if (n0 == 2) {
|
|
A[3] = (HALF)(e >> BASEB);
|
|
A[2] = (HALF)e;
|
|
A[1] = (HALF)(f >> BASEB);
|
|
A[0] = (HALF)f;
|
|
m = 2;
|
|
goto done;
|
|
}
|
|
u = (HALF)(s ? BASE1 : 0);
|
|
if (k < BASEB) {
|
|
A[m1 - 1] = (HALF)(e >> (BASEB - 1));
|
|
A[m1 - 2] = ((HALF)(e << 1) | (HALF)(s > 0));
|
|
A[m1 - 3] = (HALF)(f >> BASEB);
|
|
A[m1 - 4] = (HALF)f;
|
|
m = m1 - 2;
|
|
k1 = k + 1;
|
|
} else {
|
|
A[m1 - 1] = 1;
|
|
A[m1 - 2] = (HALF)(e >> (BASEB - 1));
|
|
A[m1 - 3] = ((HALF)(e << 1) | (HALF)(s > 0));
|
|
A[m1 - 4] = u;
|
|
A[m1 - 5] = (HALF)(f >> BASEB);
|
|
A[m1 - 6] = (HALF)f;
|
|
m = m1 - 3;
|
|
k1 = 1;
|
|
}
|
|
h = e >> k;
|
|
onebit = ((e & ((FULL)1 << (k - 1))) ? true : false);
|
|
j2 = BASEB - k1;
|
|
j1 = BASEB + j2;
|
|
while (m > n0) {
|
|
a = A + m - 1;
|
|
if (j2)
|
|
f = (FULL) *a << j1 | (FULL) a[-1] << j2 | a[-2] >> k1;
|
|
else
|
|
f = (FULL) *a << BASEB | a[-1];
|
|
if (u)
|
|
f = ~f;
|
|
x = f / h;
|
|
if (x) {
|
|
if (onebit && x > 2 * (f % h) + 2)
|
|
x--;
|
|
b = a + 1;
|
|
i = m1 - m;
|
|
a -= i + 1;
|
|
if (u) {
|
|
f = *a + x * (BASE - x);
|
|
*a++ = (HALF)f;
|
|
u = (HALF)(f >> BASEB);
|
|
while (i--) {
|
|
f = *a + x * *b++ + u;
|
|
*a++ = (HALF)f;
|
|
u = (HALF)(f >> BASEB);
|
|
}
|
|
u += *a;
|
|
x = ~x + !u;
|
|
if (!(x & TOPHALF))
|
|
a[1] -= 1;
|
|
} else {
|
|
f = *a - x * x;
|
|
*a++ = (HALF)f;
|
|
u = -(HALF)(f >> BASEB);
|
|
while (i--) {
|
|
f = *a - x * *b++ - u;
|
|
*a++ = (HALF)f;
|
|
u = -(HALF)(f >> BASEB);
|
|
}
|
|
u = *a - u;
|
|
x = x + u;
|
|
if (x & TOPHALF)
|
|
a[1] |= 1;
|
|
}
|
|
*a = ((HALF)(x << 1) | (HALF)(u > 0));
|
|
} else {
|
|
*a = u;
|
|
}
|
|
m--;
|
|
if (*--a == u) {
|
|
while (m > 1 && *--a == u)
|
|
m--;
|
|
}
|
|
}
|
|
i = n;
|
|
a = a0;
|
|
while (i--) {
|
|
*a >>= 1;
|
|
if (a[1] & 1) *a |= TOPHALF;
|
|
a++;
|
|
}
|
|
s = u;
|
|
done: if (s == 0) {
|
|
while (m > 0 && A[m - 1] == 0)
|
|
m--;
|
|
if (m == 0) {
|
|
remsign = 0;
|
|
sqrt.v = alloc(n);
|
|
sqrt.len = n;
|
|
sqrt.sign = 0;
|
|
memcpy(sqrt.v, a0, n * sizeof(HALF));
|
|
freeh(A);
|
|
*dest = sqrt;
|
|
return remsign;
|
|
}
|
|
}
|
|
if (rnd & 16) {
|
|
if (s == 0) {
|
|
if (m != n) {
|
|
up = (m > n);
|
|
} else {
|
|
i = n;
|
|
b = a0 + n;
|
|
a = A + n;
|
|
while (i > 0 && *--a == *--b)
|
|
i--;
|
|
up = (i > 0 && *a > *b);
|
|
}
|
|
} else {
|
|
while (m > 1 && A[m - 1] == BASE1)
|
|
m--;
|
|
if (m != n) {
|
|
up = (m < n);
|
|
} else {
|
|
i = n;
|
|
b = a0 + n;
|
|
a = A + n;
|
|
while (i > 0 && *--a + *--b == BASE1)
|
|
i--;
|
|
up = ((FULL) *a + *b >= BASE);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
if (rnd & 8)
|
|
up = (((rnd ^ *a0) & 1) ? true : false);
|
|
else
|
|
up = ((rnd & 1) ? true : false);
|
|
if (up) {
|
|
remsign = -1;
|
|
i = n;
|
|
a = a0;
|
|
while (i-- && *a == BASE1)
|
|
*a++ = 0;
|
|
if (i >= 0) {
|
|
(*a)++;
|
|
} else {
|
|
n++;
|
|
*a = 1;
|
|
}
|
|
} else {
|
|
remsign = 1;
|
|
}
|
|
sqrt.v = alloc(n);
|
|
sqrt.len = n;
|
|
sqrt.sign = 0;
|
|
memcpy(sqrt.v, a0, n * sizeof(HALF));
|
|
freeh(A);
|
|
*dest = sqrt;
|
|
return remsign;
|
|
|
|
}
|
|
|
|
/*
|
|
* Take an arbitrary root of a number (to the greatest integer).
|
|
* This uses the following iteration to get the K-th root of N:
|
|
* x = ((K-1) * x + N / x^(K-1)) / K
|
|
*/
|
|
void
|
|
zroot(ZVALUE z1, ZVALUE z2, ZVALUE *dest)
|
|
{
|
|
ZVALUE ztry, quo, old, temp, temp2;
|
|
ZVALUE k1; /* holds k - 1 */
|
|
int sign;
|
|
long i;
|
|
LEN highbit, k;
|
|
SIUNION sival;
|
|
|
|
/* firewall */
|
|
if (dest == NULL) {
|
|
math_error("%s: dest NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
sign = z1.sign;
|
|
if (sign && ziseven(z2)) {
|
|
math_error("Even root of negative number");
|
|
not_reached();
|
|
}
|
|
if (ziszero(z2) || zisneg(z2)) {
|
|
math_error("Non-positive root");
|
|
not_reached();
|
|
}
|
|
if (ziszero(z1)) { /* root of zero */
|
|
*dest = _zero_;
|
|
return;
|
|
}
|
|
if (zisunit(z2)) { /* first root */
|
|
zcopy(z1, dest);
|
|
return;
|
|
}
|
|
if (zge31b(z2)) { /* humongous root */
|
|
*dest = _one_;
|
|
dest->sign = (bool)((HALF)sign);
|
|
return;
|
|
}
|
|
k = (LEN)ztolong(z2);
|
|
highbit = zhighbit(z1);
|
|
if (highbit < k) { /* too high a root */
|
|
*dest = _one_;
|
|
dest->sign = (bool)((HALF)sign);
|
|
return;
|
|
}
|
|
sival.ivalue = k - 1;
|
|
k1.v = &sival.silow;
|
|
/* ignore Saber-C warning #112 - get ushort from uint */
|
|
/* OK to ignore on name zroot`sival */
|
|
k1.len = 1 + (sival.sihigh != 0);
|
|
k1.sign = 0;
|
|
z1.sign = 0;
|
|
/*
|
|
* Allocate the numbers to use for the main loop.
|
|
* The size and high bits of the final result are correctly set here.
|
|
* Notice that the remainder of the test value is rubbish, but this
|
|
* is unimportant.
|
|
*/
|
|
highbit = (highbit + k - 1) / k;
|
|
ztry.len = (highbit / BASEB) + 1;
|
|
ztry.v = alloc(ztry.len);
|
|
zclearval(ztry);
|
|
ztry.v[ztry.len-1] = ((HALF)1 << (highbit % BASEB));
|
|
ztry.sign = 0;
|
|
old.v = alloc(ztry.len);
|
|
old.len = 1;
|
|
zclearval(old);
|
|
old.sign = 0;
|
|
/*
|
|
* Main divide and average loop
|
|
*/
|
|
for (;;) {
|
|
zpowi(ztry, k1, &temp);
|
|
zquo(z1, temp, &quo, 0);
|
|
zfree(temp);
|
|
i = zrel(ztry, quo);
|
|
if (i <= 0) {
|
|
/*
|
|
* Current try is less than or equal to the root since
|
|
* it is less than the quotient. If the quotient is
|
|
* equal to the try, we are all done. Also, if the
|
|
* try is equal to the old value, we are done since
|
|
* no improvement occurred. If not, save the improved
|
|
* value and loop some more.
|
|
*/
|
|
if ((i == 0) || (zcmp(old, ztry) == 0)) {
|
|
zfree(quo);
|
|
zfree(old);
|
|
ztry.sign = (bool)((HALF)sign);
|
|
zquicktrim(ztry);
|
|
*dest = ztry;
|
|
return;
|
|
}
|
|
old.len = ztry.len;
|
|
zcopyval(ztry, old);
|
|
}
|
|
/* average current try and quotient for the new try */
|
|
zmul(ztry, k1, &temp);
|
|
zfree(ztry);
|
|
zadd(quo, temp, &temp2);
|
|
zfree(temp);
|
|
zfree(quo);
|
|
zquo(temp2, z2, &ztry, 0);
|
|
zfree(temp2);
|
|
}
|
|
}
|
|
|
|
|
|
/*
|
|
* Test to see if a number is an exact square or not.
|
|
*/
|
|
bool
|
|
zissquare(ZVALUE z)
|
|
{
|
|
long n;
|
|
ZVALUE tmp;
|
|
|
|
/* negative values are never perfect squares */
|
|
if (zisneg(z)) {
|
|
return false;
|
|
}
|
|
|
|
/* ignore trailing zero words */
|
|
while ((z.len > 1) && (*z.v == 0)) {
|
|
z.len--;
|
|
z.v++;
|
|
}
|
|
|
|
/* zero or one is a perfect square */
|
|
if (zisabsleone(z)) {
|
|
return true;
|
|
}
|
|
|
|
/* check mod 4096 values */
|
|
if (issq_mod4k[(int)(*z.v & 0xfff)] == 0) {
|
|
return false;
|
|
}
|
|
|
|
/* must do full square root test now */
|
|
n = !zsqrt(z, &tmp, 0);
|
|
zfree(tmp);
|
|
return (n ? true : false);
|
|
}
|
|
|
|
|
|
/*
|
|
* test if a number is an integer power of 2
|
|
*
|
|
* given:
|
|
* z value to check if it is a power of 2
|
|
* log2 when z is an integer power of 2 (true return), *log2 is set to log base 2 of z
|
|
* when z is NOT an integer power of 2 (false return), *log2 is not touched
|
|
*
|
|
* returns:
|
|
* true z is a power of 2
|
|
* false z is not a power of 2
|
|
*/
|
|
bool
|
|
zispowerof2(ZVALUE z, FULL *log2)
|
|
{
|
|
FULL ilogz; /* potential log base 2 return value or -1 */
|
|
HALF tophalf; /* most significant HALF in z */
|
|
LEN len; /* length of z in HALFs */
|
|
int i;
|
|
|
|
/* firewall */
|
|
if (log2 == NULL) {
|
|
math_error("%s: log2 NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
/* zero and negative values are never integer powers of 2 */
|
|
if (ziszero(z) || zisneg(z)) {
|
|
return false;
|
|
}
|
|
|
|
/*
|
|
* trim z just in case
|
|
*
|
|
* An untrimmed z will give incorrect results.
|
|
*/
|
|
ztrim(&z);
|
|
|
|
/*
|
|
* all HALFs below the top HALF must be zero
|
|
*/
|
|
len = z.len;
|
|
for (i=0, ilogz=0; i < len-1; ++i, ilogz+=BASEB) {
|
|
if (z.v[i] != 0) {
|
|
return false;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* top HALF must be a power of 2
|
|
*
|
|
* For non-zero values of tophalf,
|
|
* (tophalf & (tophalf-1)) == 0 ==> tophalf is a power of 2,
|
|
* (tophalf & (tophalf-1)) != 0 ==> tophalf is NOT a power of 2.
|
|
*/
|
|
tophalf = z.v[len-1];
|
|
if ((tophalf == 0) || ((tophalf & (tophalf-1)) != 0)) {
|
|
return false;
|
|
}
|
|
|
|
/*
|
|
* count the bits in the top HALF which is a power of 2
|
|
*/
|
|
for (i=0; i < BASEB && tophalf != ((HALF)1 << i); ++i) {
|
|
++ilogz;
|
|
}
|
|
|
|
/*
|
|
* return power of 2
|
|
*/
|
|
*log2 = ilogz;
|
|
return true;
|
|
}
|