lib/number/sqrt.cc

Go to the documentation of this file.

//
// canola - canon canola 1614p emulator
// Copyright (C) 2011, 2012 Peter Miller
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License, version 3, as
// published by the Free Software Foundation.
//
// This program 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
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program. If not, see <http://www.gnu.org/licenses/>.
//

#include <lib/ac/assert.h>

#include <lib/number.h>


number
number::sqrt(void)
    const
{
    if (is_zero())
        return number();

    enum { blue_shift = 15 };
    assert(blue_shift < 16);

    //
    // Instruction Manual, p. 8
    //
    // "Square Root Key
    //
    // "Used for finding square roots.  When this key is depressed the
    // indicated value is square rooted.  The result of the square root
    // operation is indicated by a positive value regardless of positive or
    // negative value.  It is not necessary to depress the [+=] key after
    // this key has been depressed."
    //
    number_z a(value.shift_left(2 * blue_shift - decimal));

    //
    // We are going to use a Newton-Raphson approximation
    // with quadratic convergence.
    //
    // FIXME: How about Halley's Method, with cubic convergence?
    //        The 1614P almost certainly did N-R in hardware.
    //
    number_z epsilon(1000u);
    number_z x = number_z::one().shift_left(blue_shift + decimal / 2);
    for (;;)
    {
        number_z x2 = (x + a / x) / 2u;
        if ((x2 < x ? x - x2 : x2 - x) < epsilon)
            break;
        x = x2;
    }

    //
    // Because it is a square root, the number of digits to the left of
    // the decimal point always reduces or stays the same.
    // This means we will never get an overflow.
    //
    assert(x < number_z::one().shift_left(16 + blue_shift));

    unsigned x_dec = blue_shift;
    truncate_floating_operation(x, x_dec);

    // <question>
    // Does sqrt strip redundant zero decimal places to the right of the
    // decimal point?  Does [9] [sqrt] display "3." or "3.0000000000000000"?
    // </question>
    strip_redundant_decimal_places(x, x_dec);

    return number(x, x_dec);
}


// vim: set ts=8 sw=4 et :