contrib/bc/scripts/sqrt_random.bc

   1 #! /usr/bin/bc
   2 #
   3 # SPDX-License-Identifier: BSD-2-Clause
   4 #
   5 # Copyright (c) 2018-2023 Gavin D. Howard and contributors.
   6 #
   7 # Redistribution and use in source and binary forms, with or without
   8 # modification, are permitted provided that the following conditions are met:
   9 #
  10 # * Redistributions of source code must retain the above copyright notice, this
  11 #   list of conditions and the following disclaimer.
  12 #
  13 # * Redistributions in binary form must reproduce the above copyright notice,
  14 #   this list of conditions and the following disclaimer in the documentation
  15 #   and/or other materials provided with the distribution.
  16 #
  17 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18 # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19 # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20 # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  21 # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22 # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23 # SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24 # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25 # CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26 # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27 # POSSIBILITY OF SUCH DAMAGE.
  28 #
  29
  30 scale = 0
  31
  32 bits = rand()
  33
  34 # This extracts a bit and takes it out of the original value.
  35 #
  36 # Here, I am getting a bit to say whether we should have a value that is less
  37 # than 1.
  38 bits = divmod(bits, 2, negpow[])
  39
  40 # Get a bit that will say whether the value should be an exact square.
  41 bits = divmod(bits, 2, square[])
  42
  43 # See below. This is to help bias toward small numbers.
  44 pow = 4
  45
  46 # I want to bias toward small numbers, so let's give a 50 percent chance to
  47 # values below 16 or so.
  48 bits = divmod(bits, 2, small[])
  49
  50 # Let's keep raising the power limit by 2^4 when the bit is zero.
  51 while (!small[0])
  52 {
  53         pow += 4
  54         bits = divmod(bits, 2, small[])
  55 }
  56
  57 limit = 2^pow
  58
  59 # Okay, this is the starting number.
  60 num = irand(limit) + 1
  61
  62 # Figure out if we should have (more) fractional digits.
  63 bits = divmod(bits, 2, extra_digits[])
  64
  65 if (square[0])
  66 {
  67         # Okay, I lied. If we need a perfect square, square now.
  68         num *= num
  69
  70         # If we need extra digits, we need to multiply by an even power of 10.
  71         if (extra_digits[0])
  72         {
  73                 extra = (irand(8) + 1) * 2
  74         }
  75         else
  76         {
  77                 extra = 0
  78         }
  79
  80         # If we need a number less than 1, just take the inverse, which will still
  81         # be a perfect square.
  82         if (negpow[0])
  83         {
  84                 scale = length(num) + 5
  85                 num = 1/num
  86                 scale = 0
  87
  88                 num >>= extra
  89         }
  90         else
  91         {
  92                 num <<= extra
  93         }
  94 }
  95 else
  96 {
  97         # Get this for later.
  98         l = length(num)
  99
 100         # If we need extra digits.
 101         if (extra_digits[0])
 102         {
 103                 # Add up to 32 decimal places.
 104                 num += frand(irand(32) + 1)
 105         }
 106
 107         # If we need a value less than 1...
 108         if (negpow[0])
 109         {
 110                 # Move right until the number is
 111                 num >>= l
 112         }
 113 }
 114
 115 bits = divmod(bits, 2, zero_scale[])
 116
 117 # Do we want a zero scale?
 118 if (zero_scale[0])
 119 {
 120         print "scale = 0\n"
 121 }
 122 else
 123 {
 124         print "scale = 20\n"
 125 }
 126
 127 print "sqrt(", num, ")\n"
 128
 129 halt