Bug Summary

File:src/lib/libm/src/k_rem_pio2.c
Warning:line 360, column 19
The left operand of '==' is a garbage value due to array index out of bounds

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple amd64-unknown-openbsd7.0 -analyze -disable-free -disable-llvm-verifier -discard-value-names -main-file-name k_rem_pio2.c -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 1 -fhalf-no-semantic-interposition -mframe-pointer=all -relaxed-aliasing -fno-rounding-math -mconstructor-aliases -munwind-tables -target-cpu x86-64 -target-feature +retpoline-indirect-calls -target-feature +retpoline-indirect-branches -tune-cpu generic -debugger-tuning=gdb -fcoverage-compilation-dir=/usr/src/lib/libm/obj -resource-dir /usr/local/lib/clang/13.0.0 -include namespace.h -I /usr/src/lib/libm/arch/amd64 -I /usr/src/lib/libm/src -I /usr/src/lib/libm/src/ld80 -I /usr/src/lib/libm/hidden -D PIC -internal-isystem /usr/local/lib/clang/13.0.0/include -internal-externc-isystem /usr/include -O2 -fdebug-compilation-dir=/usr/src/lib/libm/obj -ferror-limit 19 -fwrapv -D_RET_PROTECTOR -ret-protector -fgnuc-version=4.2.1 -vectorize-loops -vectorize-slp -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-valloc -fno-builtin-free -fno-builtin-strdup -fno-builtin-strndup -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /home/ben/Projects/vmm/scan-build/2022-01-12-194120-40624-1 -x c /usr/src/lib/libm/src/k_rem_pio2.c
1/* @(#)k_rem_pio2.c 5.1 93/09/24 */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * __kernel_rem_pio2(x,y,e0,nx,prec)
15 * double x[],y[]; int e0,nx,prec;
16 *
17 * __kernel_rem_pio2 return the last three digits of N with
18 * y = x - N*pi/2
19 * so that |y| < pi/2.
20 *
21 * The method is to compute the integer (mod 8) and fraction parts of
22 * (2/pi)*x without doing the full multiplication. In general we
23 * skip the part of the product that are known to be a huge integer (
24 * more accurately, = 0 mod 8 ). Thus the number of operations are
25 * independent of the exponent of the input.
26 *
27 * (2/pi) is represented by an array of 24-bit integers in ipio2[].
28 *
29 * Input parameters:
30 * x[] The input value (must be positive) is broken into nx
31 * pieces of 24-bit integers in double precision format.
32 * x[i] will be the i-th 24 bit of x. The scaled exponent
33 * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
34 * match x's up to 24 bits.
35 *
36 * Example of breaking a double positive z into x[0]+x[1]+x[2]:
37 * e0 = ilogb(z)-23
38 * z = scalbn(z,-e0)
39 * for i = 0,1,2
40 * x[i] = floor(z)
41 * z = (z-x[i])*2**24
42 *
43 *
44 * y[] output result in an array of double precision numbers.
45 * The dimension of y[] is:
46 * 24-bit precision 1
47 * 53-bit precision 2
48 * 64-bit precision 2
49 * 113-bit precision 3
50 * The actual value is the sum of them. Thus for 113-bit
51 * precison, one may have to do something like:
52 *
53 * long double t,w,r_head, r_tail;
54 * t = (long double)y[2] + (long double)y[1];
55 * w = (long double)y[0];
56 * r_head = t+w;
57 * r_tail = w - (r_head - t);
58 *
59 * e0 The exponent of x[0]. Must be <= 16360 or you need to
60 * expand the ipio2 table.
61 *
62 * nx dimension of x[]
63 *
64 * prec an integer indicating the precision:
65 * 0 24 bits (single)
66 * 1 53 bits (double)
67 * 2 64 bits (extended)
68 * 3 113 bits (quad)
69 *
70 * External function:
71 * double scalbn(), floor();
72 *
73 *
74 * Here is the description of some local variables:
75 *
76 * jk jk+1 is the initial number of terms of ipio2[] needed
77 * in the computation. The recommended value is 2,3,4,
78 * 6 for single, double, extended,and quad.
79 *
80 * jz local integer variable indicating the number of
81 * terms of ipio2[] used.
82 *
83 * jx nx - 1
84 *
85 * jv index for pointing to the suitable ipio2[] for the
86 * computation. In general, we want
87 * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
88 * is an integer. Thus
89 * e0-3-24*jv >= 0 or (e0-3)/24 >= jv
90 * Hence jv = max(0,(e0-3)/24).
91 *
92 * jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
93 *
94 * q[] double array with integral value, representing the
95 * 24-bits chunk of the product of x and 2/pi.
96 *
97 * q0 the corresponding exponent of q[0]. Note that the
98 * exponent for q[i] would be q0-24*i.
99 *
100 * PIo2[] double precision array, obtained by cutting pi/2
101 * into 24 bits chunks.
102 *
103 * f[] ipio2[] in floating point
104 *
105 * iq[] integer array by breaking up q[] in 24-bits chunk.
106 *
107 * fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
108 *
109 * ih integer. If >0 it indicates q[] is >= 0.5, hence
110 * it also indicates the *sign* of the result.
111 *
112 */
113
114
115/*
116 * Constants:
117 * The hexadecimal values are the intended ones for the following
118 * constants. The decimal values may be used, provided that the
119 * compiler will convert from decimal to binary accurately enough
120 * to produce the hexadecimal values shown.
121 */
122
123#include <float.h>
124#include <math.h>
125
126#include "math_private.h"
127
128static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
129
130/*
131 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
132 *
133 * integer array, contains the (24*i)-th to (24*i+23)-th
134 * bit of 2/pi after binary point. The corresponding
135 * floating value is
136 *
137 * ipio2[i] * 2^(-24(i+1)).
138 *
139 * NB: This table must have at least (e0-3)/24 + jk terms.
140 * For quad precision (e0 <= 16360, jk = 6), this is 686.
141 */
142static const int32_t ipio2[] = {
1430xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
1440x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
1450x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
1460xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
1470x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
1480x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
1490x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
1500xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
1510x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
1520x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
1530x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
154
155#if LDBL_MAX_EXP16384 > 1024
156#if LDBL_MAX_EXP16384 > 16384
157#error "ipio2 table needs to be expanded"
158#endif
1590x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
1600xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
1610xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
1620xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
1630x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
1640x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
1650xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
1660xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
1670xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
1680xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
1690x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
1700xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
1710x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
1720x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
1730xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
1740xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
1750xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
1760x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
1770xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
1780x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
1790xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
1800x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
1810x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
1820x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
1830xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
1840x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
1850x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
1860xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
1870x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
1880x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
1890x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
1900x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
1910x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
1920x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
1930xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
1940x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
1950xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
1960xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
1970xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
1980xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
1990x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
2000x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
2010x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
2020xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
2030x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
2040x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
2050x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
2060xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
2070x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
2080xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
2090xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
2100x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
2110x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
2120x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
2130xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
2140x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
2150x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
2160xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
2170x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
2180xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
2190xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
2200x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
2210xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
2220x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
2230xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
2240x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
2250x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
2260x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
2270xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
2280x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
2290xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
2300x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
2310xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
2320x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
2330x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
2340xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
2350x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
2360xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
2370x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
2380x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
2390x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
2400x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
2410xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
2420xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
2430x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
2440xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
2450x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
2460xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
2470xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
2480x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
2490xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
2500x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
2510x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
2520x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
2530xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
2540xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
2550x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
2560x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
2570xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
2580x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
2590x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
2600x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
2610x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
2620x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
263#endif
264
265};
266
267static const double PIo2[] = {
268 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
269 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
270 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
271 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
272 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
273 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
274 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
275 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
276};
277
278static const double
279zero = 0.0,
280one = 1.0,
281two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
282twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
283
284int
285__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec)
286{
287 int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
288 double z,fw,f[20],fq[20],q[20];
289
290 /* initialize jk*/
291 jk = init_jk[prec];
292 jp = jk;
293
294 /* determine jx,jv,q0, note that 3>q0 */
295 jx = nx-1;
296 jv = (e0-3)/24; if(jv<0) jv=0;
1
Assuming 'jv' is >= 0
→
2
←
Taking false branch
→
297 q0 = e0-24*(jv+1);
298
299 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
300 j = jv-jx; m = jx+jk;
301 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
3
←
Assuming 'i' is > 'm'
→
4
←
Loop condition is false. Execution continues on line 304
→
302
303 /* compute q[0],q[1],...q[jk] */
304 for (i=0;i<=jk;i++) {
5
←
Assuming 'i' is <= 'jk'
→
6
←
Loop condition is true. Entering loop body
→
9
←
Assuming 'i' is > 'jk'
→
10
←
Loop condition is false. Execution continues on line 308
→
305 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
7
←
Assuming 'j' is > 'jx'
→
8
←
Loop condition is false. Execution continues on line 305
→
306 }
307
308 jz = jk;
309recompute:
310 /* distill q[] into iq[] reversingly */
311 for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
11
←
Loop condition is false. Execution continues on line 318
→
312 fw = (double)((int32_t)(twon24* z));
313 iq[i] = (int32_t)(z-two24*fw);
314 z = q[j-1]+fw;
315 }
316
317 /* compute n */
318 z = scalbn(z,q0); /* actual value of z */
319 z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
320 n = (int32_t) z;
321 z -= (double)n;
322 ih = 0;
323 if(q0>0) { /* need iq[jz-1] to determine n */
12
←
Assuming 'q0' is <= 0
→
13
←
Taking false branch
→
324 i = (iq[jz-1]>>(24-q0)); n += i;
325 iq[jz-1] -= i<<(24-q0);
326 ih = iq[jz-1]>>(23-q0);
327 }
328 else if(q0==0) ih = iq[jz-1]>>23;
14
←
Assuming 'q0' is not equal to 0
→
15
←
Taking false branch
→
329 else if(z>=0.5) ih=2;
16
←
Assuming the condition is false
→
17
←
Taking false branch
→
330
331 if(ih
17.1
'ih' is <= 0
>0) { /* q > 0.5 */
18
←
Taking false branch
→
332 n += 1; carry = 0;
333 for(i=0;i<jz ;i++) { /* compute 1-q */
334 j = iq[i];
335 if(carry==0) {
336 if(j!=0) {
337 carry = 1; iq[i] = 0x1000000- j;
338 }
339 } else iq[i] = 0xffffff - j;
340 }
341 if(q0>0) { /* rare case: chance is 1 in 12 */
342 switch(q0) {
343 case 1:
344 iq[jz-1] &= 0x7fffff; break;
345 case 2:
346 iq[jz-1] &= 0x3fffff; break;
347 }
348 }
349 if(ih==2) {
350 z = one - z;
351 if(carry!=0) z -= scalbn(one,q0);
352 }
353 }
354
355 /* check if recomputation is needed */
356 if(z==zero) {
19
←
Assuming 'z' is equal to 'zero'
→
20
←
Taking true branch
→
357 j = 0;
358 for (i=jz-1;i>=jk;i--) j |= iq[i];
21
←
Loop condition is false. Execution continues on line 359
→
359 if(j
21.1
'j' is equal to 0
==0) { /* need recomputation */
22
←
Taking true branch
→
360 for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */
23
←
The left operand of '==' is a garbage value due to array index out of bounds
361
362 for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */
363 f[jx+i] = (double) ipio2[jv+i];
364 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
365 q[i] = fw;
366 }
367 jz += k;
368 goto recompute;
369 }
370 }
371
372 /* chop off zero terms */
373 if(z==0.0) {
374 jz -= 1; q0 -= 24;
375 while(iq[jz]==0) { jz--; q0-=24;}
376 } else { /* break z into 24-bit if necessary */
377 z = scalbn(z,-q0);
378 if(z>=two24) {
379 fw = (double)((int32_t)(twon24*z));
380 iq[jz] = (int32_t)(z-two24*fw);
381 jz += 1; q0 += 24;
382 iq[jz] = (int32_t) fw;
383 } else iq[jz] = (int32_t) z ;
384 }
385
386 /* convert integer "bit" chunk to floating-point value */
387 fw = scalbn(one,q0);
388 for(i=jz;i>=0;i--) {
389 q[i] = fw*(double)iq[i]; fw*=twon24;
390 }
391
392 /* compute PIo2[0,...,jp]*q[jz,...,0] */
393 for(i=jz;i>=0;i--) {
394 for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
395 fq[jz-i] = fw;
396 }
397
398 /* compress fq[] into y[] */
399 switch(prec) {
400 case 0:
401 fw = 0.0;
402 for (i=jz;i>=0;i--) fw += fq[i];
403 y[0] = (ih==0)? fw: -fw;
404 break;
405 case 1:
406 case 2:
407 fw = 0.0;
408 for (i=jz;i>=0;i--) fw += fq[i];
409 STRICT_ASSIGN(double,fw,fw)((fw) = (fw));
410 y[0] = (ih==0)? fw: -fw;
411 fw = fq[0]-fw;
412 for (i=1;i<=jz;i++) fw += fq[i];
413 y[1] = (ih==0)? fw: -fw;
414 break;
415 case 3: /* painful */
416 for (i=jz;i>0;i--) {
417 fw = fq[i-1]+fq[i];
418 fq[i] += fq[i-1]-fw;
419 fq[i-1] = fw;
420 }
421 for (i=jz;i>1;i--) {
422 fw = fq[i-1]+fq[i];
423 fq[i] += fq[i-1]-fw;
424 fq[i-1] = fw;
425 }
426 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
427 if(ih==0) {
428 y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
429 } else {
430 y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
431 }
432 }
433 return n&7;
434}