CROTG/ZROTG return incorrect rotations for scaled and unbalanced finite inputs

Start with the shared implementation in interface/zrotg.c, especially the scaled branch and small-cosine branch identified in the report. Use the two supplied C testbenches to reproduce the failures and compare the results with the documented CROTG/ZROTG identities and Reference BLAS control. Done means both Fortran and CBLAS entry points pass the scaled and unbalanced cases in single and double precision.

CROTG and ZROTG return rotations with incorrect magnitudes or sine phase for small, large, and strongly unbalanced finite operands. The returned values fail the documented normalization or second-coordinate annihilation identities.

Versions and environment

• OpenBLAS 0.3.34.dev at 33746423e74d816281b7a13196d8be55cab95bca. Fresh source build with make -j3 DYNAMIC_ARCH=1 NO_LAPACK=1.
• Reference BLAS/CBLAS 3.12.0 at 04b044e020a3560ccfa9988c8a80a1fb7083fc2e. Fresh standalone Fortran controls and a full BLAS/CBLAS build.
• Linux x86-64 on AMD EPYC 9J14, Zen 4, with AVX-512. GCC and GNU Fortran 13.3.0, Ubuntu package 13.3.0-6ubuntu2~24.04.1. GNU Make 4.3, CMake 3.28.3, glibc 2.39.
• Detected Cooperlake and forced Haswell/SkylakeX dispatch. Configured thread limits 1 and 3. These small tests do not establish execution by three workers. Round-to-nearest, without fast-math.
• A read-only check on 2026-10-11 found upstream develop still at the tested commit. No later revision was built.

All four entry points use the shared scalar interface implementation. CPU selection and thread limits do not select another ROTG algorithm.

Expected and observed

The CROTG definition and ZROTG definition require c**2 + conjg(s)*s = 1 and -conj(s)*a + c*b = 0. For nonzero a, with d=sqrt(|a|^2+|b|^2), the formulas are c=|a|/d, s=(a/|a|)*conj(b)/d, and r=(a/|a|)*d.

Routine and input	Expected	Observed in OpenBLAS
CROTG, a=3*2^-100, b=4*2^-100	c≈0.6, s≈0.8, r=5*2^-100	c=0.654653668, s=0.654653668, r=3.61501489e-30
ZROTG, a=3*2^-700, b=4*2^-700	c≈0.6, s≈0.8, r=5*2^-700	c=0.65465367070797709, s=0.6546536707079772, r=8.7118960055897298e-211
CROTG, a=b=i*2^70	positive real s≈1/sqrt(2)	real s=-0.707106769
ZROTG, a=b=i*2^600	positive real s≈1/sqrt(2)	real s=-0.70710678118654746
CROTG, a=2^-30, b=2^50	r=2^50, c=2^-80, s=1 after rounding	r=2^50, c=2^-80, s=2^-30
ZROTG, a=2^-200, b=2^400	r=2^400, c=2^-600, s=1 after rounding	r=2^400, c=2^-600, s=2^-200

The real 3-4-5 cases have c^2+|s|^2≈6/7. The equal-imaginary cases have the wrong sine sign. The strongly unbalanced cases have a near-zero sine where the documented formula rounds to 1. All listed inputs are finite normal values and all expected outputs are in range.

Reproduction

On Linux x86-64, with GCC and GNU Fortran available, build the pin:

git clone https://github.com/OpenMathLib/OpenBLAS.git OpenBLAS
git -C OpenBLAS checkout --detach 33746423e74d816281b7a13196d8be55cab95bca
make -C OpenBLAS -j3 DYNAMIC_ARCH=1 NO_LAPACK=1

Save the following unchanged testbench as rotg-scaled.c:

/* Expected values come from the published ROTG formulas (3-4-5 triangle
   and equal purely imaginary operands), evaluated in long double. */
#include <float.h>
#include <math.h>
#include <stdio.h>

extern void crotg_(float *, float *, float *, float *);
extern void zrotg_(double *, double *, double *, double *);
#ifndef FORTRAN_ONLY
extern void cblas_crotg(void *, void *, float *, void *);
extern void cblas_zrotg(void *, void *, double *, void *);
#endif

