1#! /usr/bin/bc 2# 3# SPDX-License-Identifier: BSD-2-Clause 4# 5# Copyright (c) 2018-2024 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 30scale = 0 31 32bits = 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. 38bits = divmod(bits, 2, negpow[]) 39 40# Get a bit that will say whether the value should be an exact square. 41bits = divmod(bits, 2, square[]) 42 43# See below. This is to help bias toward small numbers. 44pow = 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. 48bits = divmod(bits, 2, small[]) 49 50# Let's keep raising the power limit by 2^4 when the bit is zero. 51while (!small[0]) 52{ 53 pow += 4 54 bits = divmod(bits, 2, small[]) 55} 56 57limit = 2^pow 58 59# Okay, this is the starting number. 60num = irand(limit) + 1 61 62# Figure out if we should have (more) fractional digits. 63bits = divmod(bits, 2, extra_digits[]) 64 65if (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} 95else 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 115bits = divmod(bits, 2, zero_scale[]) 116 117# Do we want a zero scale? 118if (zero_scale[0]) 119{ 120 print "scale = 0\n" 121} 122else 123{ 124 print "scale = 20\n" 125} 126 127print "sqrt(", num, ")\n" 128 129halt 130