c_cpp/etc/calc/zfunc.c

   1 /*
   2  * zfunc - extended precision integral arithmetic non-primitive routines
   3  *
   4  * Copyright (C) 1999-2007  David I. Bell, Landon Curt Noll and Ernest Bowen
   5  *
   6  * Primary author:  David I. Bell
   7  *
   8  * Calc is open software; you can redistribute it and/or modify it under
   9  * the terms of the version 2.1 of the GNU Lesser General Public License
  10  * as published by the Free Software Foundation.
  11  *
  12  * Calc is distributed in the hope that it will be useful, but WITHOUT
  13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  15  * Public License for more details.
  16  *
  17  * A copy of version 2.1 of the GNU Lesser General Public License is
  18  * distributed with calc under the filename COPYING-LGPL.  You should have
  19  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  20  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  21  *
  22  * @(#) $Revision: 30.3 $
  23  * @(#) $Id: zfunc.c,v 30.3 2013/08/11 08:41:38 chongo Exp $
  24  * @(#) $Source: /usr/local/src/bin/calc/RCS/zfunc.c,v $
  25  *
  26  * Under source code control:   1990/02/15 01:48:27
  27  * File existed as early as:    before 1990
  28  *
  29  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  30  */
  31
  32
  33 #include "zmath.h"
  34
  35 ZVALUE _tenpowers_[TEN_MAX+1];          /* table of 10^2^n */
  36
  37 STATIC long *power10 = NULL;
  38 STATIC int max_power10_exp = 0;
  39
  40 /*
  41  * given:
  42  *
  43  *      unsigned long x
  44  * or:  unsigned long long x
  45  * or:  long x                  and  x >= 0
  46  * or:  long long x             and  x >= 0
  47  *
  48  * If issq_mod4k[x & 0xfff] == 0, then x cannot be a perfect square
  49  * else x might be a perfect square.
  50  */
  51 STATIC USB8 issq_mod4k[1<<12] = {
  52         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,
  53         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,
  54         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,
  55         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,
  56         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,
  57         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,
  58         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,
  59         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,
  60         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,
  61         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,
  62         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,
  63         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,
  64         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,
  65         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,
  66         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,
  67         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,
  68         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,
  69         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,
  70         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,
  71         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,
  72         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,
  73         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,
  74         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,
  75         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,
  76         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,
  77         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,
  78         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,
  79         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,
  80         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,
  81         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,
  82         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,
  83         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,
  84         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,
  85         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,
  86         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,
  87         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,
  88         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,
  89         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,
  90         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,
  91         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,
  92         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,
  93         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,
  94         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,
  95         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,
  96         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,
  97         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,
  98         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,
  99         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,
 100         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,
 101         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,
 102         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,
 103         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,
 104         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,
 105         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,
 106         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,
 107         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,
 108         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,
 109         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,
 110         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,
 111         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,
 112         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,
 113         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,
 114         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,
 115         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,
 116         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,
 117         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,
 118         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,
 119         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,
 120         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,
 121         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,
 122         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,
 123         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,
 124         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,
 125         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,
 126         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,
 127         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,
 128         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,
 129         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,
 130         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,
 131         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,
 132         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,
 133         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,
 134         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,
 135         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,
 136         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,
 137         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,
 138         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,
 139         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,
 140         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,
 141         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,
 142         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,
 143         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,
 144         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,
 145         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,
 146         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,
 147         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,
 148         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,
 149         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,
 150         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,
 151         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,
 152         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,
 153         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,
 154         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,
 155         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,
 156         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,
 157         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,
 158         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,
 159         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,
 160         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,
 161         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,
 162         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,
 163         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,
 164         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,
 165         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,
 166         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,
 167         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,
 168         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,
 169         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,
 170         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,
 171         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,
 172         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,
 173         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,
 174         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,
 175         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,
 176         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,
 177         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,
 178         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,
 179         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,
 180 };
 181
 182
 183 /*
 184  * Compute the factorial of a number.
 185  */
 186 void
 187 zfact(ZVALUE z, ZVALUE *dest)
 188 {
 189         long ptwo;              /* count of powers of two */
 190         long n;                 /* current multiplication value */
 191         long m;                 /* reduced multiplication value */
 192         long mul;               /* collected value to multiply by */
 193         ZVALUE res, temp;
 194
 195         if (zisneg(z)) {
 196                 math_error("Negative argument for factorial");
 197                 /*NOTREACHED*/
 198         }
 199         if (zge31b(z)) {
 200                 math_error("Very large factorial");
 201                 /*NOTREACHED*/
 202         }
 203         n = ztolong(z);
 204         ptwo = 0;
 205         mul = 1;
 206         res = _one_;
 207         /*
 208          * Multiply numbers together, but squeeze out all powers of two.
 209          * We will put them back in at the end.  Also collect multiple
 210          * numbers together until there is a risk of overflow.
 211          */
 212         for (; n > 1; n--) {
 213                 for (m = n; ((m & 0x1) == 0); m >>= 1)
 214                         ptwo++;
 215                 if (mul <= MAXLONG/m) {
 216                         mul *= m;
 217                         continue;
 218                 }
 219                 zmuli(res, mul, &temp);
 220                 zfree(res);
 221                 res = temp;
 222                 mul = m;
 223         }
 224         /*
 225          * Multiply by the remaining value, then scale result by
 226          * the proper power of two.
 227          */
 228         if (mul > 1) {
 229                 zmuli(res, mul, &temp);
 230                 zfree(res);
 231                 res = temp;
 232         }
 233         zshift(res, ptwo, &temp);
 234         zfree(res);
 235         *dest = temp;
 236 }
 237
 238
 239 /*
 240  * Compute the permutation function  M! / (M - N)!.
 241  */
 242 void
 243 zperm(ZVALUE z1, ZVALUE z2, ZVALUE *res)
 244 {
 245         SFULL count;
 246         ZVALUE cur, tmp, ans;
 247
 248         if (zisneg(z1) || zisneg(z2)) {
 249                 math_error("Negative argument for permutation");
 250                 /*NOTREACHED*/
 251         }
 252         if (zrel(z1, z2) < 0) {
 253                 math_error("Second arg larger than first in permutation");
 254                 /*NOTREACHED*/
 255         }
 256         if (zge31b(z2)) {
 257                 math_error("Very large permutation");
 258                 /*NOTREACHED*/
 259         }
 260         count = ztolong(z2);
 261         zcopy(z1, &ans);
 262         zsub(z1, _one_, &cur);
 263         while (--count > 0) {
 264                 zmul(ans, cur, &tmp);
 265                 zfree(ans);
 266                 ans = tmp;
 267                 zsub(cur, _one_, &tmp);
 268                 zfree(cur);
 269                 cur = tmp;
 270         }
 271         zfree(cur);
 272         *res = ans;
 273 }
 274
 275 /*
 276  * docomb evaluates binomial coefficient when z1 >= 0, z2 >= 0
 277  */
 278 S_FUNC int
 279 docomb(ZVALUE z1, ZVALUE z2, ZVALUE *res)
 280 {
 281         ZVALUE ans;
 282         ZVALUE mul, div, temp;
 283         FULL count, i;
 284 #if BASEB == 16
 285         HALF dh[2];
 286 #else
 287         HALF dh[1];
 288 #endif
 289
 290         if (zrel(z2, z1) > 0)
 291                 return 0;
 292         zsub(z1, z2, &temp);
 293
 294         if (zge31b(z2) && zge31b(temp)) {
 295                 zfree(temp);
 296                 return -2;
 297         }
 298         if (zrel(temp, z2) < 0)
 299                 count = ztofull(temp);
 300         else
 301                 count = ztofull(z2);
 302         zfree(temp);
 303         if (count == 0)
 304                 return 1;
 305         if (count == 1)
 306                 return 2;
 307         div.sign = 0;
 308         div.v = dh;
 309         div.len = 1;
 310         zcopy(z1, &mul);
 311         zcopy(z1, &ans);
 312         for (i = 2; i <= count; i++) {
 313 #if BASEB == 16
 314                 dh[0] = (HALF)(i & BASE1);
 315                 dh[1] = (HALF)(i >> BASEB);
 316                 div.len = 1 + (dh[1] != 0);
 317 #else
 318                 dh[0] = (HALF) i;
 319 #endif
 320                 zsub(mul, _one_, &temp);
 321                 zfree(mul);
 322                 mul = temp;
 323                 zmul(ans, mul, &temp);
 324                 zfree(ans);
 325                 zquo(temp, div, &ans, 0);
 326                 zfree(temp);
 327         }
 328         zfree(mul);
 329         *res = ans;
 330         return 3;
 331 }
 332
 333 /*
 334  * Compute the combinatorial function  M! / ( N! * (M - N)! ).
 335  * Returns 0 if result is 0
 336 *          1      1
 337 *          2      z1
 338 *         -1      -1
 339 *         -2    if too complicated
 340 *          3    result stored at res
 341  */
 342 int
 343 zcomb(ZVALUE z1, ZVALUE z2, ZVALUE *res)
 344 {
 345         ZVALUE z3, z4;
 346         int r;
 347
 348         if (z2.sign || (!z1.sign && zrel(z2, z1) > 0))
 349                 return 0;
 350         if (zisone(z2))
 351                 return 2;
 352         if (z1.sign) {
 353                 z1.sign = 0;
 354                 zsub(z1, _one_, &z3);
 355                 zadd(z3, z2, &z4);
 356                 zfree(z3);
 357                 r = docomb(z4, z2, res);
 358                 if (r == 2) {
 359                         *res = z4;
 360                         r = 3;
 361                 }
 362                 else
 363                         zfree(z4);
 364                 if (z2.v[0] & 1) {
 365                         if (r == 1)
 366                                 r = -1;
 367                         if (r == 3)
 368                                 res->sign = 1;
 369                 }
 370                 return r;
 371         }
 372         return docomb(z1, z2, res);
 373 }
 374
 375
 376 /*
 377  * Compute the Jacobi function (m / n) for odd n.
 378  *
 379  * The property of the Jacobi function is: If n>2 is prime then
 380  *
 381  *      the result is 1 if m == x^2 (mod n) for some x.
 382  *      otherwise the result is -1.
 383  *
 384  * If n is not prime, then the result does not prove that n is not prine
 385  * when the value of the Jacobi is 1.
 386  *
 387  * Jacobi evaluation of (m / n), where n > 0 is odd AND m > 0 is odd:
 388  *
 389  *      rule 0: (0 / n) == 0
 390  *      rule 1: (1 / n) == 1
 391  *      rule 2: (m / n) == (a / n)      if m == a % n
 392  *      rule 3: (m / n) == (2*m / n)    if n == 1 % 8 OR n == 7 % 8
 393  *      rule 4: (m / n) == -(2*m / n)   if n != 1 & 8 AND n != 7 % 8
 394  *      rule 5: (m / n) == (n / m)      if m == 3 % 4 AND n == 3 % 4
 395  *      rule 6: (m / n) == -(n / m)     if m != 3 % 4 OR n != 3 % 4
 396  *
 397  * NOTE: This function returns 0 in invalid Jacobi parameters:
 398  *       m < 0 OR n is even OR n < 1.
 399  */
 400 FLAG
 401 zjacobi(ZVALUE z1, ZVALUE z2)
 402 {
 403         ZVALUE p, q, tmp;
 404         long lowbit;
 405         int val;
 406
 407         /* firewall */
 408         if (ziszero(z1) || zisneg(z1))
 409                 return 0;
 410         if (ziseven(z2) || zisneg(z2))
 411                 return 0;
 412
 413         /* assume a value of 1 unless we find otherwise */
 414         if (zisone(z1))
 415                 return 1;
 416         val = 1;
 417         zcopy(z1, &p);
 418         zcopy(z2, &q);
 419         for (;;) {
 420                 zmod(p, q, &tmp, 0);
 421                 zfree(p);
 422                 p = tmp;
 423                 if (ziszero(p)) {
 424                         zfree(p);
 425                         zfree(q);
 426                         return 0;
 427                 }
 428                 if (ziseven(p)) {
 429                         lowbit = zlowbit(p);
 430                         zshift(p, -lowbit, &tmp);
 431                         zfree(p);
 432                         p = tmp;
 433                         if ((lowbit & 1) && (((*q.v & 0x7) == 3) ||
 434                             ((*q.v & 0x7) == 5)))
 435                                 val = -val;
 436                 }
 437                 if (zisunit(p)) {
 438                         zfree(p);
 439                         zfree(q);
 440                         return val;
 441                 }
 442                 if ((*p.v & *q.v & 0x3) == 3)
 443                         val = -val;
 444                 tmp = q;
 445                 q = p;
 446                 p = tmp;
 447         }
 448 }
 449
 450
 451 /*
 452  * Return the Fibonacci number F(n).
 453  * This is evaluated by recursively using the formulas:
 454  *      F(2N+1) = F(N+1)^2 + F(N)^2
 455  * and
 456  *      F(2N) = F(N+1)^2 - F(N-1)^2
 457  */
 458 void
 459 zfib(ZVALUE z, ZVALUE *res)
 460 {
 461         long n;
 462         int sign;
 463         ZVALUE fnm1, fn, fnp1;          /* consecutive fibonacci values */
 464         ZVALUE t1, t2, t3;
 465         FULL i;
 466
 467         if (zge31b(z)) {
 468                 math_error("Very large Fibonacci number");
 469                 /*NOTREACHED*/
 470         }
 471         n = ztolong(z);
 472         if (n == 0) {
 473                 *res = _zero_;
 474                 return;
 475         }
 476         sign = z.sign && ((n & 0x1) == 0);
 477         if (n <= 2) {
 478                 *res = _one_;
 479                 res->sign = (BOOL)sign;
 480                 return;
 481         }
 482         i = TOPFULL;
 483         while ((i & n) == 0)
 484                 i >>= (FULL)1;
 485         i >>= (FULL)1;
 486         fnm1 = _zero_;
 487         fn = _one_;
 488         fnp1 = _one_;
 489         while (i) {
 490                 zsquare(fnm1, &t1);
 491                 zsquare(fn, &t2);
 492                 zsquare(fnp1, &t3);
 493                 zfree(fnm1);
 494                 zfree(fn);
 495                 zfree(fnp1);
 496                 zadd(t2, t3, &fnp1);
 497                 zsub(t3, t1, &fn);
 498                 zfree(t1);
 499                 zfree(t2);
 500                 zfree(t3);
 501                 if (i & n) {
 502                         fnm1 = fn;
 503                         fn = fnp1;
 504                         zadd(fnm1, fn, &fnp1);
 505                 } else {
 506                         zsub(fnp1, fn, &fnm1);
 507                 }
 508                 i >>= (FULL)1;
 509         }
 510         zfree(fnm1);
 511         zfree(fnp1);
 512         *res = fn;
 513         res->sign = (BOOL)sign;
 514 }
 515
 516
 517 /*
 518  * Compute the result of raising one number to the power of another
 519  * The second number is assumed to be non-negative.
 520  * It cannot be too large except for trivial cases.
 521  */
 522 void
 523 zpowi(ZVALUE z1, ZVALUE z2, ZVALUE *res)
 524 {
 525         int sign;               /* final sign of number */
 526         unsigned long power;    /* power to raise to */
 527         FULL bit;               /* current bit value */
 528         long twos;              /* count of times 2 is in result */
 529         ZVALUE ans, temp;
 530
 531         sign = (z1.sign && zisodd(z2));
 532         z1.sign = 0;
 533         z2.sign = 0;
 534         if (ziszero(z2) && !ziszero(z1)) {      /* number raised to power 0 */
 535                 *res = _one_;
 536                 return;
 537         }
 538         if (zisabsleone(z1)) {  /* 0, 1, or -1 raised to a power */
 539                 ans = _one_;
 540                 ans.sign = (BOOL)sign;
 541                 if (*z1.v == 0)
 542                         ans = _zero_;
 543                 *res = ans;
 544                 return;
 545         }
 546         if (zge31b(z2)) {
 547                 math_error("Raising to very large power");
 548                 /*NOTREACHED*/
 549         }
 550         power = ztoulong(z2);
 551         if (zistwo(z1)) {       /* two raised to a power */
 552                 zbitvalue((long) power, res);
 553                 return;
 554         }
 555         /*
 556          * See if this is a power of ten
 557          */
 558         if (zistiny(z1) && (*z1.v == 10)) {
 559                 ztenpow((long) power, res);
 560                 res->sign = (BOOL)sign;
 561                 return;
 562         }
 563         /*
 564          * Handle low powers specially
 565          */
 566         if (power <= 4) {
 567                 switch ((int) power) {
 568                 case 1:
 569                         ans.len = z1.len;
 570                         ans.v = alloc(ans.len);
 571                         zcopyval(z1, ans);
 572                         ans.sign = (BOOL)sign;
 573                         *res = ans;
 574                         return;
 575                 case 2:
 576                         zsquare(z1, res);
 577                         return;
 578                 case 3:
 579                         zsquare(z1, &temp);
 580                         zmul(z1, temp, res);
 581                         zfree(temp);
 582                         res->sign = (BOOL)sign;
 583                         return;
 584                 case 4:
 585                         zsquare(z1, &temp);
 586                         zsquare(temp, res);
 587                         zfree(temp);
 588                         return;
 589                 }
 590         }
 591         /*
 592          * Shift out all powers of twos so the multiplies are smaller.
 593          * We will shift back the right amount when done.
 594          */
 595         twos = 0;
 596         if (ziseven(z1)) {
 597                 twos = zlowbit(z1);
 598                 ans.v = alloc(z1.len);
 599                 ans.len = z1.len;
 600                 ans.sign = z1.sign;
 601                 zcopyval(z1, ans);
 602                 zshiftr(ans, twos);
 603                 ztrim(&ans);
 604                 z1 = ans;
 605                 twos *= power;
 606         }
 607         /*
 608          * Compute the power by squaring and multiplying.
 609          * This uses the left to right method of power raising.
 610          */
 611         bit = TOPFULL;
 612         while ((bit & power) == 0)
 613                 bit >>= 1;
 614         bit >>= 1;
 615         zsquare(z1, &ans);
 616         if (bit & power) {
 617                 zmul(ans, z1, &temp);
 618                 zfree(ans);
 619                 ans = temp;
 620         }
 621         bit >>= 1;
 622         while (bit) {
 623                 zsquare(ans, &temp);
 624                 zfree(ans);
 625                 ans = temp;
 626                 if (bit & power) {
 627                         zmul(ans, z1, &temp);
 628                         zfree(ans);
 629                         ans = temp;
 630                 }
 631                 bit >>= 1;
 632         }
 633         /*
 634          * Scale back up by proper power of two
 635          */
 636         if (twos) {
 637                 zshift(ans, twos, &temp);
 638                 zfree(ans);
 639                 ans = temp;
 640                 zfree(z1);
 641         }
 642         ans.sign = (BOOL)sign;
 643         *res = ans;
 644 }
 645
 646
 647 /*
 648  * Compute ten to the specified power
 649  * This saves some work since the squares of ten are saved.
 650  */
 651 void
 652 ztenpow(long power, ZVALUE *res)
 653 {
 654         long i;
 655         ZVALUE ans;
 656         ZVALUE temp;
 657
 658         if (power <= 0) {
 659                 *res = _one_;
 660                 return;
 661         }
 662         ans = _one_;
 663         _tenpowers_[0] = _ten_;
 664         for (i = 0; power; i++) {
 665                 if (_tenpowers_[i].len == 0) {
 666                         if (i <= TEN_MAX) {
 667                                 zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
 668                         } else {
 669                                 math_error("cannot compute 10^2^(TEN_MAX+1)");
 670                                 /*NOTREACHED*/
 671                         }
 672                 }
 673                 if (power & 0x1) {
 674                         zmul(ans, _tenpowers_[i], &temp);
 675                         zfree(ans);
 676                         ans = temp;
 677                 }
 678                 power /= 2;
 679         }
 680         *res = ans;
 681 }
 682
 683
 684 /*
 685  * Calculate modular inverse suppressing unnecessary divisions.
 686  * This is based on the Euclidean algorithm for large numbers.
 687  * (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
 688  * Returns TRUE if there is no solution because the numbers
 689  * are not relatively prime.
 690  */
 691 BOOL
 692 zmodinv(ZVALUE u, ZVALUE v, ZVALUE *res)
 693 {
 694         FULL    q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
 695         ZVALUE  u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;
 696
 697         v.sign = 0;
 698         if (zisneg(u) || (zrel(u, v) >= 0))
 699                 zmod(u, v, &v3, 0);
 700         else
 701                 zcopy(u, &v3);
 702         zcopy(v, &u3);
 703         u2 = _zero_;
 704         v2 = _one_;
 705
 706         /*
 707          * Loop here while the size of the numbers remain above
 708          * the size of a HALF.  Throughout this loop u3 >= v3.
 709          */
 710         while ((u3.len > 1) && !ziszero(v3)) {
 711                 vh = 0;
 712 #if LONG_BITS == BASEB
 713                 uh =  u3.v[u3.len - 1];
 714                 if (v3.len == u3.len)
 715                         vh = v3.v[v3.len - 1];
 716 #else
 717                 uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
 718                 if ((v3.len + 1) >= u3.len)
 719                         vh = v3.v[v3.len - 1];
 720                 if (v3.len == u3.len)
 721                         vh = (vh << BASEB) + v3.v[v3.len - 2];
 722 #endif
 723                 A = 1;
 724                 B = 0;
 725                 C = 0;
 726                 D = 1;
 727
 728                 /*
 729                  * Calculate successive quotients of the continued fraction
 730                  * expansion using only single precision arithmetic until
 731                  * greater precision is required.
 732                  */
 733                 while ((vh + C) && (vh + D)) {
 734                         q1 = (uh + A) / (vh + C);
 735                         q2 = (uh + B) / (vh + D);
 736                         if (q1 != q2)
 737                                 break;
 738                         T = A - q1 * C;
 739                         A = C;
 740                         C = T;
 741                         T = B - q1 * D;
 742                         B = D;
 743                         D = T;
 744                         T = uh - q1 * vh;
 745                         uh = vh;
 746                         vh = T;
 747                 }
 748
 749                 /*
 750                  * If B is zero, then we made no progress because
 751                  * the calculation requires a very large quotient.
 752                  * So we must do this step of the calculation in
 753                  * full precision
 754                  */
 755                 if (B == 0) {
 756                         zquo(u3, v3, &qz, 0);
 757                         zmul(qz, v2, &tmp1);
 758                         zsub(u2, tmp1, &tmp2);
 759                         zfree(tmp1);
 760                         zfree(u2);
 761                         u2 = v2;
 762                         v2 = tmp2;
 763                         zmul(qz, v3, &tmp1);
 764                         zsub(u3, tmp1, &tmp2);
 765                         zfree(tmp1);
 766                         zfree(u3);
 767                         u3 = v3;
 768                         v3 = tmp2;
 769                         zfree(qz);
 770                         continue;
 771                 }
 772                 /*
 773                  * Apply the calculated A,B,C,D numbers to the current
 774                  * values to update them as if the full precision
 775                  * calculations had been carried out.
 776                  */
 777                 zmuli(u2, (long) A, &tmp1);
 778                 zmuli(v2, (long) B, &tmp2);
 779                 zadd(tmp1, tmp2, &tmp3);
 780                 zfree(tmp1);
 781                 zfree(tmp2);
 782                 zmuli(u2, (long) C, &tmp1);
 783                 zmuli(v2, (long) D, &tmp2);
 784                 zfree(u2);
 785                 zfree(v2);
 786                 u2 = tmp3;
 787                 zadd(tmp1, tmp2, &v2);
 788                 zfree(tmp1);
 789                 zfree(tmp2);
 790                 zmuli(u3, (long) A, &tmp1);
 791                 zmuli(v3, (long) B, &tmp2);
 792                 zadd(tmp1, tmp2, &tmp3);
 793                 zfree(tmp1);
 794                 zfree(tmp2);
 795                 zmuli(u3, (long) C, &tmp1);
 796                 zmuli(v3, (long) D, &tmp2);
 797                 zfree(u3);
 798                 zfree(v3);
 799                 u3 = tmp3;
 800                 zadd(tmp1, tmp2, &v3);
 801                 zfree(tmp1);
 802                 zfree(tmp2);
 803         }
 804
 805         /*
 806          * Here when the remaining numbers become single precision in size.
 807          * Finish the procedure using single precision calculations.
 808          */
 809         if (ziszero(v3) && !zisone(u3)) {
 810                 zfree(u3);
 811                 zfree(v3);
 812                 zfree(u2);
 813                 zfree(v2);
 814                 return TRUE;
 815         }
 816         ui3 = ztofull(u3);
 817         vi3 = ztofull(v3);
 818         zfree(u3);
 819         zfree(v3);
 820         while (vi3) {
 821                 q1 = ui3 / vi3;
 822                 zmuli(v2, (long) q1, &tmp1);
 823                 zsub(u2, tmp1, &tmp2);
 824                 zfree(tmp1);
 825                 zfree(u2);
 826                 u2 = v2;
 827                 v2 = tmp2;
 828                 q2 = ui3 - q1 * vi3;
 829                 ui3 = vi3;
 830                 vi3 = q2;
 831         }
 832         zfree(v2);
 833         if (ui3 != 1) {
 834                 zfree(u2);
 835                 return TRUE;
 836         }
 837         if (zisneg(u2)) {
 838                 zadd(v, u2, res);
 839                 zfree(u2);
 840                 return FALSE;
 841         }
 842         *res = u2;
 843         return FALSE;
 844 }
 845
 846
 847 /*
 848  * Compute the greatest common divisor of a pair of integers.
 849  */
 850 void
 851 zgcd(ZVALUE z1, ZVALUE z2, ZVALUE *res)
 852 {
 853         int h, i, j, k;
 854         LEN len, l, m, n, o, p, q;
 855         HALF u, v, w, x;
 856         HALF *a, *a0, *A, *b, *b0, *B, *c, *d;
 857         FULL f, g;
 858         ZVALUE gcd;
 859         BOOL needw;
 860
 861         if (zisunit(z1) || zisunit(z2)) {
 862                 *res = _one_;
 863                 return;
 864         }
 865         z1.sign = 0;
 866         z2.sign = 0;
 867         if (ziszero(z1) || !zcmp(z1, z2)) {
 868                 zcopy(z2, res);
 869                 return;
 870         }
 871         if (ziszero(z2)) {
 872                 zcopy(z1, res);
 873                 return;
 874         }
 875
 876         o = 0;
 877         while (!(z1.v[o] | z2.v[o])) o++;       /* Count common zero digits */
 878
 879         c = z1.v + o;
 880         d = z2.v + o;
 881
 882         m = z1.len - o;
 883         n = z2.len - o;
 884         u = *c | *d;                    /* Count common zero bits */
 885         v = 1;
 886         p = 0;
 887         while (!(u & v)) {
 888                 v <<= 1;
 889                 p++;
 890         }
 891
 892         while (!*c) {                           /* Removing zero digits */
 893                 c++;
 894                 m--;
 895         }
 896
 897         while (!*d) {
 898                 d++;
 899                 n--;
 900         }
 901
 902
 903         u = *d;                 /* Count zero bits for *d */
 904         v = 1;
 905         q = 0;
 906         while (!(u & v)) {
 907                 v <<= 1;
 908                 q++;
 909         }
 910
 911         a0 = A = alloc(m);
 912         b0 = B = alloc(n);
 913
 914         memcpy(A, c, m * sizeof(HALF));                 /* Copy c[] to A[] */
 915
 916                                 /* Copy d[] to B[], shifting if necessary */
 917         if (q) {
 918                 i = n;
 919                 b = B + n;
 920                 d += n;
 921                 f = 0;
 922                 while (i--) {
 923                         f = f << BASEB | *--d;
 924                         *--b = (HALF) (f >> q);
 925                 }
 926                 if (B[n-1] == 0) n--;
 927         }
 928         else memcpy(B, d, n * sizeof(HALF));
 929
 930         if (n == 1) {           /* One digit case; use Euclid's algorithm */
 931                 n = m;
 932                 b0 = A;
 933                 m = 1;
 934                 a0 = B;
 935                 if (m == 1) {                           /* a has one digit */
 936                         v = *a0;
 937                         if (v > 1) {                    /* Euclid's algorithm */
 938                                 b = b0 + n;
 939                                 i = n;
 940                                 u = 0;
 941                                 while (i--) {
 942                                         f = (FULL) u << BASEB | *--b;
 943                                         u = (HALF) (f % v);
 944                                 }
 945                                 while (u) { w = v % u; v = u; u = w; }
 946                         }
 947                         *b0 = v;
 948                         n = 1;
 949                 }
 950                 len = n + o;
 951                 gcd.v = alloc(len + 1);
 952                 /* Common zero digits */
 953                 if (o) memset(gcd.v, 0, o * sizeof(HALF));
 954                 /* Left shift for common zero bits */
 955                 if (p) {
 956                         i = n;
 957                         f = 0;
 958                         b = b0;
 959                         a = gcd.v + o;
 960                         while (i--) {
 961                                 f = f >> BASEB | (FULL) *b++ << p;
 962                                 *a++ = (HALF) f;
 963                         }
 964                         if (f >>= BASEB) {len++; *a = (HALF) f;}
 965                 } else {
 966                         memcpy(gcd.v + o, b0, n * sizeof(HALF));
 967                 }
 968                 gcd.len = len;
 969                 gcd.sign = 0;
 970                 freeh(A);
 971                 freeh(B);
 972                 *res = gcd;
 973                 return;
 974         }
 975
 976         u = B[n-1];                             /* Bit count for b */
 977         k = (n  - 1) * BASEB;
 978         while (u >>= 1) k++;
 979
 980         needw = TRUE;
 981
 982         w = 0;
 983         j = 0;
 984         while (m) {                             /* START OF MAIN LOOP */
 985            if (m - n < 2 || needw) {
 986                 q = 0;
 987                 u = *a0;
 988                 v = 1;
 989                 while (!(u & v)) {              /* count zero bits for *a0 */
 990                         q++;
 991                         v <<= 1;
 992                 }
 993
 994                 if (q) {                        /* right-justify a */
 995                         a = a0 + m;
 996                         i = m;
 997                         f = 0;
 998                         while (i--) {
 999                                 f = f << BASEB | *--a;
1000                                 *a = (HALF) (f >> q);
1001                         }
1002                         if (!a0[m-1]) m--;      /* top digit vanishes */
1003                 }
1004
1005                 if (m == 1) break;
1006
1007                 u = a0[m-1];
1008                 j = (m - 1) * BASEB;
1009                 while (u >>= 1) j++;            /* counting bits for a */
1010                 h = j - k;
1011                 if (h < 0) {                    /* swapping to get h > 0 */
1012                         l = m;
1013                         m = n;
1014                         n = l;
1015                         a = a0;
1016                         a0 = b0;
1017                         b0 = a;
1018                         k = j;
1019                         h = -h;
1020                         needw = TRUE;
1021                 }
1022                 if (h > 1) {
1023                         if (needw) {            /* find w = minv(*b0, h0) */
1024                                 u = 1;
1025                                 v = *b0;
1026                                 w = 0;
1027                                 x = 1;
1028                                 i = h;
1029                                 while (i-- && x) {
1030                                         if (u & x) { u -= v * x; w |= x;}
1031                                         x <<= 1;
1032                                 }
1033                                 needw = FALSE;
1034                         }
1035                         g = *a0 * w;
1036                         if (h < BASEB) g &= (1 << h) - 1;
1037                         else g &= BASE1;
1038                 }
1039                 else g = 1;
1040            } else
1041                 g = (HALF) *a0 * w;
1042                 a = a0;
1043                 b = b0;
1044                 i = n;
1045                 if (g > 1) {                    /* a - g * b case */
1046                         f = 0;
1047                         while (i--) {
1048                                 f = (FULL) *a - g * *b++ - f;
1049                                 *a++ = (HALF) f;
1050                                 f >>= BASEB;
1051                                 f = -f & BASE1;
1052                         }
1053                         if (f) {
1054                                 i = m - n;
1055                                 while (i-- && f) {
1056                                         f = *a - f;
1057                                         *a++ = (HALF) f;
1058                                         f >>= BASEB;
1059                                         f = -f & BASE1;
1060                                 }
1061                         }
1062                         while (m && !*a0) { /* Removing trailing zeros */
1063                                 m--;
1064                                 a0++;
1065                         }
1066                         if (f) {                /* a - g * b < 0 */
1067                                 while (m > 1 && a0[m-1] == BASE1) m--;
1068                                 *a0 = - *a0;
1069                                 a = a0;
1070                                 i = m;
1071                                 while (--i) {
1072                                         a++;
1073                                         *a = ~*a;
1074                                 }
1075                         }
1076                 } else {                                /* abs(a - b) case */
1077                         while (i && *a++ == *b++) i--;
1078                         q = n - i;
1079                         if (m == n) {           /* a and b same length */
1080                                 if (i) {         /* a not equal to b */
1081                                         while (m && a0[m-1] == b0[m-1]) m--;
1082                                         if (a0[m-1] < b0[m-1]) {
1083                                                 /* Swapping since a < b */
1084                                                 a = a0;
1085                                                 a0 = b0;
1086                                                 b0 = a;
1087                                                 k = j;
1088                                         }
1089                                         a = a0 + q;
1090                                         b = b0 + q;
1091                                         i = m - q;
1092                                         f = 0;
1093                                         while (i--) {
1094                                                 f = (FULL) *a - *b++ - f;
1095                                                 *a++ = (HALF) f;
1096                                                 f >>= BASEB;
1097                                                 f = -f & BASE1;
1098                                         }
1099                                 }
1100                         } else {                /* a has more digits than b */
1101                                 a = a0 + q;
1102                                 b = b0 + q;
1103                                 i = n - q;
1104                                 f = 0;
1105                                 while (i--) {
1106                                         f = (FULL) *a - *b++ - f;
1107                                         *a++ = (HALF) f;
1108                                         f >>= BASEB;
1109                                         f = -f & BASE1;
1110                                 }
1111                                 if (f) { while (!*a) *a++ = BASE1;
1112                                         (*a)--;
1113                                 }
1114                         }
1115                         a0 += q;
1116                         m -= q;
1117                         while (m && !*a0) { /* Removing trailing zeros */
1118                                 m--;
1119                                 a0++;
1120                         }
1121                 }
1122                 while (m && !a0[m-1]) m--;      /* Removing leading zeros */
1123         }
1124         if (m == 1) {                           /* a has one digit */
1125                 v = *a0;
1126                 if (v > 1) {                    /* Euclid's algorithm */
1127                         b = b0 + n;
1128                         i = n;
1129                         u = 0;
1130                         while (i--) {
1131                                 f = (FULL) u << BASEB | *--b;
1132                                 u = (HALF) (f % v);
1133                         }
1134                         while (u) { w = v % u; v = u; u = w; }
1135                 }
1136                 *b0 = v;
1137                 n = 1;
1138         }
1139         len = n + o;
1140         gcd.v = alloc(len + 1);
1141         if (o) memset(gcd.v, 0, o * sizeof(HALF));      /* Common zero digits */
1142         if (p) {                        /* Left shift for common zero bits */
1143                 i = n;
1144                 f = 0;
1145                 b = b0;
1146                 a = gcd.v + o;
1147                 while (i--) {
1148                         f = (FULL) *b++ << p | f;
1149                         *a++ = (HALF) f;
1150                         f >>= BASEB;
1151                 }
1152                 if (f) {
1153                         len++; *a = (HALF) f;
1154                 }
1155         } else {
1156                 memcpy(gcd.v + o, b0, n * sizeof(HALF));
1157         }
1158         gcd.len = len;
1159         gcd.sign = 0;
1160         freeh(A);
1161         freeh(B);
1162         *res = gcd;
1163         return;
1164 }
1165
1166 /*
1167  * Compute the lcm of two integers (least common multiple).
1168  * This is done using the formula:  gcd(a,b) * lcm(a,b) = a * b.
1169  */
1170 void
1171 zlcm(ZVALUE z1, ZVALUE z2, ZVALUE *res)
1172 {
1173         ZVALUE temp1, temp2;
1174
1175         zgcd(z1, z2, &temp1);
1176         zequo(z1, temp1, &temp2);
1177         zfree(temp1);
1178         zmul(temp2, z2, res);
1179         zfree(temp2);
1180 }
1181
1182
1183 /*
1184  * Return whether or not two numbers are relatively prime to each other.
1185  */
1186 BOOL
1187 zrelprime(ZVALUE z1, ZVALUE z2)
1188 {
1189         FULL rem1, rem2;                /* remainders */
1190         ZVALUE rem;
1191         BOOL result;
1192
1193         z1.sign = 0;
1194         z2.sign = 0;
1195         if (ziseven(z1) && ziseven(z2)) /* false if both even */
1196                 return FALSE;
1197         if (zisunit(z1) || zisunit(z2)) /* true if either is a unit */
1198                 return TRUE;
1199         if (ziszero(z1) || ziszero(z2)) /* false if either is zero */
1200                 return FALSE;
1201         if (zistwo(z1) || zistwo(z2))   /* true if either is two */
1202                 return TRUE;
1203         /*
1204          * Try reducing each number by the product of the first few odd primes
1205          * to see if any of them are a common factor.
1206          */
1207         rem1 = zmodi(z1, (FULL)3 * 5 * 7 * 11 * 13);
1208         rem2 = zmodi(z2, (FULL)3 * 5 * 7 * 11 * 13);
1209         if (((rem1 % 3) == 0) && ((rem2 % 3) == 0))
1210                 return FALSE;
1211         if (((rem1 % 5) == 0) && ((rem2 % 5) == 0))
1212                 return FALSE;
1213         if (((rem1 % 7) == 0) && ((rem2 % 7) == 0))
1214                 return FALSE;
1215         if (((rem1 % 11) == 0) && ((rem2 % 11) == 0))
1216                 return FALSE;
1217         if (((rem1 % 13) == 0) && ((rem2 % 13) == 0))
1218                 return FALSE;
1219         /*
1220          * Try a new batch of primes now
1221          */
1222         rem1 = zmodi(z1, (FULL)17 * 19 * 23);
1223         rem2 = zmodi(z2, (FULL)17 * 19 * 23);
1224         if (((rem1 % 17) == 0) && ((rem2 % 17) == 0))
1225                 return FALSE;
1226         if (((rem1 % 19) == 0) && ((rem2 % 19) == 0))
1227                 return FALSE;
1228         if (((rem1 % 23) == 0) && ((rem2 % 23) == 0))
1229                 return FALSE;
1230         /*
1231          * Yuk, we must actually compute the gcd to know the answer
1232          */
1233         zgcd(z1, z2, &rem);
1234         result = zisunit(rem);
1235         zfree(rem);
1236         return result;
1237 }
1238
1239
1240 /*
1241  * Compute the integer floor of the log of an integer to a specified base.
1242  * The signs of the integers and base are ignored.
1243  * Example:  zlog(123456, 10) = 5.
1244  */
1245 long
1246 zlog(ZVALUE z, ZVALUE base)
1247 {
1248         ZVALUE *zp;                     /* current square */
1249         long power;                     /* current power */
1250         ZVALUE temp;                    /* temporary */
1251         ZVALUE squares[32];             /* table of squares of base */
1252
1253         /* ignore signs */
1254
1255         z.sign = 0;
1256         base.sign = 0;
1257
1258         /*
1259          * Make sure that the numbers are nonzero and the base is > 1
1260          */
1261         if (ziszero(z) || ziszero(base) || zisone(base)) {
1262                 math_error("Zero or too small argument argument for zlog!!!");
1263                 /*NOTREACHED*/
1264         }
1265
1266         /*
1267          * Some trivial cases.
1268          */
1269         power = zrel(z, base);
1270         if (power <= 0)
1271                 return (power + 1);
1272
1273         /* base - power of two */
1274         if (zisonebit(base))
1275                 return (zhighbit(z) / zlowbit(base));
1276
1277         /* base = 10 */
1278         if (base.len == 1 && base.v[0] == 10)
1279                 return zlog10(z, NULL);
1280         /*
1281          * Now loop by squaring the base each time, and see whether or
1282          * not each successive square is still smaller than the number.
1283          */
1284         zp = &squares[0];
1285         *zp = base;
1286         while (zp->len * 2 - 1 <= z.len  && zrel(z, *zp) > 0) {
1287                 /* while square not too large */
1288                 zsquare(*zp, zp + 1);
1289                 zp++;
1290         }
1291
1292         /*
1293          * Now back down the squares,
1294          */
1295         power = 0;
1296         for (; zp > squares; zp--) {
1297                 if (zrel(z, *zp) >= 0) {
1298                         zquo(z, *zp, &temp, 0);
1299                         if (power)
1300                                 zfree(z);
1301                         z = temp;
1302                         power++;
1303                 }
1304                 zfree(*zp);
1305                 power <<= 1;
1306         }
1307         if (zrel(z, *zp) >= 0)
1308                 power++;
1309         if (power > 1)
1310                 zfree(z);
1311         return power;
1312 }
1313
1314
1315 /*
1316  * Return the integral log base 10 of a number.
1317  *
1318  * If was_10_power != NULL, then this flag is set to TRUE if the
1319  * value was a power of 10, FALSE otherwise.
1320  */
1321 long
1322 zlog10(ZVALUE z, BOOL *was_10_power)
1323 {
1324         ZVALUE *zp;                     /* current square */
1325         long power;                     /* current power */
1326         ZVALUE temp;                    /* temporary */
1327         ZVALUE pow10;                   /* power of 10 */
1328         FLAG rel;                       /* relationship */
1329         int i;
1330
1331         if (ziszero(z)) {
1332                 math_error("Zero argument argument for zlog10");
1333                 /*NOTREACHED*/
1334         }
1335
1336         /* Ignore sign of z */
1337         z.sign = 0;
1338
1339         /* preload power10 table if missing */
1340         if (power10 == NULL) {
1341                 long v;
1342
1343                 /* determine power10 table size */
1344                 for (v=1, max_power10_exp=0;
1345                      v <= (long)(MAXLONG/10L);
1346                      v *= 10L, ++max_power10_exp) {
1347                 }
1348
1349                 /* create power10 table */
1350                 power10 = malloc(sizeof(long) * (max_power10_exp+1));
1351                 if (power10 == NULL) {
1352                         math_error("cannot malloc power10 table");
1353                         /*NOTREACHED*/
1354                 }
1355
1356                 /* load power10 table */
1357                 for (i=0, v = 1L; i <= max_power10_exp; ++i, v *= 10L) {
1358                         power10[i] = v;
1359                 }
1360         }
1361
1362         /* assume not a power of ten unless we find out otherwise */
1363         if (was_10_power != NULL) {
1364                 *was_10_power = FALSE;
1365         }
1366
1367         /* quick exit for small values */
1368         if (! zgtmaxlong(z)) {
1369                 long value = ztolong(z);
1370
1371                 for (i=0; i <= max_power10_exp; ++i) {
1372                         if (value == power10[i]) {
1373                                 if (was_10_power != NULL) {
1374                                         *was_10_power = TRUE;
1375                                 }
1376                                 return i;
1377                         } else if (value < power10[i]) {
1378                                 return i-1;
1379                         }
1380                 }
1381         }
1382
1383         /*
1384          * Loop by squaring the base each time, and see whether or
1385          * not each successive square is still smaller than the number.
1386          */
1387         zp = &_tenpowers_[0];
1388         *zp = _ten_;
1389         while (((zp->len * 2) - 1) <= z.len) {  /* while square not too large */
1390                 if (zp >= &_tenpowers_[TEN_MAX]) {
1391                         math_error("Maximum storable power of 10 reached!");
1392                         /*NOTREACHED*/
1393                 }
1394                 if (zp[1].len == 0)
1395                         zsquare(*zp, zp + 1);
1396                 zp++;
1397         }
1398
1399         /*
1400          * Now back down the squares, and multiply them together to see
1401          * exactly how many times the base can be raised by.
1402          */
1403         /* find the tenpower table entry < z */
1404         do {
1405                 rel = zrel(*zp, z);
1406                 if (rel == 0) {
1407                         /* quick exit - we match a tenpower entry */
1408                         if (was_10_power != NULL) {
1409                                 *was_10_power = TRUE;
1410                         }
1411                         return (1L << (zp - _tenpowers_));
1412                 }
1413         } while (rel > 0 && --zp >= _tenpowers_);
1414         if (zp < _tenpowers_) {
1415                 math_error("fell off bottom of tenpower table!");
1416                 /*NOTREACHED*/
1417         }
1418
1419         /* the tenpower value is now our starting comparison value */
1420         zcopy(*zp, &pow10);
1421         power = (1L << (zp - _tenpowers_));
1422
1423         /* try to build up a power of 10 from tenpower table entries */
1424         while (--zp >= _tenpowers_) {
1425
1426                 /* try the next lower tenpower value */
1427                 zmul(pow10, *zp, &temp);
1428                 rel = zrel(temp, z);
1429                 if (rel == 0) {
1430                         /* exact power of 10 match */
1431                         power += (1L << (zp - _tenpowers_));
1432                         if (was_10_power != NULL) {
1433                                 *was_10_power = TRUE;
1434                         }
1435                         return power;
1436
1437                 /* ignore this entry if we went too large */
1438                 } else if (rel > 0) {
1439                         zfree(temp);
1440
1441                 /* otherwise increase power and keep going */
1442                 } else {
1443                         power += (1L << (zp - _tenpowers_));
1444                         zfree(pow10);
1445                         pow10 = temp;
1446                 }
1447         }
1448         return power;
1449 }
1450
1451
1452 /*
1453  * Return the number of times that one number will divide another.
1454  * This works similarly to zlog, except that divisions must be exact.
1455  * For example, zdivcount(540, 3) = 3, since 3^3 divides 540 but 3^4 won't.
1456  */
1457 long
1458 zdivcount(ZVALUE z1, ZVALUE z2)
1459 {
1460         long count;             /* number of factors removed */
1461         ZVALUE tmp;             /* ignored return value */
1462
1463         if (ziszero(z1) || ziszero(z2) || zisunit(z2))
1464                 return 0;
1465         count = zfacrem(z1, z2, &tmp);
1466         zfree(tmp);
1467         return count;
1468 }
1469
1470
1471 /*
1472  * Remove all occurrences of the specified factor from a number.
1473  * Also returns the number of factors removed as a function return value.
1474  * Example:  zfacrem(540, 3, &x) returns 3 and sets x to 20.
1475  */
1476 long
1477 zfacrem(ZVALUE z1, ZVALUE z2, ZVALUE *rem)
1478 {
1479         register ZVALUE *zp;            /* current square */
1480         long count;                     /* total count of divisions */
1481         long worth;                     /* worth of current square */
1482         long lowbit;                    /* for zlowbit(z2) */
1483         ZVALUE temp1, temp2, temp3;     /* temporaries */
1484         ZVALUE squares[32];             /* table of squares of factor */
1485
1486         z1.sign = 0;
1487         z2.sign = 0;
1488         /*
1489          * Reject trivial cases.
1490          */
1491         if ((z1.len < z2.len) || (zisodd(z1) && ziseven(z2)) ||
1492                 ziszero(z2) || zisone(z2) ||
1493                 ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))) {
1494                 rem->v = alloc(z1.len);
1495                 rem->len = z1.len;
1496                 rem->sign = 0;
1497                 zcopyval(z1, *rem);
1498                 return 0;
1499         }
1500         /*
1501          * Handle any power of two special.
1502          */
1503         if (zisonebit(z2)) {
1504                 lowbit = zlowbit(z2);
1505                 count = zlowbit(z1) / lowbit;
1506                 rem->v = alloc(z1.len);
1507                 rem->len = z1.len;
1508                 rem->sign = 0;
1509                 zcopyval(z1, *rem);
1510                 zshiftr(*rem, count * lowbit);
1511                 ztrim(rem);
1512                 return count;
1513         }
1514         /*
1515          * See if the factor goes in even once.
1516          */
1517         zdiv(z1, z2, &temp1, &temp2, 0);
1518         if (!ziszero(temp2)) {
1519                 zfree(temp1);
1520                 zfree(temp2);
1521                 rem->v = alloc(z1.len);
1522                 rem->len = z1.len;
1523                 rem->sign = 0;
1524                 zcopyval(z1, *rem);
1525                 return 0;
1526         }
1527         zfree(temp2);
1528         z1 = temp1;
1529         /*
1530          * Now loop by squaring the factor each time, and see whether
1531          * or not each successive square will still divide the number.
1532          */
1533         count = 1;
1534         worth = 1;
1535         zp = &squares[0];
1536         *zp = z2;
1537         while (((zp->len * 2) - 1) <= z1.len) { /* while square not too large */
1538                 zsquare(*zp, &temp1);
1539                 zdiv(z1, temp1, &temp2, &temp3, 0);
1540                 if (!ziszero(temp3)) {
1541                         zfree(temp1);
1542                         zfree(temp2);
1543                         zfree(temp3);
1544                         break;
1545                 }
1546                 zfree(temp3);
1547                 zfree(z1);
1548                 z1 = temp2;
1549                 *++zp = temp1;
1550                 worth *= 2;
1551                 count += worth;
1552         }
1553
1554         /*
1555          * Now back down the list of squares, and see if the lower powers
1556          * will divide any more times.
1557          */
1558         /*
1559          * We prevent the zp pointer from walking behind squares
1560          * by stopping one short of the end and running the loop one
1561          * more time.
1562          *
1563          * We could stop the loop with just zp >= squares, but stopping
1564          * short and running the loop one last time manually helps make
1565          * code checkers such as insure happy.
1566          */
1567         for (; zp > squares; zp--, worth /= 2) {
1568                 if (zp->len <= z1.len) {
1569                         zdiv(z1, *zp, &temp1, &temp2, 0);
1570                         if (ziszero(temp2)) {
1571                                 temp3 = z1;
1572                                 z1 = temp1;
1573                                 temp1 = temp3;
1574                                 count += worth;
1575                         }
1576                         zfree(temp1);
1577                         zfree(temp2);
1578                 }
1579                 if (zp != squares)
1580                         zfree(*zp);
1581         }
1582         /* run the loop manually one last time */
1583         if (zp == squares) {
1584                 if (zp->len <= z1.len) {
1585                         zdiv(z1, *zp, &temp1, &temp2, 0);
1586                         if (ziszero(temp2)) {
1587                                 temp3 = z1;
1588                                 z1 = temp1;
1589                                 temp1 = temp3;
1590                                 count += worth;
1591                         }
1592                         zfree(temp1);
1593                         zfree(temp2);
1594                 }
1595                 if (zp != squares)
1596                         zfree(*zp);
1597         }
1598
1599         *rem = z1;
1600         return count;
1601 }
1602
1603
1604 /*
1605  * Keep dividing a number by the gcd of it with another number until the
1606  * result is relatively prime to the second number.  Returns the number
1607  * of divisions made, and if this is positive, stores result at res.
1608  */
1609 long
1610 zgcdrem(ZVALUE z1, ZVALUE z2, ZVALUE *res)
1611 {
1612         ZVALUE tmp1, tmp2;
1613         long count, onecount;
1614         long sh;
1615
1616         if (ziszero(z1) || ziszero(z2)) {
1617                 math_error("Zero argument in call to zgcdrem!!!");
1618                 /*NOTREACHED*/
1619         }
1620         /*
1621          * Begin by taking the gcd for the first time.
1622          * If the number is already relatively prime, then we are done.
1623          */
1624         z1.sign = 0;
1625         z2.sign = 0;
1626         if (zisone(z2))
1627                 return 0;
1628         if (zisonebit(z2)) {
1629                 sh = zlowbit(z1);
1630                 if (sh == 0)
1631                         return 0;
1632                 zshift(z1, -sh, res);
1633                 return 1 + (sh - 1)/zlowbit(z2);
1634         }
1635         if (zisonebit(z1)) {
1636                 if (zisodd(z2))
1637                         return 0;
1638                 *res = _one_;
1639                 return zlowbit(z1);
1640         }
1641
1642         zgcd(z1, z2, &tmp1);
1643         if (zisunit(tmp1) || ziszero(tmp1))
1644                 return 0;
1645         zequo(z1, tmp1, &tmp2);
1646         count = 1;
1647         z1 = tmp2;
1648         z2 = tmp1;
1649         /*
1650          * Now keep alternately taking the gcd and removing factors until
1651          * the gcd becomes one.
1652          */
1653         while (!zisunit(z2)) {
1654                 onecount = zfacrem(z1, z2, &tmp1);
1655                 if (onecount) {
1656                         count += onecount;
1657                         zfree(z1);
1658                         z1 = tmp1;
1659                 }
1660                 zgcd(z1, z2, &tmp1);
1661                 zfree(z2);
1662                 z2 = tmp1;
1663         }
1664         *res = z1;
1665         return count;
1666 }
1667
1668
1669 /*
1670  * Return the number of digits (base 10) in a number, ignoring the sign.
1671  */
1672 long
1673 zdigits(ZVALUE z1)
1674 {
1675         long count, val;
1676
1677         z1.sign = 0;
1678         if (!zge16b(z1)) {      /* do small numbers ourself */
1679                 count = 1;
1680                 val = 10;
1681                 while (*z1.v >= (HALF)val) {
1682                         count++;
1683                         val *= 10;
1684                 }
1685                 return count;
1686         }
1687         return (zlog10(z1, NULL) + 1);
1688 }
1689
1690
1691 /*
1692  * Return the single digit at the specified decimal place of a number,
1693  * where 0 means the rightmost digit.  Example:  zdigit(1234, 1) = 3.
1694  */
1695 long
1696 zdigit(ZVALUE z1, long n)
1697 {
1698         ZVALUE tmp1, tmp2;
1699         long res;
1700
1701         z1.sign = 0;
1702         if (ziszero(z1) || (n < 0) || (n / BASEDIG >= z1.len))
1703                 return 0;
1704         if (n == 0)
1705                 return zmodi(z1, 10L);
1706         if (n == 1)
1707                 return zmodi(z1, 100L) / 10;
1708         if (n == 2)
1709                 return zmodi(z1, 1000L) / 100;
1710         if (n == 3)
1711                 return zmodi(z1, 10000L) / 1000;
1712         ztenpow(n, &tmp1);
1713         zquo(z1, tmp1, &tmp2, 0);
1714         res = zmodi(tmp2, 10L);
1715         zfree(tmp1);
1716         zfree(tmp2);
1717         return res;
1718 }
1719
1720
1721 /*
1722  * z is to be a nonnegative integer
1723  * If z is the square of a integer stores at dest the square root of z;
1724  *      otherwise stores at z an integer differing from the square root
1725  *      by less than 1.  Returns the sign of the true square root minus
1726  *      the calculated integer.  Type of rounding is determined by
1727  *      rnd as follows: rnd = 0 gives round down, rnd = 1
1728  *      rounds up, rnd = 8 rounds to even integer, rnd = 9 rounds to odd
1729  *      integer, rnd = 16 rounds to nearest integer.
1730  */
1731 FLAG
1732 zsqrt(ZVALUE z, ZVALUE *dest, long rnd)
1733 {
1734         HALF *a, *A, *b, *a0, u;
1735         int i, j, j1, j2, k, k1, m, m0, m1, n, n0, o;
1736         FULL d, e, f, g, h, s, t, x, topbit;
1737         int remsign;
1738         BOOL up, onebit;
1739         ZVALUE sqrt;
1740
1741         if (z.sign) {
1742                 math_error("Square root of negative number");
1743                 /*NOTREACHED*/
1744         }
1745         if (ziszero(z)) {
1746                 *dest = _zero_;
1747                 return 0;
1748         }
1749         m0 = z.len;
1750         o = m0 & 1;
1751         m = m0 + o;             /* m is smallest even number >= z.len */
1752         n0 = n = m / 2;
1753         f = z.v[z.len - 1];
1754         k = 1;
1755         while (f >>= 2)
1756                 k++;
1757         if (!o)
1758                 k += BASEB/2;
1759         j = BASEB - k;
1760         m1 = m;
1761         if (k == BASEB) {
1762                 m1 += 2;
1763                 n0++;
1764         }
1765         A = alloc(m1);
1766         A[m1] = 0;
1767         a0 = A + n0;
1768         memcpy(A, z.v, m0 * sizeof(HALF));
1769         if (o)
1770                 A[m - 1] = 0;
1771         if (n == 1) {
1772                 if (j)
1773                         f = (FULL) A[1] << j | A[0] >> k;
1774                 else
1775                         f = A[1];
1776                 g = (FULL) A[0] << (j + BASEB);
1777                 d = e = topbit = (FULL)1 << (k - 1);
1778         } else {
1779                 if (j)
1780                         f = (FULL) A[m-1] << (j + BASEB) | (FULL) A[m-2] << j |
1781                                 A[m-3] >> k;
1782                 else
1783                         f = (FULL) A[m-1] << BASEB | A[m-2];
1784                 g = (FULL) A[m-3] << (j + BASEB) | (FULL) A[m-4] << j;
1785                 d = e = topbit = (FULL)1 << (BASEB + k - 1);
1786         }
1787
1788         s = (f & topbit);
1789         f <<= 1;
1790         if (g & TOPFULL)
1791                 f++;
1792         g <<= 1;
1793         if (s) {
1794                 f -= 4 * d;
1795                 e = 2 * d - 1;
1796         }
1797         else
1798                 f -= d;
1799         while (d >>= 1) {
1800                 if (!(s | f | g))
1801                         break;
1802                 while (d && (f & topbit) == s) {
1803                         d >>= 1;
1804                         f <<= 1;
1805                         if (g & TOPFULL)
1806                                 f++;
1807                         g <<= 1;
1808                 }
1809                 if (d == 0)
1810                         break;
1811                 if (s)
1812                         f += e + 1;
1813                 else
1814                         f -= e;
1815                 t = f & topbit;
1816                 f <<= 1;
1817                 if (g & TOPFULL)
1818                         f++;
1819                 g <<= 1;
1820                 if (t == 0 && f < d)
1821                         t = topbit;
1822                 f -= d;
1823                 if (s)
1824                         e -= d - !t;
1825                 else
1826                         e += d - (t > 0);
1827                 s = t;
1828         }
1829         if (n0 == 1) {
1830                 A[1] = (HALF)e;
1831                 A[0] = (HALF)f;
1832                 m = 1;
1833                 goto done;
1834         }
1835         if (n0 == 2) {
1836                 A[3] = (HALF)(e >> BASEB);
1837                 A[2] = (HALF)e;
1838                 A[1] = (HALF)(f >> BASEB);
1839                 A[0] = (HALF)f;
1840                 m = 2;
1841                 goto done;
1842         }
1843         u = (HALF)(s ? BASE1 : 0);
1844         if (k < BASEB) {
1845                 A[m1 - 1] = (HALF)(e >> (BASEB - 1));
1846                 A[m1 - 2] = ((HALF)(e << 1) | (HALF)(s > 0));
1847                 A[m1 - 3] = (HALF)(f >> BASEB);
1848                 A[m1 - 4] = (HALF)f;
1849                 m = m1 - 2;
1850                 k1 = k + 1;
1851         } else {
1852                 A[m1 - 1] = 1;
1853                 A[m1 - 2] = (HALF)(e >> (BASEB - 1));
1854                 A[m1 - 3] = ((HALF)(e << 1) | (HALF)(s > 0));
1855                 A[m1 - 4] = u;
1856                 A[m1 - 5] = (HALF)(f >> BASEB);
1857                 A[m1 - 6] = (HALF)f;
1858                 m = m1 - 3;
1859                 k1 = 1;
1860         }
1861         h = e >> k;
1862         onebit = ((e & ((FULL)1 << (k - 1))) ? TRUE : FALSE);
1863         j2 = BASEB - k1;
1864         j1 = BASEB + j2;
1865         while (m > n0) {
1866                 a = A + m - 1;
1867                 if (j2)
1868                         f = (FULL) *a << j1 | (FULL) a[-1] << j2 | a[-2] >> k1;
1869                 else
1870                         f = (FULL) *a << BASEB | a[-1];
1871                 if (u)
1872                         f = ~f;
1873                 x = f / h;
1874                 if (x) {
1875                         if (onebit && x > 2 * (f % h) + 2)
1876                                 x--;
1877                         b = a + 1;
1878                         i = m1 - m;
1879                         a -= i + 1;
1880                         if (u) {
1881                                 f = *a + x * (BASE - x);
1882                                 *a++ = (HALF)f;
1883                                 u = (HALF)(f >> BASEB);
1884                                 while (i--) {
1885                                         f = *a + x * *b++ + u;
1886                                         *a++ = (HALF)f;
1887                                         u = (HALF)(f >> BASEB);
1888                                 }
1889                                 u += *a;
1890                                 x = ~x + !u;
1891                                 if (!(x & TOPHALF))
1892                                         a[1] -= 1;
1893                         } else {
1894                                 f = *a - x * x;
1895                                 *a++ = (HALF)f;
1896                                 u = -(HALF)(f >> BASEB);
1897                                 while (i--) {
1898                                         f = *a - x * *b++ - u;
1899                                         *a++ = (HALF)f;
1900                                         u = -(HALF)(f >> BASEB);
1901                                 }
1902                                 u = *a - u;
1903                                 x = x + u;
1904                                 if (x & TOPHALF)
1905                                         a[1] |= 1;
1906                         }
1907                         *a = ((HALF)(x << 1) | (HALF)(u > 0));
1908                 } else {
1909                         *a = u;
1910                 }
1911                 m--;
1912                 if (*--a == u) {
1913                         while (m > 1 && *--a == u)
1914                                 m--;
1915                 }
1916         }
1917         i = n;
1918         a = a0;
1919         while (i--) {
1920                 *a >>= 1;
1921                 if (a[1] & 1) *a |= TOPHALF;
1922                 a++;
1923         }
1924         s = u;
1925 done:   if (s == 0) {
1926                 while (m > 0 && A[m - 1] == 0)
1927                         m--;
1928                 if (m == 0) {
1929                         remsign = 0;
1930                         sqrt.v = alloc(n);
1931                         sqrt.len = n;
1932                         sqrt.sign = 0;
1933                         memcpy(sqrt.v, a0, n * sizeof(HALF));
1934                         freeh(A);
1935                         *dest = sqrt;
1936                         return remsign;
1937                 }
1938         }
1939         if (rnd & 16) {
1940                 if (s == 0) {
1941                         if (m != n) {
1942                                 up = (m > n);
1943                         } else {
1944                                 i = n;
1945                                 b = a0 + n;
1946                                 a = A + n;
1947                                 while (i > 0 && *--a == *--b)
1948                                         i--;
1949                                 up = (i > 0 && *a > *b);
1950                         }
1951                 } else {
1952                         while (m > 1 && A[m - 1] == BASE1)
1953                                 m--;
1954                         if (m != n) {
1955                                 up = (m < n);
1956                         } else {
1957                                 i = n;
1958                                 b = a0 + n;
1959                                 a = A + n;
1960                                 while (i > 0 && *--a + *--b == BASE1)
1961                                         i--;
1962                                 up = ((FULL) *a + *b >= BASE);
1963                         }
1964                 }
1965         }
1966         else
1967         if (rnd & 8)
1968                 up = (((rnd ^ *a0) & 1) ? TRUE : FALSE);
1969         else
1970                 up = ((rnd & 1) ? TRUE : FALSE);
1971         if (up) {
1972                 remsign = -1;
1973                 i = n;
1974                 a = a0;
1975                 while (i-- && *a == BASE1)
1976                         *a++ = 0;
1977                 if (i >= 0) {
1978                         (*a)++;
1979                 } else {
1980                         n++;
1981                         *a = 1;
1982                 }
1983         } else {
1984                 remsign = 1;
1985         }
1986         sqrt.v = alloc(n);
1987         sqrt.len = n;
1988         sqrt.sign = 0;
1989         memcpy(sqrt.v, a0, n * sizeof(HALF));
1990         freeh(A);
1991         *dest = sqrt;
1992         return remsign;
1993
1994 }
1995
1996 /*
1997  * Take an arbitrary root of a number (to the greatest integer).
1998  * This uses the following iteration to get the Kth root of N:
1999  *      x = ((K-1) * x + N / x^(K-1)) / K
2000  */
2001 void
2002 zroot(ZVALUE z1, ZVALUE z2, ZVALUE *dest)
2003 {
2004         ZVALUE ztry, quo, old, temp, temp2;
2005         ZVALUE k1;                      /* holds k - 1 */
2006         int sign;
2007         long i;
2008         LEN highbit, k;
2009         SIUNION sival;
2010
2011         sign = z1.sign;
2012         if (sign && ziseven(z2)) {
2013                 math_error("Even root of negative number");
2014                 /*NOTREACHED*/
2015         }
2016         if (ziszero(z2) || zisneg(z2)) {
2017                 math_error("Non-positive root");
2018                 /*NOTREACHED*/
2019         }
2020         if (ziszero(z1)) {      /* root of zero */
2021                 *dest = _zero_;
2022                 return;
2023         }
2024         if (zisunit(z2)) {      /* first root */
2025                 zcopy(z1, dest);
2026                 return;
2027         }
2028         if (zge31b(z2)) {       /* humongous root */
2029                 *dest = _one_;
2030                 dest->sign = (BOOL)((HALF)sign);
2031                 return;
2032         }
2033         k = (LEN)ztolong(z2);
2034         highbit = zhighbit(z1);
2035         if (highbit < k) {      /* too high a root */
2036                 *dest = _one_;
2037                 dest->sign = (BOOL)((HALF)sign);
2038                 return;
2039         }
2040         sival.ivalue = k - 1;
2041         k1.v = &sival.silow;
2042         /* ignore Saber-C warning #112 - get ushort from uint */
2043         /*        ok to ignore on name zroot`sival */
2044         k1.len = 1 + (sival.sihigh != 0);
2045         k1.sign = 0;
2046         z1.sign = 0;
2047         /*
2048          * Allocate the numbers to use for the main loop.
2049          * The size and high bits of the final result are correctly set here.
2050          * Notice that the remainder of the test value is rubbish, but this
2051          * is unimportant.
2052          */
2053         highbit = (highbit + k - 1) / k;
2054         ztry.len = (highbit / BASEB) + 1;
2055         ztry.v = alloc(ztry.len);
2056         zclearval(ztry);
2057         ztry.v[ztry.len-1] = ((HALF)1 << (highbit % BASEB));
2058         ztry.sign = 0;
2059         old.v = alloc(ztry.len);
2060         old.len = 1;
2061         zclearval(old);
2062         old.sign = 0;
2063         /*
2064          * Main divide and average loop
2065          */
2066         for (;;) {
2067                 zpowi(ztry, k1, &temp);
2068                 zquo(z1, temp, &quo, 0);
2069                 zfree(temp);
2070                 i = zrel(ztry, quo);
2071                 if (i <= 0) {
2072                         /*
2073                          * Current try is less than or equal to the root since
2074                          * it is less than the quotient. If the quotient is
2075                          * equal to the try, we are all done.  Also, if the
2076                          * try is equal to the old value, we are done since
2077                          * no improvement occurred.  If not, save the improved
2078                          * value and loop some more.
2079                          */
2080                         if ((i == 0) || (zcmp(old, ztry) == 0)) {
2081                                 zfree(quo);
2082                                 zfree(old);
2083                                 ztry.sign = (BOOL)((HALF)sign);
2084                                 zquicktrim(ztry);
2085                                 *dest = ztry;
2086                                 return;
2087                         }
2088                         old.len = ztry.len;
2089                         zcopyval(ztry, old);
2090                 }
2091                 /* average current try and quotient for the new try */
2092                 zmul(ztry, k1, &temp);
2093                 zfree(ztry);
2094                 zadd(quo, temp, &temp2);
2095                 zfree(temp);
2096                 zfree(quo);
2097                 zquo(temp2, z2, &ztry, 0);
2098                 zfree(temp2);
2099         }
2100 }
2101
2102
2103 /*
2104  * Test to see if a number is an exact square or not.
2105  */
2106 BOOL
2107 zissquare(ZVALUE z)
2108 {
2109         long n;
2110         ZVALUE tmp;
2111
2112         /* negative values are never perfect squares */
2113         if (zisneg(z)) {
2114                 return FALSE;
2115         }
2116
2117         /* ignore trailing zero words */
2118         while ((z.len > 1) && (*z.v == 0)) {
2119                 z.len--;
2120                 z.v++;
2121         }
2122
2123         /* zero or one is a perfect square */
2124         if (zisabsleone(z)) {
2125                 return TRUE;
2126         }
2127
2128         /* check mod 4096 values */
2129         if (issq_mod4k[(int)(*z.v & 0xfff)] == 0) {
2130                 return FALSE;
2131         }
2132
2133         /* must do full square root test now */
2134         n = !zsqrt(z, &tmp, 0);
2135         zfree(tmp);
2136         return (n ? TRUE : FALSE);
2137 }