static int run_float(int api, int phase) {
    float scale = ldexpf(1, phase ? 70 : -100);
    float a[2] = {phase ? 0 : 3*scale, phase ? scale : 0};
    float b[2] = {phase ? 0 : 4*scale, phase ? scale : 0};
    float c, s[2];
    long double er = phase ? sqrtl(2)*scale : 5.L*scale;
    long double ec = phase ? 1/sqrtl(2) : 3.L/5;
    long double es = phase ? 1/sqrtl(2) : 4.L/5;
    if (api == 0) crotg_(a,b,&c,s);
#ifndef FORTRAN_ONLY
    else cblas_crotg(a,b,&c,s);
#endif
    int ok = fabsl(c-ec) <= 8*FLT_EPSILON*ec &&
             fabsl(s[0]-es) <= 8*FLT_EPSILON*es && s[1] == 0 &&
             (phase ? fabsl(a[1]-er) <= 8*FLT_EPSILON*er : a[0] == (float)er) &&
             (phase ? a[0] : a[1]) == 0;
    printf("%s CROTG %s: c=%.9g s=(%.9g,%.9g) r=(%.9g,%.9g); "
           "expected c=%.9Lg s=(%.9Lg,0) |r|=%.9Lg: %s\n",
           api ? "CBLAS" : "Fortran", phase ? "phase" : "3-4-5",
           c,s[0],s[1],a[0],a[1],ec,es,er,ok ? "PASS" : "FAIL");
    return !ok;
}
static int run_double(int api, int phase) {
    double scale = ldexp(1, phase ? 600 : -700);
    double a[2] = {phase ? 0 : 3*scale, phase ? scale : 0};
    double b[2] = {phase ? 0 : 4*scale, phase ? scale : 0};
    double c, s[2];
    long double er = phase ? sqrtl(2)*scale : 5.L*scale;
    long double ec = phase ? 1/sqrtl(2) : 3.L/5;
    long double es = phase ? 1/sqrtl(2) : 4.L/5;
    if (api == 0) zrotg_(a,b,&c,s);
#ifndef FORTRAN_ONLY
    else cblas_zrotg(a,b,&c,s);
#endif
    int ok = fabsl(c-ec) <= 8*DBL_EPSILON*ec &&
             fabsl(s[0]-es) <= 8*DBL_EPSILON*es && s[1] == 0 &&
             (phase ? fabsl(a[1]-er) <= 8*DBL_EPSILON*er : a[0] == (double)er) &&
             (phase ? a[0] : a[1]) == 0;
    printf("%s ZROTG %s: c=%.17g s=(%.17g,%.17g) r=(%.17g,%.17g); "
           "expected c=%.17Lg s=(%.17Lg,0) |r|=%.17Lg: %s\n",
           api ? "CBLAS" : "Fortran", phase ? "phase" : "3-4-5",
           c,s[0],s[1],a[0],a[1],ec,es,er,ok ? "PASS" : "FAIL");
    return !ok;
}
int main(void) {
    int failures=0;
    int apis=1;
#ifndef FORTRAN_ONLY
    apis=2;
#endif
    for (int api=0;api<apis;api++)
        for (int phase=0;phase<2;phase++)
            failures += run_float(api,phase)+run_double(api,phase);
    return failures ? 1 : 0;
}

Save this second unchanged testbench as rotg-unbalanced.c:

/* The exact ROTG formulas round to r=b, c=a/b and s=1 for these powers
   of two: the relative correction from (a/b)^2 is far below one ulp. */
