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
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