#include <math.h>
#include <stdio.h>
extern void crotg_(float *, float *, float *, float *);
extern void zrotg_(double *, double *, double *, double *);
#ifndef FORTRAN_ONLY
extern void cblas_crotg(void *, void *, float *, void *);
extern void cblas_zrotg(void *, void *, double *, void *);
#endif
int main(void) {
    int failures=0,apis=1;
#ifndef FORTRAN_ONLY
    apis=2;
#endif
    for (int api=0;api<apis;api++) {
        float a[2]={0x1p-30f,0},b[2]={0x1p50f,0},c,s[2];
        double da[2]={0x1p-200,0},db[2]={0x1p400,0},dc,ds[2];
        if (!api) { crotg_(a,b,&c,s); zrotg_(da,db,&dc,ds); }
#ifndef FORTRAN_ONLY
        else { cblas_crotg(a,b,&c,s); cblas_zrotg(da,db,&dc,ds); }
#endif
        int ok=a[0]==0x1p50f && a[1]==0 && c==0x1p-80f && s[0]==1 && s[1]==0;
        int dok=da[0]==0x1p400 && da[1]==0 && dc==0x1p-600 && ds[0]==1 && ds[1]==0;
        printf("%s CROTG: r=%a c=%a s=(%a,%a); expected r=%a c=%a s=(1,0): %s\n",
               api?"CBLAS":"Fortran",a[0],c,s[0],s[1],0x1p50,0x1p-80,ok?"PASS":"FAIL");
        printf("%s ZROTG: r=%a c=%a s=(%a,%a); expected r=%a c=%a s=(1,0): %s\n",
               api?"CBLAS":"Fortran",da[0],dc,ds[0],ds[1],0x1p400,0x1p-600,dok?"PASS":"FAIL");
        failures+=!ok+!dok;
    }
    return failures?1:0;
}

Compile and run both from the directory containing the source files and OpenBLAS/:

for name in rotg-scaled rotg-unbalanced; do
  gcc -std=c11 -O2 -Wall -Wextra "$name.c" \
    -L"$PWD/OpenBLAS" -Wl,-rpath,"$PWD/OpenBLAS" -lopenblas -lm -o "$name"
  env -u LD_LIBRARY_PATH -u OPENBLAS_CORETYPE \
    OPENBLAS_NUM_THREADS=1 OPENBLAS_VERBOSE=2 "./$name"
done

rotg-scaled prints eight FAIL results and exits 1. rotg-unbalanced prints four FAIL results and exits 1. Both testbenches call the Fortran and CBLAS entry points in single and double precision.

For the standalone Reference BLAS control:

git clone --branch v3.12.0 --depth 1 https://github.com/Reference-LAPACK/lapack.git lapack-3.12.0
git -C lapack-3.12.0 rev-parse HEAD
# Expected commit: 04b044e020a3560ccfa9988c8a80a1fb7083fc2e
mkdir refcheck
gfortran -shared -fPIC -O0 \
  lapack-3.12.0/BLAS/SRC/crotg.f90 \
  lapack-3.12.0/BLAS/SRC/zrotg.f90 \
  -o refcheck/libreference.so
gcc -std=c11 -O2 -Wall -Wextra -DFORTRAN_ONLY rotg-scaled.c \
  -L"$PWD/refcheck" -Wl,-rpath,"$PWD/refcheck" -lreference -lm -o reference-rotg-scaled
env -u LD_LIBRARY_PATH ./reference-rotg-scaled
gcc -std=c11 -O2 -Wall -Wextra -DFORTRAN_ONLY rotg-unbalanced.c \
  -L"$PWD/refcheck" -Wl,-rpath,"$PWD/refcheck" -lreference -lm -o reference-rotg-unbalanced
env -u LD_LIBRARY_PATH ./reference-rotg-unbalanced

The reference programs exit 0 with four and two PASS results respectively. The separately built full Reference BLAS/CBLAS control also passes all eight and four checks.

Source and scope

The shared scaled branch takes square roots when forming f2 and g2, then uses them as squared magnitudes. Its real sine expression subtracts the imaginary-component product where a*conj(b) requires addition. The small-cosine branch uses conj(b)*|a|^2/adahsq instead of the complex a in that product. This source reading explains the examples. No fix was tested.

Both original testbenches fail under each of the six dispatch/thread-limit combinations above. An independent 500-digit formula check reproduces these cases and additional imaginary unbalanced inputs. Ordinary (a,b)=(3,4), zero-A, and zero-B controls pass. Other CPU architectures, subnormal inputs, NaN/Inf inputs, and complete application workloads were not tested for this report.

Related history: #4130 requests reference safe scaling, originally for real DROTG. #4909 covers an uninitialized denominator fixed before this pin. #4289 changes workspace allocation and contains the small-cosine expression without reporting its arithmetic error. Read-only searches on 2026-10-11 found no exact match for the magnitude, phase, or small-cosine examples above.

The example testbench code can be used under CC0.

Found by ChipLab's AI verification agents and independently reproduced.

ChipLab ref: CLAB-2026-219